Algorithme de Savitzky-Golay

Un article de Wikipédia, l'encyclopédie libre.
Aller à : navigation, rechercher

L'algorithme de Savitzky-Golay est une méthode utilisée en traitement du signal pour lisser une courbe et en extraire les dérivées successives. Il a été décrit en 1964 par Abraham Savitzky et Marcel Golay[1].

Description de l'algorithme[modifier | modifier le code]

Algorithme de Savitzky-Golay appliqué avec un polynôme de degré 3 et une largeur de 9 points
Lissage point par point (mêmes conditions que l'image précédente).

Considérons une courbe y = ƒ(x), et présentant des « aspérités », des oscillations de faibles amplitudes ; on parle de signal bruité. Il s'agit d'une courbe discrète, c'est-à-dire définie par un nuage de points (x(i), y(i))1 ≤ in.

L'algorithme de lissage le plus simple est la méthode des moyennes glissantes :

  • on considère une fenêtre, un intervalle, de « demie largeur » l (nombre de points) ;
  • on calcule la moyenne yl + 1 de la fonction sur l'intervalle [1 ; 2×l + 1] (l'intervalle a donc une largeur de 2l + 1, l n'est pas exactement la demie largeur) ;
  • on définit un nouveau point (x(l + 1), yliss(l + 1)) avec yliss(l + 1) = yl + 1 ;
  • on fait glisser l'intervalle d'un point et l'on recommence.

Cet algorithme est simple, mais « violent » : il atténue énormément les fortes variations, il écrête les pics.

Pour la suite, nous désignons l'intervalle i comme étant [i - l ; i + l], l'intervalle de milieu i et de largeur 2l + 1.

Une manière plus fine consiste à considérer un polynôme de degré d, avec d < 2l + 1. Pour chaque intervalle i, on effectue une régression pour déterminer le polynôme Pi minimisant l'erreur au sens des moindres carrés. On définit la valeur lissée

yliss = Pi(x(i)).

Par ailleurs, si le polynôme est au moins de degré 1, on peut déterminer la dérivée

y' = P'i(x(i))

et de manière générale, on peut déterminer la d-ième dérivée.

On voit que la méthode des moyennes glissantes est un lissage de Savitzky-Golay avec un polynôme de degré 0.

Choix des paramètres[modifier | modifier le code]

Dans la pratique :

  • un polynôme de degré 2 permet de prendre en compte la courbure ; un polynôme de degré 3 permet de prendre en compte des points d'inflexion ; il est rarement nécessaire d'aller au-delà ;
  • le nombre de points d'un intervalle doit être suffisamment grand devant le degré du polynôme pour que le lissage soit effectif ; s'ils sont égaux, le polynôme passe exactement par tous les points de l'intervalle, il n'y a donc pas de lissage.

On prend donc en général un polynôme de degré 3 et une fenêtre glissante d'au moins 5 points.

Plus la fenêtre est large, plus la courbe est lissée, mais plus on atténue les variations de grande longueur d'onde. La largeur limite est de l'ordre de la longueur d'onde des détails que l'on veut conserver. Si par exemple on s'intéresse à des pics — cas très fréquents en spectrométrie et en diffractométrie —, on prend en général pour largeur de fenêtre la largeur à mi-hauteur du pic le plus fin.

Cas d'un pas constant[modifier | modifier le code]

Concrètement, les coefficients de la régression polynomiale sont calculés par une combinaison linéaire des points de l'intervalle glissant. La valeur lissée, et les valeurs des dérivées, sont donc elles-mêmes des combinaisons linéaires :

y_{\mathrm{liss}\ i} = b_0 \cdot y_{i - l} + b_1 \cdot y_{i - l + 1} + \ldots + b_{l} \cdot y_i + \ldots +  b_{2l} \cdot y_{i + l} = \sum_{k = 0}^{2l} b_k \cdot y_{i - l + k} ;
y'_{\mathrm{liss}\ i} = b'_0 \cdot y_{i - l} + b'_1 \cdot y_{i - l + 1} + \ldots + b'_{l} \cdot y_i + \ldots +  b'_{2l} \cdot y_{i + l}  = \sum_{k = 0}^{2l} b'_k \cdot y_{i - l + k} ;

Si le pas h = xi - xi - 1 est constant[2], les coefficients bk, b'k, sont les mêmes partout ; on se retrouve dans le cas d'un filtre à réponse impulsionnelle finie (filtre RIF). On peut donc déterminer au départ les coefficients bk, b'k, …, ce qui fournit un algorithme rapide et facile à mettre en œuvre — par exemple avec un tableur ; on parle alors de coefficients de convolution.

Dans ces conditions, on trouve que :

  • le lissage est identique que l'on utilise un polynôme de degré 2 (quadratique) ou 3 (cubique) ;
  • la dérivation est identique que l'on utilise un polynôme de degré 3 (cubique) ou 4 (« quartique ») ;
  • la dérivée seconde est identique que l'on utilise un polynôme de degré 4 (« quatrique ») ou 5 (« quintique ») ;
  • les coefficients sont symétriques ou anti-symétriques, selon l'ordre de dérivation :
bk = b2l - k ; b'k = -b'2l - k ; …
  • on a toujours b'l = 0 : le point xi ne contribue pas au calcul de y'liss i.

Les valeurs des coefficients sont tabulées pour des tailles de fenêtre (intervalle glissant) données. Par exemple, pour une fenêtre de 5 point et un polynôme de degré 3 :

\textstyle y_{\mathrm{liss}\ i} = \frac{1}{35} (-3 \cdot y_{i - 2} + 12 \cdot y_{i - 1} + 17 \cdot y_i + 12 \cdot y_{i + 1} -3 \cdot y_{i + 2}) ;
\textstyle y'_{\mathrm{liss}\ i} = \frac{1}{10h} (-2 \cdot y_{i - 2} - 1  \cdot y_{i - 1} + 0 \cdot y_i + 1 \cdot y_{i + 1} + 2 \cdot y_{i + 2}) ;
\textstyle y^{\prime \prime}_{\mathrm{liss}\ i} = \frac{1}{7h^2} (2 \cdot y_{i - 2} - 1  \cdot y_{i - 1} - 2 \cdot y_i - 1 \cdot y_{i + 1} + 2 \cdot y_{i + 2}).

Les tableaux ci-dessous donnent les coefficients de convolution pour des degrés de polynôme et des largeurs de fenêtre données.

Il ne faut pas oublier d'appliquer un coefficient 1/h pour la dérivée et 1/h2 pour une dérivée seconde.
Coefficients de convolution pour le lissage
Degré 2/3 (quadratique/cubique) 4/5 (quartique/quintique)
Largeur 5 7 9 7 9
-4 -21 15
-3 -2 14 5 -55
-2 -3 3 39 -30 30
-1 12 6 54 75 135
0 17 7 59 131 179
1 12 6 54 75 135
2 -3 3 39 -30 30
3 -2 14 5 -55
4 -21 15
Normalisation 35 21 231 231 429
Coefficients de convolution pour la dérivée
Degré 2 (quadratique) 3/4 (cubique/quartique)
Largeur 5 7 9 5 7 9
-4 -4 86
-3 -3 -3 22 -142
-2 -2 -2 -2 1 -67 -193
-1 -1 -1 -1 -8 -58 -126
0 0 0 0 0 0 0
1 1 1 1 8 58 126
2 2 2 2 -1 67 193
3 3 3 -22 142
4 4 -86
Normalisation 10 28 60 12 252 1 188
Coefficients de convolution pour la dérivée seconde
Degré 2/3 (quadratique/cubique) 4/5 (quartique/quintique)
Largeur 5 7 9 5 7 9
-4 28 -4 158
-3 5 7 -117 12 243
-2 2 0 -8 -3 603 4 983
-1 -1 -3 -17 48 -171 -6 963
0 -2 -4 -20 -90 -630 -12 210
1 -1 -3 -17 48 -171 -6 963
2 2 0 -8 -3 603 4 983
3 5 7 -117 12 243
4 28 -4 158
Normalisation 7 42 462 36 1 188 56 628

Calcul des coefficients de convolution[modifier | modifier le code]

On se place en un point k donné. Les coordonnées du point expériemntal sont (xk, yk).

Pour déterminer les coefficients de corrélation, on effectue un changement de variable afin d'avoir une abscisse centrée réduite — c'est-à-dire que l'intervalle sur lequel on effectue la régression devient [-l ; l] : on pose

z = \frac{x - \bar{x}}{h}

  • x est la moyenne des valeurs de x sur l'intervalle [xk - l ; xk + l] considéré, c'est le centre de la fenêtre ;
  • h est le pas d'échantillonnage (distance entre deux points voisins).

L'abscisse z prend les valeurs (-l ; -l + 1 ; … ; -1 ; 0 ; 1 ; … l). Par exemple, pour une fenêtre glissante de 5 points, on a (zi)1 ≤ i ≤ 5 = (-2 ; -1 ; 0 ; 1 ; 2).

Le polynôme s'exprime alors sous la forme :

y = \mathrm{P}(z) = a_0 + a_1 \cdot z + a_2 \cdot z^2 + \ldots

Considérons la fonction vectorielle

\mathrm{F} : \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_{2l + 1} \end{pmatrix}
\longmapsto \begin{pmatrix} \mathrm{P}(z_1) \\ \mathrm{P}(z_2) \\
\vdots \\ \mathrm{P}(z_{2l + 1}) \end{pmatrix}

avec, rappelons-le,

\begin{pmatrix} \mathrm{P}(z_1) \\ \mathrm{P}(z_2) \\
\vdots \\ \mathrm{P}(z_{2l + 1}) \end{pmatrix}
\simeq \begin{pmatrix} y_{k-l} \\ y_{k - l + 1} \\ \vdots \\ y_{k + l} \end{pmatrix}

Sa matrice jacobienne s'écrit

\mathrm{J} = \left ( \frac{\partial \mathrm{P}}{\partial a_j}(z_i) \right )
= \begin{pmatrix} 1 & z_1 & z_1^2 & \cdots & z_1^d \\
1 & z_2 & z_2^2 & \cdots & z_2^d \\
\vdots & \vdots & \vdots & & \vdots \\
1 & z_{2l + 1} & z_{2l + 1}^2 & \cdots & z_{2l + 1}^d
\end{pmatrix}

La solution des équations normales

{}^\mathrm{t}\mathrm{J} \mathrm{J} (a_i) = {}^\mathrm{t}\mathrm{J}(y_j)

est donc

(a_i) = \left ( {}^\mathrm{t}\mathrm{J} \mathrm{J} \right )^{-1}{}^\mathrm{t}\mathrm{J}(y_j)

La valeur de la courbe lissée est la valeur du polynôme au centre de l'intervalle :

yliss = P(0) = a0 ;

Pour les dérivées successives :

y'_\mathrm{liss} = \frac{\partial \mathrm{P}}{\partial x}(x_k) = \frac{\partial z}{\partial x}(0)\frac{\partial \mathrm{P}}{\partial z}(0) = \frac{1}{h}a_1 ;
y^{\prime\prime}_\mathrm{liss} = \frac{\partial^2 \mathrm{P}}{\partial x^2}(x_k) = \frac{\partial^2 z}{\partial x^2}(0)\frac{\partial \mathrm{P}}{\partial z}(0) + \frac{\partial z}{\partial x}\frac{\partial^2 \mathrm{P}}{\partial z^2}(0) = \frac{2}{h^2}a_2.

Par exemple, pour un polynôme de degré d = 3 et une fenêtre de 5 points, donc l = 2 :

\mathrm{J} = \begin{pmatrix} 1 & -2 & 4 & -8 \\
1 & -1 & 1 & -1 \\
1 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 \\
1 & 2 & 4 & 8
\end{pmatrix} \text{ ; }
{}^\mathrm{t}\mathrm{J} = \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\
-2 & -1 & 0 & 1 & -1 \\
4 & 1 & 0 & 1 & 4 \\
-8 & -1 & 0 & 1 & 8
\end{pmatrix}
\begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} =
\left ( {}^\mathrm{t}\mathrm{J} \mathrm{J} \right )^{-1}{}^\mathrm{t}\mathrm{J}
\begin{pmatrix}
y_{k-2} \\ y_{k-1} \\ y_k \\ y_{k+1} \\ y_{k+2}
\end{pmatrix}
= \begin{pmatrix} - 3/35 & 12/35 & 17/35 & 12/35 & -3/35 \\
1/12 & -8/12 & 0 & 8/12 & -1/12 \\  
2/14 & -1/14 & -2/14 & -1/14 & 2/14 \\ 
-1/12 & 2/12 & 0 & -2/12 & 1/12
\end{pmatrix}
\begin{pmatrix}
y_{k-2} \\ y_{k-1} \\ y_k \\ y_{k+1} \\ y_{k+2}
\end{pmatrix}.

On a donc :

y_{\mathrm{liss} \ k} = a_0 = \frac{1}{35}(-3y_{k-2} + 12y_{k-1} + 17y_k + 12y_{k+1} - 3y_{k+2}) ;
y'_{\mathrm{liss} \ k} = a_1/h = \frac{1}{12h}(y_{k-2} - 8y_{k-1} + 0y_k + 8_{k+1} - y_{k+2}) ;
y^{\prime\prime}_{\mathrm{liss} \ k} = 2a_2/h^2 = \frac{1}{7h^2}(2y_{k-2} - y_{k-1} -2y_k -y_{k+1} + 2y_{k+2}).

Notes et références[modifier | modifier le code]

  1. Abraham Savitzky et Marcel J. E. Golay, « Smoothing and Differentiation of Data by Simplified Least Squares Procedures », Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI 10.1021/ac60214a047)
  2. on n'a pas toujours un pas constant, on peut avoir un pas variable
    • pour des raisons d'optimisation : on veut mesurer plus finement une zone et plus grossièrement une autre (pour gagner du temps) ;
    • pour des raisons matérielles :
      • le système générant la valeur d'entrée x n'a pas une réponse linéaire,
      • pour des mesures très précises, la boucle de régulation peut prendre beaucoup de temps pour stabiliser la valeur d'entrée à la valeur de consigne xth, on peut alors gagner du temps en mesurant un point alors que la consigne n'est pas totalement atteinte, la valeur de x étant relevée de manière très précise.