Savitzky-Golayfilter

Een Savitzky-Golayfilter is een wiskundig filter uit de signaalverwerking voor toepassing op equidistante getalsmatige data. De methode is gebaseerd op bestaande wiskundige technieken en werd voor het eerst in 1964 beschreven door Abraham Savitzky en Marcel J. E. Golay.^[1] Het filter maakt het oorspronkelijke signaal "gladder", dat wil zeggen vlakt snelle wisselingen af, waardoor ruis wordt gefilterd. Tevens kan het filter gebruikt worden om van het gefilterde signaal afgeleiden te bepalen. De afvlakking gebeurt door in het midden van een venster de gefilterde waarde te bepalen als gewogen som van de oorspronkelijke data in het venster. De wegingsfactoren worden bepaald met een polynoom die in een omgeving van het punt goed bij de data past.

Principe

Een signaal is op equidistante afstanden gegeven door de rij:

${\displaystyle \ldots ,x_{-2},x_{-1},x_{0},x_{1},x_{2},x_{3},\ldots }$

Rondom het ${\displaystyle k}$-de punt in de rij wordt een venster

${\displaystyle m=-n,\ldots ,-2,-1,0,1,2,\ldots ,n}$

beschouwd, dat dus zicht biedt op de data

${\displaystyle x_{k-n},\ldots ,x_{k-2},x_{k-1},x_{k},x_{k+1},x_{k+2},\ldots ,x_{k+n}}$

Op het venster wordt een polynoom van graad ${\displaystyle r}$ gedefinieerd door

${\displaystyle P_{k}(m)=a_{k0}+a_{k1}m+a_{k2}m^{2}+\ldots +a_{kr}m^{r}}$,

die zo goed mogelijk bij de data in het venster past in de zin van kleinste kwadraten, dus zo dat

${\displaystyle \sum _{m}\left(P_{k}(m)-x_{k+m}\right)^{2}}$ minimaal is.

De kleinste-kwadratenoplossing ${\displaystyle a_{k}=(a_{k0},a_{k1},a_{k2},\ldots ,a_{kr})^{\top }}$ voldoet aan de normaalvergelijkingen:

${\displaystyle {\begin{alignedat}{7}a_{k0}\sum _{m}m^{0}&+&a_{k1}\sum _{m}m^{1}&+&a_{k2}\sum _{m}m^{2}&+\ldots +&a_{kr}\sum _{m}m^{r+0}=\sum _{m}m^{0}x_{k-m}\\a_{k0}\sum _{m}m^{1}&+&a_{k1}\sum _{m}m^{2}&+&a_{k2}\sum _{m}m^{3}&+\ldots +&a_{kr}\sum _{m}m^{r+1}=\sum _{m}m^{1}x_{k-m}\\a_{k0}\sum _{m}m^{2}&+&a_{k1}\sum _{m}m^{3}&+&a_{k2}\sum _{m}m^{4}&+\ldots +&a_{kr}\sum _{m}m^{r+2}=\sum _{m}m^{2}x_{k-m}\\a_{k0}\sum _{m}m^{3}&+&a_{k1}\sum _{m}m^{4}&+&a_{k2}\sum _{m}m^{5}&+\ldots +&a_{kr}\sum _{m}m^{r+3}=\sum _{m}m^{3}x_{k-m}\\\ldots \\a_{k0}\sum _{m}m^{r}&+&a_{k1}\sum _{m}m^{r+1}&+&a_{k2}\sum _{m}m^{r+2}&+\ldots +&a_{kr}\sum _{m}m^{r+r}=\sum _{m}m^{r}x_{k-m}\\\end{alignedat}}}$

Deze vergelijkingen hebben de vorm:

${\displaystyle X^{\top }Xa_{k}=X^{\top }b_{k}}$,

met als oplossing

${\displaystyle a_{k}=(X^{\top }X)^{-1}X^{\top }b_{k}=Wb_{k}}$,

waarin de ${\displaystyle (2n+1)\times (r+1)}$-matrix ${\displaystyle X=[X_{0},X_{1},\ldots ,X_{r}]}$ vanwege de equidistante afstanden een eenvoudige vorm heeft, die voor alle punten dezelfde is. De kolommen zijn:

${\displaystyle X_{0}=(1,\ldots ,1,1,1,1,1,\ldots ,1)^{\top }}$

${\displaystyle X_{1}=(-n,\ldots ,-2,-1,0,1,2,\ldots ,n)^{\top }}$

${\displaystyle X_{2}=((-n)^{2},\ldots ,(-2)^{2},(-1)^{2},0^{2},1^{2},2^{2},\ldots ,n^{2})^{\top }}$

${\displaystyle X_{3}=((-n)^{3},\ldots ,(-2)^{3},(-1)^{3},0^{3},1^{3},2^{3},\ldots ,n^{3})^{\top }}$

${\displaystyle \vdots }$

${\displaystyle X_{r}=((-n)^{r},\ldots ,(-2)^{r},(-1)^{r},0^{r},1^{r},2^{r},\ldots ,n^{r})^{\top }}$

en

${\displaystyle b_{k}=(x_{k-n},\ldots ,x_{k-2},x_{k-1},x_{k},x_{k+1},x_{k+2},\ldots ,x_{k+n})^{\top }}$

is de kolomvector met de waarden van de data in het venster.

De matrix ${\displaystyle X^{\top }X}$ is een matrix met constante antidiagonalen (diagonalen dwars op de hoofddiagonaal), dus bij snijden van de hooddiagonaal met elementen gelijk aan het element op de hoofddiagonaal:

${\displaystyle (X^{\top }X)_{i-j,i+j}=\sum _{m=-n}^{n}m^{i}}$,

en bij kruisen van de hoofddiagonaal de waarden 0.

De matrix ${\displaystyle W=(X^{\top }X)^{-1}X^{\top }}$ is alleen afhankelijk van de halve vensterbreedte ${\displaystyle n}$ en de graad ${\displaystyle r}$ van de polynoom, en kan dus berekend worden onafhankelijk van de data.

Voor bijvoorbeeld ${\displaystyle n=2}$ en ${\displaystyle r=2}$ is:

${\displaystyle X^{\top }X={\begin{bmatrix}\sum m^{0}&\sum m^{1}&\sum m^{2}\\\sum m^{1}&\sum m^{2}&\sum m^{3}\\\sum m^{2}&\sum m^{3}&\sum m^{4}\end{bmatrix}}={\begin{bmatrix}5&0&10\\0&10&0\\10&0&34\end{bmatrix}}}$.

Voor ${\displaystyle n=3}$ en ${\displaystyle r=2}$ is:

${\displaystyle X^{\top }X={\begin{bmatrix}\sum m^{0}&\sum m^{1}&\sum m^{2}\\\sum m^{1}&\sum m^{2}&\sum m^{3}\\\sum m^{2}&\sum m^{3}&\sum m^{4}\end{bmatrix}}={\begin{bmatrix}7&0&28\\0&28&0\\28&0&196\end{bmatrix}}}$

en voor ${\displaystyle n=3}$ en ${\displaystyle r=3}$ is:

${\displaystyle X^{\top }X={\begin{bmatrix}\sum m^{0}&\sum m^{1}&\sum m^{2}&\sum m^{3}\\\sum m^{1}&\sum m^{2}&\sum m^{3}&\sum m^{4}\\\sum m^{2}&\sum m^{3}&\sum m^{4}&\sum m^{5}\\\sum m^{3}&\sum m^{4}&\sum m^{5}&\sum m^{6}\end{bmatrix}}={\begin{bmatrix}7&0&28&0\\0&28&0&196\\28&0&196&0\\0&196&0&1588\end{bmatrix}}}$

Het gefilterde signaal bestaat uit de rij:

${\displaystyle \ldots ,a_{-2,0},a_{-1,0},a_{0,0},a_{1,0},a_{2,0},a_{3,0},\ldots }$,

met

${\displaystyle a_{k,0}=\sum _{i=-n}^{n}W_{0,i}x_{k+i}}$

Het gefilterde signaal is dus een vorm van convolutie van de eerste rij (met rij-index 0) van de matrix ${\displaystyle W}$ en het oorspronkelijke signaal, en wordt gegeven door de gewogen som van de waarden van de data in het venster.

Afgeleiden

De matrix ${\displaystyle W}$ kan ook gebruikt worden om afgeleiden van het gefilterde signaal te berekenen. In het ${\displaystyle k}$-de punt wordt het gefilterde signaal voorgesteld door de polynoom

${\displaystyle P_{k}(m)=a_{k0}+a_{k1}m+a_{k2}m^{2}+\ldots +a_{kr}m^{r}}$

Deze heeft de afgeleiden:

${\displaystyle P'_{k}(m)=a_{k1}+2a_{k2}m+\ldots +ra_{kr}m^{r-1}}$,

${\displaystyle P_{k}^{(2)}(m)=2a_{k2}+\ldots +r(r-1)a_{kr}m^{r-2}}$,

met in het ${\displaystyle k}$-de punt de waarde

${\displaystyle P'_{k}(0)=a_{k1}=\sum _{i=-n}^{n}W_{1,i}x_{k+i}}$

en

${\displaystyle P_{k}^{(2)}(0)=2a_{k2}=\sum _{i=-n}^{n}W_{2,i}x_{k+i}}$

Analoog geldt voor andere afgeleiden:

${\displaystyle P_{k}^{(s)}(0)=s!\,a_{ks}=\sum _{i=-n}^{n}W_{s,i}x_{k+i}}$

De rij van ${\displaystyle W}$ met rij-index ${\displaystyle s}$ levert dus de wegingsfactoren voor het berekenen van de afgeleide van orde ${\displaystyle s}$.

Voorbeeld

Voor ${\displaystyle n=3}$ en ${\displaystyle r=3}$ is de matrix ${\displaystyle W}$:

Gewichtsfactor	Waarde gewichtsfactor
w0	−0,095238095	0,142857143	0,285714286	0,333333333	0,285714286	0,142857143	−0,095238095
w1	0,087301587	−0,265873016	−0,23015873	0	0,23015873	0,265873016	−0,087301587
w2	0,05952381	0	−0,035714286	−0,047619048	−0,035714286	0	0,05952381
w3	−0,027777778	0,027777778	0,027777778	0	−0,027777778	−0,027777778	0,027777778

Eigenschappen en beperkingen

Het Savitzky-Golayfilter heeft ook zijn beperkingen. Als ruisfilter zijn de eigenschappen vaak eerder cosmetisch dan van praktisch belang omdat het filter het frequentiespectrum van de ruis nadelig beïnvloedt. Het onderdrukt een deel van de hoogfrequentie ruis, maar laat de langzamere componenten die vaak meer de analyse van de gegevens in de war sturen ongemoeid. Het is immers de laagfrequente ruis die het moeilijk maakt een piek van de signaalachtergrond te onderscheiden. Na filtering door middel van Savitzky-Golay is de ruis niet langer 'wit' en de meetfout niet langer normaal verdeeld.

Het aantal punten en het aantal termen in het filter met zorg gekozen moeten worden. De eigenschappen van een S-G filter hangen in belangrijke mate af van het gekozen venster en de gekozen orde van de polynoom. Als het model niet past is de uitkomst van de regressie niet te vertrouwen en komen er soms volkomen verkeerde getallen uit de berekening. Dit geldt eveneens in sterke mate voor de aanwezigheid van uitbijters in het signaal: een enkel volkomen fout punt (bijvoorbeeld ten gevolge van een schakelpuls in de stroom). Hoe groter het venster des te meer middeling en des te meer de ruis onderdrukt wordt. Hoe meer termen in de polynoom, hoe meer het oorspronkelijke signaal (inclusief ruis) benaderd wordt. Het benodigde venster voor het berekenen van de gekozen polynoom mag eveneens niet kleiner zijn dan het meetsignaal.

Bronnen, noten en/of referenties ↑ A. Savitkzy and M.J.E. Golay (1964). Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 36(8): 1627-1964