IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO.

Y, ZZZ 20WW 1

Variational Mode Decomposition

Konstantin Dragomiretskiy and Dominique Zosso∗ , Member, IEEE

Abstract—In the late nineties, Huang introduced the Hilbert- signal into principal modes, however the resulting decomposi-
Huang transform, also known as Empirical Mode Decomposition. tion is highly dependent on methods of extremal point finding,
The goal is to recursively decompose a signal into different modes interpolation of extremal points into carrier envelopes, and the
of separate spectral bands, which are unknown beforehand. The
HHT/EMD algorithm is widely used today, although there is no stopping criteria imposed. The lack of mathematical theory and
exact mathematical model corresponding to this algorithm, and, the aformentioned degrees of freedom reducing the algorithm’s
consequently, the exact properties and limits are widely unknown. robustness all leave room for theoretical development and
A few limitations are quite apparent, though: the algorithm is improvement on the robustness of the decomposition [2], [3].
sensitive to noise and sampling. Therefore, EMD for example has In some experiments it has been shown that EMD shares
difficulties separating tones of similar frequencies. Several more
mathematical attempts to this decomposition problem have been important similarities with wavelets and (adaptive) filter banks
made, like synchrosqueezing, empirical wavelets or recursive [4].
variational decomposition into smooth signals and residuals. Despite the limited mathematical understanding and some
Here, we propose an entirely non-recursive variational mode obvious shortcomings, the EMD method, also known as the
decomposition model, where the modes are extracted concur-
Hilbert-Huang transform (HHT), has had significant impact
rently. The model looks for a number of modes and their
respective center frequencies, such that the modes reproduce and is widely used in a broad variety of time-frequency anal-
the input signal, while being smooth after demodulation into ysis applications. Applications involve signal decomposition in
baseband. In Fourier domain, this corresponds to a narrow-band audio engineering [5], climate analysis [6], and various flux,
prior. We show important relations to Wiener filter denoising. respiratory, and neuromuscular signals found in medicine and
Indeed, the proposed method is a generalization of the classic
biology [7], [8], [9], [10], to name just a few examples.
Wiener filter into adaptive, multiple bands. Our model provides
a solution to the decomposition problem that is theoretically well With EMD, and in all of the previous signals, the core
founded and still easy to understand. The variational model is assumption on the individual modes is that they have compact
efficiently optimized using an alternating direction method of Fourier support. In the original description, in such a mode
multipliers approach. Preliminary results show excellent perfor- the number of local extrema and zero-crossings differ at most
mance with respect to existing mode decomposition models. In
particular, single harmonics can be reconstructed independently by one [1]. In most related works, the definition is slightly
of their frequency and with precision controlled by a simple changed into so-called Intrinsic Mode Functions (IMF).
convergence tolerance criterion. Further, in contrast to EMD,
the proposed VMD model is able to precisely separate any pair Definition Intrinsic Mode Functions are amplitude-
of harmonics, largely irrespective of their relative amplitudes modulated-frequency-modulated (AM-FM) signals, written
and how close their frequencies are. Finally, we show promising as:
practical decomposition results on a series of artificial and real
data. uk (t) = Ak (t) cos(φk (t)), (1)
Index Terms—Mode decomposition, variational problem,
Wiener filter, AM-FM, spectral decomposition, Hilbert trans- where the phase φk (t) is a non-decreasing function, φ0k (t) ≥ 0,
form, Fourier transform, augmented Lagrangian. the envelope is non-negative Ak (t) ≥ 0, and, very importantly,
both the envelope Ak (t) and the instantaneous frequency
ωk (t) := φ0k (t) vary much slower than the phase φk (t) [11],
I. I NTRODUCTION [12].
Empirical Mode Decompositon (EMD) proposed by Huang
In other words, on a sufficiently long interval [t − δ, t + δ],
et al. [1] is an algorithmic method to detect and decompose a
δ ≈ 2π/φ0k (t), the mode uk (t) can be considered to be a
signal into principal “modes” - a signal with mostly compact
pure harmonic signal with amplitude Ak (t) and instantaneous
supported Fourier spectrum. This algorithm recursively detects
frequency φ0k (t) [11]. Note that the newer definition of signal
local minima/maxima in a signal, estimates lower/upper en-
components is slightly more restrictive than the original one.
velopes by spline-interpolation of these extrema, removes the
The immediate consequence of the IMF assumption is limited
average of the envelopes as “low-pass” centerline, thus isolat-
bandwidth.
ing the high-frequency oscillations as “mode” of a signal, and
Indeed, if ωk is the mean frequency of a mode, then its
continues recursively on the extracted “low-pass” centerline.
practical bandwidth increases both, with the maximum devia-
In some cases, this sifting algorithm does indeed decompose a
tion of the instantaneous frequency, ∆f ∼ max(|ωk (t) − ωk |),
∗ To whom correspondence should be addressed. and with the rate of change of the instantaneous frequency,
This work is supported by the Swiss National Science Foundation (SNF) under fFM ∼ ω 0 (t), according to Carson’s rule: BW = 2(∆f + fFM )
grant PBELP2 137727. [13]. In addition to this comes the bandwidth of the envelope
The authors are with the Department of Mathematics, University of California,
Los Angeles (UCLA), Box 951555, Los Angeles, CA 90095-1555, USA, Ak (t) modulating the amplitude of the FM signal, given by its
{konstantin,zosso}@[Link]. highest frequency fAM . Hence we estimate the total bandwidth
2 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO. Y, ZZZ 20WW

103 dubbed synchrosqueezing, was proposed by Daubechies et al.

1 [11], [16]. They remove unimportant wavelet coefficients (both
102
0
101 in time and scale) by thresholding of the respective signal en-
−1 ergy in that portion. Conversely, locally relevant wavelets are
100 1 selected as local maxima of the continuous wavelet transform,
0 0.2 0.4 0.6 0.8 1 10 102 103
that are shown to be tuned with the local signals, and from
a) AM signal and spectrum. (∆f = 0)
which the current instantaneous frequency of each mode can
1 103 be recovered.
0.5 102 Other recent work pursuing the same goal is the Empirical
0 Wavelet Transform (EWT) to explicitly build an adaptive
−0.5 101
wavelet basis to decompose a given signal into adaptive
−1 100 1
0 0.2 0.4 0.6 0.8 1 10 102 103 subbands [12]. This model relies on robust preprocessing for
b) FM signal and spectrum. (fFM  ∆f ) peak detection, then performs spectrum segmentation based on
detected maxima, and constructs a corresponding wavelet filter
1 103 bank. The filter bank includes flexibility for some mollification
0.5 102 (spectral overlap), but explicit construction of frequency bands
0
101 still appears slightly strict.
−0.5
−1 100 1 In this paper, we propose a new, fully intrinsic and adaptive,
0 0.2 0.4 0.6 0.8 1 10 102 103 variational method, the minimization of which leads to a
c) FM signal and spectrum. (fFM  ∆f ) decomposition of a signal into its principal modes. Indeed, the
current decomposition models are mostly limited by 1) their al-
103
1 gorithmic ad-hoc nature lacking mathematical theory (EMD),
102 2) the recursive sifting in most methods, which does not allow
0
101 for backward error correction, 3) the inability to properly cope
−1
100 1 with noise, 4) the hard band-limits of wavelet approaches,
0 0.2 0.4 0.6 0.8 1 10 102 103 and 5) the requirement of predefining filter bank boundaries
d) AM-FM signal and spectrum. (fAM ∼ fFM ∼ ∆f ) in EWT. In contrast, we propose a variational model that
Fig. 1. AM-FM signals with limited bandwidth. Here, we use a signal f (t) =
determines the relevant bands adaptively, and estimates the
(1 + 0.5 cos(2πfAM t)) · cos(2πfc t + ∆f /fFM cos(2πfFM t)). a) Pure AM corresponding modes concurrently, thus properly balancing
signal. b) Pure FM signal with little but rapid frequency deviations. c) Pure errors between them. Motivated by the narrow-band properties
FM signal with slow but important frequency oscillations. d) Combined AM-
FM signal. The solid vertical line in the spectrum shows the carrier frequency
corresponding to the current common IMF definition, we
fc , the dotted lines correspond to the estimated band limits at fc ± BW/2, look for an ensemble of modes that reconstruct the given
based on (2). input signal optimally (either exactly, or in a least-squares
sense), while each being band-limited about a center frequency
estimated on-line. Here, our variational model specifically can
of an IMF as address the presence of noise in the input signal. Indeed,
BW = 2(∆f + fFM + fAM ). (2) the tight relations to the Wiener filter actually suggest that
our approach has some optimality in dealing with noise. The
Depending on the actual IMF, either of these terms may be variational model assesses the bandwidth of the modes as
dominant. An illustration of four typical cases is provided in H1-norm, after shifting the Hilbert-complemented, analytic
figure 1, where the last example is rather extreme in terms of signal down into baseband by complex harmonic mixing. The
required bandwidth (for illustrational purposes). resulting optimization scheme is very simple and fast: each
Some recent works create a partially variational approach to mode is iteratively updated directly in Fourier domain, as
EMD where the signal is explicitly modeled as an IMF [14]. the narrow-band Wiener filter corresponding to the current
This method still relies on interpolation, selection of a Fourier estimate of the mode’s center-frequency being applied to the
low-pass filter, and sifting of high-frequency components. signal estimation residual of all other modes; then the center
Here, the candidate modes are extracted variationally. The frequency is re-estimated as the center-of-gravity of the mode’s
signal is recursively decomposed into an IMF with TV3- power spectrum. Our quantitative results on tone detection
smooth envelope, and a TV3-smooth residual. The resulting and separation show excellent performance irrespective of
algorithm is very similar to EMD in structure, but somewhat harmonic frequencies, in particular when compared to the
more robust to noise. apparent limits of EMD in this respect. Further, qualitative
A slightly more variational, but still recursive decomposition results on synthetic and real test signals are convincing, also
scheme has been proposed in [15], for the analysis of time- regarding robustness to signal noise.
varying vibration. Here, the dominant vibration is extracted by The rest of this paper is organized as follows: Section
estimating its instantaneous frequency as average frequency II introduces the notions of the Wiener filter, the Hilbert
after the Hilbert transform. Again, this process is repeated transform, and the analytic signal. Also, we briefly review the
recursively on the residual signal. concept of frequency shifting through harmonic mixing. These
An approach based on selecting appropriate wavelet scales, concepts are the very building blocks of our variational mode
DRAGOMIRETSKIY AND ZOSSO: VARIATIONAL MODE DECOMPOSITION 3

decomposition model. Section III presents and explains our B. Hilbert transform and analytic signal
variational model in detail, our algorithm to minimize it, and Here, we cite the definition of the Hilbert transform given
finer technicalities on boundaries, periodicity, and windowing. in [22]:
Section IV contains our experiments and results, namely some
simple quantitative performance evaluations, and comparisons Definition The 1-D Hilbert transform is the linear, shift-
to EMD, and various synthetic multi-mode signals and our invariant operator H that maps all 1-D cosine functions into
method’s decomposition of them. Specifically, tone detection their corresponding sine functions. It is an all-pass filter that
and separation will be analyzed and compared to that of is characterized by the transfer function ĥ(ω) = −j sgn(ω) =
EMD. Additionally, real signals will be considered. Section −jω/|ω|.
V concludes on our proposed variational mode decomposi-
Thus, the Hilbert transform is a multiplier operator in the spec-
tion method, including some future directions and expected
tral domain. The corresponding impulse response is h(t) =
improvements.
1/(πt). Because h(t) is not integrable the integrals defining
the convolution do not converge. Instead, the Hilbert transform
II. T OOLS FROM S IGNAL P ROCESSING Hf (t) of a signal f (t) is therefore obtained as the Cauchy
principal value (denoted p.v.) of :
In this section we briefly review a few concepts and tools 1
Z
f (v)
from signal processing that will constitute the building blocks Hf (t) = p.v. dv. (7)
of our variatonal mode decomposition model. First, we present
π R t−v

a classical case of Wiener filtering for image denoising. Finally, the inverse Hilbert transform is given by its negative,
Next, we describe the Hilbert transform and its use in the H−1 = −H, thus:
construction of a single-side band analytic signal. Finally, we
show how multiplication with pure complex harmonics is used H2 f (t) = −f (t). (8)
to shift the frequencies in a signal. For further properties and analysis of the Hilbert transform,
we refer e.g. to [23]. The most prominent use of the Hilbert
transform is in the construction of an analytic signal from a
A. Gaussian regularizer and Wiener filtering purely real signal, as proposed by Gabor [24].
Let us start with a simple denoising problem. Consider the Definition Let f (t) be a purely real signal. The complex
observed signal f0 (t) to be a copy of the original signal f (t) to analytic signal is now defined as:
be recovered, affected by additive zero-mean Gaussian noise:
fA (t) = f (t) + jHf (t) = A(t)ejφ(t) . (9)
f0 = f + η (3)
This analytic signal has the following important properties.
Recovering the unknown signal f is a typical ill-posed inverse The complex exponential term ejφ(t) is a phasor describing the
problem [17]. If the original signal is known to vary smoothly, rotation of the complex signal in time, φ(t) being the phase,
one would typically write the following Tikhonov regularized while the amplitude is governed by the real envelope A(t).
minimization problem in order to estimate the noise-free signal This representation is particularly useful in the analysis of
[18], [19]: time-varying amplitude and instantaneous frequency, defined
as ω(t) = dφ(t)/dt. The second property is the unilateral
min kf − f0 k22 + αk∂t f k22 (4)

f spectrum of the analytic signal, consisting only of non-negative
frequencies, hence its use in single-sideband modulation. Fi-
This is a standard, Gaussian regularized minimum mean nally, we note that from such an analytical signal, the original
squares, i.e. “L2-H1” problem, of which the Euler-Lagrange real signal is easily retrieved as the real part:
equations are easily obtained as
f (t) = <{fA (t)}. (10)
f − f0 = α∂t2 f. (5)
It is worthwhile highlighting the simple relations between the
These EL equations are typically solved in Fourier domain: Fourier spectra of the real signal and its analytic counterpart,
as defined by (9). First, we recall that the (Fourier) spectrum
fˆ0 of a real signal is a Hermitian function:
fˆ(ω) = , (6)
1 + αω 2
fˆ(−ω) = fˆ(ω). (11)
√ R
where fˆ(ω) := F{f (·)}(ω) := 1/ 2π R f (t)e−jωt dt, with
In contrast, the spectrum of the analytic signal has only non-
j 2 = −1, is the Fourier transform of the signal f (t). Clearly,
negative frequencies. In particular:
the recovered signal f is a low-pass narrow-band selection
of the input signal f0 around ω = 0. Indeed, the solution

0
 ω<0
corresponds to convolution with a Wiener filter, where α fˆA (ω) = fˆ(0) (12)
ω=0
represents the variance of the white noise, and the signal has  ˆ

a lowpass 1/ω 2 power spectrum prior [20], [21]. 2f (ω) ω > 0.
4 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO. Y, ZZZ 20WW

C. Frequency mixing and heterodyne demodulation The reconstruction constraint can be addressed in different
The last concept that we wish to recall before introducing ways. Here, we suggest making use of both a quadratic penalty
the proposed variational mode decomposition, is the principle term and Lagrangian multipliers in order to render the prob-
of frequency mixing. Mixing is the process of combining two lem unconstrained. Therefore, we introduce the augmented
signals non-linearily, thus introducing cross-frequency terms in Lagrangian L as follows [25], [26]:
the output. The simplest mixer is multiplication. Multiplying
2
two real signals with frequencies f1 and f2 , respectively,
  
X j
L(uk , ωk , λ) =α ∂t δ(t) + ∗ uk (t) e−jωk t
creates mixed frequencies in the output at f1 − f2 and f1 + f2 , πt 2
k
which is easily illustrated by the following trigonometric X 2 D X E
identity: + f− uk + λ, f − uk . (17)
2
2 cos(2πf1 t) cos(2πf2 t) = cos(2π(f1 +f2 )t)+cos(2π(f1 −f2 )t).
(13) The solution to the original minimization problem (16) is now
Typical applications are the heterodyne downmixing of the found as the saddle point of the augmented Lagrangian L
modulated high-frequency carrier signal with a local (het- in a sequence of iterative sub-optimizations called alternate
erodyne) oscillator in a radio receiver. In such devices, the direction method of multipliers (ADMM), see algorithm 1. In
selection of either of the two output terms is achieved by the next paragraphs, we detail how the respective sub-problems
filtering. Here, instead of filtering the output, we mix the two can be solved.
respective analytic signals:
Algorithm 1 ADMM optimization concept for VMD
ej2πf1 t ej2πf2 t = ej2π(f1 +f2 )t , (14)
Initialize u1k , ωk1 , λ1 , n ← 0
i.e., the mixed signal is automatically “mono-tone” (consti- repeat
tuted of a single frequency only). In Fourier terms, this is n←n+1
well known as the following transform pair: for k = 1 : K do
F Update uk :
fA (t)e−jω0 t ←→ fˆA (ω) ∗ δ(ω + ω0 ) = fˆA (ω + ω0 ), (15)
n+1
uk ← arg min L(un+1 1 , . . . , un+1 n n
k−1 , uk , uk+1 , . . . , uK ,
where δ is the Dirac distribution and ∗ denotes convolution. uk

Thus, multiplying an analytic signal with a pure exponential ω1n , . . . , ωK n

, λn ) (18)
results in simple frequency shifting.
end for
for k = 1 : K do
III. VARIATIONAL M ODE D ECOMPOSITION Update ωk :
In this section we introduce our proposed model for vari-
ωkn+1 ← arg min L(un+1 , . . . , un+1
K , ω1
n+1 n+1
, . . . , ωk−1 ,
ational mode decomposition, essentially based on the three ωk
1

concepts outlined in the previous section. ωk , ωk+1n n

, . . . , ωK , λn ) (19)
The goal of VMD is to decompose an input signal into
a discrete number of sub-signals (modes), that have specific end for
sparsity properties while reproducing the input. Here, the Dual ascent:
sparsity prior of each mode is chosen to be its bandwidth
!
X
in spectral domain. In other words, we require each mode k λn+1 n
←λ +τ f − ukn+1
(20)
to be mostly compact around a center pulsation ωk , which is k
to be determined along with the decomposition. until convergence: k kun+1
P
− unk k22 /kunk k22 < .
k
In order to assess the bandwidth of a mode, we propose
the following scheme: 1) for each mode uk , compute the
associated analytic signal by means of the Hilbert transform
in order to obtain a unilateral frequency spectrum. 2) for each
mode, shift the mode’s frequency spectrum to “baseband”, by
mixing with an exponential tuned to the respective estimated A. Minimization w.r.t. uk
center frequency. 3) The bandwidth is now estimated through
the H1 Gaussian smoothness of the demodulated signal, i.e. To update the modes uk , we first rewrite the subproblem
the squared L2-norm of the gradient. The resulting constrained (18) as the following equivalent minimization problem:
variational problem is the following: (    2
(
2
) j
n+1
∗ uk (t) e−jωk t
  
X j −jωk t uk = arg min α ∂t δ(t) +
min ∂t δ(t) + ∗ uk (t) e uk ∈R πt 2
uk ,ωk πt 2
k 2
)
X X λ
s.t. uk = f (16) + f− ui + . (21)
k
2 2
DRAGOMIRETSKIY AND ZOSSO: VARIATIONAL MODE DECOMPOSITION 5

Making use of the Parseval/Plancherel Fourier isometry under which puts the new ωk at the center of gravity of the corre-
the L2 norm, this problem can be solved in spectral domain: sponding mode’s power spectrum. This mean carrier frequency
( is the frequency of a least squares linear regression to the
ûn+1 = arg min α kjω [(1 + sgn(ω + ωk ))ûk (ω + ωk )]k2
2 instantaneous phase observed in the mode.
k
ûk ,ûk =û∗
k Plugging the solutions of the sub-optimizations into the
ADMM algorithm 1, and directly optimizing in Fourier do-

2
λ̂ 
main where appropriate, we get the complete algorithm for
+ fˆ −
X
ûi + . (22)
2  variational mode decomposition, summarized in algorithm 2.
2

We now perform a change of variables ω → ω + ωk in the Algorithm 2 Complete optimization of VMD

first term:
( Initialize û1k , ωk1 , λ̂1 , n ← 0
2 repeat
ûn+1
k = arg min α kj(ω − ωk ) [(1 + sgn(ω))ûk (ω)]k2
ûk ,ûk =û∗
k
n←n+1
2
 for k = 1 : K do
λ̂ Update ûk for all ω ≥ 0:

+ fˆ −
X
ûi + . (23)
2 λ̂n
fˆ − i<k ûn+1 − i>k ûni +
 P P
2
i 2
ûn+1
k ← (29)
Exploiting the Hermitian symmetry of the real signals in the 1 + 2α(ω − ωkn )2
reconstruction fidelity term, we can write both terms as half- Update ωk :
space integrals over the non-negative frequencies: R∞
ω|ûn+1
k (ω)|2 dω
ωkn+1 ← R0 ∞ n+1 (30)
Z ∞
ûkn+1 = arg min 4α(ω − ωk )2 |ûk (ω)|2 (24) 0
|ûk (ω)|2 dω
ûk ,ûk =û∗
k 0
!2  end for
λ̂ Dual ascent for all ω ≥ 0:

+ 2 fˆ −
X
ûi + dω .
2  !
fˆ −
X
λ̂n+1 n
← λ̂ + τ ûn+1
k (31)
The solution of this quadratic optimization problem is readily k
found by letting the first variation vanish for the positive
until convergence: kûn+1 − ûnk k22 /kûnk k22 < .
P
frequencies: k k
 
λ̂ 1
= fˆ −
X
ûn+1
k ûi +  , (25)
2 1 + 2α(ω − ωk )2 C. Inexact reconstruction and denoising
i6=k

which is clearly identified as a Wiener filtering of the current Here, the role of the Lagrangian multiplier is to enforce the
residual, with signal prior 1/(ω − ωk )2 . The full spectrum of constraint, while the quadratic penalty improves convergence.
the real mode is then simply obtained by Hermitian symmetric If exact reconstruction is not required, but some slack is to be
completion. Conversely, the mode in time domain is obtained allowed, using the quadratic penalty only while dropping the
as the real part of the inverse Fourier transform of this filtered Lagrangian multiplier would be the appropriate choice. Indeed,
analytic signal. the quadratic penalty on its own represents the least-squares
fidelity prior associated with additive Gaussian noise.

B. Minimization w.r.t. ωk
D. On boundaries, periodicity, and windowing
The center frequencies ωk do not appear in the reconstruc-
tion fidelity term, but only in the bandwidth prior. The relevant Up until now, the signals f and the modes uk have been
problem thus writes: considered continuous over the whole axis t ∈ R. However,
( ) in signal processing we are much more likely to be working
2
with signals that are both finite in time and resolution. Let
  
n+1 j −jωk t
ωk = arg min ∂t δ(t) + ∗ uk (t) e . us say we restrict the time window to t ∈ [0, 1]. Luckily
ωk πt 2
(26) the results presented so far equally hold for discrete, finite
As before, the optimization can take place in Fourier domain, time signals, where simply the continuous Fourier transform
and we end up optimizing: is replaced by its discrete counterpart. The only problems arise
Z ∞  at the boundaries of the signal.
n+1
ωk = arg min 2 2
(ω − ωk ) |ûk (ω)| dω , (27) Indeed, when considering short-time signals, the implicit
ωk 0 assumption here is that the signal considered is just a one-
This quadratic problem is easily solved as: period extract of an infinitely long, periodic signal. Conse-
R∞ quently, the spectrum of a seemingly simple “general trend”-
ω|ûk (ω)|2 dω function on a short interval, say f : [0, 1] 7→ R : f (t) = t,
ωkn+1 = R0 ∞ , (28)
0
|ûk (ω)|2 dω contains an important amount of high-frequency harmonics,
6 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO. Y, ZZZ 20WW

since we are effectively looking at the spectrum of the periodic

sawtooth function. Conversely in time domain, we realize that 100 100
at the endpoints of the domain, the periodized function is
discontinuous, thus severely affecting the H1 smoothing term. 10−3 10−3
There are two remedies to this. Ideally, one should exclude
10−6 −2 10−6 −2
the boundaries of the domain in the evaluation of the smooth- 10 10−1 10 10−1
ness, i.e. restrict its evaluation to the open interval (0, 1). a)  = 1e − 3 b)  = 1e − 5
However, this clearly breaks the Parseval/Plancherel Fourier 0 0
isometry and the whole beauty of the spectral solution is 10 10
lost. Therefore, we suggest a less far-reaching remedy, that
10−3 10−3
is classically used in short-term Fourier analysis: smooth win-
dowing. This approach is particularly useful in cases, where 10−6 −2 10−6 −2
the variational mode decomposition is anyway performed on 10 10−1 10 10−1
short chunks of a much longer time series signal. c)  = 1e − 7 d)  = 1e − 9
For simplicity, in the following examples, we will use
Fig. 2. Mode decomposition of a pure harmonic: Relative error for a range
a Gaussian window. This window is simply applied to the of 257 frequencies, for different convergence tolerance levels . The relative
input signal f prior to performing the VMD algorithm. After error does not correlate with the tone frequency. Further, reconstruction error
decomposition, the individual modes can be “unwindowed” by can be controlled by decreasing the stopping criterion’s convergence tolerance,
except for frequencies very close to the Nyquist frequency. In contrast, EMD’s
simple division. This, however, will largely affect reconstruc- relative tone reconstruction error is bounded by a quadratic increase with
tion fidelity close to the window borders. This is particularly frequency (dotted line) [2].
apparent in the single frame decomposition. In a sliding win-
dow short-time analysis of a larger time series signal, however,
instead of window division, the modes can be stitched together
by simple addition without error amplification. ν2 < ν1

ν2 < ν1
IV. E XPERIMENTS AND R ESULTS
In this section, we apply the proposed VMD algorithm
to a series of test signals in order to assess the validity of
our approach. First, we focus on a few problems that have
been successfully employed for highlighting the strengths and
ν1 ν1
shortcomings of the EMD / Hilbert-Huang-Transform, namely
tones versus sampling, and tones separation [2]. Then we a) EMD, ρ = 1/4 b) VMD, ρ = 1/4
briefly investigate noise robustness of VMD. Finally, we shift
our attention to more complex signals, which have already
been used in [14] and [12].
ν2 < ν1

ν2 < ν1

A. Tones and sampling

When the input signal f = fν (t) = cos(2πνt) is composed
of a pure harmonic, then the mode decomposition is expected
ν1 ν1
to output exactly this harmonic. As reported in [2], this does
not happen to be the case with EMD, since the local extrema c) EMD, ρ = 1 d) VMD, ρ = 1
can suffer from important jittering with increasing frequency.
In [2], the relative error
ν2 < ν1

ν2 < ν1

e(ν) = kfν (t) − u1 (t)k2 /kfν (t)k2 (32)

was introduced, and a quadratic increase with frequency of
an upper bound to this relative error was reported for EMD.
Further, EMD has pronounced spikes of near-perfect recon-
struction when the sampling frequency is an even multiple of ν1 ν1
the tone’s frequency. e) EMD, ρ = 4 f) VMD, ρ = 4
Here, we perform this analysis for the proposed VMD
Fig. 3. Tones separation. In a superposition of two tones of frequencies
model. We refrain from windowing and consider exactly ν2 < ν1 < fs /2 and equal amplitudes, the mode decompositions between
the same signals as for the EMD analysis. The results for EMD and VMD vary significantly. The plot indicates relative error, with
different convergence tolerance levels  are shown in figure 2. values between 0 (white) and 0.5 (black). a,c,e) EMD has important areas
of confusion (dark), where the tones cannot be separated correctly [2]. b,d,f)
It can be clearly seen that the relative reconstruction error is In contrast, VMD achieves good tones separation almost everywhere but for
largely independent of the harmonic’s frequency. Moreover, ν1 too close to the Nyquist frequency.
the relative error is nicely controlled by the tolerance level .
DRAGOMIRETSKIY AND ZOSSO: VARIATIONAL MODE DECOMPOSITION 7

1 1 1 0.2
0 0 0 0
−1 −1
−1 −0.2
−2 −2
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
a) fn (t) b) u1 (t) a) fn (t) b) u1 (t)
0.2 0.1 0.2 0.4
0 0 0.1 0.2
−0.2 −0.1 0 0
−0.4 −0.2 −0.1 −0.2
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −0.2 −0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
c) u2 (t) d) u3 (t)
c) u2 (t) d) u3 (t)
104 0.4 0.5
0.2
10 1
0 0
−0.2 −0.5
10−2 1 2 3 −0.4
10 10 10 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
e) |fˆn |(ω) e) u4 (t) f) u5 (t)
1
104 0.5 0.2
0 0.1
101 −0.5 0
−1 −0.1
10−2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
101 102 103
g) u6 (t) h) u7 (t)
e) |ûk |(ω)
Fig. 5. EMD decomposition of noisy tri-harmonic. (a) noisy input signal.
Fig. 4. VMD decomposition of noisy tri-harmonic. (a) The noisy input signal. (b)-(h) The seven modes extracted by EMD. None of the modes corresponds
(b)-(d) The three modes extracted by denoising VMD, and the theoretical to a pure harmonic.
mode (dotted). (e) The spectrum of the input signal, and (f) its distribution
over the three modes.

where η ∼ N (0, 1) represents the Gaussian additive noise. The

B. Tones separation noise level is quite important with respect to the amplitude of
The next slightly more complicated challenge is the separa- the highest harmonic. We perform variational modes decom-
tion of two different superimposed tones [2]. Here, the input position into three modes, without Lagrangian multipliers in
signal is composed of two different, pure harmonics: order to remove the noise. The signal, and the three compo-
nents estimated using VMD are shown in Fig. 4. The strong,
fν1 ,ν2 (t) = a1 cos(2πν1 t) + a2 cos(2πν2 t), (33) lowest frequency signal is recovered almost flawlessly. The
with ν2 < ν1 < fs /2, and a1,2 two possibly different medium-strength medium-frequency signal is still detected at
amplitudes. As a function of the amplitude ratio ρ = a1 /a2 , acceptable quality. The weak, high-frequency signal, however,
EMD exhibits different, important regions of confusion, where is difficult. The VMD algorithm correctly tunes the third
the two signals are too close in frequency to be separated center-frequency on this harmonic, but the recovered mode is
correctly, as reported in [2] and illustrated in figure 3. highly affected by the noise. Here, decreasing the bandwidth
Again, we apply the same analysis to the proposed VMD by increasing α comes at the risk of not properly capturing
model, and again we do not employ any windowing. The the correct center frequency, while too low an α includes
results for varying amplitude ratios ρ ∈ {1/4, 1, 4} are shown more noise in the estimated mode. The mode could, however,
in figure 3 along with the corresponding EMD results. As can be cleaned further in post-processing. For reference, we note
be clearly seen, the proposed VMD achieves good tones sepa- that the estimated VMD center frequencies are off by 0.27%,
ration over the whole domain except at the Nyquist frequency. 1.11% and 0.18% only.
In particular, the decomposition quality is not significantly We provide a comparison with EMD1 based on exactly the
worse for close harmonics. same signal in Fig. 5. The EMD produces 7 estimated modes.
The first two modes contain the highest-frequency harmonic,
C. Noise robustness and important amounts of noise. The forth mode comes closest
to the middle harmonic, however important features have been
To illustrate the VMD robustness with respect to noise in the
attributed to the third and fifth mode. The sixth mode picks up
input signal, we test using the following tri-harmonic signal,
most of the low frequency harmonic, but is severely distorted.
affected by noise:
1 1 1 Implementation by Gabriel Rilling, available at [Link]
fn (t) = cos(4πt)+ cos(48πt)+ cos(576πt)+0.1η, (34) [Link]/[Link]
4 16
8 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO. Y, ZZZ 20WW

8 6 1
6 4 4
4 0.8
2 2 0.6 2
0 0 0.4 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
a) fSig1 (t) b) 6t a) w(t) b) fSig1
w
(t)
1 0.4
0.5 0.2 104
0 0
−0.5 −0.2 101
−1 −0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 10−2
c) cos(8πt) d) 0.5 cos(40πt) 101 102 103
104 c) |fˆSig1
w
|(ω)
800

Iterations
101 600
400
10−2 200
101 102 103 0
e) |fˆSig1 |(ω) 101 102 103

Fig. 6. a) fSig1 (t), b–d) its constituent modes. e) The signal’s spectrum.
d) ωk

104

101
D. Complex multimode signals
10−2
Now we look at slightly more complex signals to be 101 102 103
decomposed. In particular, we consider the same test signals e) |ûw
k |(ω)
that were previously suggested in [14] and also used in [12], 4
with the purpose of increased comparability. 3 1
1) Example 1: The first signal is a composition of three 2 0
simple components, namely a general linear trend and two 1
different harmonics: −1
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
f) uw
1 (t) g) uw
2 (t)
fSig1 (t) = 6t + cos(8πt) + 0.5 cos(40πt). (35) 0.5
0
The signal, its three constituent modes, and the composite −0.5
Fourier spectrum are shown in figure 6. The main challenge 0 0.2 0.4 0.6 0.8 1
of this signal is the linear growth term. Without windowing, h) uw
3 (t)
the higher order harmonics of the periodized sawtooth signal
6
spread over the whole spectrum. 2
4
In order to reduce the effects of periodization, we apply 0
Gaussian windowing. The corresponding windowed signal, 2
and the respective VMD results are illustrated in detail in 0 −2
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
figure 7. In particular, we show how the two non-zero center
frequencies ω2 and ω3 quickly converge towards the exact har- i) u1 (t) j) u2 (t)
monics. The corresponding modes constitute a nice partition 2
of the input spectrum, with each mode being clearly dominant 1
around its respective center frequency. The three modes in 0
time domain show nice separation into three distinct signals of −1
−2
characteristic oscillations. After “unwindowing” by pointwise 0 0.2 0.4 0.6 0.8 1
division of the estimated modes by the Gaussian window, we k) u3 (t)
recover good estimates of the true underlying modes (dotted
Fig. 7. VMD decomposition of fSig1 . a) The applied window, b) the
lines), valid on the central 60% of the signal. windowed signal, and c) its spectrum. d) Evolution of the detected center
2) Example 2: The second example uses a quadratic trend, frequencies, and e) the corresponding spectrum decomposition. f–h) the
reconstructed modes prior to, and i–k) after Gaussian window removal.
a chirp signal, and a third mode with sharp transition between
DRAGOMIRETSKIY AND ZOSSO: VARIATIONAL MODE DECOMPOSITION 9

8 6 1 4
6 0.8
4 4 0.6 2
2 2 0.4
0 0.2 0
−2 0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
a) fSig2 (t) b) 6t 2
a) w(t) b) fSig2
w
(t)
1 1
0.5 0.5 104
0 0
−0.5 −0.5 101
−1 −1
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 10−2
101 102 103
c) cos(10πt + 10πt2 ) d)
c) |fˆSig2
w
|(ω)
(
cos(60πt) t ≤ 0.5
cos(80πt − 10π) t > 0.5 200

Iterations
104 150
100
101 50
101 102 103
10−2 1 2 3
10 10 10 d) ωk
e) |fˆSig2 |(ω) 104
Fig. 8. a) fSig2 , b–d) its constituent modes. e) The signal’s spectrum.
101

10−2
two constant frequencies2 : 101 102 103
( e) |ûw
k |(ω)
cos(60πt) t ≤ 0.5 2
fSig2 (t) = 6t2 +cos(10πt+10πt2 )+ 1
cos(80πt − 10π) t > 0.5 1
(36) 0
The signal, its three constituent modes, and the composite 0
−1
Fourier spectrum are shown in figure 8. The instantaneous 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
frequency of the chirp is given by the time derivative of its f) uw g) uw
1 (t) 2 (t)
phase:
1 1
ω(t) := ∂t φ(t) = 10π + 20πt. (37) 0.5 0.5
0 0
−0.5 −0.5
Thus, for t ∈ [0, 1] the instantaneous frequency varies linearly −1 −1
between 10π an 30π. Consequently, the theoretical center 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
frequency of the mode is located at 20π. The piecewise- h) uw
3 (t) i) uw
4 (t)
constant bi-harmonic has spectral peaks expected at 60π and 6
80π. 1
4
Here, too, we employ Gaussian windowing to alleviate 2 0
periodization artifacts. Indeed, the windowed signal has a 0 −1
much cleaner spectrum, and the expected peaks of the signal’s 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
components become more prominent, as illustrated in figure 9. j) u1 (t) k) u2 (t)
Again, the estimated center frequencies ωk converge to the
2
expected frequencies precisely. Here, we chose to decompose 1
1
into four modes, thus assigning each half of the piecewise- 0
constant frequency signal to a separate mode. The spectral 0
−1 −1
partitioning can be nicely appreciated in the spectral plot of
the different modes. The unwindowed mode estimates fit well 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
the theoretical signals, except again for boundary issues. l) u3 (t) m) u4 (t)
Fig. 9. Results of VMD on fSig2 . a) The applied window, b) the windowed
2 Here,
signal, and c) its spectrum. d) Evolution of the detected center frequencies, and
we changed the phase shift in the third component, with piecewise- e) the corresponding spectrum decomposition. f–i) the reconstructed modes
constant frequency, from 15π to 10π, in order to have a continuous signal. prior to, and j–m) after Gaussian window removal.
10 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO. Y, ZZZ 20WW

6 6 1 6
4 4 0.8 4
2 0.6 2
0 2
0.4 0
−2 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
a) fSig3 (t) b) (1.2 + cos(2πt)) −1
a) w(t) b) fSig3
w
(t)
2
1 104
0
−1 101
−2
0 0.2 0.4 0.6 0.8 1 10−2
101 102 103
c) (1.5 + sin(2πt))−1 cos(32πt + 0.2 cos(64πt))
c) |fˆSig3
w
|(ω)
104
300

Iterations
101 200
10−2 100
101 102 103
101 102 103
d) |fˆSig2 |(ω)
d) ωk
Fig. 10. a) fSig3 , b–c) its constituent modes. d) The signal’s spectrum.
104

3) Example 3: The third synthetic signal has intrawave 101

frequency modulation:
10−2
1 cos(32πt + 0.2 cos(64πt)) 101 102 103
fSig3 (t) = + .
1.2 + cos(2πt) 1.5 + sin(2πt) e) |ûw
k |(ω)
(38) 2
The signal, its three constituent modes, and the composite 4 1
Fourier spectrum are shown in figure 10. While the first, bell- 0
2
shaped component has mostly low-pass content, the second −1
mode’s main peak is clearly identified at 32π. However, 0 −2
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
due to the non-linear intrawave frequency modulation, an
important amount of higher-order harmonics are also observed f) uw
1 (t) g) uw
2 (t)
at 32π +64π = 96π, 32π +2·64π = 160π and 32π +3·64π = 6 2
224π, respectively. This second component obviously violates 4 1
the narrowband assumption, and one would naturally expect 0
2 −1
some difficulties recovering this mode using VMD. Indeed,
0 −2
by Carson’s rule, the mode’s bandwidth here is dominantly 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
controlled by the relatively high frequency of the modulating h) u1 (t) i) u2 (t)
term cos(64πt), essentially spreading the mode over the whole
practical spectrum. Fig. 11. Results of VMD on fSig3 . a) The applied window, b) the windowed
signal, and c) its spectrum. d) Evolution of the detected center frequencies, and
The slightly windowed signal and the corresponding VMD e) the corresponding spectrum decomposition. f–g) the reconstructed modes
results are illustrated in figure 11. The non-zero ω2 quickly prior to, and h–i) after Gaussian window removal.
converges to the correct main frequency 32π. The higher order
harmonics are not uniquely attributed to the second mode,
but shared between both modes. Consequently, the intrawave the expected spikes-train driven by the rhythm of the heartbeat,
frequency modulation is shared by both modes, creating some one can clearly see an oscillating low-frequency pattern. At
ripples in the otherwise low-frequency mode. Nonetheless, the the other end of the spectrum, there is distinct high-frequency
reconstructed estimated modes fit well the constituent signals noise at a single high-pitch harmonic, most likely the electric
(dotted lines). Most of the error occurs at the boundaries, and power-line frequency. The distinct spikes of the ECG signal
at the very center of the signal, where the low-frequency mode create important higher-order harmonics.
has a sharp peak, involving some higher frequency features The spectrum after slight Gaussian windowing, and the
wrongly attributed to the higher-frequency mode. results of VMD are depicted in figure 13. We chose a high-
4) Example 4: The forth example is a real signal from number of 10 modes to be detected, to accommodate the
an electrocardiogram (ECG), data shared by [12]. These data numerous higher-order harmonics of the spikes. The respective
present numerous components, as seen in figure 12. Beyond center frequencies nicely converge to these spectral peaks. The
DRAGOMIRETSKIY AND ZOSSO: VARIATIONAL MODE DECOMPOSITION 11

6 4
4 2 1
2 4
0 0.8 2
0
−2 0.6 0
−2
0 0.2 0.4 0.6 0.8 1 0.5 0.55 0.6 0.4 −2
a) fECG (t) b) fECG (t) (Detail) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
104 a) w(t) b) fECG
w
(t)

101 104

10−2 101
1 2 3
10 10 10
10−2
c) |fˆECG |(ω) 101 102 103
Fig. 12. a) ECG signal 7. b) Detail. c) The signal’s spectrum. c) |fˆECG
w
|(ω)
400

Iterations
300
first, low-frequency mode captures the low-frequency oscil-
200
lation of the baseline. The highest frequency mode contains 100
the most noise. The first actual ECG specific mode oscillates
precisely at the frequency of the heartbeat. The higher ECG 101 102 103
modes then contain the higher-order wave-packages around d) ωk
the highly non-sinusoidal spikes. A “clean” ECG signal can 1 0.5
be reconstructed by summing all but the first and last VMD 0
0
modes, thus discarding the low-frequency baseline oscillation −0.5
and most of the high-frequency noise. −1
−1
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
e) uw
1 (t) f) uw
2 (t)
V. C ONCLUSIONS AND O UTLOOK
0.5 0.5
In this paper, we have presented a novel variational method 0
for decomposing a signal into an ensemble of band-limited 0
−0.5
intrinsic mode functions, that we call Variational Mode De- −0.5
composition, (VMD). In contrast to existing decomposition 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
models, like the empirical mode decomposition (EMD), we g) uw
3 (t) h) uw
4 (t)
refrain from modeling the individual modes as signals with 0.5 0.5
explicit IMFs. Instead, we replace the most recent definition
0 0
of IMFs, namely their characteristic description as AM-FM
signals, by the corresponding narrow-band property. Indeed, −0.5 −0.5
we provide a formula that relates the parameters of the explicit 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
AM-FM descriptors to the estimated signal bandwidth. i) uw j) uw
5 (t) 6 (t)
Our decomposition model solves the inverse problem as fol-
0.5 0.4
lows: decompose a signal into a given number of modes, either 0.2
exactly or in a least squares sense, such that each individual 0 0
mode has limited bandwidth. We assess the mode’s bandwidth −0.5 −0.2
−0.4
as the squared H1 norm of its Hilbert complemented analytic
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
signal with only positive frequencies, shifted to baseband
by mixing with a complex exponential of the current center k) uw
7 (t) l) uw
8 (t)
frequency estimate. The variational problem is solved very 4
4
efficiently in a classical ADMM approach: The modes are 2 2
updated by simple Wiener filtering, directly in Fourier domain
0 0
with a filter tuned to the current center frequency, then the
−2
center frequencies are updated as the center of gravity of the 0 0.2 0.4 0.6 0.8 1 0.5 0.55 0.6
mode’s power spectrum, and finally the Lagrangian multiplier P8 P8
m) k=2 uk (t) n) uk (t) (Detail)
enforcing exact signal reconstruction is updated as dual ascent. k=2
Fig. 13. Results of ECG signal 7. a) The applied window, b) the windowed
In our experiments, we show that the proposed VMD signal, and c) its spectrum. d) Evolution of the detected center frequencies.
scheme clearly outperforms EMD with regards to tone de- e–l) The reconstructed modes prior to Gaussian window removal. m) Cleaned
tection, tone separation, and noise robustness. Further, we ECG, and n) detail.
apply our model to more complicated signals for comparison
12 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. X, NO. Y, ZZZ 20WW

with other state-of-the-art methods, and can show successful [13] J. Carson, “Notes on the Theory of Modulation,” Proceedings of the
decomposition. IRE, vol. 10, no. 1, pp. 57–64, Feb. 1922.
The most important limitation of the proposed VMD is [14] T. Y. Hou and Z. Shi, “Adaptive Data Analysis via Sparse Time-
Frequency Representation,” Advances in Adaptive Data Analysis,
with boundary effects, and sudden signal onset in general. vol. 03, no. 1 & 2, pp. 1–28, Apr. 2011.
This is strongly related to the use of an L2-based smoothness [15] M. Feldman, “Time-varying vibration decomposition and analysis based
term, that overly penalizes jumps at the domain borders and on the Hilbert transform,” Journal of Sound and Vibration, vol. 295, no.
3-5, pp. 518–530, Aug. 2006.
within; conversely, this is also reflected by implicit periodicity [16] H.-T. Wu, P. Flandrin, and I. Daubechies, “One or Two Frequencies?
assumptions when optimizing in Fourier domain, and by the the Synchrosqueezing Answers,” Advances in Adaptive Data Analysis,
narrow-band violation caused by discontinuous envelopes in vol. 03, no. 01n02, pp. 29–39, Apr. 2011.
[17] M. Bertero, T. A. Poggio, and V. Torre, “Ill-Posed Problems in Early
such AM-FM signals. Another point that critics might high- Vision,” in Proceedings of the IEEE, vol. 76, no. 8, 1988, pp. 869–889.
light, is the required explicit (manual) selection of the number [18] A. N. Tichonov, “Solution of incorrectly formulated problems and the
of active modes in the decomposition, like in EWT but as regularization method,” Soviet Mathematics, vol. 4, pp. 1035–1038,
1963.
opposed to EMD. Current work addresses these shortcomings, [19] V. A. Morozov, “Linear and nonlinear ill-posed problems,” Journal of
and we are also working on suitable extensions to signals on Mathematical Sciences, vol. II, no. 6, pp. 706–736, 1975.
domains of dimension greater than one. [20] N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary
Time Series, 1949.
[21] R. C. Gonzalez and R. E. Woods, Digital Image Processing. Addison
ACKNOWLEDGMENTS Wesley, 1992.
[22] M. Unser, D. Sage, and D. Van De Ville, “Multiresolution Monogenic
The authors gratefully acknowledge valuable discussions Signal Analysis Using the Riesz Laplace Wavelet Transform,” IEEE
with and input from: Prof. A. L. Bertozzi, Dr. J. Gilles, and Transactions on Image Processing, vol. 18, no. 11, pp. 2402–2418, 2009.
Dr. Y. van Gennip. [23] S. L. Hahn, Hilbert transforms in signal processing. Artech House,
Inc., 1996.
[24] D. Gabor, “Theory of Communication,” Journal of the Institution of
R EFERENCES Electrical Engineers - Part III: Radio and Communication Engineering,
vol. 93, no. 26, pp. 429–457, 1946.
[1] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-
C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition [25] D. P. Bertsekas, “Multiplier methods: A survey,” Automatica, vol. 12,
and the Hilbert spectrum for nonlinear and non-stationary time series no. 2, pp. 133–145, 1976.
[26] J. Nocedal and S. J. Wright, Numerical optimization, 2nd ed. Springer,
analysis,” Proceedings of the Royal Society A: Mathematical, Physical
and Engineering Sciences, vol. 454, no. 1971, pp. 903–995, Mar. 1998. Berlin, 2006.
[2] G. Rilling, P. Flandrin, and P. Gonçalvès, “On empirical mode decom-
position and its algorithms,” in IEEE-EURASIP workshop on Nonlinear
Signal and Image Processing, 2003.
[3] G. Rilling and P. Flandrin, “One or Two Frequencies? The Empirical
Mode Decomposition Answers,” IEEE Transactions on Signal Process- Konstantin Dragomiretskiy received the
ing, vol. 56, no. 1, pp. 85–95, Jan. 2008. B.S. (Hons.) degree in mathematics and the
[4] P. Flandrin, P. Gonçalvès, and G. Rilling, “EMD equivalent filter banks, B.A. degree in economics from the University of
from interpretation to applications,” in Hilbert-Huang Transform and Its California, San Diego, CA, USA, in 2010.
Applications, 2005, pp. 57–74. He was an Intern Researcher at Sun
[5] N. Klügel, “Practical Empirical Mode Decomposition for Audio Syn- Microsystems, Menlo Park, CA, USA, in 2006,
thesis,” in Int. Conference on Digital Audio Effects (DAFx-12), no. 2, and an Undergraduate Researcher at California
2012, pp. 15–18. State University, Chico, CA, USA, in 2008. He is
[6] B. Barnhart and W. Eichinger, “Empirical Mode Decomposition applied currently a Research Assistant at the Department
to solar irradiance, global temperature, sunspot number, and CO2 con- of Mathematics at the University of California,
centration data,” Journal of Atmospheric and Solar-Terrestrial Physics, Los Angeles, CA, USA, working towards his Ph.D.
vol. 73, no. 13, pp. 1771–1779, Aug. 2011. degree with Prof. Andrea L. Bertozzi. His current research interests include
[7] S. Assous, A. Humeau, and J.-P. L’huillier, “Empirical mode decompo- variational and PDE based methods applied to signal and image processing
sition applied to laser Doppler flowmetry signals : diagnosis approach.” problems.
IEEE Engineering in Medicine and Biology Conference (EMBC), vol. 2,
pp. 1232–5, Jan. 2005.
[8] A. O. Andrade, S. Nasuto, P. Kyberd, C. M. Sweeney-Reed, and F. Van
Kanijn, “EMG signal filtering based on Empirical Mode Decomposi-
tion,” Biomedical Signal Processing and Control, vol. 1, no. 1, pp. 44– Dominique Zosso (S’06–M’11) received the [Link].
55, Jan. 2006. degree in electrical and electronics engineering and
[9] S. Liu, Q. He, R. X. Gao, and P. Freedson, “Empirical mode decompo- the Ph.D. degree from École Polytechnique Fédérale
sition applied to tissue artifact removal from respiratory signal.” IEEE de Lausanne (EPFL), Lausanne, Switzerland, in
Engineering in Medicine and Biology Conference (EMBC), pp. 3624– 2006 and 2011, respectively.
3627, Jan. 2008. He was a Researcher with the Structural Bioin-
[10] I. Mostafanezhad, O. Boric-Lubecke, V. Lubecke, and D. P. Mandic, formatics Group at the Swiss Institute of Bioinfor-
“Application of empirical mode decomposition in removing fidgeting matics and Biozentrum, University of Basel, Basel,
interference in doppler radar life signs monitoring devices,” IEEE Switzerland, from 2006 to 2007. He was Research
Engineering in Medicine and Biology Conference (EMBC), pp. 340– and Teaching Assistant at the Signal Processing Lab-
343, Jan. 2009. oratory, EPFL, from 2007 to 2012. He is currently
[11] I. Daubechies, J. Lu, and H.-T. Wu, “Synchrosqueezed wavelet trans- a Post-Doctoral Scholar with the Department of Mathematics, University
forms: An empirical mode decomposition-like tool,” Applied and Com- of California, Los Angeles, CA, USA, with Luminita A. Vese, Andrea L.
putational Harmonic Analysis, vol. 30, no. 2, pp. 243–261, Mar. 2011. Bertozzi and Stanley J. Osher. His current research interests include PDE and
[12] J. Gilles, “Empirical Wavelet Transform,” IEEE Transactions on Signal variational models for inverse problems in computer vision, signal and image
Processing. processing, and efficient algorithms to solve them.