## 2 Variational Mode Decomposition  
This section summarizes the principle of VMD based on [21]. VMD is an entirely non-recursive decomposition method, and it extracts the constituent AM-FM components of a complex multi-component signal adaptively and concurrently. It defines IMFs as explicit AM-FM models, and relates the parameters of AM-FM models to the bandwidth of IMFs. According to the narrow-band property of IMFs, the AM-FM parameters can be found by minimizing the bandwidth, thus obtaining IMFs. This method has good merits over other available mode decomposition methods, such as theoretical rationale and robustness to noise and sampling. In essence, IMFs are AM-FM signals  
\[u_{k}(t)=A_{k}(t)\cos[\phi_{k}(t)], \tag{1}\]  
where \(A_{k}(t)\) is the amplitude envelope, \(\phi_{k}(t)\) is the instantaneous phase, and both the amplitude envelope \(A_{k}(t)\) and the instantaneous frequency \(\omega_{k}(t)=\frac{d\phi_{k}(t)}{dt}\) vary much slower than the instantaneous phase \(\phi_{k}(t)\). Such IMFs of AM-FM nature have a limited bandwidth. VMD decomposes a signal \(x(t)\) into an ensemble of IMFs \(u_{k}(t)\) that are band-limited about their respective center frequency \(\omega_{k}\), while reconstructing the signal optimally. It iteratively updates each IMF \(u_{k}(t)\) in the frequency domain, and then estimates the center frequency \(\omega_{k}\) as the center of gravity of the IMF power spectrum. Motivated by the narrow-band properties of the AM-FM IMF definition, each IMF is assumed to be mostly compact around a center frequency \(\omega_{k}\), i.e. it has specific sparsity properties. The sparsity prior of each IMF is described by its bandwidth. For each IMF \(u_{k}(t)\), in order to assess its bandwidth, its analytic signal is firstly computed by means of Hilbert transform to obtain a spectrum of unilateral non-negative frequency. Then, its spectrum is shifted to baseband by multiplying with an exponential harmonic tuned to the respective center frequency. Finally, the bandwidth can be estimated through the squared \(l_{2}\) norm of the gradient. The resulting constrained variational optimization problem is  
\[\min_{\{u_{k}(t)\},\{\omega_{k}\}}\sum_{k}\left\|\frac{\partial}{\partial t} \{[\beta(t)+j\frac{1}{\pi t}]*u_{k}(t)\}\exp(-j\omega_{k}t)\right\|_{2}^{2}, \text{ s.t.}\sum_{k=1}^{K}u_{k}(t)=x(t), \tag{2}\]where \(\delta(\cdot)\) is the Dirac delta function, the symbol \(*\) denotes the convolution operator, and \(K\) is the number of IMFs to be extracted. To render the optimization problem, Equation (2), into an unconstrained form, a quadratic penalty term and a Lagrangian multiplier are introduced, for quick convergence and strict enforcement of the constraint. Then, the objective function to be minimized becomes an augmented Lagrangian  
\[ \begin{array}{l}L[\{u_{k}(t)\},\{\omega_{k}\},\lambda(t)]=\alpha\underset{k} {\sum}\|\frac{\partial}{\partial t}[\delta(t)+j\frac{1}{\pi t}]*u_{k}(t)\}\exp( -j\omega_{k}t)\|_{2}^{2}\\ \qquad\qquad\qquad\qquad\qquad+\|x(t)-\underset{k}{\sum}u_{k}(t)\|_{2}^{2}+ \left\langle\lambda(t),x(t)-\underset{k}{\sum}u_{k}(t)\right\rangle,\end{array}  \tag{3}\]  
where \(\lambda(t)\) is the Lagrange multiplier, \(\alpha\) is the balancing parameter of the data-fidelity constraint, and \(\left\langle\cdot,\cdot\right\rangle\) stands for inner product. The solution to the minimization problem, Equation (2), can now be found as the saddle point of the augmented Lagrangian, Equation (3), in a sequence of iterative sub-optimizations [21]. Each IMF \(u_{k}(t)\) can be updated as a solution to a minimization problem equivalent to Equation (3)  
\[ \begin{array}{l}u_{k}(t)\ \ =\underset{u_{k}}{\operatorname{argmin}}L(\{u_{k}(t)\},\{ \omega_{k}\},\lambda(t)]\\ \ \ =\underset{u_{k}}{\operatorname{argmin}}\bigg{(}\alpha\underset{k}{ \sum}\|\frac{\partial}{\partial t}[\delta(t)+j\frac{1}{\pi t}]*u_{k}(t)\}\exp (-j\omega_{k}t)\|_{2}^{2}+\ \|x(t)-\underset{k}{\sum}u_{k}(t)+\frac{\lambda(t)}{2}\|_{2}^{2} \bigg{)}.\end{array}