Multivariate Variational Mode Decomposition

TL;DR

This paper introduces a multivariate extension of variational mode decomposition (VMD) that effectively extracts joint oscillatory modes from multichannel data, demonstrated on synthetic and real-world biomedical signals.

Contribution

It develops a novel multivariate VMD model based on a joint frequency component, optimized via ADMM, to extract multichannel modes with common frequencies, advancing multivariate signal analysis.

Findings

Successfully separates joint oscillatory modes in multivariate data

Effective in real-world EEG and cardiotocographic signal analysis

Outperforms traditional methods in extracting multichannel oscillations

Abstract

In this paper, a generic extension of variational mode decomposition (VMD) algorithm for multivariate or multichannel data sets is presented. We first define a model for multivariate modulated oscillations that is based on the presence of a joint or common frequency component among all channels of input data. Using that model for multivariate oscillations, we construct a variational optimization problem that aims to extract an ensemble of band-limited modes containing inherent multivariate modulated oscillations present in multivariate input signal. The cost function to be minimized is the sum of bandwidths of all signal modes across all input data channels, which is a generic extension of the cost function used in standard VMD to multivariate data. Minimization of the resulting variational model is achieved through the alternate direction method of multipliers (ADMM) approach. That yields an optimal set of multivariate modes in terms of narrow bandwidth and corresponding center frequencies that are assumed to be commonly present among all channels of a multivariate modulated oscillation. We demonstrate the effectiveness of the proposed method through results obtained from extensive simulations involving test (synthetic) and real world multivariate data sets. Specifically, we focus on the ability of the proposed method to yield joint oscillatory modes in multivariate data which is a prerequisite in many real world applications involving nonstationary multivariate data. We also highlight the utility of the proposed method in two real world applications which include the separation of alpha rhythms in multivariate electroencephalogram (EEG) data and the decomposition of bivariate cardiotocographic signals that consist of fetal heart rate and maternal uterine contraction (FHR-UC) as its two channels.

Equations

$$x(t)=\sum _{k=1}^K{c}_k(t)+r(t),$$

$$c_k(t)=a_k(t)\cos ({\varphi }_k(t)).$$

$$u(t)=a(t)\cos (\varphi (t)),$$

$$u_{+}(t)=u(t)+j\mathcal{H}u(t)=a(t)e^{j\varphi (t)},$$

$$\mathcal{H}u(t)=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{u(t^{\prime })}{t-t^{\prime }}d{t}^{\prime },$$

$${\hat{u}}_{+}(\omega )={\int }_{-\infty }^{\infty }{u}_{+}(t)e^{-j\omega t}dt=1+\text{sgn}(\omega )\hat{u}(\omega ).$$

$$u(t)=\mathcal{R}\{u_{+}(t)\}.$$

$$\mathbf{u}(t)=\left[\begin{array}{c}{u}_1(t)\\ u_2(t)\\ ⋮\\ u_c(t)\end{array}\right]=\left[\begin{array}{c}{a}_1(t)\cos ({\varphi }_1(t))\\ a_2(t)\cos ({\varphi }_2(t))\\ ⋮\\ a_c(t)\cos ({\varphi }_c(t))\end{array}\right].$$

$${\mathbf{u}}_{+}(t)$$

$$=$$

$$\mathbf{u}(t)=\mathcal{R}\{{\mathbf{u}}_{+}(t)\}.$$

$${\mathbf{u}}_{+}(t)=||{\mathbf{u}}_{+}(t)||e^{j\varphi (t)}.$$

$$x(t)=\sum _{k=1}^K{u}_k(t).$$

$$\underset{\{u_k\},\{{\omega }_k\}}{\mathrm{minimize}}$$

$$\sum _k{u}_k(t)=x(t),$$

$$\displaystyle\mathcal{L}\left({\{u_{k}\},\{\omega_{k}\},\lambda}\right)=\alpha\ {\sum}_{k}\Bigg{\|}\partial_{t}\Big{[}u_{+}^{k}(t)e^{-j\omega_{k}t}\Big{]}\Bigg{\|}^{2}_{2}$$

$$\displaystyle+\Bigg{\|}x(t)-\sum_{k}u_{k}(t)\Bigg{\|}_{2}^{2}+\Bigg{\langle}\lambda(t),x(t)-\sum_{k}u_{k}(t)\Bigg{\rangle}.$$

$$u_k^{n+1}\leftarrow \underset{u_k}{\text{arg min }}\mathcal{L}\left(\{u_{i<k}^{n+1}\},\{u_{i\ge k}^n\},\{{\omega }_i^n\},{\lambda }^n\right).$$

$${\hat{u}}_k^{n+1}(\omega )=\frac{\hat{x}(\omega )-{\sum }_{i\ne k}{\hat{u}}_i(\omega )+\frac{\hat{\lambda }(\omega )}{2}}{1+2\alpha (\omega -{\omega }_k{)}^2}.$$

$${\omega }_k^{n+1}\leftarrow \underset{{\omega }_k}{\text{arg min }}\mathcal{L}\left(\{u_i^{n+1}\},\{{\omega }_{i<k}^{n+1}\}\{{\omega }_{i\ge k}^n\},{\lambda }^n\right),$$

$${\omega }_k^{n+1}=\frac{{\int }_0^{\infty }\omega |{\hat{u}}_k(\omega ){|}^{2}d\omega }{{\int }_0^{\infty }|{\hat{u}}_k(\omega ){|}^{2}d\omega },$$

$$\mathbf{x}(t)=\sum _{k=1}^K{\mathbf{u}}_k(t),$$

$$f={\sum}_{k}\Bigg{\|}\partial_{t}\Big{[}e^{-j\omega_{k}t}\mathbf{u}_{+}^{k}(t)\Big{]}\Bigg{\|}^{2}_{2}$$

$$f={\sum_{k}\sum_{c}}\Bigg{\|}{\partial_{t}\Big{[}u_{+}^{k,c}(t)e^{-j\omega_{k}t}\Big{]}\Bigg{\|}^{2}_{2}}$$

$$\underset{\{u_{k,c}\},\{{\omega }_k\}}{\mathrm{minimize}}$$

$$\sum _k{u}_{k,c}(t)=x_c(t),c=1,2,\dots ,C.$$

$$\displaystyle\mathcal{L}\left({\{u_{k,c}\},\{\omega_{k}\},\lambda_{c}}\right)=\alpha\ {\sum_{k}\sum_{c}}\Bigg{\|}\partial_{t}\Big{[}u_{+}^{k,c}(t)e^{-j\omega_{k}t}\Big{]}\Bigg{\|}^{2}_{2}$$

$$\displaystyle+{\sum_{c}}\Bigg{\|}x_{c}(t)-\sum_{k}u_{k,c}(t)\Bigg{\|}_{2}^{2}+{\sum_{c}}\Bigg{\langle}\lambda_{c}(t),x_{c}(t)-\sum_{k}u_{k,c}(t)\Bigg{\rangle}.$$

$$\{u_{k,c}^1\}\text{mbox},\{{\omega }_k^1\}\text{mbox},{\lambda }^1\text{mbox},n\leftarrow 0$$

$$u_{k,c}^{n+1}\leftarrow \underset{u_{k,c}}{\text{arg min }}\mathcal{L}\left(\{u_{i<k,c}^{n+1}\},\{u_{i\ge k,c}^n\},\{{\omega }_i^n\},{\lambda }_c^n\right)$$

$${\omega }_k^{n+1}\leftarrow \underset{{\omega }_k}{\text{arg min }}\mathcal{L}\left(\{u_{i,c}^{n+1}\},\{{\omega }_{i<k}^{n+1}\}\{{\omega }_{i\ge k}^n\},{\lambda }_c^n\right)$$

$$\lambda_{c}^{n+1}=\lambda_{c}^{n}+\tau\Big{(}x_{c}-\sum_{k}u_{k,c}^{n+1}\Big{)}$$

$$\displaystyle u_{k,c}^{n+1}=\underset{u_{k,c}}{\text{arg min}}\Bigg{\{}\alpha\Bigg{\|}\partial_{t}\Big{[}u_{+}^{k,c}(t)e^{-j\omega_{k}t}\Big{]}\Bigg{\|}^{2}_{2}+$$

Full text

I Introduction

Nonstationary signal processing methods have attracted a lot of interest over the last few decades owing to their relevance and applicability to a large class of real world signals. Here we are often interested in tracking time variations of packets of close frequencies in input data. One way of achieving that is through linear time-frequency transforms that correlate or project input data with a family of signals, termed as time frequency atoms, that are well concentrated in both time and frequency. Short Time Fourier Transform (STFT) [1] and wavelet transform [2] are two popular examples of linear transform methods. In addition to providing a graphical display in the form of time-frequency (T-F) spectrum (spectrogram or scalogram) for the purpose of exploratory data analysis, these methods enable signal separation at multiple resolution thus enabling signal processing at each scale separately and subsequent reconstruction.

A downside of linear transform methods is that the time frequency atoms are fixed and predefined. They exhibit limited joint time and frequency resolution which blurs the resulting T-F representation due to the famous Heisenberg uncertainty principle. That may be problematic in many real world applications that involve intermittent signal patterns i.e. nonstationarity. To that end, there has been a surge of interest in data-driven signal decomposition and T-F methods that make little or no a priori assumptions on input data.

The trend started in the late nineties when Huang et al. proposed a recursive algorithm, termed empirical mode decomposition (EMD) [3]. EMD decomposes input data into its inherent oscillatory modes, called as intrinsic mode functions (IMFs), through a recursive sifting process that makes use of signal extrema. Specifically, signal extrema are interpolated to yield upper and lower envelopes which are averaged to obtain a local mean of the signal. Local mean, which are low-frequency estimates of the data, are then recursively removed from input data to yield a high-frequency (or fast oscillating) mode in the signal. The process is repeated until all principal oscillatory modes present in the data are recovered. Owing to its fully data driven nature, EMD avoids several limitations of linear T-F methods, such as limited T-F resolution (Heisenberg uncertainty theorem) and physically meaningless modes obtained due to the use of fixed T-F atoms. Consequently, EMD has made a significant impact in scientific community and is routinely employed in wide-ranging engineering applications such as biomedical signal classification [4, 5], climate analysis [6] and off-line machine monitoring [7].

In essence, the goal of EMD is to extract multicomponent signals (IMFs) $c_{k}(t)$ and trend $r(t)$, from input data $x(t)$, as follows

[TABLE]

where $c_{k}(t)$ are defined such that they are amplitude- and frequency-modulated (AM-FM) signals

[TABLE]

Here, $\phi_{k}(t)$ denotes instantaneous phase which satisfies $\phi^{{}^{\prime}}(t)\geq 0\mbox{, }\forall t$ and the instantaneous amplitude $a_{k}(t)$ is nonnegative i.e., $a_{k}(t)\geq 0\mbox{, }\forall t$. The instantaneous frequency (IF), denoted by $\omega_{k}(t)$, is then defined by $\omega_{k}(t)=\frac{d\phi(t)}{dt}$ and is assumed to vary much slower than the phase $\phi_{k}(t)$.

Despite its advantages, EMD remains a recursive process whose output depends heavily on the choice of interpolation scheme and stopping criteria to stop the sifting process. Due to the lack of mathematical theory, there is a little or no margin to obtain necessary theoretical guarantees for EMD decomposition.

To overcome that difficulty, earlier attempts included replacing envelope and local mean estimation in EMD with robust constrained optimization based methods [8, 9]. Another class of data driven methods that aim to extract EMD-like decomposition include synchrosqueezed transform (SST) [10] and empirical wavelet transform (EWT) [11]. SST belongs to a class of methods that obtains AM-FM signal components as given in (2), albeit through a sound and convenient mathematical framework in contrast to EMD. It can operate both in the STFT and wavelet domain and works by first sharpening the signals’ STFT spectrogram (or wavelet scalogram) through a frequency reassignment operator, similar to the one proposed in [12]. Then, assuming that the total number of signal components are known a priori, ridge extraction techniques are employed to estimate ridges of IFs i.e., $\omega_{k}(t)$. Finally, signal components are recovered by integrating the reassigned STFT (or reassigned wavelet transform) in the vicinity of the corresponding ridges of $\omega_{k}(t)$. EWT, on the other hand, adopts an empirical approach that relies on robust peak detection mechanism and spectrum segmentation techniques based on the detected peaks to develop a wavelet filter bank for effective data decomposition.

Recently, a variational approach to data decomposition, termed variational mode decomposition (VMD), has been proposed that obtains principal modes of a signal, similar in form to (2), adaptively and concurrently by solving a convex optimization problem [13]. In this approach, an optimization problem is formulated in order minimize the collective bandwidth of modes, subject to the constraint that the modes fully reconstruct the input signal. The bandwidth is measured as $\mathcal{H}^{1}$ norm of corresponding analytic signal shifted to the baseband via complex harmonic mixing. The resulting variational model is minimized through alternate direction method of multipliers (ADMM) approach to find an ensemble of modes that are band-limited around a center frequency $\omega_{k}$; $\omega_{k}$ is determined online during the optimization process.

While efforts to develop fully data driven yet mathematically sound algorithms for signal decomposition and T-F analysis of nonstationary data continue to grow, there has been a lot of interest in extending existing data driven approaches to process nonstationary multidimensional and multivariate (multichannel) data sets. Owing to advances in sensor and computing technologies, such data types are routinely acquired and need to be processed in modern engineering and scientific applications e.g., classification based on multichannel electroencephalogram (EEG) [14] and ECG signals [15] and denoising [16]. The main challenges involved in multivariate T-F algorithms for nonstationary data include: i) alignment of frequency information across multiple channels in each mode, termed as mode-alignment; ii) incorporating any correlation between multiple data channels. Typically, developing extensions that operate directly in multidimensional spaces where input signal resides are able to fulfill the above requirements.

Multivariate signal processing using standard linear transform methods is straightforward. Fixed basis functions (filters) associated with these transforms can be applied to each channel of a multivariate signal separately. For fully data driven methods like EMD, however, specialized extensions are employed for multivariate data that operate directly in multidimensional space where signal resides. A generic multivariate extension of EMD, termed as multivariate EMD (MEMD) [17] belongs to this class of methods that can operate on multichannel data with any number of input channels. While univariate EMD works on the principle of separating faster oscillations from slower oscillations, MEMD aims to separate faster multivariate oscillations from slower ones. To accomplish that, local mean of a multivariate signal is directly estimated in multidimensional space using a uniform projection based approach. MEMD and other extensions of EMD, especially for bivariate [18] and trivariate [19] signals, have been utilized in wide ranging applications involving multivariate data e.g., in biomedical engineering applications involving EEG and ECG, acoustics, fault detection and machine monitoring based applications to name a few. However, these extensions inherit all the limitations of standard EMD such as sensitivity to sampling rate, lack of robustness to noise and issues related to the empirical and algorithmic nature of the EMD algorithm.

Multivariate extensions for wavelet based methods such as SST [20] and EWT [21] have also emerged recently. The multivariate extension of SST operates by first applying the standard SST algorithm to each channel separately. That is followed by an adaptive partitioning of the T-F domain in order to separate monocomponent multivariate oscillations in the input data. Finally, multivariate instantaneous frequency and amplitude are estimated and a multivariate synchrosqueezed transform operator is computed based on those estimates. That results in a sharp T-F spectrum of multivariate data which is useful for exploratory data analysis. However, individual modes cannot be reconstructed using this approach which limits its applications. The multivariate extension of EWT [21] employs a mode estimation procedure to obtain an optimum signal in the multivariate data set and then segments its corresponding spectra to recover relevant modes across all the channels of input data. The method suffers from inherent flaws of EWT in that it requires explicit construction of adaptive wavelet filters based on spectrum segmentation.

Recently, an extension of VMD to process complex-valued signals has been proposed [22]. The paper investigates filterbank structure of complex VMD for wGn and also proposes bi-directional T-F spectrum for complex-valued data based on complex VMD. The method, by design, is however not suitable for extracting matched frequency components across its channels [30, 17]. In addition, it is applicable only for bivariate data containing up to two channels only.

In this paper, we present a novel extension of the VMD algorithm, termed multivariate VMD (MVMD), to process multivariate data containing any number of channels. This is the first fully multivariate extension of VMD and is important due to the following limitations of existing multivariate data driven approaches: i) multichannel extension of EMD, namely MEMD, suffers from all the problems of standard EMD, including absence of mathematical framework and lack of robustness to low sampling rates and noise. mode-mixing among different channels of MEMD [23] is also problematic; ii) multivariate extension of SST [20] does not have mode separation capability. It only gives a graphical representation of T-F spectrum thereby limiting its scope of applications to exploratory data analysis only; iii) multivariate extension of EWT [21] is simplistic and relies on predefined boundaries of wavelet filterbank which may be difficult to define in real applications.

In the proposed variational model of MVMD, we seek a collection of common multivariate modulated oscillations residing in input data that exhibit minimum collective bandwidth while fully reconstructing all channels of input data. Multivariate oscillations are modeled using canonical analytic signal representation based on Hilbert transform with a constraint that a joint oscillation (or joint frequency component) exists among all data channels, i.e., multivariate instantaneous frequency. The proposed multivariate variational model is then constructed as a generic extension of standard univariate VMD model to multivariate data. That leads to our proposed method inheriting a lot of desirable properties of standard VMD. Like standard VMD, we solve the optimization problem directly in the Fourier domain through the application of ADMM optimization approach [24]. The effectiveness of the proposed method is demonstrated through extensive simulations performed on both test (synthetic) signals and real world data sets. We specifically focus on the mode-alignment property of the proposed method for multivariate data that is critical in many real world applications including fusion [25] and denoising [16].

The remaining part of the paper is organized as follows: Section II introduces the concept of univariate and multivariate modulated oscillations which are utilized in the construction of VMD and the proposed multivariate VMD models respectively. In Section III, we review the standard VMD algorithm. The proposed multivariate VMD method is illustrated in detail in Section IV. To demonstrate the effectiveness of our method, Section V presents detailed experimental results on both synthetic and real world signals and also includes related discussion. The section is divided into two parts: firstly, we show the ability of the proposed method to extract multivariate modulated oscillations (or common rotational modes) in multivariate data, its ability to align common modes across multiple channels and its noise robustness. Finally, we demonstrate the utility of the proposed method through its decomposition capability on real world data sets including multivariate EEG data and cardiotocographic traces.

II Univariate and Multivariate Modulated Oscillations

In this section, we describe mathematical representations of modulated oscillations in single and multidimensional spaces that will be utilized in the construction of standard single-channel VMD and the proposed multivariate VMD models, respectively. We only illustrate important and relevant concepts here; interested readers are referred to [26, 27, 28] for a detailed treatment on the topic.

II-A Univariate Oscillations

We start with a univariate AM-FM signal $u(t)$ which is given by

[TABLE]

where $a(t)$ is amplitude function while $\phi(t)$ represents time varying phase function due to frequency modulation. The above representation (3) is not unique since more than one functions for $a(t)$ and $\phi(t)$ may be associated with real valued $u(t)$.

By using an analytic representation of $u(t)$ by means of Hilbert transform, however, a unique canonic set of amplitude and phase functions can be defined. That analytic representation $u_{+}(t)$ is defined by pairing the signal with its own Hilbert transform to yield the following complex signal

[TABLE]

where the Hilbert transform of $u(t)$ is given by

[TABLE]

where $\int$ is the Cauchy principal integral value. The Fourier transform of $u_{+}(t)$ yields a unilateral signal spectrum and is defined as

[TABLE]

The original real-valued signal $u(t)$ can finally be recovered from the corresponding complex-valued analytic signal $u_{+}(t)$ by taking its real part

[TABLE]

II-B Multivariate Oscillations

To obtain a convenient mathematical form for multivariate oscillations, we first represent a set of $C$ real valued AM-FM signals, denoted by ${\{}u_{i}(t){\}}_{i=1}^{C}$, in a vector form

[TABLE]

where $a_{i}(t)$ and $\phi_{i}(t)$ denote amplitude and phase function corresponding to the $i$th signal component respectively.

Next, we obtain an analytic representation of the vector signal $\mathbf{u}(t)$ by employing the Hilbert transform operator, as defined in (5), to each component of $\mathbf{u}(t)$ separately. That results in the following analytic signal vector

[TABLE]

where the canonical (and unique) set of amplitude and phase functions for the $i$th component are given by $a_{i}(t)=|u_{+}^{i}(t)|$ and $\phi_{i}(t)=\text{arg}(u_{+}^{i}(t))$. Here, “arg” denotes the complex argument or the phase function. Finally, the original real-valued signal vector can be obtained as the real part $\mathcal{R}$ of $\mathbf{u}_{+}(t)$ i.e.,

[TABLE]

Note that the above mathematical description of $\mathbf{u}(t)$ considers the $C$ number of AM-FM components in isolation from each other. For a multivariate modulated oscillation, there may be one or more common frequency components present in $\mathbf{u}(t)$. In this work, we consider a simplified multivariate analytic signal model [28] for $\mathbf{u}_{+}(t)$ that assumes a single common component, i.e., $\omega=\frac{d\phi(t)}{dt}$, among all data channels

[TABLE]

The above model will be used in the formulation of the proposed variational decomposition model for multivariate data in Section IV.

III Variational Mode Decomposition

In this section, we describe the variational mode decomposition (VMD) algorithm for univariate data [13]. The goal of VMD is to decompose an input data into a finite and predefined $K$ number of intrinsic principal modes $u_{k}(t)$, that are similar in form to (2), i.e.,

[TABLE]

The modes are chosen based on a variational model that aims to minimize the sum of bandwidths of all modes ${\{}u_{k}(t){\}}_{k=1}^{K}$ while also reconstructing the input signal fully or in least squares sense through the addition of modes. In VMD, the bandwidth of each mode $u_{k}(t)$ is estimated by: i) computing the analytic signal representation $u_{+}^{k}(t)$ that exhibits unilateral frequency spectrum; ii) shifting the resulting unilateral spectrum to baseband via harmonic mixing with a complex exponential of frequency $\omega_{k}$; iii) taking the squared $L^{2}$ norm of the gradient of the harmonically mixed signal obtained in the previous step. The associated constrained variational optimization problem therefore becomes

[TABLE]

where $u_{+}^{k}(t)$ denotes analytic signal corresponding to $u_{k}(t)$; the symbol $\partial_{t}$ represents partial derivative operation with respect to time; $\{u_{k}\}$ and $\{\omega_{k}\}$ respectively denote sets of all $K$ number of modes and their center frequencies.

Next, an augmented Lagrangian is constructed that renders the problem given in (14) unconstrained. Two penalty terms are added in the Lagrangian: a quadratic term to enforce reconstruction fidelity and a term with Lagrangian multipliers $\lambda$ to ensure that the constraints are fulfilled strictly. The resulting augmented Lagrangian function is given below

[TABLE]

The solution to the original optimization problem as given in (14) can now be found as the saddle point of the Lagrangian function (15). In VMD, (15) is solved using the alternate direction method of multipliers (ADMM) approach. Using ADMM, the complete optimization problem is solved as a sequence of iterative sub-optimization problems. These sub-problems are convenient to handle since they seek to minimize cost function iteratively for a single parameter/function of interest rather than optimizing cost function for all optimization variables simultaneously. For instance, to update the mode $u_{k}(t)$, the following sub-optimization problem is considered at the $n$th iteration

[TABLE]

The above minimization problem is solved in the spectral domain within VMD \big{(}see (19) - (22) in [13]\big{)}, resulting in the following update in the frequency domain

[TABLE]

In the next stage, in order to get an update on center frequency $\omega_{k}$, the following sub-optimization problem is solved iteratively

[TABLE]

to get an update for $\omega_{k}^{n+1}$ in the dual frequency domain

[TABLE]

which estimates the new frequency as the center of gravity of the associated modes’ power spectrum.

The mode and the center frequency update relations in the spectral domain, given in (17) and (19) respectively, form the crux of the VMD algorithm. We refer readers to Algorithm 2 in [13] for details.

IV Multivariate Variational Mode Decomposition

A novel multivariate extension of VMD, namely multivariate VMD (MVMD) is presented in this section. As a generalized extension of the original VMD algorithm for multivariate data residing in multidimensional spaces, the main goal of MVMD is to extract predefined $K$ number of multivariate modulated oscillations $\mathbf{u}_{k}(t)$ from input data $\mathbf{x}(t)$ containing $C$ number of data channels, i.e., $\mathbf{x}(t)=[x_{1}(t),x_{2}(t),\cdots x_{C}(t)]$

[TABLE]

where $\mathbf{u}_{k}(t)=[u_{1}(t),u_{2}(t),\cdots u_{C}(t)]$.

The goal is to extract an ensemble of multivariate modulated oscillations $\{\mathbf{u}_{k}(t)\}_{k=1}^{K}$ in input data such that: i) the sum of bandwidths of the extracted modes is minimum and ii) the sum of the extracted modes exactly recover the original signal $\mathbf{u}_{k}(t)$. To accomplish that, let us denote the vector analytic representation of $\mathbf{u}_{k}(t)$ by $\mathbf{u}_{+}^{k}(t)$ (see (9)). The bandwidth of $\mathbf{u}_{k}(t)$ can then be estimated by taking the $L_{2}$ norm of the gradient function of the harmonically shifted $\mathbf{u}_{+}^{k}(t)$. The resulting cost function $f$ for MVMD then becomes a multivariate extension of the cost function used in the corresponding VMD optimization problem in (14) and is given by

[TABLE]

Another important point to highlight in (21) is that a single frequency component $\omega_{k}$ is used in the harmonic mixing of the whole vector $\mathbf{u}_{+}^{k}(t)$. By design, therefore, we are looking to find multivariate oscillations in $\mathbf{u}_{k}(t)$ that have a single common frequency component $\omega_{k}$ in all channels. That is exactly how we defined our model for multivariate modulated oscillation in (12). The bandwidth of modulated multivariate oscillation is therefore estimated by shifting the unilateral frequency spectrum of all channels of $\mathbf{u}_{+}^{k}(t)$ by $\omega_{k}$ and taking the Frobenius norm 111The Frobenius norm is taken as a direct extension of the $L_{2}$ norm used in original VMD to the matrices which appear due to multiple channel signal representation in MVMD. of the resulting matrix. That leads to the following convenient representation of $f$

[TABLE]

where $u_{+}^{k,c}(t)$ denotes the analytic modulated signal corresponding to channel number $c$ and mode number $k$. In contrast to the vector signal in (21), $u_{+}^{k,c}(t)$ is a complex-valued signal with single component in (22). We now present the constrained optimization problem for MVMD

[TABLE]

Note that there are multiple linear equality constraints in the above model corresponding to the total number of channels. The corresponding augmented Lagrangian function then becomes

[TABLE]

The above unconstrained optimization problem is solved using ADMM approach as illustrated in Algorithm 1. Note from Algorithm 1 that how the ADMM approach converts a complex optimization problem (23) into multiple more simpler suboptimization problems, as given by the update equations (25)-(27).

IV-A Mode update

We first focus on the minimization problem related to the mode update (25). Its equivalent optimization problem is given below

[TABLE]

This subproblem is similar in form to the mode update problem of original VMD (16). There, the problem was solved in the Fourier domain to yield a convenient update relation in the frequency domain (17). We use that result to give the following mode update relation in our case

[TABLE]

IV-B Center frequency update

Now we turn our attention to the optimization problem corresponding to center frequency update (26). Considering that the last two terms of the Lagrangian (24) do not depend on $\omega_{k}$, the relevant problem simplifies to

[TABLE]

The optimization can be performed in the frequency domain more conveniently. Using the Plancherel theorem that relates the inner product of a function in time and frequency domain, we obtain an equivalent problem of (30) in the Fourier domain which is given below

[TABLE]

Minimizing the above sum of quadratic functions by setting its first derivative to zero, followed by a few simple algebraic manipulations yields the following relation

[TABLE]

By comparing the above frequency update relation of MVMD with the corresponding relation for VMD (19), we note that while updating $\omega_{k}$ for each mode in the proposed MVMD method, contributions from the power spectrum of all $C$ number of channels is taken into account.

Having the mode (29) and the center frequency update relations (32) in the frequency domain at our disposal, we can perform the optimization in frequency domain using the ADMM approach; the steps are listed in Algorithm 2.

V Experimental Results

To evaluate the performance of the proposed method, we conducted experiments and simulations on a wide range of multivariate signals; their detailed results are presented in this section.

Being a direct extension of the original VMD algorithm for multivariate data, MVMD inherits all the important advantages of the original VMD such as its robustness to noise and low sampling rates and the ability to separate tones very effectively. All these properties of VMD have been explored in detail in [13] and we will not verify them for MVMD, in this section. Instead, we will focus on the ability of the proposed method to align common frequency scales across multiple data channels. That is a critical requirement in many scientific and engineering applications involving multivariate data such as image fusion [29], biomedical signal based classification [14] and denoising [16], to name a few.

We first show the ability of the proposed method to decompose a real world bivariate data into its inherent bivariate modulated oscillations, which are also referred to as rotational modes. Secondly, we demonstrate the mode-alignment property of MVMD on synthetic test signals containing a combination of tones across its different channels and white Gaussian noise (wGn). We also compare the results with those obtained from applying original VMD to each channel separately. We next illustrate the effect of noise in input data on mode-alignment property of MVMD. We show comparison of our method against multivariate EMD (MEMD) algorithm. Finally, the effectiveness of the proposed method is shown on a couple of real world data sets including multivariate EEG and bivariate cardiotocographic (CTG) data.

V-A Separation of multivariate modulated oscillations

We demonstrate the ability of the proposed method to separate multivariate modulated oscillations or rotational modes from a real bivariate data set. The data was taken during the Eastern Basin experiment and consists of position record of a subsurface oceanographic float that was deployed in North Atlantic ocean to track the trajectory of salty water flowing from the Mediterranean Sea. The data can be downloaded from the World Ocean Circulation Experiment Subsurface Float Data Assembly Center (WFDAC) at http://wfdac.whoi.edu.

Separate input channels of the data are shown in Figure 1(a) (top row) along with the illustration of the data in 2D space (see Figure 1(b) (top left)) to visualize multivariate oscillations. From the graphical 2D representation of the input data, bivariate modulated oscillations (or 2D rotations) are quite evident. Looking at the corresponding time plots, joint or common frequency scales across data channels can be seen that explains the presence of rotations or multivariate modulated oscillations in the 2D plot of the data (See Section II.B for detail).

We applied the proposed MVMD to the bivariate data with an aim to separate its principal multivariate modulated oscillations. We decomposed the data into $K=4$ modes which are shown in Figure 1, both as time plots (Figure 1(a)) and 2D graphical plots (Figure 1(b)). Note from the 2D plot that the first two modes consist of multivariate modulated oscillations or rotating modes which may correspond to the presumed coherent vortex. The non-rotating modes are characterized by the absence of joint frequency scale in one of the channels. This example demonstrates the ability of MVMD to effectively separate multivariate oscillations from input data while also verifying the mode-alignment property of the algorithm on a practical data set.

V-B Mode-alignment property

Mode-alignment of multivariate data refers to the alignment of common or joint oscillations (containing similar frequency content) across multiple channels of a single mode. This property is of fundamental importance to ensure coherent multivariate T-F analysis in many practical applications involving multivariate data such as fusion, denoising and classification. Typically, in such applications, signal processing methods are applied to each mode separately followed by their reconstruction. In case a mode contains oscillations that exhibit different frequency content across multiple channels, application of signal processing methods will result in physically meaningless estimates or decisions due to misalignment of similar information content across channels. Readers are referred to [30] for detailed discussion on the importance of mode-alignment in fusion data related application.

We illustrate the ability of MVMD to identify and align principal modulated oscillations present in the data across multiple channels and modes in Figure 2(a). The input data was a bivariate signal whose individual components were a mixture of a 36-Hz sinusoid that was common to both the data channels; a 2-Hz tone in the channel-1 and a 24-Hz tone in the channel-2. The $K=3$ number of modes obtained by applying MVMD to the above dataset are shown. Observe that all modes are aligned in terms of their frequency content: the 36-Hz tone present in all data channels is localized in a single mode $u_{3}$. The 2-Hz signal is located in the channel-1 of $u_{1}$ while the 24-Hz tone is localized in the channel-2 of $u_{2}$. In Figure 2(b), we illustrate that similar mode-alignment across channels is not possible by applying VMD to each channel separately. Notice that frequency content across channels in mode $u_{2}$ is not aligned i.e., 36-Hz in channel-1 and 24-Hz in channel-2.

V-C MVMD vs MEMD filterbanks for wGn

In [31, 32], EMD and MEMD had respectively been shown to exhibit quasi-dyadic filterbank structure for wGn, which is similar to the wavelet filterbank. That has led to diverse range of applications of (M)EMD in science and engineering. We show here that while MVMD also exhibits filterbank structure for wGn, it appears to be different from the quasi-dyadic structure observed within EMD based algorithms.

To verify that, we applied MVMD on 50 realizations of a four-channel wGn process of length $L=1000$. We also applied MEMD to the same data set. The power spectral density (PSD) plots averaged over 100 realizations are plotted for both MVMD and MEMD in Figure 3. By comparing the two set of plots, we can observe that MEMD exhibits quasi-dyadic filterbank structure for wGn as evident by almost similar bandwidths of different IMFs in the log-frequency domain. On the other hand, MVMD modes exhibit different bandwidths in the log-frequency domain. A thorough investigation of the MVMD filterbanks is beyond the scope of this paper and would be a good topic for future research.

Also observe from the Figure 3 that there is alignment of frequency content among different channels within similar modes of both MVMD and MEMD: alignment within MVMD appears slightly more prominent as compared to MEMD, resulting in well defined subband filters when compared against MEMD.

V-D Quasi-orthogonality of MVMD modes

In linear signal decomposition and T-F methods such as STFT and wavelet transform, predefined basis functions are by construction chosen to be orthogonal. That ensures that there is no ‘leakage’ of information across different modes (and channels). Since basis functions within data driven decomposition and T-F methods are adaptive in nature, it is important to demonstrate quasi-orthogonality empirically.

In this section, we establish the quasi-orthogonality of MVMD modes for multiple realizations of the four-channel wGn process, which were used in the previous section to demonstrate the filterbank structure for MVMD and MEMD. Here, we use correlation coefficient as a means to quantify the dependence between different modes that were generated from MVMD and MEMD. Specifically, to compute the correlation coefficient between the modes $u_{i}$ and $u_{j}$, whose variances are denoted by $\sigma_{i}$ and $\sigma_{j}$ respectively, we have

[TABLE]

where $\text{cov}(.)$ denotes the covariance between a set of signals.

If the correlation coefficient is close to zero for a set of modes, those modes can be considered to be quasi-orthogonal; a value close to unity suggests very strong correlation. We obtained the correlation coefficient matrix for the set of $K=8$ modes obtained from MVMD and MEMD for 4-channel wGn, averaged them over 100 realizations, and plotted the resulting matrix as gray-scale images in Figure 4. Notice the almost diagonal structure of the correlation matrix in the case of MVMD, suggesting a strong quasi-orthogonality among MVMD modes. The correlation matrix corresponding to MEMD, on the other hand, shows some ‘leakage’ between adjacent modes as apparent from its tridiagonal nature.

V-E Noise robustness

One of the highlights of the VMD related algorithms is its robustness to noise. The origin of this robustness is the inherent flexibility of the method to forgo strict reconstruction property in the presence of noise. We demonstrate here that MVMD also inherits that property from standard VMD as a result of being a natural extension of the method. In the case of MVMD, we achieve that by making the Lagrangian term in (24) to vanish by forcing the parameter $\tau$ in (35) to be equal to zero. In this section, our interest specifically is in showing that the mode-alignment property of MVMD is preserved by the presence of noise in input channels. We also compare the performance of our method with MEMD which is known to be very sensitive to noise in input data channels.

For our experiments, we obtained a bivariate test signal that consisted of a mixture of three tones in its two channels. The 2-Hz and 288-Hz tone were present in both the channels whereas 24-Hz tone was only present in the channel-2. We added an unbalanced wGn with zero mean and variances of $\sigma^{2}=0.1$ and $\sigma^{2}=0.05$ that resulted in the SNR of $-dB$ and $-dB$ in channel-1 and channel-2 respectively. The resulting decomposed modes obtained from MVMD are shown in Figure 5(a). Notice that the MVMD recovers the common modes of 2-Hz and 288-Hz tones in both channels in mode 1 and mode 3 respectively. The 24-Hz tone was only present in channel-2 of mode 2. We also show the decomposed intrinsic mode functions (IMFs) obtained by applying the MEMD algorithm to the same data set. A visual comparison between the two sets of modes reveal that MVMD is more robust to noise in that the mode-alignment is more prominent across MVMD modes as compared to those obtained from MEMD. Some leakage of information can be seen in IMF 4 across channels and there is definite mode-mixing within the first three IMFs: $c_{3}$ does not fully represent the 288-Hz tone since the first two IMFs also share that information i.e., mode-mixing which may be problematic in practical applications. Similar results were reported in [13] while comparing single channel VMD vs EMD in terms of robustness to noise.

By design, the VMD and MVMD algorithms are inherently robust to noise due to control over Lagrangian multipliers via the parameter $\tau$. In the case of noisy input data, it may not be desirable to obtain a decomposition that fully reconstructs the noisy input data. That can be controlled within VMD and MVMD by forcing the parameter $\tau$ to be equal to zero. In MEMD, on the other hand, there is no mechanism to stop the noise entering into the decomposition process. As a result, in many practical application involving noisy data, a common preprocessing step in EMD involves rejecting noisy IMFs manually which may be prone to errors.

V-F Performance evaluation for increasing number of channels

The performance of MVMD is evaluated for varying number of data channels. In our experiments, the input were 200 realizations of wGn process with $C=2$, $4$ ,$8$ and $16$ noise channels. After obtaining the MVMD modes in each case, we investigated how the mode-alignment property, filterbank structure and quasi-orthogonality of the modes were affected with increasing number of input channels. We show the resulting filterbank plots and the correlation coefficient matrices in Figure 6. It can be observed from the top row of the Figure 6 that the filterbank and mode-alignment properties within MVMD are perfectly retained as the number of channels are increased. Similarly, quasi-orthogonality of MVMD decomposition is also not affected by increasing number of channels as highlighted by the diagonal nature of all correlation plots shown in the second row of Figure 6.

V-G Real World Examples

The utility of the proposed MVMD method in signal decomposition and T-F analysis of real world data is illustrated via two examples. Those include separation of $\alpha$-rhythms in multivariate EEG data and decomposition of fetal heart rate (FHR) and maternal uterine contraction (UC) recordings from cardiotocographic (CTG) traces [34]. Both these examples involve multiscale signal decomposition of multivariate data and require alignment of similar frequency content across multiple channels of each decomposed mode. MVMD, by its design, achieves such a decomposition and is therefore a suitable candidate for such applications.

V-G1 Separation of $\alpha$ rhythms in EEG

An example of multichannel Electroencephalogram (EEG) signal processing for brain computer interface (BCI) using MVMD is described. A 4-channel EEG data was obtained in an experiment that involved a subject who remained in the relaxed state with eyes closed for a certain period of time. It is well established that $\alpha$-rhythms, corresponding to frequency range of 8-12 Hz, are detected in EEG data during the relaxed and eyes-closed state. The aim here is to investigate the frequency localization property of MVMD to detect $\alpha$-rhythms at multiple channels. The EEG data were obtained through the OpenBCI Cyton board at a sampling frequency of 250 samples per second.

Figure 7 shows the time plots of the original data and selected modes ($u_{2}$, $u_{3}$, $u_{5}$) obtained from MVMD (a) along with their estimated power spectral density (PSD) plots (b). Note that the $\alpha$-rhythm within EEG is localized in the second mode $u_{2}$; it was also observed that all channels corresponding to $u_{2}$ contained $\alpha$-rhythm thus highlighting the mode-alignment property of MVMD across multiple channels of the data. Similar mode-alignment was observed across multiple channels of all the extracted modes. $u_{5}$ corresponds to artifacts due to the AC mains power line. The smoothed PSD plots of the three modes in Figure 7(b), with each showing distinct peak along a certain frequency band, and the alignment of PSD estimates from different channels highlight the mode-mixing and mode-alignment property of MVMD respectively.

V-G2 Decomposition of fetal heart rate (FHR) and maternal uterine contraction (UC) recordings from cardiotocographic (CTG) signal

Fetal heart rate (FHR) signal monitoring provides a useful means to assess the fetal health in obstetric practice. Precise monitoring of FHR can be achieved through Cardiotocography (CTG) that also captures maternal uterine contractions (UCs), making CTG an attractive technique in obstetrics [33]. A recent study utilized the mode-alignment property of bivariate extension of EMD (BEMD) to specify robust features and measures for classification of vaginal vs cesarean delivery [34]. The data was taken from a freely available CTU-UHB intrapartum cardiotocography database available at Physionet: http://www.physionet.org/physiobank/database/ctu-uhb-ctgdb/ [35]. We highlight here the potential of MVMD to avoid mode-mixing and achieve mode-alignment in a single record of FHR-UC data and compare the results from those obtained from BEMD.

Figure 8 shows the plots of selected decomposed modes obtained by applying the BEMD and MVMD to one of the patients’ CTG record i.e., FHR-UC bivariate data; the total number of modes (or IMFs) obtained in both cases were $M=7$. Looking at the BEMD decomposition of the data, multiple instances of mode-mixing can be observed, for instance, at the time=20 sec in IMF1, IMF3 and IMF4; and at the time=100 sec in IMF3, IMF5 and IMF6 and at the time instant of around 65 seconds in IMF2 and IMF3. MVMD modes, on the other hand, are free of any artifacts arising due to mode-mixing in that each mode consists of narrow band frequencies. MVMD modes can also be seen to exhibit mode-alignment across all channels of the selected modes in the Figure 8. Given that mode-mixing and misalignment of frequency information across channels can severely limit the performance of data driven multiscale algorithms, MVMD is expected to perform better than multivariate extensions of EMD across a range of practical applications.

VI Conclusion

We have proposed a novel extension of the variational mode decomposition (VMD) algorithm to multivariate data. VMD, in its original construction, aims to decompose a multicomponent single-channel signal into univariate modulated oscillations exhibiting limited bandwidth across a center frequency. That is achieved by formulating and solving a convex optimization problem that minimizes sum of bandwidths of all modes; bandwidth of a mode is estimated in VMD through squared $H_{1}$ norm of corresponding analytic signal, shifted to baseband by harmonic mixing with a complex exponential of center frequency estimate. We have proposed a generic extension of the VMD’s variational model to be applied for multivariate data. The crux of the proposed method is to decompose multivariate data into its inherent multivariate modulated oscillations. We have first defined a model for multivariate modulated oscillations that is based on constructing a Hilbert transformed multivariate analytic signal that has a common or joint frequency component across all its channels. Based on that model we formulate a convex optimization problem that aims to minimize the sum of bandwidths of all modes across all channels while fulfilling multiple constraints of exact signal reconstruction across all data channels. The resulting estimates of multivariate signal modes along with their center frequencies have been shown to fulfill the narrow band requirements typically associated with EMD-inspired methods like VMD and synchrosqueezed transform.

All multiscale data driven approaches for multivariate data require alignment of similar frequency content in a single scale or mode across all data channels, the so called mode-alignment property. The proposed method has been shown to exhibit this property on a variety of input multivariate signals ranging from test synthetic signal consisting of mixture of tones to multivariate wGn data and finally on data obtained from real world applications such as electroencephalogram (EEG) and cardiotocographic traces. As a result, we have demonstrated initial potential of the proposed method in data fusion application.

Being a generic extension of the VMD method, our proposed approach inherits all the advantages and flaws of the original VMD. Some of the challenges in the original VMD approach include: i) need to give predefined number of modes $K$ a priori; ii) inability to separate dc component in original data; and iii) problem to deal with signals exhibiting high level of nonstationarity such as sudden signal onset. We observed similar limitations in our proposed method. However, given the tremendous interest shown in VMD over the last few years since its inception, novel solutions to these shortcomings are expected which will improve both VMD and our proposed multivariate extension.