# 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. IEEE Signal Processing Letters, 2007, 14 (12), pp.936-939. 10.1109/LSP.2007.904710. inria-00570563v2

HAL Id: inria-00570563

https://inria.hal.science/inria-00570563v2

Submitted on 6 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific 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 diffusion de documents
scientifiques 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bivariate Empirical Mode Decomposition

Gabriel Rilling, Patrick Flandrin, Fellow, IEEE, Paulo Gonçalves, and Jonathan M. Lilly

**Abstract**—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, 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.

**Index Terms**—Bivariate time series, complex-valued signals, empirical mode decomposition.

## 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 letter 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 coordinates of a point moving in a two-dimensional space. In particular, the meaning of the signal should not depend on the choice of such Cartesian coordinates. It is worth noticing that two other bivariate extensions have been introduced very recently [2], [3]. The first one is very different from the one we propose here since it cleverly uses the original EMD to decompose bivariate time series, based on properties of the complex field, while the second one and ours directly adapt the rationale underlying the EMD to the bivariate framework and might be extensible to more dimensions. As will be seen later, the algorithm in [3] turns out to be a simplified version of one of the algorithms proposed in this letter. 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(t)$, supporting a local “high frequency” component as a zero-mean oscillation or local detail $d_1(t)$, so that we can express $x(t)$ as

$$x(t) = m_1(t) + d_1(t). \quad (1)$$

By construction, $d_1(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(t)$ is that it locally oscillates more slowly than $d_1(t)$. We can then apply the same decomposition to it, leading to $m_1(t) = m_2(t) + d_2(t)$, and, recursively applying this on the $m_k(t)$, we get a representation of $x(t)$ of the form

$$x(t) = m_K(t) + \sum_{k=1}^{K} d_k(t). \quad (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 stopping criterion is met. Given a signal $x(t)$, the operator $S$ is defined by the following procedure.

1) Identify all extrema of $x(t)$.

2) Interpolate (using a cubic spline) between minima (resp. maxima), ending up with some “envelope” $e_{min}(t)$ (resp. $e_{max}(t)$).

3) Compute the mean $m(t) = (e_{min}(t) + e_{max}(t))/2$.

4) Subtract from the signal to obtain $S[x](t) = x(t) - m(t)$.

If the stopping criterion is met after $n$ iterations of the sifting process, the local detail and the local trend are defined as $d_1(t) = S^n[x](t)$ and $m_1(t) = x(t) - d_1(t)$.

### B. Envelopes in Three Dimensions

The EMD is based on the intuitive notion of “oscillation” that naturally relates to local extrema. However, the notion of oscillation is much more confusing when the analyzed data are intrinsically bivariate, and it is unclear how to define and interpret

<footnote>G. Rilling and P. Flandrin are with the Physics Department (UMR 5672 CNRS), Ecole Normale Supérieure de Lyon, 69364 Lyon Cedex 07, France (e-mail: grilling@ens-lyon.fr; flandrin@ens-lyon.fr).</footnote>

<footnote>P. Gonçalves is with the Laboratoire de l'Informatique du Parallélisme (UMR CNRS/INRIA 5668), Ecole Normale Supérieure de Lyon, 69364 Lyon Cedex 07, France (e-mail: paulo.goncalves@inria.fr).</footnote>

<footnote>J. M. Lilly is with Earth and Space Research, Seattle, WA 98102 USA (e-mail: Jonathanlilly@gmail.com).</footnote>
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, 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 three-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 four 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) barycenter of the four points, considering each to have unit mass [see Fig. 2(a)];

2) intersection of two straight lines, one being halfway between the two horizontal tangents and 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 overlooked since the original EMD has been shown to be very sensitive to sampling [4].

The desired goal concerning the interpolation is the same as for the classical EMD: a smooth interpolation 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 [5].

In the preceding discussion, we have limited ourselves to four directions for the sake of simplicity, but there is of course no such limitation in practice. Moreover, a large number of directions may be interesting insofar as it reduces the dependence of

<figure><img src="image_2.png" /><figcaption>Fig. 1. Principle of the bivariate extensions. (a) Composite rotating signal. (b) Signal enclosed in its 3-D 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).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>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).</figcaption></figure>

<caption>**Algorithm 1:** EMD bivariate extension: scheme 1</caption>

1 **for** 1 ≤ k ≤ N **do**
2 &nbsp;&nbsp;&nbsp;&nbsp;Project the complex-valued signal $x(t)$ on direction $\varphi_k$:
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p_{\varphi_k}(t) = \text{Re}(e^{-i\varphi_k}x(t))$
3 &nbsp;&nbsp;&nbsp;&nbsp;Extract the locations $\{t_j^k\}$ of the maxima of $p_{\varphi_k}(t)$.
4 &nbsp;&nbsp;&nbsp;&nbsp;Interpolate the set $\{(t_j^k, x(t_j^k))\}$ to obtain the envelope curve in direction $\varphi_k: e_{\varphi_k}(t)$.
5 **Compute the mean of all envelope curves:**
&nbsp;&nbsp;&nbsp;&nbsp;$m(t) = \frac{1}{N} \sum_k e_{\varphi_k}(t)$
6 **Subtract the mean to obtain** $S^{B1}[x](t) = x(t) - m(t)$

the final decomposition with respect to rotations of the spatial coordinates. For convenience of the presentation, the bivariate time series in the following will be treated as complex-valued time series. Given a set of directions $\varphi_k = 2k\pi/N$, $1 \le k \le N$, the proposed bivariate extensions are defined by the same algorithm as the basic EMD, only with new sifting elementary operators $S^{B1}$ and $S^{B2}$ corresponding, respectively, to the algorithms Algo. 1 and Algo. 2. The algorithm proposed in [3] is the same as Algo. 1 with $N = 2$ directions.
<caption>**Algorithm 2:** EMD bivariate extension: scheme 2</caption>

1 **for** 1 ≤ k ≤ N **do**
2 &nbsp;&nbsp;&nbsp;&nbsp;Project the complex-valued signal $x(t)$ on direction $\varphi_k$:
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p_{\varphi_k}(t) = \operatorname{Re}(e^{-i\varphi_k}x(t))$
3 &nbsp;&nbsp;&nbsp;&nbsp;Extract the maxima of $p_{\varphi_k}(t)$: $\{t_j^k, p_j^k\}$.
4 &nbsp;&nbsp;&nbsp;&nbsp;Interpolate the set $\{(t_j^k, e^{i\varphi_k}p_j^k)\}$ to obtain the
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;"tangent" in direction $\varphi_k$: $e'_{\varphi_k}(t)$.
5 Compute the mean of all tangents: $m(t) = \frac{2}{N} \sum_k e'_{\varphi_k}(t)$
6 Subtract the mean to obtain $S^{B2}[x](t) = x(t) - m(t)$

<caption>**Algorithm 3:** Reformulation of scheme 2</caption>

1 **for** 1 ≤ k ≤ N/2 **do**
2 &nbsp;&nbsp;&nbsp;&nbsp;Project the complex-valued signal $x(t)$ on direction $\varphi_k$:
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p_{\varphi_k}(t) = \operatorname{Re}(e^{-i\varphi_k}x(t))$
3 &nbsp;&nbsp;&nbsp;&nbsp;Compute the partial estimate in direction $\varphi_k$:
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$s_{\varphi_k}(t) = e^{i\varphi_k}S[p_{\varphi_k}](t)$
4 Compute the final estimate: $S^{B2}[x](t) = \frac{4}{N} \sum_{k=1}^{N/2} s_{\varphi_k}(t)$

Furthermore, the second algorithm can be greatly simplified if we notice that the set that is interpolated at step 4 is in fact included in the plane containing the time axis and direction $\varphi_k$. Thus, the interpolation that is performed is very similar to that in the original EMD sifting. Hence, if the number of considered directions is an even number, the second algorithm can be expressed in terms of the sifting operator for the univariate EMD (see Algo. 3). From a theoretical viewpoint, this is interesting as it allows to study the behavior of the algorithm in light of what is already known about the classical EMD.

### III. BIVARIATE INTRINSIC MODE FUNCTIONS

The proposed EMD bivariate extensions have been designed so that signals rotating around zero are admissible outputs. As this is a rather vague notion, the purpose of this section is to clarify what signals the algorithms actually consider admissible outputs, that is, we ask what sort of signal $x(t)$ is nearly a fixed point of the sifting operator

$$S^{B1}[x](t) \approx x(t) \quad \text{or} \quad S^{B2}[x](t) \approx x(t). \quad (3)$$

For the sake of simplicity, we will only address the case of simple periodic solutions, which are convenient to describe while still being rather general. 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. Moreover, signals will here also be considered infinitely continuously differentiable, since the purpose of this section is to provide insight into some possible solutions rather than to make an exhaustive study.

What we find from various simulations is that both algorithms generally accept two types of solutions, the first one corresponding to rotating signals as intended and the second one to cases where the method actually fails to extract rotating components and therefore outputs signals that wander around zero in

a more complicated way. 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 complex-valued white Gaussian noise signal. Moreover, it seems that these non-rotating solutions almost never occur when the analyzed signal is analytic (or anti-analytic). On the other hand, solutions from the first type are signals in which the local sense of rotation never changes. The sense of rotation can be defined, e.g., from the sense of the vector product of the velocity and acceleration vectors

$$\operatorname{sign} \left\{ \operatorname{Im} \left\{ \frac{dx}{dt} \cdot \left( \frac{d^2 x}{dt^2} \right)^* \right\} \right\}. \quad (4)$$

Alternatively, if the signal rotates in a counterclockwise fashion (and $|dx/dt| > 0, \forall t$), this exactly means that the signal's derivative can be expressed as

$$\frac{dx}{dt}(t) = r(t)e^{i\psi(t)}, \quad r(t) > 0, \quad \frac{d\psi(t)}{dt} > 0. \quad (5)$$

Finally, periodic signals satisfying (5) are not necessarily fixed points of the sifting operators (3). One of the reasons for this is that there is no condition specifying that they are rotating around zero yet. In order to clarify what is meant by “around zero,” we can consider the very simple case in which the signal performs only one full rotation around zero per period ($\psi(t+T) = \psi(t)+2\pi$, where T is the period). This property simplifies greatly the study, as one can easily show that it implies there is only one maximum per period at step 2 in the algorithms. Therefore, all envelope curves are in fact constants with respect to time, which allows one to derive the mean analytically. The envelope curve associated with the direction $\varphi_k = 2k\pi/N$ is then equal to the maximum signal value in that direction, where the phase of the signal's derivative is $\psi(t) = 2k\pi/N + \pi/2$. Thus, the mean for the first algorithm reads

$$m^{B1}(t) = \frac{1}{N} \sum_k x \left( \psi^{-1} \left( \frac{2k\pi}{N} + \frac{\pi}{2} \right) \right). \quad (6)$$

In the limit where the number of directions tends to infinity, this results in

$$\begin{aligned} m^{B1}(t) &= \frac{1}{2\pi} \int_{0}^{2\pi} x(\psi^{-1}(\varphi)) d\varphi \\ &= \frac{1}{2\pi} \int_{0}^{T} x(t) \frac{d\psi}{dt} dt \end{aligned} \quad (7)$$

which is simply the mean of the signal over a period weighted by $d\psi(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

$$\begin{aligned} m^{B2}(t) &= \frac{1}{\pi} \int_{0}^{T} e^{i\psi(t)} - \frac{i\pi}{2} \operatorname{Re} \left( e^{-i\psi(t)} + \frac{i\pi}{2} x(t) \right) \frac{d\psi}{dt} dt \\ &= m^{B1}(t) - \frac{1}{2\pi} \int_{0}^{T} e^{2i\psi(t)} x^*(t) \frac{d\psi}{dt} dt \end{aligned} \quad (8)$$
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

and hence, $m^{B2}(t) = m^{B1}(t)$ since

$$ 
\begin{aligned}
& \int_{0}^{T} e^{2i\psi(t)} x^{*}(t) \frac{d\psi}{dt} dt \\
&= \frac{i}{2} \left( \left[ -e^{2i\psi(t)} x^{*}(t) \right]_{0}^{T} + \int_{0}^{T} \frac{dx^{*}}{dt} e^{2i\psi(t)} dt \right) \\
&= \frac{i}{2} \int_{0}^{T} r(t) e^{i\psi(t)} dt \\
&= \frac{i}{2} \int_{0}^{T} \frac{dx}{dt} dt = 0.
\end{aligned} 
\quad (9) $$

Thus, in this very simple case, the mean is in fact the same for both algorithms, and therefore, such a simple signal is a fixed point of both sifting operators iff the integral (7) is close to zero.

More generally, the outputs of the two proposed algorithms are very similar when the data clearly contain rotating components, but they may differ significantly when they fail to extract rotating components. Notice, though, that rarely one method succeeds in retrieving rotating components when the other one fails.

## IV. ILLUSTRATION

A typical application of the proposed algorithms is proposed in Fig. 3. The data are 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 [6]. The data are available online from the World Ocean Circulation Experiment Subsurface Float Data Assembly Center (WFDAC) at http://wfdac.whoi.edu. Looping trajectories are indicative of intense swirling currents around an isolated packet of Mediterranean Sea water. Such structures, called “coherent vortices,” are frequently observed in the ocean [7] and are more generally a ubiquitous feature of rotating turbulent fluids [8]. Applied to such signals that 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. As mentioned earlier, though, not all components are rotating but primarily the first two ones, which correspond to the presumed coherent vortex. In this example, however, we do not precisely know what information can be extracted from the decomposition yet, but we expect the large-scale non-rotating components to reveal useful information regarding the background fluctuations determining the vortex position while the rotating components can a priori be used to extract finer informations, such as amplitude, angular frequency, etc. Such advanced study of the rotating components has already been performed using wavelet ridges to extract the coherent vortex signal [9]. A comparative study of the EMD bivariate extensions for similar purposes is currently under investigation.

## REFERENCES

[1] N. E. Huang, Z. Shen, S. R. Long, M. L. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proc. R. Soc. London A., vol. 454, pp. 903–995, 1998.

[2] T. Tanaka and D. P. Mandic, “Complex empirical mode decomposition,” IEEE Signal Process. Lett., vol. 14, no. 2, pp. 101–104, Feb. 2007.

[3] M. U. Altaf, T. Gautama, T. Tanaka, and D. P. Mandic, “Rotation invariant empirical mode decomposition,” in IEEE Int. Conf. Acoustics, Speech, Signal Processing, Honolulu, HI, 2007.

[4] G. Rilling and P. Flandrin, “On the influence of sampling on the empirical mode decomposition,” in IEEE Int. Conf. Acoustics, Speech, Signal Processing, Toulouse, France, 2006.

[5] N. E. Huang and S. S. P. Shen, Eds., Hilbert-Huang Transform and Its Applications. Singapore: World Scientific, 2005.

[6] P. Richardson, D. Walsh, L. Armi, M. Schroder, and J. F. Price, “Tracking three Meddies with SOFAR floats,” J. Phys. Oceanogr., vol. 19, pp. 371–383, 1989.

[7] J. C. McWilliams, “Submesoscale coherent vortices in the ocean,” Rev. Geophys., vol. 23, no. 2, pp. 165–182, 1985.

[8] J. C. McWilliams, “The vortices of geostrophic turbulence,” J. Fluid Mech., vol. 219, pp. 387–404, 1990.

[9] J. M. Lilly and J.-C. Gascard, “Wavelet ridge diagnosis of time-varying elliptical signals with application to an oceanic eddy,” Nonlin. Processes Geophys., vol. 13, pp. 467–483, 2006.