Bivariate Empirical Mode Decomposition

TL;DR: The empirical mode decomposition is extended to bivariate time series that generalizes the rationale underlying the EMD to the bivariate framework and is designed to extract zero-mean rotating components.

Abstract: The empirical mode decomposition (EMD) has been introduced quite recently to adaptively decompose nonstationary and/or nonlinear time series. The method being initially limited to real-valued time series, we propose here an extension to bivariate (or complex-valued) time series that generalizes the rationale underlying the EMD to the bivariate framework. Where the EMD extracts zero-mean oscillating components, the proposed bivariate extension is designed to extract zero-mean rotating components. The method is illustrated on a real-world signal, and properties of the output components are discussed. Free Matlab/C codes are available at http://perso.ens-lyon.fr/patrick.flandrin.

Summary

A. Classical EMD

• Basically, the EMD considers a signal at the scale of its local oscillations.
• The main idea of EMD is then to formalize the idea that, locally: "signal = fast oscillations superimposed on slow oscillations".
• Practically, this primarily implies that all its maxima are positive and all its minima are negative.
• The discrimination between "fast" and "slow" oscillations is obtained through an algorithm referred to as the sifting process [1] which iterates a nonlinear elementary operator S on the signal until some stopping criterion is met.

B. Envelopes in 3 dimensions

• The EMD is based on the intuitive notion of "oscillation" which naturally relates to local extrema.
• In order to separate the more rapidly rotating component from slower ones, the idea is once again to define the slowly rotating component as the mean of some "envelope".
• In practice, the top point, for example, is uniquely defined only when the signal reaches a local maximum in the vertical direction and is therefore tangent to the top of the tube.
• In practice, however, the second scheme may be preferred because it is naturally more robust to sampling errors.
• The desired goal concerning the interpolation is the same as for the classical EMD: a smooth interpolation with as few "spurious bumps" as possible.

Time (days)

• Fig. 3 . A signal and its Bivariate Empirical Mode Decomposition.
• Indeed, the method operating at a local scale, its behavior on signals whose properties evolve slowly with respect to the local period is very similar to its operation on exactly periodic signals of constant amplitude.
• It is worth noting that examples of the second type of solutions are generally encountered when the analyzed signal does not clearly contain rotating components, as in e.g. a complexvalued white gaussian noise signal.
• In the limit where the number of directions tends to infinity, this results in EQUATION which is simply the mean of the signal over a period weighted by dψ(t)/dt > 0, where the weighting conveys the fact that the distribution of sampling points on the tube section is denser where the curvature is larger.
• Likewise, the same reasoning for the second algorithm results in EQUATION and hence m EQUATION.

IV. ILLUSTRATION

• The data is a position record from an acoustically tracked, neutrally buoyant subsurface oceanographic float, one of a number deployed in the eastern subtropical North Atlantic Ocean in order to track the motion of dense salty water flowing out from the Mediterranean Sea during the "Eastern Basin" experiment [5] .
• The data is available online from the World Ocean Circulation Experiment Subsurface Float Data Assembly Center at http://wfdac.whoi.edu.
• Looping trajectories are indicative of intense swirling currents around an isolated packet of Mediterranean Sea water.
• Applied to such signals which a priori contain meaningful rotating components, the output of the bivariate extensions typically provide the given decomposition, where the rotations that were apparent in the original signal have been isolated in separate components.
• Such advanced study of the rotating components has already been performed using wavelet ridges to extract the coherent vortex signal [8] .

Figures

Fig. 2. Definition of the mean of the envelope for each algorithm. The accuracies of the estimation of the sampling points and of their mean are represented schematically as the thick lines of variable length under each point. (the provided script allows to easily test other configurations)

Fig. 3. A signal and its Bivariate Empirical Mode Decomposition. The real parts are plotted as solid blue lines and the imaginary parts as dashed black lines. The decomposition has been obtained with the second algorithm using 64 directions and 10 iterations to extract each component. Portions where components are rotating can be identified by a constant phase shift between the real and imaginary parts: if these are in quadrature, then the signal is rotating circularly; other phase shifts correspond to elliptic rotations.

Fig. 1. The principle of the bivariate extensions. (a) A composite rotating signal. (b) The signal enclosed in its 3D envelope. The black thick lines stand for the envelope curves that are used to derive the mean. (c) Rapidly rotating component. (d) More slowly rotating component corresponding to the mean of the tube in (b).

Full text

HAL Id: ensl-00137611
https://hal-ens-lyon.archives-ouvertes.fr/ensl-00137611
Preprint submitted on 20 Mar 2007
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Bivariate Empirical Mode Decomposition
Gabriel Rilling, Patrick Flandrin, Paulo Gonçalves, Jonathan M. Lilly
To cite this version:
Gabriel Rilling, Patrick Flandrin, Paulo Gonçalves, Jonathan M. Lilly. Bivariate Empirical Mode
Decomposition. 2007. �ensl-00137611�

IEEE SIGNAL PROCESSING LETTERS 1
Bivariate Empirical Mode Decomposition
Gabriel Rilling, Patrick Flandrin, Fellow, IEEE, Paulo Gonçalves, Jonathan M. Lilly
Abstract
The Empirical Mode Decomposition (EMD) has been introduced quite recently to adaptively decom-
pose nonstationary and/or nonlinear time series [1]. The method being initially limited to real-valued time
series, we propose here an extension to bivariate (or complex-valued) time series which generalizes the
rationale underlying the EMD to the bivariate framework. Where the EMD extracts zero-mean oscillating
components, the proposed bivariate extension is designed to extract zero-mean rotating components. The
method is illustrated on a real-world signal and properties of the output components are discussed. Free
Matlab/C codes are available at http://perso.ens-lyon.fr/patrick.flandrin.
Index Terms
Empirical Mode Decomposition, complex-valued signals, bivariate time series
EDICS Category: DSP-TFSR
I. INTRODUCTION
In its original formulation [1], the Empirical Mode Decomposition (EMD) can only be applied to
real-valued time series. The purpose of this paper is to introduce a new extension of the EMD destined
to handle bivariate (or complex-valued) time series. Note however that not all bivariate time series can be
processed by this new method but only those where the two components can be assimilated to Cartesian
Manuscript submitted March 20, 2007. G. Rilling and P. Flandrin are with the Physics Department (UMR 5672 CNRS), Ecole
Normale Supérieure de Lyon, 46 allée d’Italie, 69364 Lyon Cedex 07 France. Phone: +33 (0)4 72 72 81 60; Fax: +33 (0)4 72
72 80 80; E-mail: {grilling,flandrin}@ens-lyon.fr.
P. Gonçalves is with the Laboratoire de l’Informatique du Parallélisme (UMR CNRS/INRIA 5668), Ecole Normale Supérieure
de Lyon, 46 allée d’Italie, 69364 Lyon Cedex 07 France. Phone: +33 (0)4 72 72 83 89; E-mail: paulo.goncalves@inria.fr.
J. M. Lilly is with the Earth and Space Research, 1910 Fairview Ave E, Suite 210, Seattle WA 98102 USA. E-mail:
Jonathanlilly@gmail.com

IEEE SIGNAL PROCESSING LETTERS 2
coordinates of a point moving in a 2-dimensional space. In particular, the meaning of the signal should not
depend on the choice of such Cartesian coodinates. It is worth noticing that another bivariate extension
has been introduced very recently [2]. The difference with the one we propose here is significant, but in
a nutshell that other method cleverly uses the original EMD to decompose bivariate time series, while
ours is a new algorithm that adapts the rationale underlying the EMD to the bivariate framework. Further
comparison of the approaches however is out of the scope of this paper. The communication is organized
as follows. The bivariate extension is introduced in Section II. Section III is about the components of
the resulting decomposition and an illustration is proposed in Section IV. Additionally, free Matlab/C
codes corresponding to the proposed algorithms are made available at http://perso.ens-lyon.
fr/patrick.flandrin along with small scripts aimed at reproducing the figures and other EMD-
related software.
II. FROM UNIVARIATE EMD TO BIVARIATE EMD
A. Classical EMD
Basically, the EMD considers a signal at the scale of its local oscillations. The main idea of EMD is
then to formalize the idea that, locally: “signal = fast oscillations superimposed on slow oscillations”.
Looking at a single oscillation (defined, e.g., as the signal between two consecutive local minima), the
EMD is designed to define a local “low frequency” component as the local trend m
1
[x](t), supporting
a local “high frequency” component as a zero-mean oscillation or local detail d
1
[x](t), so that we can
express x(t) as
x(t) = m
1
[x](t) + d
1
[x](t). (1)
By construction, d
1
[x](t) is an oscillatory signal and, if it is furthermore required to be locally zero-mean
everywhere, it corresponds to what is referred to as an Intrinsic Mode Function (IMF) [1]. Practically,
this primarily implies that all its maxima are positive and all its minima are negative. On the other hand,
all we know about m
1
[x](t) is that it locally oscillates more slowly than d
1
[x](t). We can then apply the
same decomposition to it, leading to m
1
[x](t) = m
2
[x](t) + d
2
[x](t) and, recursively applying this on
the m
k
[x](t), we get a representation of x(t) of the form
x(t) = m
K
[x](t) +
K
X
k=1
d
k
[x](t). (2)
The discrimination between “fast” and “slow” oscillations is obtained through an algorithm referred
to as the sifting process [1] which iterates a nonlinear elementary operator S on the signal until some

IEEE SIGNAL PROCESSING LETTERS 3
Fig. 1. The principle of the bivariate extensions. (a) A composite rotating signal. (b) The signal enclosed in its 3D envelope.
The black thick lines stand for the envelope curves that are used to derive the mean. (c) Rapidly rotating component. (d) More
slowly rotating component corresponding to the mean of the tube in (b).
stopping criterion is met. Given a signal x(t), the operator S is defined by the following procedure:
Identify all extrema of x(t)1
Interpolate (using a cubic spline) between minima (resp.2
maxima), ending up with some “envelope” e
min
(t) (resp. e
max
(t))
Compute the mean m(t) =
(e
min
(t)+e
max
(t))
2
3
Subtract [to] from the signal to obtain S[x](t) = x(t) − m(t)4
If the convergence criterion is met after n iterations, the local detail and the local trend are defined
as d
1
[x](t) = S
n
[x](t) and m
1
[x](t) = x(t) − d
1
[x](t).
B. Envelopes in 3 dimensions
The EMD is based on the intuitive notion of “oscillation” which naturally relates to local extrema. But
the notion of oscillation is much more confusing when the analyzed data is intrinsically bivariate and it
is unclear how to define and interpret local extrema. What is rather clear on the other hand is the notion
of rotation, which moreover is arguably a two-dimensional extension of the usual notion of a univariate
oscillation. Therefore, the basic idea underlying the proposed bivariate EMD is to formalize the following
idea: “bivariate signal = fast rotations superimposed on slower rotations”. As with the classical EMD,

IEEE SIGNAL PROCESSING LETTERS 4
it is clear that the adopted viewpoint is a priori rather restrictive as, e.g., a white noise signal is not
meaningfully treated as a sum of oscillations (or rotations). Still this does not prevent the algorithm from
producing a decomposition for any signal, as with the univariate EMD.
In order to separate the more rapidly rotating component from slower ones, the idea is once again
to define the slowly rotating component as the mean of some “envelope”. Yet the envelope is now a
3-dimensional tube that tightly encloses the signal (see
Fig. 1 (b)). Given this, the slowly rotating portion
of the signal at any point in time can then be defined as the center of the enclosing tube. To this end,
only a given number of points on the tube’s periphery are considered, each one being associated with a
specific direction. If only 4 points are used, these can be the extreme points in the directions top, bottom,
left and right (see Fig. 2). In practice, the top point, for example, is uniquely defined only when the signal
reaches a local maximum in the vertical direction and is therefore tangent to the top of the tube. Between
such characteristic moments in time, the top point is then simply defined using interpolation, ending up
with the “envelope” associated with the upwards direction (cf the black thick lines in Fig. 1 (b)). Now,
given some set of points on the tube periphery at a given instant in time, there are at least two ways to
define their mean:
1) define the mean as the barycenter of the 4 points, considering each to have unit mass (see Fig. 2 (a)).
2) define the mean as the intersection of two straight lines, one being halfway between the two
horizontal tangents, the other one halfway between the vertical ones (see Fig. 2 (b)).
In practice, however, the second scheme may be preferred because it is naturally more robust to sampling
errors. More precisely, the reason for this is that the envelope points are defined up to an uncertainty that
is not isotropic. Indeed, the order of magnitude of this uncertainty can be estimated through a Taylor
expansion, which results in an uncertainty jointly proportional to dx/dt and to the sampling period. Thus,
the uncertainty is greatest in the direction locally tangent to the signal, and much smaller (of second
order) in the orthogonal direction. As the second scheme only uses information from the orthogonal
direction, it is naturally more accurate, especially when the signal is sampled sparsely with respect to its
period. Note that sampling effects shall not be taken too lightly as the original EMD has been shown to
be very sensitive to sampling[3].
The desired goal concerning the interpolation is the same as for the classical EMD: a smooth interpo-
lation with as few “spurious bumps” as possible. Among common interpolation schemes, this calls for
cubic spline as it is well known for its minimum curvature property and, in practice, it is still considered
the best interpolation scheme for the EMD[4].
In the preceding discussion, we have limited ourselves to 4 directions for the sake of simplicity,

Frequently asked questions

Q1. What have the authors contributed in "Bivariate empirical mode decomposition" ?

The Empirical Mode Decomposition ( EMD ) has been introduced quite recently to adaptively decompose nonstationary and/or nonlinear time series [ 1 ]. The method being initially limited to real-valued time series, the authors propose here an extension to bivariate ( or complex-valued ) time series which generalizes the rationale underlying the EMD to the bivariate framework. The method is illustrated on a real-world signal and properties of the output components are discussed. 

Q2. What is the envelope curve associated with the direction k = 2k/N?

The envelope curve associated with the direction ϕk = 2kπ/N is then equal to the maximum signal value in that direction, where the phase of the signal’s derivative is ψ(t) = 

Q3. What is the main idea of the EMD?

Looking at a single oscillation (defined, e.g., as the signal between two consecutive local minima), the EMD is designed to define a local “low frequency” component as the local trend m1[x](t), supporting a local “high frequency” component as a zero-mean oscillation or local detail d1[x](t), so that the authors can express x(t) asx(t) = m1[x](t) + d1[x](t). 

Q4. What is the main idea of EMD?

The main idea of EMD is then to formalize the idea that, locally: “signal = fast oscillations superimposed on slow oscillations”. 

Q5. What is the advantage of a large number of directions?

Moreover a large number of directions may be interesting insofar as it reduces the dependance of the final decomposition with respect to rotations of the spatial coordinates. 

Q6. What is the meaning of the mean of the two straight lines?

2) define the mean as the intersection of two straight lines, one being halfway between the twohorizontal tangents, the other one halfway between the vertical ones (see Fig. 2 (b)). 

Q7. What is the EMD for a linear time series?

2.Algorithm 1: EMD bivariate extension: scheme 1for 1 ≤ k ≤ N do1 Project the complex-valued signal x(t) on direction ϕk:2pϕk(t) = Re ( e−iϕkx(t) ) 

Q8. What is the basic idea underlying the proposed bivariate EMD?

the basic idea underlying the proposed bivariate EMD is to formalize the following idea: “bivariate signal = fast rotations superimposed on slower rotations”.