Sur la décomposition modale empirique

Patrick F LANDRIN1 , Paulo G ONÇALV ÈS2

1 CNRS — Laboratoire de Physique (UMR 5672)

École Normale Supérieure de Lyon
46 allée d’Italie, 69364 Lyon Cedex 07
2 INRIA
Rhône-Alpes — Projet IS2
655 avenue de l’Europe, 38334 St Ismier Cedex
flandrin@[Link], [Link]@[Link]

Résumé – On cherche ici à mieux comprendre une technique récente de décomposition locale et auto-adaptative de signaux appelée
Décomposition Modale Empirique (ou EMD, pour “Empirical Mode Decomposition”). On se base sur des simulations numériques extensives
permettant de mettre en évidence des comportements typiques de la méthode et d’en quantifier les performances dans des cas simples et contrôlés.
On donne une interprétation spectrale de la méthode et on souligne ses ressemblances et différences avec des techniques de décomposition mul-
tirésolution à base d’ondelettes.

Abstract – This paper is intended to give new insights into a recent technique of local and data-dependent signal decomposition, referred to as
“Empirical Mode Decomposition” (EMD). The study relies heavily on numerical simulations which allow for evidencing typical characteristics
of the method and for quantifying performance in simple well-controlled situations. A spectral interpretation of EMD is given, and emphasis is
put on its similarities and differences with wavelet-based multiresolution decomposition techniques.

1 Introduction pour t− ≤ t ≤ t+ . Si ce point de vue est adopté pour l’ensem-

ble des oscillations constituant le signal, la procédure peut alors
N.E. Huang et ses collaborateurs ont récemment introduit être appliquée sur le résidu m1 [x](t) formé par l’ensemble des
une technique de décomposition de signaux non stationnaires, tendances locales et considéré comme un nouveau signal, con-
appelée Décomposition Modale Empirique (ou EMD, pour “Em- duisant à un nouveau détail d2 [x](t) et à un nouveau résidu
pirical Mode Decomposition”) [1]. Si celle-ci possède un car- m2 [x](t).
actère très intuitif, elle est restée jusqu’à présent rêtive à une On obtient ainsi une décomposition dont les différents modes
véritable évaluation théorique dans la mesure où sa définition (ou IMF, pour “Intrinsic Mode Functions”) dk [x](t) sont ex-
n’est en fait donnée que par l’algorithme de décomposition pro- traits itérativement, conduisant à une représentation du type :
posé.
L’EMD a néanmoins été appliquée à un certain nombre de

K
x(t) = mK [x](t) + dk [x](t)
situations réelles (cf., par exemple, [1, 2, 3]), motivant des études
k=1
systématiques de la méthode dans des conditions contrôlées
afin de mieux cerner ses possibilités et limitations. pour une profondeur de décomposition K. On pourra trouver
dans [1, 4] davantage de détails sur la méthode, et dans [5] des
codes M ATLAB pour la mettre en œuvre.
2 Principe de l’EMD
Dans son principe, l’EMD considère les signaux à l’échelle 3 L’EMD comme filtrage
de leurs oscillations locales, sans que celles-ci soient néces-
sairement harmoniques au sens de Fourier. D’une manière plus Pour la simplicité de la présentation, les différents modes et
précise, si l’on cherche à décrire un signal x(t) entre deux ex- résidus ont été introduits sur des arguments “spectraux”, mais il
trema consécutifs (par exemple, deux minima situés aux temps est bon de noter que leur discrimination “hautes fréquences” vs.
t− et t+ ), on peut définir de façon heuristique une contribu- “basses fréquences” est essentiellement qualitative (les oscilla-
tion “hautes fréquences” locale {d1 [x](t), t− ≤ t ≤ t+ }, ou tions ne sont pas nécessairement sinusoı̈dales), qu’elle ne s’ap-
détail local, qui correspond à l’oscillation se terminant aux plique que localement et qu’elle ne correspond en aucune façon
deux minima considérés et passant par le maximum qui ex- à un filtrage prédéterminé. La sélection des modes s’apparente
iste nécessairement entre eux. Pour que la description du com- davantage à un filtrage non linéaire dépendant du temps, au-
portement local soit complète, il suffit d’identifier la contribu- tomatique et auto-adaptatif (c’est-à-dire localement dépendant
tion “basses fréquences” locale correspondante m1 [x](t), ou du signal). Néanmoins, un des objectifs de cette communica-
tendance locale, de telle sorte que l’on ait tion est de justifier que, dans des cas idéalisés et bien contrôlés,
l’EMD réalise de manière automatique une décomposition en
x(t) = m1 [x](t) + d1 [x](t) sous-bandes très proche de ce que donnerait une analyse mul-
 H = 0.2 H = 0.1 H = 0.2 H = 0.3
10
0 IMF # : 7 6 5 4 3 2 1

-10
dB

-20
-30
-40
-9 -8 -7 -6 -5 -4 -3 -2 -1
H = 0.5
10
H = 0.4 H = 0.5 H = 0.6
0 IMF # : 7 6 5 4 3 2 1

-10
dB

-20
-30
-40
-9 -8 -7 -6 -5 -4 -3 -2 -1
H = 0.8
10
0 IMF # : 7 6 5 4 3 2 1
H = 0.7 H = 0.8 H = 0.9
-10
dB

-20
-30
-40
-9 -8 -7 -6 -5 -4 -3 -2 -1
log (fréquence)
2

F IG . 1: Filtres équivalents — Dans le cas du bruit gaussien

fractionnaire, l’EMD peut s’interpréter comme un banc de fil- F IG . 2: Variance-covariance — Dans le cas du bruit gaussien
tres passe-bande pour les modes d’indices k ≥ 2, le mode fractionnaire, la variance des IMF suit une loi de puissance liée
#1 correspondant essentiellement à un filtre passe-haut. Pour simplement à l’exposant de Hurst, les IMF étant de plus ap-
chaque valeur de l’exposant de Hurst (H = 0.2, 0.5 et 0.8), on proximativement décorrélés entre eux. Pour chaque valeur de
a simulé 5000 séries indépendantes de 512 points chacune, et l’exposant de Hurst (H = 0.1 à 0.9), on a simulé 100 séries
on a tracé dans un diagramme log-log le spectre moyen des 7 indépendantes de 2048 points chacune, et on a tracé la matrice
premières IMF en fonction de la fréquence normalisée. de variance-covariance moyenne relative aux 7 premières IMF
(codage logarithmique en niveaux de gris, dynamique de 35
tirésolution de type ondelettes. Pour ce faire, on s’appuie ici sur dB).
la conduite de deux types d’expériences numériques.
Banc de filtres — Le résultat central de ces simulations est
3.1 Approche spectrale donné à la Figure 1 qui représente, mode par mode, les estimées
des densités spectrales
Dans le premier cas, on aborde la question de la décom-  +∞
position d’un bruit large bande en s’intéressant au comporte- Γk [gH ](f ) = E {dk [gH ](t)dk [gH ](t + τ )} e−i2πf τ dτ,
ment statistique de l’EMD lorsqu’on l’applique au modèle très −∞
général du bruit gaussien fractionnaire [8]. mettant en évidence une organisation spontanée de type banc
Rappelons que le bruit gaussien fractionnaire (fGn) est défini de filtres à surtension constante (des résultats de même nature
comme le processus des accroissements du mouvement brown- ont depuis été obtenus indépendamment [9], quoique dans le
ien fractionnaire (fBm) [6]. À temps discret, le fGn est une série cas du seul bruit blanc).
{gH [n], n = . . . , −1, 0, 1, . . . }, indexée par un paramètre réel
0 < H < 1 (son exposant de Hurst), et telle que sa séquence Variance-covariance — La nature quasi-dyadique de ces bancs
d’autocorrélation de filtres a été justifiée dans [8], où l’on a de plus montré com-
rgH [k] := E{gH [n]gH [n + k]} ment la puissance spectrale se répartissait exponentiellement
dans les différents IMF, de telle sorte que cette répartition per-
s’écrit: met en théorie une estimation de H. Afin de préciser plus avant
σ2   la pertinence de cette approche, en particulier lorsqu’on ne dis-
rgH [k] = |k − 1|2H − 2|k|2H + |k + 1|2H .
2 pose que d’une réalisation, nous avons complété l’étude pré-
Le cas particulier H = 1/2 se réduit à celui du bruit blanc sentée en [8] en évaluant le lien statistique pouvant exister en-
(donc, décorrélé), alors que les autres valeurs de H introduisent tre modes. À cette fin nous avons estimé la matrice de variance-
des corrélations non nulles pouvant être, soit négatives quand covariance
0 < H < 1/2, soit positives quand 1/2 < H < 1 (longue   
1 T
mémoire). CH [m, n] = E dm [gH ](t)dn [gH ](t) dt
Partant de ce modèle, nous avons conduit des simulations T 0
numériques extensives pour des valeurs de H allant de 0.1 à sur la base de 100 réalisations indépendantes, pour une taille
0.9. La longueur typique des séries étudiées était de N = 512 de données de N = 2048 points et m et n variant de 1 à 7.
points et, pour chaque valeur de H, 5000 réalisations indé- La Figure 2 présente le résultat normalisé CH [m, n]/CH [1, 1].
pendantes de fGn ont été générées par l’algorithme de Wood Comme attendu, on y observe sur la diagonale principale (soit,
et Chan [7] (on ne donnera pas davantage de détails ici sur m = n) le comportement
la mise en œuvre effective de l’EMD, mais les codes Matlab
2(H−1)m
utilisés pour les simulations sont disponibles [5]). CH [m, m] = var dm [gH ](t) ∝ ρH ,
 H = 0.2 H = 0.5 H = 0.8

signal
imf1
pk(x)

p (x)

pk(x)
k

imf2
amplitude amplitude amplitude

imf3
H = 0.2 H = 0.5 H = 0.8
10 10 10

imf4
0 0 0
log pk(x) (dB)

log p (x) (dB)

log pk(x) (dB)

-10 -10 -10

imf5
-20 -20 -20
k

-30 -30 -30

-40 -40 -40

imf6
-50 -50 -50
amplitude amplitude amplitude
temps

F IG . 3: Statistiques marginales des IMF — Hormis pour l’IMF F IG . 5: Réponse impulsionnelle — La réponse impulsionnelle
1 qui est bimodale, les densités de probabilité marginales équivalente de l’EMD peut être obtenue en moyennant mode
{pk (x), k = 2, . . . 5} des IMF du fGn (calculées ici sur 500 par mode un grand nombre de décompositions calculées sur
réalisations de N = 512 points) sont approximativement une impulsion faiblement bruitée. Dans le cas présent, on a
gaussiennes et réductibles à un gabarit unique moyennant une simulé 5000 réalisations indépendantes de 256 points chacune
renormalisation bien choisie (diagrammes du bas : superposi- et on a tracé (sous le signal moyen) la moyenne d’ensemble,
tion des log pk (x) renormalisées). normalisée en amplitude mode par mode, des 6 premières IMF.

H = 0.5
L’IMF 1 s’écarte du cas gaussien en possédant une densité
de probabilité bimodale, cohérente avec sa nature de bruit filtré
passe-haut (cf. Figure 1). En effet, dans le cas (H = 0.5) du
bruit blanc gaussien à temps discret, un calcul simple montre
p1(x)

que les densités de probabilité p+ (x) et p− (x) des maxima et

des minima s’écrivent :
  2
1 1
amplitude p± (x) = ± − erf(x) γ(x),
2 2
avec
F IG . 4: Densité de probabilité de l’IMF 1 — Dans le cas H =
1
γ(x) = √ e−x /2
2
0.5, la densité de probabilité p1 (x) de l’IMF 1 (calculée ici sur
500 réalisations de N = 512 points et représentée par les croix) 2π
est très proche de celle prédite pour les maxima positifs et les et
minima négatifs d’un bruit blanc gaussien à temps discret (trait  +∞
continu). erf(x) := γ(θ) dθ.
x

avec ρH ≈ 2 [8]. Ce que l’on observe de façon additionnelle, L’EMD ne retenant, par construction, que les maxima posi-
c’est un effet global de décorrélation effective entre modes, à tifs et les minima négatifs, on constate alors (cf. Figure 4) que
l’exception cependant du mode 1 dont la contribution spectrale la densité de probabilité p∗ (x) := max(p− (x), p+ (x)) de ces
non négligeable dans la moitié inférieure du spectre (cf. Figure derniers est très proche de celle observée pour l’IMF 1, en
1) maintient des liens forts avec les autres modes. cohérence avec son interpétation en termes de bruit blanc filtré
passe-haut.
Statistiques marginales — Afin d’aller au-delà de l’ordre 2, il
est intéressant d’évaluer les densités de probabilité marginales
des IMF. La Figure 3 révêle à cet égard le caractère approxi- 3.2 Approche temporelle
mativement gaussien des IMF de rang au moins 2. Ces densités Une vue duale du comportement qui vient d’être décrit peut
de probabilité gaussiennes pk (x) sont de plus renormalisables être obtenue par une approche de type réponse impulsionnelle
au sens où elles se superposent toutes sur un gabarit unique en évaluant la collection des détails moyens
moyennant la transformation
  dk [δ](t) = lim+ E {dk [δ + εb](t)} ,
1 x
pk (x) → pk ε→0
β β
où δ(t) est une impulsion et b(t) un bruit additif dont le rôle
lorsque l’on choisit β := 2−(H−1)k . est de garantir l’existence d’extrema permettant d’appliquer
 12

11
1 problématique. Attractive, car intuitive, flexible, non paramé-
10 0.5
trique et auto-adaptative; mais problématique aussi car forte-
log2(maxima)

ment non linéaire, d’interprétation délicate et définie seulement

réponses
9

8 0 par un algorithme difficilement réductible à l’analyse.

7
-0.5
À mi-chemin entre l’application d’une méthode encore mal
6 comprise à des données réelles et l’attente d’une théorie bien
5
0 2 4 6 8
-1
temps établie, il nous a semblé nécessaire de recourir à des simula-
IMF
tions extensives dans des situations simples et bien contrôlées.
5 Les résultats présentés ici complètent d’autres études conduites
1
0 - - - : IMF 6 dans le même esprit [4, 8], apparentant l’EMD à une méthode
---- : spline
-5 0.5 cubique multirésolution que l’on pourrait qualifier de “spontanée” dans
spectres(dB)

réponses

-10 la mesure où elle repose sur une hiérarchie d’échelles “na-
0
-15
turelles” extraites itérativement en fonction du signal analysé.
-20
-0.5 À la différence de méthodes plus classiques à base d’ondelettes
-25
(et variations), cette extraction est faite localement et de façon
-30 -1
log(fréquence) temps fondamentalement auto-adaptative, aucun filtre pré-déterminé
ne présidant a priori à la décomposition.
F IG . 6: Auto-similarité — Le logarithme de l’amplitude maxi- Ces premiers éléments nécessitent maintenant de disposer
male des modes 2 à 5 de la Figure 5 varie linéairement en fonc- d’une véritable compréhension théorique qui permettrait une
tion de l’indice des modes (en haut, à gauche). Renormaliser de utilisation raisonnée de la méthode comme analyse exploratoire
façon couplée ces modes, soit en temps (en haut, à droite), soit de signaux réels.
en fréquence (en bas, à gauche) fournit un gabarit unique. La
Remerciements — Les auteurs ont plaisir à remercier G. Rilling
“réponse impulsionnelle” ainsi identifiée peut être comparée
(ENS Lyon) pour l’aide qu’il leur a apportée dans le dévelop-
avec profit à une ondelette de type spline cubique (en bas, à
pement des algorithmes utilisés dans cette étude.
droite).

l’EMD. On obtient là encore (cf. Figure 5) une signature car- Références
actéristique directement évocatrice de ce que donnerait une anal-
yse en ondelettes dans une situation analogue. [1] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q.
Pour préciser ce point plus avant, on peut — tout comme Zheng, N.C. Yen, C.C. Tung and H.H. Liu, “The empi-
on l’avait fait dans l’approche spectrale [8] — chercher dans la rical mode decomposition and Hilbert spectrum for non-
décomposition obtenue une signature éventuelle d’auto-similarité linear and non-stationary time series analysis,” Proc. Roy.
qui raménerait la collection des modes à un motif unique moyen- Soc. London A, Vol. 454, pp. 903–995, 1998.
nant une renormalisation judicieusement choisie. [2] E.P. Souza Neto et al., “Assessment of cardiovascular au-
La figure 6 montre qu’il en est bien ainsi. Dans un premier tonomic control by the Empirical Mode Decomposition,”
temps, le tracé de l’amplitude maximale des modes en fonction 4th Int. Workshop on Biosignal Interpretation, Como (I),
de leur indice permet d’identifier une loi exponentielle du type : pp. 123-126, 2002.
[3] R. Fournier, “Analyse stochastique modale du signal sta-
log2 dk [δ](0) = C − pk, bilométrique. Application à l’étude de l’équilibre chez
avec p ≈ 0.85. On montre alors que le choix concomitant d’un l’Homme,” Thèse de Doctorat, Univ. Paris XII Val de
facteur de dilatation α = 2p ≈ 1.80 permet de superposer les Marne, 2002.
gabarits temporels et spectraux des différents IMF, justifiant [4] G. Rilling, P. Flandrin and P. Gonçalvès, “On Empir-
une relation de la forme : ical Mode Decomposition and its algorithms,” IEEE-
EURASIP Workshop on Nonlinear Signal and Image Pro-
  cessing NSIP-03, Grado (I), 2003.
1 t
dk [δ](t) = k ψ , [5] [Link]/˜flandrin/[Link].
α αk
où ψ(t) est une forme d’onde de référence (analogue à une [6] B.B. Mandelbrot and J.W. van Ness, “Fractional Brow-
ondelette-mère dans une analyse multirésolution). Dans le cas nian motions, fractional noises and applications,” SIAM
considéré ici, la procédure d’interpolation utilisée par l’EMD Rev., Vol. 10, pp. 422–437, 1968.
repose sur des splines cubiques [1, 4, 5], et il est intéressant [7] A.T. Wood and G. Chan, “Simulation of stationary pro-
de constater que le motif ψ(t) obtenu est très proche de l’on- cesses in [0, 1]d ,” J. Comp. Graph. Stat., Vol. 3, pp. 409–
delette spline cubique obtenue comme troisème puissance de 432, 1994.
convolution de l’ondelette de Haar. [8] P. Flandrin, G. Rilling and P. Gonçalvès, “Empirical
Mode Decomposition as a filterbank,” IEEE Sig. Proc.
Lett., à paraı̂tre 2003.
4 Conclusion [9] Z. Wu and N.E. Huang, “A study of the characteristics
of white noise using the Empirical Mode Decomposition
La décomposition modale empirique (ou “EMD”) est une
method,” soumis à Proc. Roy. Soc. London A, 2002.
méthode d’analyse des signaux qui est à la fois attractive et