Thomas Oberlin - Différentes méthodes pour la décomposition modale empirique - Résultats

De Ensiwiki
Aller à : navigation, rechercher


Différentes méthodes pour la décomposition modale empirique

Labo LJK
Equipe MGMI
Encadrants valerie.perrier@imag.fr

Etudiant

OBERLIN, Thomas

Introduction

La décomposition modale empirique (ou EMD pour Empirical Mode Decomposition) est une méthode algorithmique de décomposition spectrale adaptative : au lieu d’analyser le signal dans une base fixe comme avec Fourier, on construit au fur et à mesure les fonctions de base, appelées IMF pour Intrinsic Mode Function. Elle a été introduite pour traiter de manière relativement souple des données quelconques, pouvant être non-stationnaires et non-linéaires (cf. [1], 1998). Cette méthode, relativement jeune, a fait depuis son introduction l’objet de nombreux travaux, la plupart essayant de la formaliser et de la justifier mathématiquement. En effet, le gros défaut de cette méthode est l’absence de théorie, notamment pour le sous-programme de lissage, qui extrait les IMF, appelé Sifting Process (SP), ou procédé de tamisage. Ainsi, de nombreuses études ont modifié le fonctionnement de ce sifting process, ou bien l’ont formalisé à l’aide d’outils mathématiques solides, afin de justifier la méthode. Parmi ces travaux, nous nous sommes intéressé une méthode basée sur de l’optimisation sous contraintes, introduite par V. Perrier et S. Meignen (cf. [2]), ainsi qu’à une toute récente tentative pour modéliser le SP par une équation différentielle (cf. [3]).

Le TER a consisté en une étude de la décomposition modale empirique et de ces deux variantes. On a ensuite implémenté la méthode 3 sur Matlab (§ 1.4.2), et réalisé de nombreux tests de comparaison (§ 1.5).

EMD : méthode originale

Principe

Considérons un signal s, l'EMD va le décomposer en une somme finie de modes oscillants. On pourra alors écrire :

s = \sum_{k=1}^{K}d_k+r

Les modes oscillants d_k sont appelés IMF (Intrinsic Mode Function), sont des fonctions oscillantes autour de 0, et de moyenne locale nulle. Le résidu r est non-oscillant, c'est-à-dire qu'il contient au plus 3 extrema.

La méthode originale définit l'enveloppe supérieure d'un signal (resp. inférieure) comme l'interpolation par splines cubiques des maxima (resp. minima) du signal. On peut alors définir l'enveloppe moyenne comme moyenne de ces deux enveloppes.

L'algorithme consiste donc à extraire successivement les modes en soustrayant au signal son enveloppe moyenne. Cependant, les modes ainsi extraits ne sont pas des IMF : bien qu'ils oscillent autour de 0, ils ont une moyenne locale non nulle. Le rôle du sifting process est de tamiser ces modes pour obtenir de véritables IMF. Pour cela, on soustrait au mode son enveloppe moyenne, plusieurs fois, jusqu'à obtenir une moyenne locale (ou enveloppe moyenne quasi-nulle).

Algorithme

Détaillons précisément l'algorithme de l'EMD :

  • 1. Initialisation : r = s, k = 1
  • 2. Calculer la moyenne locale e de r
  • 3. Extraire le mode oscillant  p_i = r - e , poser r = e
  • 4. Répéter, jusqu’à satisfaction d’un critère  stop(p_i)  :

Calculer l’enveloppe moyenne e_i de p_i

 p_{i+1} = p_i - e_i , i=i+1

  • 5. d_k = p_i,  r = r - d_k
  • 6. Si r a plus de trois extrema, aller à l’étape 2 en posant k = k + 1. Sinon, la décomposition est achevée

Un premier exemple

On donne ici la décomposition obtenue pour une somme de cosinus (IMF 1 à 3 et résidu) : l'EMD sépare les trois composantes du signal, et donne un résidu non oscillant de faible amplitude.

Signal original : s(t)=cos(7t)+2cos(3t)+0.3cos(t)


Défauts de l'EMD, présentation de deux variantes

Si la méthode originale est appliquée avec succès dans de nombreux domaines, il lui manque toujours des justifications théoriques : il n'y a pas de définition claire pour la "moyenne locale", pas plus que pour les IMF, ce qui exclue toute preuve de la convergence du SP. Nous allons présenter ci-dessous deux variantes, moins empiriques mais justifiées mathématiquement.

EMD par optimisation sous contraintes

Cette méthode a été introduite par S. Meignen et V. Perrier. Elle consiste à déterminer l'enveloppe moyenne du signal en résolvant un problème d'optimisation sous contraintes. On remplace l’interpolation spline par une interpolation de Hermite, avec raccordements C1 aux noeuds (au lieu de C2 pour les splines). On gagne donc plusieurs degrés de liberté. Cela nous permet d'introduire des contraintes d'égalité et d'inégalité, notamment pour préserver la monotonicité des enveloppes. L'enveloppe moyenne sera la solution d'énergie minimale.


Modélisation du SP par une EDP

Cette méthode est très récente, et n'a pas fait l'objet de beaucoup d'études. Elle consiste à simuler le sifting process par une équation aux dérivées partielles.

On se propose ici de changer la définition de l'enveloppe moyenne. Soit \delta un réel strictement positif fixé, définissons l'enveloppe supérieure d'un signal s par :

S_{\delta}s(x)=\sup_{|y|<\delta}s(x+y)

De même, on définit l'enveloppe inférieure :

I_{\delta}s(x)=\inf_{|y|<\delta}s(x+y)

Puis enfin l'enveloppe moyenne e, qui est la moyenne de ces deux enveloppes. D'après cette définition, on voit que le choix de \delta est fondamental :

  • Si \delta est trop petit, on ne filtrera pas les hautes fréquences. On peut d'ailleurs remarquer, en faisant l'hypothèse de la continuité de s, que :
\forall x, \lim_{\delta \rightarrow 0} E_{\delta}s(x) = s(x)
  • Si \delta est trop grand, on introduira des paliers artificiels dans la décomposition. Le cas limite est donné par :
\forall x,\lim_{\delta \rightarrow +\infty} E_{\delta}s(x) = \frac{\sup(s)+\inf(s)}{2}

Observons à présent le comportement de ces enveloppes, selon la valeur de delta. Ici, la valeur convenable de \delta est \delta_0=0.19. On observe également l'allure des enveloppes pour \delta=\frac{\delta_0}{3} et \delta=3\delta_0.

Forme des enveloppes pour s(t)=sin(3*t)+sin(15t), selon la valeur de \delta

Introduisons à présent une fonction h(x,t) telle que, \forall n \in \mathbb{N}, h(x,n\delta)=p_{n}(x). Cela signifie que h(.,n\delta) est égale à la n-ième itération du sifting process du signal. Les auteurs de la méthode montrent que, sous certaines hypothèses (pour un signal de classe C^2 notamment), la fonction h vérifie :


\begin{cases}
\frac{\partial h}{\partial t}+\frac{1}{\delta^2}h+\frac{1}{2}\frac{\partial^2h}{\partial x^2} & =0 \\
h(x,0) & =s(x) \forall x \in \Omega
\end{cases}

Cette équation est celle de la chaleur, mais avec un coefficient de diffusion négatif. Cela pose problème lors de la résolution numérique, car il n'y a pas de condition de stabilité explicite. Dans ce TER, nous avons utilisé d'une part une méthode aux différences finies, comme les créteurs de la méthode, mais aussi une résolution dans le domaine de Fourier, comme expliqué ci-après.


Notons H(x,t)=h exp \left(\frac{1}{\delta^2}\right). L'équation s'écrit alors :


\begin{cases}
\frac{\partial H}{\partial t}+\frac{1}{2}\frac{\partial^2H}{\partial x^2} & =0 \\
H(x,0) &=s(x)exp \left(\frac{1}{\delta^2}\right) \forall x \in \Omega
\end{cases}

Dans le domaine de Fourier, cette équation donne :



\begin{cases}
\frac{\partial \hat{H} }{\partial t}+ \left(1-2\pi^2\nu^2\right)\hat{H} & = 0 \\
\hat{H}(\nu,0) & =\hat{s}(\nu)exp \left(\frac{1}{\delta^2}\right) \forall x \in \Omega
\end{cases}

On en déduit l'expression explicite de \hat{h} :

\hat{h}(\nu,t)=\hat{s}(\nu) e^{-\frac{t}{\delta^2}}e^{2\pi^2\nu^2t}

Si le signal s est à fréquences bornées, \hat{s} est à support compact, on peut alors écrire la transformée de Fourier inverse :

h(x,t)=\int_{\R}\hat{s}(\nu) e^{-\frac{t}{\delta^2}}e^{2\pi^2\nu^2t+2i\pi\nu x}d\nu


Comparaison des méthodes

Les tests et comparaisons des trois méthode représente une partie conséquente de ce TER. Nous allons ici nous limiter à deux critères de comparaison : la capacité des méthodes à séparer deux composantes de fréquences proches, et celle à débruiter un signal quelconque. Remarquons qu'il est également intéressant d'observer le comportement des méthodes pour des signaux non linéaires, mais nous ne présentons pas ici de tels tests, par manque de place.

Séparation de composantes

Le signal est une simple somme de sinusoïdes. On donne ici, pour chaque méthode, la première IMF obtenue. On peut observer ici que la méthode 1 extrait bien le sinus de plus haute fréquence, malgré de faibles effets de bord. La méthode 2, quant à elle, donne un mode de même fréquence, mais modulé en amplitude. En fait, cela est dû aux contraintes d'égalité qui sont relachées lorsque deux maxima successifs sont égaux. Cela ne s'observerait pas sur des signaux réels. Enfin, la méthode 3 réagi convenablement : on a obervé sur d'autres tests que c'est celle qui sépare le mieux deux sinus de fréquence très proche. Néanmoins, à cause des paramètres \delta et n, on obtient un mode d'amplitude aberrante.


Signal original : s(t)=sin(4\pi t)+0.6sin(12\pi t)
Première IMF - Méthode 1
Première IMF - Méthode 2
Première IMF - Méthode 3































Débruitage de signaux

Ce test évalue le comportement des IMF pour des signaux bruités : après avoir décomposé le signal en différents modes, on lui soustrait les IMF correspondant au bruit (les 2 à 4 premières, cela dépend). Ici encore, la méthode originale s'avère très efficace. La méthode 2 conserve le signal mais laisse un bruit résiduel, et on retrouve pour la méthode 3 les problèmes précédents de normalisation. Par ailleurs, le choix du paramètre \delta est difficile lorsqu'on ne connaît pas analytiquement le signal, ce qui est le cas ici.

Signal original : s(t)=sin(5\pi t)+0.3sin(12\pi t)
Signal bruité
Signal débruité - Méthode 1
Signal débruité - Méthode 2
Signal débruité - Méthode 3














































Conclusion

Ce TER a été l'occasion de découvrir la décomposition modale empirique, ses possibilités et ses limites, mais aussi d'étudier en détail des variantes très récentes. Ainsi, on a pu implémenter et tester la méthode 3, pour la comparer aux deux autres.

Avec plus de temps, il aurait été intéressant d'approfondir l'étude des deux variantes : observer précisément les effets des contraintes dans la méthode 2, et tenter de mieux définir les paramètres pour la méthode 3.

Références

  • [1] N. E. Huang et al. The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis. The Royal Society, 1998.
  • [2] S. Meignen and V. Perrier. A new formulation for empirical mode decomposition based on constrained optimisazion. IEEE Signal Processing Letters, 2007.
  • [3] El Hadji S. Diop, R. Alexandre, and A. O. Boudraa. A pde formulation of the sifting process.
  • [4] G. Rilling and P. Flandrin : codes Matlab
  • [5] Gabriel Rilling. Décompositions modales empiriques. PhD thesis, Ecole normale supéieure de Lyon, 2007.
  • [6] G. Rilling and P. Flandrin. Empirical mode decomposition as a filter bank. IEEE Signal Processing Letters, 2004.
  • [7] C. Damerval. Empirical mode decomposition. Master’s thesis, 2004.
  • [8] R. C. Sharpley and V. Vatchev. Analysis of the Intrinsic Mode Functions. 2006.


Documents additionnels