Méthode de Welch

En estimation spectrale, la méthode de Welch fournit un estimateur consistant de la densité spectrale de puissance. Cette méthode a été proposée par Peter D. Welch en 1967. Le biais de l'estimation est diminué en moyennant temporellement. Elle est à comparer à la méthode de Bartlett où on utilise les propriétés d'ergodicité du signal avec des moyennes statistiques. La méthode de Welch, comme la méthode de Bartlett, utilise une estimation du spectre du périodogramme ; dans les deux cas, on réduit le bruit aux dépens de la résolution en fréquence.

Principe de la méthode

La méthode de Welch s’apparente à la méthode de Bartlett, mais avec deux différences :

1. L’intervalle des M points du signal est découpé en segments qui se recouvrent partiellement : le segment original est divisé en K segments de longueur N, qui sont décalés deux à deux de D points. Par exemple :
   1. Si D = N / 3, le taux de recouvrement est de 66 %.
   2. Si D = N, le taux de recouvrement est nul, et cette situation correspond à la méthode de Bartlett.
2. Chacun des K segments est ensuite « fenêtré », ce qui consiste à introduire une pondération temporelle des données. Le fenêtrage permet d’abord d’atténuer le phénomène de Gibbs.
   1. Les fonctions de fenêtrage accordent généralement plus d’importance aux données du centre du segment qu’à celles du bord, ce qui entraîne une perte d'information. Le recouvrement des segments permet de réduire cet effet.
   2. Le fenêtrage de la méthode de Welch génère ce qui est appelé un périodogramme modifié.

Une fois les données préparées de cette manière, le périodogramme de chaque segment est défini par les carrés des modules des coefficients de la transformée de Fourier discrète, puis le résultat final est la moyenne arithmétique des périodogrammes individuels. Ce procédé permet de réduire la variance des puissances individuelles. On obtient finalement des puissances en fonction des fréquences.

Algorithme

1. Découpe du signal en segments se chevauchant partiellement.
2. Fenêtrage de chaque segment.
3. Transformée de Fourier de chaque segment.
4. Moyenne des résultats de chaque segment.

Formulation Mathématique

Un signal de longueur M est découpé en K segments de longueur N et décalés l’un par rapport à l’autre de D points. Le k-ième segments s'écrira $x_{k}(n)=x(n+(k-1)\mathrm {D} )$. Ainsi, la longueur du recouvrement sera de $\max(\mathrm {N} -\mathrm {D} ,0)$.

La densité spectrale de puissance estimée s'écrit :

${\hat {\mathrm {S} }}_{\mathrm {Welch} }(\omega )={\frac {1}{\mathrm {K} \mathrm {N} \mathrm {U} }}\sum _{k=1}^{\mathrm {K} }{\left|\sum _{n=1}^{\mathrm {N} }h(n)x(n+(k-1)\mathrm {D} )\exp(-j\omega n)\right|^{2}}$

• h(n) une fenêtre de pondération de longueur N ;
• U est la constante de normalisation des pondérations assurant que ${\hat {\mathrm {S} }}_{\mathrm {W} }(\omega )$ soit asymptotiquement non biaisé. Elle est définie par $\mathrm {U} ={\frac {1}{\mathrm {N} }}\sum _{n=1}^{\mathrm {N} }{|h(n)|}^{2}$

Biais

En partant du périodogramme modifié :

${\hat {\mathrm {S} }}_{\mathrm {Welch} }(\omega )={\frac {1}{\mathrm {K} }}\sum _{k=1}^{\mathrm {K} }{\hat {\mathrm {S} }}_{\mathrm {permod} }^{(k)}(\omega )$

On peut montrer que le biais s'écrit :

$\mathrm {E} \{{\hat {\mathrm {S} }}_{\mathrm {Welch} }(\omega )=\mathrm {E} \{{\hat {\mathrm {S} }}_{\mathrm {permod} }(\omega )\}={\frac {1}{2\pi \mathrm {N} \mathrm {U} }}\mathrm {S} _{x}(\omega )*{|\mathrm {H} (\omega )|}^{2}$

où H(ω) est la transformée de Fourrier de la fenêtre h(n).

Comme le périodogramme modifié est non biaisé, la méthode de Welch est non biaisée.

Problèmes de stationnarité

Puisque chaque segment est moyenné, cet algorithme part de l’hypothèse que les caractéristiques spectrales n'évoluent pas dans le temps, autrement dit que le signal est stationnaire.

Mise en œuvre

Le choix des paramètres de la méthode influence le résultat. Pour une première approche, lorsque l'on ne connaît pas les caractéristiques du signal, on peut réaliser une analyse avec les paramètres suivants en partant du nombre M de points du signal :

• La taille N des segments est choisie comme la puissance de 2 qui est juste inférieure à M / 2,
  N = 2ⁿ avec n ∈ ℕ et 2 N ≤ M < 4 N ; ce choix permet d’exploiter la FFT dans sa version simple.
• Le nombre de segments K est fixé en fonction d’un taux de recouvrement R proche de 50% et de manière à maximiser la taille de la réunion des segments.
• Le spectre est évalué K fois sur N points.
• Le fenêtrage se base sur la fenêtre de Hamming.

Dans la pratique :

• Taille des segments : N = 2ⁿ avec $n=\lfloor \ln(\mathrm {M} )/\ln(2)\rfloor -1$ où $\lfloor .\rfloor$ désigne la partie entière.
• Nombre de segments : $\mathrm {K} =\left\lfloor {\frac {2\mathrm {M} }{\mathrm {N} }}\right\rfloor -1$ qui peut varier entre 3 et 6.
• Décalage : $\mathrm {D} =\left\lfloor {\frac {\mathrm {M} -\mathrm {N} }{\mathrm {K} -1}}\right\rfloor$.

La réunion des segments touche N + (K – 1)D points qui est parfois inférieur à M. Cependant, ainsi définis et puisque K ≤ 6, ce sont au plus 4 points du signal d’origine qui seront ignorés.

La méthode de Welch est intégrée :

• dans la toolbox Signal Processing de Matlab, avec la fonction pwelch() ;
• dans Scilab, avec la fonction pspect().
• dans Python, avec la fonction scipy.signal.welch().

Bibliographie

• (en) Alan V. Oppenheim, Ronald W. Schafer et J. R. Buck, Discrete-time signal processing, Upper Saddle River, N.J., Prentice Hall, 1999 (ISBN 0-13-754920-2).

Voir aussi

• Fixed point per octave (FPPO)