# Méthodes de traitement du signal par décomposition en modes empiriques

Nicolas Di Palma, Alain Batailly, Mathias Legrand

► To cite this version:

Nicolas Di Palma, Alain Batailly, Mathias Legrand. Méthodes de traitement du signal par décomposition en modes empiriques. [Rapport de recherche] Université McGill. 2018. hal-01819670

HAL Id: hal-01819670

https://hal.science/hal-01819670v1

Submitted on 20 Jun 2018

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 CC BY 4.0 - Attribution - International License
# Signal analysis techniques based on empirical mode decomposition

Nicolas Di Palma<fnref n="1" />, Alain Batailly<fnref n="2" />*, Mathias Legrand<fnref n="1" />

## Abstract

This document aims at presenting, assessing and comparing several signal analysis techniques relying on the empirical mode decomposition. This adaptive method, developed by Huang, is frequently used for the analysis of nonlinear signals, i.e. signals stemming from nonlinear systems. One weakness of the method lies in its lack of theoretical foundation, though it should be mentioned that some authors attempted to make a more mathematical presentation of the method. Nonetheless, the EMD has been applied for a large variety of applications—such as industrial machines monitoring and medical signals analysis to name a few. The EMD is the cornerstone of increasingly popular signal analysis techniques such as the local mean decomposition and the Hilbert-Huang Transform. These methods are detailed in this note along with the associated Python implementation that is freely downloadable with this document.

## Keywords

Empirical mode decomposition; Local mean decomposition; Hilbert-Huang Transform; Edge effects; Signal analysis

# Méthodes de traitement du signal par décomposition en modes empiriques

Nicolas Di Palma<fnref n="1" />, Alain Batailly<fnref n="2" />*, Mathias Legrand<fnref n="1" />

## Résumé

Cet article présente, évalue et compare différentes techniques de traitement du signal reposant sur la décomposition en modes empiriques. Initialement introduite par Huang, la décomposition en modes empiriques est une méthode adaptative fréquemment utilisée pour l'analyse de réponses vibratoires de systèmes mécaniques non-linéaires. Malheureusement, l'assise théorique de la méthode est manquante même si quelques contributions abouties proposent une présentation mathématique. Elle a néanmoins été appliquée avec succès dans des domaines aussi variés que le diagnostic de machines industrielles ou l'analyse de signaux médicaux, et constitue le socle commun à plusieurs méthodes récentes de traitement du signal telles que la décomposition en moyenne locale et la méthode de Hilbert-Huang. Ces méthodes sont détaillées et leur mise en œuvre, disponible en annexe en langage Python, est évaluée pour quelques signaux académiques.

## Mots-clés

Décomposition en mode empirique; Décomposition en moyenne locale; Transformée de Hilbert-Huang; Effets de bord; Traitement du signal

## 1. Introduction

Le développement de méthodes numériques pour le traitement de signaux non-linéaires (signaux issus de systèmes mécaniques gouvernés par des équations du mouvement non-linéaires) et non-stationnaires (signaux dépendant du temps) a été l'objet de nombreuses publications, notamment au cours des vingt dernières années. De tels signaux ne peuvent pas être analysés de façon optimale par les méthodes usuelles de traitement du signal (transformée de Fourier ou ondelettes) qui ont été développées pour des signaux linéaires. En effet, lorsque la transformée de Fourier est utilisée pour des signaux non-linéaires non-stationnaires, typiquement sous forme de spectrogramme, des limites de résolution apparaissent. D'autre part, la transformée en ondelettes présente l'inconvénient de devoir adapter les familles d'ondelettes au signal (densité du contenu fréquentiel, discontinuités éventuelles...).

Les méthodes adaptatives, dont il est question ici, offrent une perspective intéressante pour l'analyse de signaux non-linéaires et non-stationnaires. L'incarnation la plus connue est probablement la méthode dite de Hilbert-Huang [1]. Elle comprend deux étapes distinctes: (1) une première décomposition du signal en une somme de signaux appelés modes empiriques; on parle de décomposition en modes empiriques<fnref n="1" /> et (2) une transformée de Hilbert des modes empiriques afin d'obtenir les fréquences instantanées associées. L'EMD est une procédure décrite dans [1] qui ne dispose pas, actuellement, de fondement mathématique clair mais repose plutôt sur des considérations intuitives. Elle a été l'objet de nombreux travaux théoriques [2, 3, 4, 5, 6, 7] et appliqués [8, 9, 10, 11, 12, 13].

<footnote n="1">1. Structural Dynamics and Vibration Laboratory, McGill University, 817 Sherbrooke West, McConnell Engineering Bldg, Room 122, H3A 0C3, Montréal, Québec, Canada</footnote>

<footnote n="2">2. Laboratory of Acoustics and Vibration Analysis, École Polytechnique de Montréal, C.P. 6079, Succ. Centre-ville, Montréal, Québec, CANADA, H3C 3A7</footnote>

<footnote>* Corresponding author: Tel.: +1 514 340 4711 - 3932. E-mail address: alain.batailly@polymtl.ca</footnote>

<footnote n="1">1. Laboratoire de Dynamique des Structures et Vibrations, Université McGill, 817 Sherbrooke Ouest, Pavillon McConnell, Salle 122, H3A-0C3, Montréal, Québec, Canada</footnote>

<footnote n="2">2. Laboratoire d'Analyse Vibratoire et Acoustique, École Polytechnique de Montréal, C.P. 6079, Succ. Centre-ville, Montréal, Québec, CANADA, H3C 3A7</footnote>

<footnote>* Corresponding author: Tel.: +1 514 340 4711 - 3932. E-mail address: alain.batailly@polymtl.ca</footnote>

<footnote n="1">¹notée EMD pour *Empirical Mode Decomposition*</footnote>
L'absence de fondement théorique a permis l'éclosion de formulations dérivées de l'EMD [14, 15, 16], la plus célèbre étant certainement la méthode de décomposition en moyenne locale<fnref n="2" /> [15]. Ce foisonnement a généré une confusion lorsqu'il s'agit de déterminer la formulation la plus adaptée à un signal. Il ne semble pas pertinent d'envisager une classification générale puisque l'EMD a été employée aussi bien dans le domaine médical que dans le diagnostic des machines tournantes, et qu'il en va de même pour ses dérivées.

## 2. Décomposition en modes empiriques

La décomposition en modes empiriques constitue la première étape de la transformée de Hilbert-Huang [1]. Par un processus itératif, cette méthode extrait d'un signal un nombre fini de composantes, de la plus oscillante au résidu final monotone.

### Conventions et notations

La décomposition d'un signal $x(t)$ quelconque, et fonction du temps $t$ qui sera une quantité discrète dans ce qui suit, est notée

$$x(t) = \sum_{i=1}^{n} y_{i}(t) + r(t), \quad \forall t \qquad (1)$$

où le résidu de la décomposition est noté $r$ et les autres composantes, nommées modes intrinsèques [1] et obtenues par itérations successives, $\gamma_i$, $i = 1, \dots, n$. D'un point de vue mathématique, il convient de vérifier que la famille des fonctions modales intrinsèques<fnref n="3" /> $\gamma_i$ est complète et orthogonale. La complétude est assurée par la définition (1) alors que des critères existent pour vérifier les propriétés d'orthogonalité [1]:

1. la différence entre le nombre d'extrema et le nombre de zéros de la fonction doit être nulle ou au plus égale à 1;

2. la valeur moyenne de la somme des enveloppes supérieure, définies par les maxima, et inférieure, définies par les minima de la fonction, doit être nulle.

<figure><img src="image_2.png" /><figcaption>Figure 1. Exemple de décomposition en modes empiriques. L'axe horizontal est le temps $t$.</figcaption></figure>

Une décomposition par EMD est affichée sur la figure 1: le signal $x(t)$ est décomposé en somme de deux fonctions oscillantes $\gamma_1(t)$ et $\gamma_2(t)$ et d'un résidu $r(t)$. La première IMF $\gamma_1(t)$, obtenue lors de la première itération, présente plus d'oscillations que la seconde composante $\gamma_2(t)$: le contenu fréquentiel d'un signal est discriminé localement et réparti sur les différentes IMF. Le résidu $r(t)$ est monotone: aucune autre composante oscillante ne peut en être extraite.

La construction de chaque IMF s'effectue à l'aide d'une seconde boucle itérative: c'est le processus de tamisage décrit en détail dans les sections suivantes.

**Notations** Dans la suite de cet article, l'indice $i$ est utilisé pour les itérations sur la boucle principale relative à l'obtention d'une IMF $\gamma_i(t)$ alors que l'indice $k$ est utilisé pour les étapes du processus de tamisage. Les fonctions obtenues à l'itération $k$ du calcul de la $i^{\text{e}}$ IMF $\gamma_i(t)$ sont ainsi indicées $i, k$.

### Processus de tamisage

Soit $x(t)$ un signal discret quelconque à décomposer. Il est seulement supposé que le signal $x(t)$ contient au moins deux extrema (un minimum et un maximum). Quelques remarques s'imposent:

<footnote n="2">²notée LMD pour Local Mean Decomposition</footnote>

<footnote n="3">³notée IMF pour Intrinsic Modal Function.</footnote>
1. l'échelle de temps caractéristique est définie à partir du temps séparant deux extrema;

2. si le signal ne comporte pas d'extrema mais plusieurs points d'inflexion, il peut être dérivé une ou plusieurs fois afin de faire apparaître les extrema. Le résultat final sera ensuite obtenu par l'intégration de ses composantes extraites.

## Calcul des modes intrinsèques

Lors de l'obtention de l'IMF $\gamma_i(t)$, la première étape du processus de tamisage consiste à déterminer les maxima et les minima du signal. Une spline cubique relie ensuite chacun des minima afin de construire l'enveloppe inférieure $e_{i,1}^{\inf}(t)$; par une procédure similaire sur les maxima, on obtient l'enveloppe supérieure $e_{i,1}^{\sup}(t)$. La moyenne $m_{i,1}(t)$ des enveloppes supérieure $e_{i,1}^{\sup}$ et inférieure $e_{i,1}^{\inf}$ est alors donnée par la moyenne

$$m_{i,1}(t) = (e_{i,1}^{\sup}(t) + e_{i,1}^{\inf}(t))/2 \qquad (2)$$

qui est soustraite du signal initial $r_i(t)$ du processus de tamisage (avec $r_1(t) = x(t)$) afin d'obtenir la première composante

$$h_{i,1} = r_i(t) - m_{i,1}. \qquad (3)$$

Les étapes associées aux équations (2) et (3) sont illustrées sur la figure 2. Une fois $h_{i,1}(t)$ obtenue, deux situations peuvent se

<figure><img src="image_3.png" /><figcaption>Figure 2. Première itération du processus de tamisage: signal initial $r_i(t)$ (—), enveloppes (—), moyenne $m_{i,1}(t)$ (—) et résultat $h_{i,1}(t)$ (—)</figcaption></figure>

présenter:

1. le signal $h_{i,1}$ vérifie les deux critères définissant une IMF (voir section 2.1) et alors $\gamma_i(t) = h_{i,1}(t)$;

2. $h_{i,1}(t)$ n'est pas une IMF et le processus est répété jusqu'à ce que $h_{i,k}(t)$ soit une IMF notée $\gamma_i(t) = h_{i,k}(t)$.

Les critères définissant une IMF peuvent être coûteux à vérifier numériquement. Ils sont généralement traduits sous forme d'un critère global basé sur l'écart type entre deux itérations successives du signal $h_{i,k}(t)$ [1]:

$$\sigma_{i,k} = \sum_{t} \left( \frac{|(h_{i,k-1}(t) - h_{i,k}(t))|^2}{h_{i,k-1}^2(t)} \right). \qquad (4)$$

L'algorithme de l'EMD est résumé par l'organigramme de la figure 3.

## Précisions sur les critères utilisés

Puisque le signal $h_{i,k}(t)$ peut s'annuler en un grand nombre de points, sa présence au dénominateur peut générer des problèmes numériques. On lui préfère souvent le critère [17]

$$\sigma_{i,k}^* = \frac{\|h_{i,k-1}(t) - h_{i,k}(t)\|^2}{\|h_{i,k-1}(t)\|^2}. \qquad (5)$$

Les deux critères globaux $\sigma_{i,k}$ et $\sigma_{i,k}^*$ peuvent, sous certaines conditions, être trop permissifs. Des approches locales, permettant de respecter plus scrupuleusement les critères d'obtention d'une IMF, ont également été développées [5].

Il est possible de montrer que $\sigma_{i,k}^*$ est un critère moins contraignant que $\sigma_{i,k}$, autrement dit $\sigma_{i,k}^* < \sigma_{i,k}$. Afin de simplifier la présentation, on note $a_j = h_{i,k-1}(t) - h_{i,k}(t)$ et $b_j = h_{i,k-1}(t)$. Les deux critères se réécrivent alors

$$\sigma_{i,k} = \sum_{j} \frac{a_j^2}{b_j^2} \quad \text{et} \quad \sigma_{i,k}^* = \frac{\sum_j a_j^2}{\sum_j b_j^2} \quad \text{avec} \quad \sum_j b_j^2 > 0. \qquad (6)$$

Il est donc possible d'écrire:

$$\sigma_{i,k} = \left( \sum_j \frac{a_j^2}{b_j^2} \sum_j b_j^2 \right) / \left( \sum_j b_j^2 \right) \quad \text{et} \quad \sigma_{i,k} = \left( \sum_k a_k^2 \frac{\sum_j b_j^2}{b_k^2} \right) / \left( \sum_j b_j^2 \right). \qquad (7)$$
<figure><img src="image_4.png" /><figcaption>Figure 3. Principe de la décomposition en modes empiriques avec, en fin de procédure, $r(t) \leftarrow r_{i+1}(t)$</figcaption></figure>

Le simple fait de supposer l'existence $\sigma_{i,k}$ implique $b_k^2 > 0$. De plus, en observant que $\sum_j b_j^2 = b_k^2 + \sum_{j \neq k} b_j^2 > b_k^2$, il est évident que $\sum_j b_j^2 > b_k^2$. Ainsi

$$ \sum_k a_k^2 \frac{\sum_j b_j^2}{b_k^2} > \sum_j a_j^2 \qquad (8) $$

et on en conclut donc $\sigma_{i,k} > \sigma_{i,k}^*$. Les valeurs usuelles des critères $\sigma_{ik}$ et $\sigma_{i,k}^*$ pour déterminer si $h_{i,k}(t)$ est l'IMF $\gamma_i(t)$ sont $\sigma_{i,k} < 0,3$ et $\sigma_{i,k}^* < 0,0003$.

## Effets de bord

Un des points sensibles de la décomposition en modes empiriques est sa grande sensibilité aux effets de bord. Ils correspondent à la difficulté de définir le comportement des enveloppes supérieure $e_{i,k}^{sup}(t)$ et inférieure $e_{i,k}^{inf}(t)$ au niveau des bords du domaine de définition du signal durant le processus de tamisage, avec un impact immédiat sur le calcul de la moyenne $m_{i,k}(t)$ définie dans l'équation (2). La figure 4 illustre cette difficulté. Les enveloppes supérieure $e_{i,k}^{sup}(t)$ et inférieure $e_{i,k}^{inf}(t)$ sont obtenues à partir de

<figure><img src="image_5.png" /><figcaption>Figure 4. Effets de bord, ici à gauche uniquement, le signal est défini pour $t > 5$ et le bord droit n'est donc pas représenté. Enveloppes supérieure et inférieure (—) moyenne $m_{1,1}(t)$ du signal (—) composante $h_{1,1}(t)$ (—).</figcaption></figure>

tous les maxima et minima. Elles doivent donc être prolongées pour être définies sur le même domaine que le signal $x(t)$. À partir du signal $x(t)$ de la figure 4, il est évident qu'il existe une infinité d'enveloppes admissibles. Selon le choix des enveloppes, la moyenne $m_{1,1}(t)$ et la composante $h_{1,1}(t)$ qui en résulte peuvent varier de façon très significative au voisinage du bord gauche, c'est-à-dire pour $t = 0$ sur la figure. Bien que non représenté, le problème symétrique émerge à $t = t_{max}$.
Sans prise en compte des effets de bord, le processus de tamisage peut ne pas converger et $h_{i,k}(t=0) \to \pm\infty$ ou $h_{i,k}(t=t_{\max}) \to \pm\infty$

Plusieurs méthodes sont à notre disposition pour minimiser les effets de bord. Il s'agit en général de prolonger le signal à gauche et à droite pour finalement repousser les effets de bord en dehors du domaine d'étude. C'est donc essentiellement la méthode prolongeant le signal qui détermine la qualité de la décomposition en modes empiriques. Parmi les méthodes évoquées dans la littérature, notons l'utilisation de machines à vecteurs de support [18, 19], de réseaux de neurones [20], de techniques de prolongation des extrema [5] et d'ajustement de forme [21].

Par défaut, c'est la méthode proposée dans [5] qui est utilisée pour toutes les décompositions: des maxima sont ajoutés sur chaque bord par symétrie axiale, l'axe de symétrie étant l'axe $x=b$ où $b$ correspond à chacune des bornes de l'intervalle de temps considéré. Cette méthode présente l'avantage d'être relativement peu coûteuse en temps de calcul, avec, pour les signaux testés, une convergence rapide.

## Validation

La validation de la décomposition en modes empiriques est conduite à la fois sur des signaux empiriques/académiques. Compte-tenu du relatif faible nombre de points des signaux académiques (entre $5 \cdot 10^3$ et $10^5$ points sur l'intervalle de temps considéré), les temps de calculs sont négligeables et ne sont pas reportés dans la première sous-section.

## Signaux académiques

La décomposition en modes empiriques est réalisée sur des signaux académiques, c'est-à-dire des signaux issus de fonctions usuelles, parmi lesquels:

**signal 1:** $x(t) = \sin(2\pi t) + \cos(6\pi t) + \sin(2\pi t)^2 \sin(6\pi t)$, $t \in [0; 10]$, $f = 1000$ Hz;

**signal 2:** $x(t) = t^2 \sin(t^2)$, $t \in [0; 5\pi]$, $f = 1000$ Hz;

**signal 3:** $x(t) = \cos(t^2) + 3 \cos(3t) + \cos(t)$, $t \in [0; 8\pi]$, $f = 1000$ Hz;

**signal 4:** $x(t) = (\exp(-0,2t)t^2 + 10) \sin(\exp(-0,1t)t^2 + 3t) + \cos(2t) + \sin(150t) + \cos(74t)$, $t \in [0; 20\pi]$, $f = 1590$ Hz;

**signal 5:** $x(t) = (\exp(-0,2t)t^2 + 10)(\sin(\exp(-0,1t)t^2 + 3t) + \cos(2t))$, $t \in [0; 30\pi]$, $f = 100$ Hz;

où $f$ est la fréquence d'échantillonnage de la discrétisation du signal. Chacun de ces signaux présente un intérêt particulier: répétition d'un motif, évolution régulière de la fréquence en fonction du temps, présence de bruit, combinaison de motifs...

La décomposition de la figure 5a illustre le comportement de l'algorithme sur un signal relativement simple, contenant peu de motifs et donc *a priori* sensible aux effets de bords. Deux IMF non négligeables sont extraites du signal et les effets de bord sont perceptibles sur un faible intervalle de temps à chaque extrémité du signal. Cet intervalle de temps est de l'ordre de grandeur de la période du signal original.

Le deuxième signal considéré est $x(t) = t^2 \sin(t^2)$, sa décomposition est donnée sur la figure 5c. Ce signal peut représenter, par exemple, un phénomène divergent: en fonction du temps, son amplitude évolue quadratiquement alors que sa fréquence croît linéairement. Il est remarquable que la décomposition en modes empiriques se fait sur une seule IMF. La deuxième et dernière IMF est d'amplitude négligeable. Pour un tel signal, les effets de bords sont négligeables aux deux extrémités de l'intervalle et attestent ainsi de la bonne tenue de la procédure qui converge rapidement (le résidu monotone est bien trouvé).

Le troisième signal est sensiblement plus complexe que les précédents puisqu'il combine une composante de haute fréquence de type *bruit* à une composante basse fréquence déterministe comme indiqué sur la figure 5b. La superposition de contenus fréquentiels distincts entraîne logiquement la détection de plusieurs IMF successives. Ce signal met en évidence le caractère local de la discrimination du contenu fréquentiel d'une IMF à une autre. En effet, sur les premiers 20 % de l'intervalle de temps, $x(t)$ n'est pas bruité: le contenu de plus haute fréquence correspond au signal de plus grande amplitude ce qui explique que l'amplitude de $\gamma_1$ sur cet intervalle est identique à celle de $x$. En dehors de cet intervalle, le bruit issu du terme $\cos(t^2)$ devient non négligeable mais son amplitude est plus faible (entre -1 et 1) et c'est pourquoi la première IMF oscille entre ces deux bornes. Le terme $3 \cos(3t)$ est capté quasi parfaitement par $\gamma_2$ et le terme $\cos(t)$, de fréquence plus faible, se retrouve dans la dernière IMF $\gamma_3$.

Le signal décomposé sur la figure 6a combine deux termes à relative haute fréquence et un terme de grande amplitude dont la fréquence et l'amplitude varient sensiblement en fonction du temps. La décomposition en modes empiriques permet de débruiter le signal très efficacement et les deux composantes à haute fréquence apparaissent respectivement dans les fonctions $\gamma_1$ et $\gamma_2$. La fonction IMF $\gamma_3$ contient logiquement le terme restant, de plus grande amplitude et de fréquence variable, les IMF successives sont d'amplitude négligeable.

Le dernier signal académique est présenté sur la figure 6b. Son amplitude et sa fréquence varient fortement en fonction du temps. En particulier, les plus hautes fréquences sont observées sur les premiers 20 % de l'intervalle de temps: une multitude d'extrema locaux se succèdent sans intersecter la ligne moyenne du signal (ici égale à 0) ce qui se traduit par la présence d'amplitudes non négligeables sur $\gamma_1$ et $\gamma_2$. Pour ce signal comme les précédents, les effets de bords sont très limités et leur influence est localisée, essentiellement sur les dernières IMF.

Il est important de noter que les décompositions en modes empiriques de tous les signaux ont été réalisées par le même code sans modifier la procédure de gestion des effets de bord ce qui atteste ainsi de son bon comportement pour différents signaux.

## Robustesse de la décomposition à...

**... la fréquence d'échantillonnage** Il s'agit de vérifier que l'algorithme de décomposition en modes empiriques est indépendant de la fréquence d'échantillonnage. Plus précisément, il doit exister une fréquence d'échantillonnage $f_{min}$ au-delà de laquelle les IMF se superposent parfaitement. À titre d'exemple, le signal de la figure 6b est ré-échantillonné à 25, 50, et 200 Hz. Les IMF correspondantes sont tracées sur la figure 7. Elles se superposent parfaitement. La même observation tient pour tous les autres
<figure><img src="image_6.png" /><figcaption>Figure 5. Exemples de décomposition en modes empiriques</figcaption></figure>

signaux étudiés et il apparaît de façon assez intuitive que la fréquence d'échantillonnage minimale est la fréquence qui permet de représenter proprement le signal $x(t)$, c'est-à-dire la fréquence maximale du spectre de ce signal. Pour plus d'informations sur l'influence de la fréquence d'échantillonnage sur l'EMD, le lecteur peut se référer à [4].

### l'intervalle de temps considéré

Ce point est particulièrement sensible pour la méthode de décomposition en modes empiriques notamment lorsqu'une application en temps réel est envisagée. Le signal des figures 6b et 7 est utilisé dans cette section. Comme mentionné dans la section 2.4.1, ce signal est considéré sur l'intervalle de temps $t \in [0; 30\pi]$ avec $f = 100$ Hz. La fréquence d'échantillonnage $f$ est constante dans cette section et plusieurs sous-blocs de ce signal sont décomposés par EMD sur les intervalles de temps $t \in [0; 20\pi]$, $[2\pi; 22\pi]$, $[4\pi; 24\pi]$, $[6\pi; 26\pi]$, $[8\pi; 28\pi]$ et $[10\pi; 30\pi]$. Les IMF affichées sur la figure 8 pour chacun des sous-blocs se superposent bien. Cela montre que la méthode est peu sensible à l'intervalle de temps d'analyse. Toutefois, l'influence des effets de bord est visible et, très localement au niveau des bornes de l'intervalle de temps, des différences non négligeables sont observées. Ceci est particulièrement vrai pour $\gamma_2$ et $\gamma_3$.

Le domaine temporel sur lequel les IMF ne sont pas parfaitement superposées est ici de l'ordre de 5 à 10 % de l'intervalle de temps sur lequel est défini un sous-bloc du signal. Il est important de noter que ce ratio est ici exagéré par le fait que, pour des raisons de lisibilité des graphiques, les signaux sont considérés sur des intervalles de temps courts par rapport à leur période caractéristique (laps de temps moyen entre deux extrema consécutifs). Dans le cas d'une application industrielle, par exemple le diagnostic et la caractérisation vibratoire de composants industriels, comme les machines tournantes par exemple, le laps de temps entre deux extrema consécutifs (de l'ordre de grandeur de la microseconde) deviendrait négligeable devant la durée totale du signal (de l'ordre de la seconde) et ce ratio serait très fortement diminué.
<figure><img src="image_7.png" /><figcaption>Figure 6. EMD de deux signaux académiques</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. EMD de $x(t) = (\exp(-0,2t)t^2 + 10)(\sin(\exp(-0,1t)t^2 + 3t) + \cos(2t))$ pour plusieurs fréquences d'échantillonnage: $f = 25$ Hz (—), $f = 50$ Hz (—), $f = 100$ Hz (—) et $f = 200$ Hz (—)</figcaption></figure>

## Comparaison à d'autres versions de l'EMD

Plusieurs codes dédiés à la décomposition en modes empiriques sont librement accessibles sur internet. Certains d'entre-eux sont répertoriés dans le tableau 1. L'un de ces codes [22] a été développé par Gabriel Rilling et Patrick Flandrin avex Matlab comme support. L'algorithme utilisé dans ce travail est comparé au code développé par Rilling et Flandrin pour les raisons suivantes:

1. Rilling et Flandrin ont une expérience significative avec la méthode de décomposition en modes empiriques comme en attestent leurs nombreuses publications sur le sujet [5, 2];

2. le code de Rilling et Flandrin n'a pas été développé pour une application industrielle précise, il est ainsi plus général que d'autres codes disponibles sur internet;

3. Rilling et Flandrin se sont intéressés de près aux effets de bords et leur code semble proposer les options les plus avancées pour réduire leur répercussion sur la décomposition.

Quelques uns des signaux étudiés dans les sections précédentes sont décomposés par le code de Rilling et Flandrin<fnref n="4" />. Les décompositions de deux signaux académiques, déjà étudiés et affichés sur les figures 5b et 6b, obtenues par chacun des deux codes sont respectivement superposées sur les figures 9a et 9b. Dans l'ensemble les IMF construites par chaque code se superposent

<footnote n="4">⁴Dans cette section, nos résultats sont tracés en vert (—) et ceux de Rilling et Flandrin sont tracés en orange (—).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** EMD de $x(t) = (\exp(-0,2t)t^2 + 10)(\sin(\exp(-0,1t)t^2 + 3t) + \cos(2t))$ sur différents intervalles de temps: $t \in [0; 20\pi]$ (—), $[2\pi; 22\pi]$ (—), $[4\pi; 24\pi]$ (—), $[6\pi; 26\pi]$ (—), $[8\pi; 28\pi]$ (—) et $[10\pi; 30\pi]$ (—)</figcaption></figure>

<table><caption>**Table 1.** Liste non-exhaustive de codes EMD disponibles sur internet</caption><thead><tr><td>auteurs</td><td>langage</td><td>contexte</td><td>téléchargement</td></tr></thead><tbody><tr><td>Rilling</td><td>Matlab</td><td>mathématique</td><td>lien [22]</td></tr><tr><td>Laszuk</td><td>Python</td><td>ingénierie</td><td>lien [23]</td></tr><tr><td>Lambert</td><td>Matlab</td><td>génie électrique</td><td>lien [24]</td></tr><tr><td>Torres</td><td>Matlab</td><td>médecine</td><td>lien [25]</td></tr><tr><td>Mandic</td><td>Matlab</td><td>traitement du signal</td><td>lien [26]</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparaison de deux versions de l'EMD</figcaption></figure>

bien mais quelques différences majeures émergent sur les bords: elles sont directement liées à la façon de les considérer dans les procédures numériques. Néanmoins, toutes ces IMF respectent les critères définis par Huang. Il est ainsi difficile de déterminer si une décomposition est *meilleure* en se limitant à l'observation des IMF. Pour porter un tel jugement, il faut exploiter ces IMF au travers d'une transformation de Hilbert afin d'obtenir leur spectre, seule quantité réellement digne d'interprétation.
## 3. Spectre de Hilbert

**Note** Dans cette section, les spectrogrammes sont représentés avec un code couleur allant du blanc au vert alors que les spectres de Hilbert sont affichés avec un code couleur allant du blanc au rouge, comme indiqué figure 10. Dans les deux cas, les couleurs claires correspondent à de faibles amplitudes.

<figure><img src="image_11.png" /><figcaption>Figure 10. Codes couleurs pour spectre de Hilbert (gauche) et spectrogramme (droite). Seuls des niveaux relatifs sont indiqués pour alléger la lecture.</figcaption></figure>

### Définition

Le spectre de Hilbert est la représentation dans un plan temps/fréquence de l'amplitude de chacune des IMF précédemment construites. Ces fonctions subissent une transformée de Hilbert d'opérateur $\mathcal{H}$ qui donne accès à leur fréquence instantanée

$$ \delta_i(t) = \mathcal{H}(\gamma_i(t)). \qquad (9) $$

Comme souligné dans [1], les IMF sont de bons candidats pour subir une transformée de Hilbert. La combinaison de $\gamma_i(t)$ et $\delta_i(t)$ permet, à chaque instant $t$ de connaître l'amplitude et la fréquence associées à une IMF. Compilées sur un même graphique pour toutes les IMF, ces informations constituent le spectre de Hilbert.

Un exemple complet est donné pour $x(t) = t^2 \sin(t^2)$ en figure 11. Seules la première IMF $\gamma_1(t)$ et la première fréquence

<figure><img src="image_12.png" /><figcaption>Figure 11. IMF $\gamma_1(t)$ et fréquence instantanée $\delta_1(t)$ de $x(t) = t^2 \sin(t^2)$ avec le spectre de Hilbert correspondant. Ce dernier prend en compte toutes les IMF, même d'amplitude négligeable, d'où les tâches jaunes claires à basse fréquence.</figcaption></figure>

instantanée $\delta_1(t)$ sont tracées sur la figure 11. Il faut bien comprendre que le spectre de Hilbert $f_{Hi}$ tient également compte de la deuxième IMF, indiquée figure 5c, ce qui explique la présence de participations d'amplitude très faible à basse fréquence. Le code couleur du spectre de Hilbert correspond à l'amplitude des IMF, du blanc (amplitude nulle) au rouge (amplitude maximale). À cause de la sensibilité numérique de la transformation de Hilbert, les fréquences instantanées $\delta_i(t)$ sont lissées en moyenne mobile, ce qui explique, sur la figure 11, pourquoi la croissance linéaire de la fréquence du signal est plus facilement identifiable sur le spectre de Hilbert.

### Comparaison aux spectrogrammes

Les spectres de Hilbert sont maintenant comparés à des spectrogrammes conventionnels. La construction de spectrogrammes requiert le calibrage de plusieurs paramètres afin d'obtenir des résultats satisfaisants.

Un compromis apparaît toutefois entre l'intervalle de temps sur lequel le spectrogramme peut fournir des résultats pertinents et sa résolution en fréquence. Ce compromis est particulièrement délicat pour les signaux considérés dans ce travail pour lesquels la période caractéristique est relativement grande en comparaison à l'intervalle de temps de définition. Deux spectrogrammes sont affichés pour chaque signal<fnref n="5" />: le premier favorise la résolution en fréquence avec un domaine de validité réduit alors que le second privilégie le domaine au détriment de la résolution en fréquence. Ces deux spectrogrammes sont significativement différents selon que la résolution en fréquence ou le domaine de validité est avantagé. Lorsque la résolution en fréquence est garantie, les résultats peuvent être très précis comme le montrent les figures 12a et 12b. Grâce aux spectrogrammes pour lesquels c'est la largeur

<footnote n="5">⁵L'étiquette $f_{Fo}$ est utilisée sur l'axe des ordonnées pour indiquer un spectrogramme alors que $f_{Hi}$ est utilisée pour un spectre de Hilbert.</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 12. Comparaison entre spectrogramme et spectre de Hilbert</figcaption></figure>

de l'intervalle de temps utilisé qui est favorisée, il est possible de voir apparaître assez précisément les variations de fréquences en fonction du temps comme sur les figures 12a et 12c. L'analyse des spectrogrammes est alors laborieuse car les variations de fréquences apparaissent comme des variations de couleurs qui sont relatives à l'amplitude du signal, il y a alors un risque de mauvaise interprétation des résultats obtenus.

Quel que soit le signal, l'allure des spectres de Hilbert est reconnaissable sur les spectrogrammes mais ces derniers ne proposent pas une résolution comparable. Aussi bien pour la détection des variations d'amplitude que pour les variations de fréquences, la résolution des spectres de Hilbert est grandement supérieure à celle des spectrogrammes. Comme mentionné précédemment, ceci est particulièrement vrai pour des signaux pour lesquels la période caractéristique est relativement longue par rapport à l'intervalle de temps du signal complet.

Toutefois, le cadre d'étude défini par de tels signaux est cohérent avec la problématique qui a motivé l'utilisation de méthodes adaptatives. En effet, pour les applications industrielles envisagées, notre intérêt se porte essentiellement sur les changements brutaux dans la réponse vibratoire. Ces changements peuvent être extrêmement localisés en temps et il est alors important, sur un intervalle de temps potentiellement petit (c'est-à-dire supérieur mais d'un ordre de grandeur comparable à la période caractéristique de la réponse analysée) de détecter des changements de fréquences ou d'amplitudes.
## 4. Décomposition en moyenne locale

De même que pour la méthode de Hilbert-Huang, la LMD est adaptative et comprend deux étapes: (1) une étape de décomposition du signal $x(t)$ en une somme de fonctions IMF suivie de (2) l'obtention du spectre de Hilbert associé. Ainsi, la LMD est très similaire à la méthode de Hilbert-Huang, elle remplace simplement la décomposition en modes empiriques par un autre type de décomposition détaillé dans cette section. Initialement introduite par Smith [15], la LMD a été appliquée, tout comme la méthode de Hilbert-Huang, dans de nombreux domaines de traitement du signal.

### Caractéristiques

L'objectif de la LMD reste une décomposition du signal telle que celle indiquée dans l'équation (1) et les IMF sont caractérisées de la même manière que pour l'EMD, les critères donnés dans la section 2.1 sont donc toujours valables pour cette méthode. La principale différence entre LMD et EMD se situe au niveau du processus de tamisage. L'interpolation par splines des minima et maxima est remplacée, à chaque itération, par le calcul d'une moyenne mobile du signal. Cette dernière permet ensuite d'obtenir les enveloppes supérieure et inférieure. Le processus de décomposition associé à la LMD est représenté schématiquement sur la figure 13.

<figure><img src="image_14.png" /><figcaption>Figure 13. Principe de la décomposition par moyenne locale</figcaption></figure>

### Construction des IMF

Les notations précédentes sont conservées et on note $r_i(t)$, le signal initial du processus de tamisage pour obtenir l'IMF $\gamma_i(t)$ ($r_i(t) = x(t)$ pour $i = 1$). Les extrema locaux de $r_i(t)$, notés $n_{i,1}$, sont déterminés. Ces extrema sont ensuite utilisés pour le calcul des moyennes locales $m_{i,1}$ et magnitudes locales $a_{i,1}$ définies par les formules [15]

$$m_{i,1} = \frac{n_{i,1} + n_{i+1,1}}{2}, \qquad (10)$$

$$a_{i,1} = \frac{|n_{i,1} - n_{i+1,1}|}{2}. \qquad (11)$$

Les moyenne et magnitude locales sont ensuite prolongées de façon à former une fonction continue en escalier $m_{esc}(t)$ et $a_{esc}(t)$. La largeur des paliers est déterminée par l'abscisse séparant deux extrema consécutifs. Ces fonctions sont lissées à l'aide d'une moyenne glissante dont la fenêtre d'étude est définie à partir de l'écart entre les deux extrema consécutifs les plus distants. Le résultat du lissage fournit la moyenne $m_{i,1}(t)$ et l'enveloppe $a_{i,1}(t)$, comme expliqué sur la figure 14. L'étape suivante consiste à soustraire au signal $x(t)$ sa moyenne $m_{i,1}(t)$ et à diviser le résultat par l'enveloppe $a_{i,1}(t)$, soit:

$$h_{i,1}(t) = r_i(t) - m_{i,1}(t) \qquad (12)$$
<figure><img src="image_15.png" /><figcaption>Figure 14. Calcul de la moyenne $m_{i,1}(t)$ et de l'enveloppe $a_{i,1}(t)$</figcaption></figure>

et:

$$s_{i,1}(t) = \frac{h_{i,1}(t)}{a_{i,1}(t)} \qquad (13)$$

où $s_{i,1}(t)$ est un signal démodulé en fréquence et en amplitude. Les étapes associées aux équations 12 et 13 sont représentées par la figure 15. Les étapes (10) à (13) sont répétées tant que le signal résultant $s_{i,k}(t)$ ne présente pas une enveloppe constante égale à 1

<figure><img src="image_16.png" /><figcaption>Figure 15. Première itération du calcul de la première IMF</figcaption></figure>

(représentée en pointillés gris sur la figure 15): $s_{i,1}(t)$ devient le signal initial pour l'itération suivante du processus de tamisage, soit $r_{i,2}(t) \leftarrow s_{i,1}(t)$. Lorsque les enveloppes du signal $s_{i,k}(t)$ sont assimilables à des droites d'équation $y = 1$ et $y = -1$, le processus de tamisage est stoppé. Par définition, $s_{i,k}(t) \in [-1; 1]$, $\forall t$ et on peut donc écrire:

$$s_{i,k}(t) = \cos(\phi_i(t)) \qquad (14)$$

c'est-à-dire

$$\phi_i(t) = \arccos(s_{i,k}(t)) \qquad (15)$$

où $\phi_i(t)$ est la phase instantanée du signal analytique $s_{i,n}(t)$. La fréquence instantanée se déduit immédiatement de l'expression (15):

$$\omega_i(t) = \frac{d\phi_i(t)}{dt} \qquad (16)$$

Enfin, l'IMF $\gamma_i(t) = s_{i,k}(t) \prod_{q=1}^{k} a_{i,q}(t)$ est soustraite à $r_i(t)$:

$$r_{i+1,1}(t) = r_i(t) - \gamma_i(t) \qquad (17)$$

et la procédure est itérée pour le calcul de $\gamma_{i+1}(t)$ jusqu'à ce que le signal résiduel soit monotone ou non oscillant.

## Application aux signaux académiques

Les cinq réponses précédentes décomposées par EMD, section 2.4.1, sont maintenant décomposées par LMD. Les résultats, indiqués sur les figures 16a, 16b, 16c, 16d et 17 sont à comparer respectivement aux figures 5a, 5c, 5b, 6a et 6b. Globalement, les IMF obtenues par LMD et EMD sont très similaires pour tous les signaux académiques considérés.
<figure><img src="image_17.png" /></figure>

$$ (a) \sin(2\pi t) + \cos(6\pi t) + \sin(2\pi t)^2 \sin(6\pi t) $$

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

$$ (c) \cos(t^2) + 3 \cos(3t) + \cos(t) $$

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

$$ (b) t^2 \sin(t^2) $$

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

$$ (d) \left( \exp(-0.2t) t^2 + 10 \right) \sin\left(\exp(-0.1t) t^2 + 3t\right) + \cos(2t) + \sin(150t) + \cos(74t) $$

<caption>Figure 16. Décomposition en moyenne locale des signaux $x(t)$ considérés pour l'EMD</caption>

## 5. Conclusion

Dans cet article, deux méthodes adaptatives de traitement du signal, la décomposition en modes empiriques d'une part et la décomposition en moyenne locale d'autre part, sont présentées en détails. L'implémentation proposée (disponible en annexe en langage Python) est validée sur plusieurs signaux académiques et est également confrontée à d'autres codes disponibles sur internet. Les résultats obtenus soulignent la similarité des décompositions en fonctions IMF par LMD ou EMD. Les différences se situent essentiellement au niveau de leur comportement numérique. Il a en effet été observé que la LMD converge plus difficilement que l'EMD. La nécessité d'obtenir une fonction comprise entre -1 et +1 au sein du processus de tamisage est une opération délicate qui peut ralentir sensiblement le temps requis pour la construction des IMF. Une plus grande sensibilité aux effets de bord a également été observée pour la LMD. Pour résumer, la rapidité de la décomposition par EMD et la qualité des résultats en font une méthode robuste et facilement applicable à différents domaines d'étude.

## References

[1] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, H. H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal
<figure><img src="image_21.png" /><figcaption>Figure 17. Décomposition en moyenne locale de $x(t) = (\exp(-0,2t)t^2 + 10)(\sin(\exp(-0,1t)t^2 + 3t) + \cos(2t))$</figcaption></figure>

Society of London A: Mathematical, Physical and Engineering Sciences 454 (1971) (1998) 903-995. doi:10.1098/rspa.1998.0193.

[2] G. Rilling, P. Flandrin, P. Gonçalvès, On empirical mode decomposition and its algorithms, in: Proceedings of the 6th IEEE/EURASIP Workshop on Nonlinear Signal and Image Processing, Grado, Italy, 2003.

[3] G. Rilling, P. Flandrin, P. Goncalves, Empirical mode decomposition as a filter bank, IEEE Transactions on signal processing 11 (2004) 112-114. doi:10.1109/LSP.2003.821662.

[4] G. Rilling, P. Flandrin, On the influence of sampling on the empirical mode decomposition, in: IEEE International Conference on Acoustics Speech and Signal Processing Proceedings, Vol. 3, 2006, pp. 444-447. doi:10.1109/ICASSP.2006.1660686.

[5] G. Rilling, Décompositions modales empiriques. contributions à la théorie, l'algorithmie et l'analyse de performances, Ph.D. thesis, Ecole Normale Supérieure de Lyon (2007). doi:tel-00442634.

[6] G. Rilling, P. Flandrin, One or two frequencies? the empirical mode decomposition answers, IEEE Transactions on signal processing 56 (2008) 85-95. doi:10.1109/TSP.2007.906771.

[7] P. Goncalves, P. Abry, G. Rilling, P. Flandrin, Fractal dimension estimation: empirical mode decomposition versus wavelets, IEEE Acoustics, Speech and Signal Processing 3 (2007) 1153-1156. doi:10.1109/ICASSP.2007.366889.

[8] Y. Zhang, C. Bingham, Z. Yang, B. W.-K. Ling, M. Gallimore, Machine fault detection by signal denoising—with application to industrial gas turbines, Measurement 58 (2014) 230-240. doi:10.1016/j.measurement.2014.08.020.

[9] A. Batailly, M. Legrand, A. Millecamps, F. Garcin, Conjectural bifurcation analysis of the contact-induced vibratory response of an aircraft engine blade, Journal of Sound and Vibration 348 (2015) 239-262. doi:10.1016/j.jsv.2015.03.005.

[10] T. Xiong, Y. Bao, Z. Hu, Interval forecasting of electricity demand: A novel bivariate EMD-based support vector regression modeling framework, International Journal of Electrical Power & Energy Systems 63 (2014) 353-362. doi:10.1016/j.ijepes.2014.06.010.

[11] Z. Shen, X. Chen, X. Zhang, Z. He, A novel intelligent gear fault diagnosis model based on EMD and multi-class TSVM, Measurement 45 (1) (2012) 30-40. doi:10.1016/j.measurement.2011.10.008.

[12] Z. Yunlong, Z. Peng, Vibration fault diagnosis method of centrifugal pump based on EMD complexity feature and least square support vector machine, Energy Procedia 17 (Part A) (2012) 939-945. doi:10.1016/j.egypro.2012.02.191.

[13] F. Wu, L. Qu, Diagnosis of subharmonic faults of large rotating machinery based on EMD, Mechanical Systems and Signal Processing 23 (2) (2009) 467-475. doi:10.1016/j.ymssp.2008.03.007.

[14] Z. Peng, P. W. Tse, F. Chu, An improved Hilbert-Huang transform and its application in vibration signal analysis, Journal of Sound and Vibration 286 (1-2) (2005) 187-205. doi:10.1016/j.jsv.2004.10.005.

[15] J. S. Smith, The local mean decomposition and its application to EEG perception data, Journal of The Royal Society Interface 2 (5) (2005) 443-454. doi:10.1098/rsif.2005.0058.

[16] Y. Yang, J. Cheng, K. Zhang, An ensemble local means decomposition method and its application to local rub-impact fault diagnosis of the rotor systems, Measurement 45 (3) (2012) 561-570. doi:10.1016/j.measurement.2011.10.010.
[17] B. Balachandran, T. Kalmar-Nagy, D. Gilsinn, Delay differential equations: Recent advances and new directions, Springer Verlag, 2009.

[18] J. Cheng, D. Yu, Y. Yang, Application of support vector regression machines to the processing of end effects of Hilbert-Huang transform, Mechanical Systems and Signal Processing 21 (3) (2007) 1197-1211. doi:10.1016/j.ymssp.2005.09.005.

[19] W. Wang, X. Li, R. Zhang, Boundary processing of HHT using support vector regression machines, in: Computational Science - ICCS 2007, Vol. 4489 of Lecture Notes in Computer Science, Springer Berlin Heidelberg, 2007, pp. 174-177. doi:10.1007/978-3-540-72588-6_31.

[20] Y. Deng, W. Wang, C. Qian, Z. Wang, D. Dai, Boundary-processing-technique in EMD method and Hilbert transform, Chinese Science Bulletin 46 (11) (2001) 954-960. doi:10.1007/BF02900475.

[21] Y. Ye, Adaptive boundary effect processing for empirical mode decomposition using template matching, Applied Mathematics & Information Sciences 7 (2013) 61-66.

[22] G. Rilling, P. Flandrin, http://perso.ens-lyon.fr/patrick.flandrin/emd.html.

[23] D. Laszuk, http://laszukdawid.com/codes/.

[24] M. Lambert, A. Engroff, M. Dyer, B. Byer, https://www.clear.rice.edu/elec301/Projects02/empiricalMode/code.html.

[25] M. Torres, M. Colominas, G. Schlotthauer, P. Flandrin, http://bioingenieria.edu.ar/grupos/ldnlys/metorres/re_inter.htm.

[26] D. Mandic, http://www.commsp.ee.ic.ac.uk/~mandic/research/emd.htm.