Savitzky-Golayfilter

Uit Wikipedia, de vrije encyclopedie
(Doorverwezen vanaf Savitsky-Golay filter)
Ga naar: navigatie, zoeken
Savitzky-golay pic gaussien bruite.svg
Lissage sg3 anim.gif

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", d.w.z. vlakt snelle wisselingen af, waardoor ruis wordt uitgefilterd. 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 oorsronkelijke data in het venster. De wegingsfactoren worden bepaald met een polynoom die in een omgeving van het punt goed bij de data past.

Principe[bewerken]

Een signaal is op equidistante afstanden gegeven door de rij:

\ldots,x_{-2},x_{-1},x_0,x_1,x_2,x_3,\ldots

Rondom het k-de punt in de rij wordt een venster

m=-n,\ldots,-2,-1,0,1,2,\ldots,n

beschouwd, dat dus zicht biedt op de data

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 r gedefinieerd door

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

\sum_m\left(P_k(m)-x_{k+m}\right)^2 minimaal is.

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


\begin{alignat}{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{alignat}

Deze vergelijkingen hebben de vorm:

X^\top X a_k=X^\top b_k,

met als oplossing

a_k=(X^\top X)^{-1}X^\top b_k=Wb_k,

waarin de (2n+1)×(r+1)-matrix 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:

X_0=(1,\ldots,1,1,1,1,1,\ldots,1)^\top
X_1=(-n,\ldots,-2,-1,0,1,2,\dots,n)^\top
X_2=((-n)^2,\ldots,(-2)^2,(-1)^2,0^2,1^2,2^2,\dots,n^2)^\top
X_3=((-n)^3,\ldots,(-2)^3,(-1)^3,0^3,1^3,2^3,\dots,n^3)^\top
\vdots
X_r=((-n)^r,\ldots,(-2)^r,(-1)^r,0^r,1^r,2^r,\dots,n^r)^\top

en

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.

A is het Ramanspectrum van cyclohexaan met veel ruis, B is hetzelfde spectrum na filtering met het Savitsky-Golayfilter

De matrix 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:

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

en bij kruisen van de hoofddiagonaal de waarden 0.

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

Voor bijvoorbeeld n=2 en r=2 is:

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 n=3 en r=2 is:

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 n=3 en r=3 is:

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:

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

met

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 W en het oorspronkelijke signaal, en wordt gegeven door de gewogen som van de waarden van de data in het venster.

Afgeleiden[bewerken]

A is het Ramanspectrum van cyclohexaan, B is de eerste afeleide en C de tweede afgeleide van het spectrum

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

P_k(m)=a_{k0}+a_{k1} m+a_{k2}m^2+\ldots+a_{kr}m^r.

Deze heeft de afgeleiden:

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

met in het k-de punt de waarde

P'_k(0)=a_{k1}=\sum_{i=-n}^n W_{1,i}x_{k+i}

en

P^{(2)}_k(0)=2a_{k2}=\sum_{i=-n}^n W_{2,i}x_{k+i}.

Analoog geldt voor andere afgeleiden:

P^{(s)}_k(0)=s!\,a_{ks}=\sum_{i=-n}^n W_{s,i}x_{k+i}.

De rij van W met rij-index s levert dus de wegingsfactoren voor het berekenen van de afgeleide van orde s.

Voorbeeld[bewerken]

Voor n=3 en r=3 is de matrix 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[bewerken]

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 minder 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
  1. A. Savitkzy and M.J.E. Golay (1964). Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 36(8): 1627-1964