# Algorithme de Savitzky-Golay

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 (en) et Marcel Golay.

## Description de l'algorithme

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≤i≤n}$.
L'algorithme de lissage le plus simple est la méthode des moyennes glissantes :
* on considère une fenêtre, un intervalle, de « demi-largeur » ℓ (nombre de points) ;
* on calcule la moyenne y$_{ℓ+1}$ de la fonction sur l'intervalle [1 ; 2 × ℓ + 1] (l'intervalle a donc une largeur de 2ℓ + 1, ℓ n'est pas exactement la demi-largeur) ;
* on définit un nouveau point (x(ℓ + 1), y$_{liss}$(ℓ + 1)) avec y$_{liss}$(ℓ + 1) = y$_{ℓ+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 – ℓ ; i + ℓ ], l'intervalle de milieu i et de largeur 2ℓ + 1.
Une manière plus fine consiste à considérer un polynôme de degré d, avec d < 2ℓ + 1. Pour chaque intervalle i, on effectue une régression pour déterminer le polynôme P$_i$ minimisant l'erreur au sens des moindres carrés. On définit la valeur lissée
    $y_{\textrm {liss}}=P_{i}(x(i)).$
Par ailleurs, si le polynôme est au moins de degré 1, on peut déterminer la dérivée
    $y'_{\textrm {liss}}=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

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 petite 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

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 = x$_i$ – x$_{i-1}$ est constant, les coefficients b$_k$, 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 b$_k$, 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 (« quartique ») ou 5 (« quintique ») ;
* les coefficients sont symétriques ou anti-symétriques, selon l'ordre de dérivation :
    b$_k$ = b$_{2ℓ–k}$ ; b'$_k$ = –b'$_{2ℓ–k}$ ; …
* on a toujours b'$_ℓ$ = 0 : le point y$_i$ ne contribue pas au calcul de y'$_{lissi}$.
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 points et un polynôme de degré 3 :
    $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})$ ;
    $y'_{\mathrm {liss} \ i}={\frac {1}{12h}}(1\cdot y_{i-2}-8\cdot y_{i-1}+0\cdot y_{i}+8\cdot y_{i+1}-1\cdot y_{i+2})$ ;
    $y_{\mathrm {liss} \ i}^{\prime \prime }={\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.
Coefficients de convolution pour le lissage :
| Degré | 2/3 (quadratique/cubique) | 2/3 (quadratique/cubique) | 2/3 (quadratique/cubique) | 4/5 (quartique/quintique) | 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) | 2 (quadratique) | 2 (quadratique) | 3/4 (cubique/quartique) | 3/4 (cubique/quartique) | 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 | 10h | 28h | 60h | 12h | 252h | 1 188h |
Remarque : la normalisation utilise le pas d'échantillonnage h.
Coefficients de convolution pour la dérivée seconde :
| Degré | 2/3 (quadratique/cubique) | 2/3 (quadratique/cubique) | 2/3 (quadratique/cubique) | 4/5 (quartique/quintique) | 4/5 (quartique/quintique) | 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 | 7h$^2$ | 42h$^2$ | 462h$^2$ | 36h$^2$ | 1 188h$^2$ | 56 628h$^2$ |

## Calcul des coefficients de convolution

On se place en un point k donné. Les coordonnées du point expérimental sont (x$_k$, y$_k$).
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 [–ℓ ; ℓ ] : on pose
    $z={\frac {x-{\bar {x}}}{h}}$
où
* x est la moyenne des valeurs de x sur l'intervalle [x$_{k–ℓ}$ ; x$_{k+ℓ}$] 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 (–ℓ ; –ℓ + 1 ; … ; –1 ; 0 ; 1 ; … ℓ ). Par exemple, pour une fenêtre glissante de 5 points, on a (z$_i$)$_{1≤i≤5}$ = (–2 ; –1 ; 0 ; 1 ; 2).
Le polynôme de degré d s'exprime alors sous la forme :
    $y=P(z)=a_{0}+a_{1}\cdot z+a_{2}\cdot z^{2}+\ldots +a_{d}\cdot z^{d}$
Considérons la fonction vectorielle
    $\mathrm {F} :{\begin{pmatrix}a_{1}\\a_{2}\\\vdots \\a_{d}\end{pmatrix}}\longmapsto {\begin{pmatrix}\mathrm {P} (z_{1})\\\mathrm {P} (z_{2})\\\vdots \\P(z_{2l+1})\end{pmatrix}}$
avec
    ${\begin{pmatrix}P(z_{1})\\P(z_{2})\\\vdots \\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 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 :
    $y_{\textrm {liss}}=P(0)=a_{0}$ ;
Pour les dérivées successives :
    $y'_{\mathrm {liss} }={\frac {\partial P}{\partial x}}(x_{k})={\frac {\partial z}{\partial x}}(0){\frac {\partial P}{\partial z}}(0)={\frac {1}{h}}a_{1}$ ;
    $y_{\mathrm {liss} }^{\prime \prime }={\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}})^{2}{\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 ℓ = 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&2\\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}={\frac {a_{1}}{h}}={\frac {1}{12h}}(y_{k-2}-8y_{k-1}+0y_{k}+8_{k+1}-y_{k+2})$ ;
    $y_{\mathrm {liss} \ k}^{\prime \prime }={\frac {2}{h^{2}}}a_{2}={\frac {1}{7h^{2}}}(2y_{k-2}-y_{k-1}-2y_{k}-y_{k+1}+2y_{k+2})$.