Savitzky-Golay-Filter

Ein Savitzky-Golay-Filter ist ein digitaler Filter , der auf einen Satz digitaler Datenpunkte angewendet werden kann , um die Daten zu glätten und somit deren Genauigkeit zu erhöhen, ohne den Signalverlauf zu verfälschen. Dies geschieht durch Faltung, indem sukzessive Teilmengen benachbarter Datenpunkte mithilfe der Methode der kleinsten Quadrate durch ein Polynom niedrigen Grades angenähert werden . Bei äquidistanten Datenpunkten lässt sich eine analytische Lösung der Ausgleichsgleichungen in Form eines einzigen Satzes von Faltungskoeffizienten finden. Dieser kann auf alle Datenteilmengen angewendet werden, um Schätzwerte des geglätteten Signals (oder seiner Ableitungen) im Mittelpunkt jeder Teilmenge zu erhalten. Die Methode, die auf etablierten mathematischen Verfahren basiert ^{[ 1 , 2 ], wurde} von Abraham Savitzky und Marcel J. E. Golay bekannt gemacht , die 1964 Tabellen mit Faltungskoeffizienten für verschiedene Polynome und Teilmengengrößen veröffentlichten ^{[ 3 , 4 ]. Einige} Fehler in den Tabellen wurden korrigiert ^{[ 5 ]} . Die Methode wurde zur Behandlung zwei- und dreidimensionaler Daten erweitert.

Die Arbeit von Savitzky und Golay gehört zu den meistzitierten Artikeln in der Zeitschrift Analytical Chemistry ^{[ 6 ]} und wird von dieser Zeitschrift als eine ihrer „10 grundlegenden Arbeiten“ eingestuft, mit der Begründung: „Man kann argumentieren, dass der Beginn des computergesteuerten Analyseinstruments auf diesen Artikel zurückgeführt werden kann.“ ^{[ 7 ]}

Die Daten bestehen aus einer Menge von Punkten , wobei eine unabhängige Variable und ein beobachteter Wert ist. Sie werden mit einem Satz von Faltungskoeffizienten, , gemäß der Formel behandelt . Ausgewählte Faltungskoeffizienten sind in den folgenden Tabellen dargestellt . Beispielsweise ist für die Glättung mit einem quadratischen Polynom 5. Ordnung der -te geglättete Datenpunkt, , gegeben durch ${\displaystyle \left\{(x_{j},y_{j})\right\}_{j=1}^{n}}$${\displaystyle x_{j}}$${\displaystyle y_{j}}$${\displaystyle m}$${\displaystyle C_{i}}$$${\displaystyle Y_{j}=\sum _{i=-s}^{s}C_{i}\,y_{j+i},\qquad s+1\leq j\leq ns,\;\;s={\frac {m-1}{2}}}$$${\displaystyle m=5,i=-2,-1,0,1,2}$${\displaystyle j}$${\displaystyle Y_{j}}$

$${\displaystyle Y_{j}={\frac {1}{35}}\left(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2}\right),}$$

Es gibt zahlreiche Anwendungsgebiete für Glättung, beispielsweise um die Ausbreitung von Rauschen durch eine Algorithmuskette zu vermeiden oder manchmal einfach nur, um die Daten weniger verrauscht erscheinen zu lassen, als sie tatsächlich sind .${\displaystyle C_{-2}=-3/35,C_{-1}=12/35}$

Nachfolgend sind Anwendungen der numerischen Differentiation von Daten aufgeführt. ^{[ 8 ]} Hinweis: Bei der Berechnung der n -ten Ableitung kann ein zusätzlicher Skalierungsfaktor von auf alle berechneten Datenpunkte angewendet werden, um Absolutwerte zu erhalten (siehe Ausdrücke für , unten, für Details). ${\displaystyle {\frac {n!}{h^{n}}}}$${\displaystyle {\frac {d^{n}Y}{dx^{n}}}}$

(1) Synthetische Lorentz-Funktion + Rauschen (blau) und 1. Ableitung (grün)

(2) Titrationskurve (blau) für Malonsäure und 2. Derivat (grün). Der Ausschnitt im hellblauen Kasten ist 10-fach vergrößert.

(3) Lorentz-Funktion auf exponentieller Basislinie (blau) und 2. Ableitung (grün)

(4) Summe zweier Lorentz-Funktionen (blau) und 2. Ableitung (grün)

(5) 4. Ableitung der Summe zweier Lorentz-Funktionen

1. Lage von Maxima und Minima in experimentellen Datenkurven. Dies war die Anwendung, die Savitzky ursprünglich motivierte. ^{[ 4 ]} Die erste Ableitung einer Funktion ist an einem Maximum oder Minimum null. Das Diagramm zeigt Datenpunkte einer synthetischen Lorentz- Kurve mit hinzugefügtem Rauschen (blaue Rauten). Die Daten sind auf einer Skala der Halbwertsbreite relativ zum Maximum bei Null dargestellt. Die geglättete Kurve (rote Linie) und die erste Ableitung (grün) wurden mit kubischen Savitzky-Golay-Filtern mit 7 Punkten berechnet. Durch lineare Interpolation der Werte der ersten Ableitung an Positionen beidseits des Nulldurchgangs ergibt sich die Position des Maximums. Auch dritte Ableitungen können zu diesem Zweck verwendet werden.
2. Lage des Endpunktes in einer Titrationskurve . Ein Endpunkt ist ein Wendepunkt, an dem die zweite Ableitung der Funktion null ist. ^{[ 9 ]} Die Titrationskurve für Malonsäure verdeutlicht die Leistungsfähigkeit der Methode. Der erste Endpunkt bei 4 ml ist kaum sichtbar, aber die zweite Ableitung ermöglicht es, seinen Wert durch lineare Interpolation zur Bestimmung des Nulldurchgangs leicht zu ermitteln .
3. Basislinienglättung. In der analytischen Chemie ist es mitunter erforderlich, die Höhe einer Absorptionsbande relativ zu einer gekrümmten Basislinie zu messen. ^{[ 10 ]} Da die Krümmung der Basislinie deutlich geringer ist als die der Absorptionsbande, glättet die zweite Ableitung die Basislinie effektiv. Drei Maße für die Ableitungshöhe, die proportional zur Höhe der Absorptionsbande ist, sind die Abstände zwischen Maximum und Minimum h1 und h2 sowie der Abstand zur Basislinie h3. ^{[ 11 ]}
4. Verbesserung der Auflösung in der Spektroskopie. Die Banden in der zweiten Ableitung einer Spektroskopiekurve sind schmaler als die Banden im Spektrum: Sie weisen eine reduzierte Halbwertsbreite auf . Dadurch lassen sich teilweise überlappende Banden in separate (negative) Peaks auflösen. ^{[ 12 ]} Die Abbildung veranschaulicht, wie dies auch für die chemische Analyse genutzt werden kann , indem die Abstände zwischen Peaks und Tälern gemessen werden. In diesem Fall sind die Täler eine Eigenschaft der zweiten Ableitung einer Lorentz-Funktion. ( Die Position auf der x -Achse bezieht sich auf die Position des Peakmaximums auf einer Skala der Halbwertsbreite bei halber Höhe .)
5. Verbesserung der Auflösung durch die 4. Ableitung (positive Peaks). Die Minima sind eine Eigenschaft der 4. Ableitung einer Lorentz-Funktion.

Der „gleitende Mittelwertfilter“ ist ein einfaches Beispiel für einen Savitzky-Golay-Filter, der häufig bei Zeitreihendaten verwendet wird , um kurzfristige Schwankungen auszugleichen und längerfristige Trends oder Zyklen hervorzuheben. Jede Teilmenge der Daten wird durch eine horizontale Gerade anstatt durch ein Polynom höherer Ordnung angenähert. Ein ungewichteter gleitender Mittelwertfilter ist der einfachste Faltungsfilter.

Der gleitende Durchschnitt wird häufig für eine schnelle technische Analyse von Finanzdaten wie Aktienkursen, Renditen oder Handelsvolumina verwendet. Er findet auch in der Volkswirtschaftslehre Anwendung, um das Bruttoinlandsprodukt , die Beschäftigung oder andere makroökonomische Zeitreihen zu untersuchen.

Es wurde in einigen Tabellen der Savitzky-Golay-Faltungskoeffizienten nicht berücksichtigt, da alle Koeffizientenwerte identisch sind, nämlich der Wert . ${\displaystyle {\frac {1}{m}}}$

When the data points are equally spaced, an analytical solution to the least-squares equations can be found.^[2] This solution forms the basis of the convolution method of numerical smoothing and differentiation. Suppose that the data consists of a set of n points (x_j, y_j) (j = 1, ..., n), where x_j is an independent variable and y_j is a datum value. A polynomial will be fitted by linear least squares to a set of m (an odd number) adjacent data points, each separated by an interval h. Firstly, a change of variable is made $${\displaystyle z={\frac {x-{\bar {x}}}{h}}}$$ where ${\displaystyle {\bar {x}}}$ is the value of the central point. z takes the values ${\displaystyle -s,\dots ,0,\dots ,s}$ where ${\displaystyle s={\tfrac {m-1}{2}}}$ (e.g. m = 5 → z = −2, −1, 0, 1, 2).^{[note 1]} The polynomial, of degree k is defined as^{[note 2]}$${\displaystyle Y=a_{0}+a_{1}z+a_{2}z^{2}+\dots +a_{k}z^{k}.}$$ The coefficients a₀, a₁ etc. are obtained by solving the normal equations (bold a represents a vector, bold J represents a matrix).

$${\displaystyle \mathbf {a} =\left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)^{-1}\mathbf {J} ^{\mathsf {T}}\mathbf {y} ,}$$ where ${\displaystyle \mathbf {J} }$ is a Vandermonde matrix, that is ${\displaystyle i}$-th row of ${\displaystyle \mathbf {J} }$ has values ${\displaystyle 1,z_{i},z_{i}^{2},\dots }$.

For example, for a cubic polynomial fitted to 5 points, z = −2, −1, 0, 1, 2 the normal equations are solved as follows.

$${\displaystyle \mathbf {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}}}$$

$${\displaystyle {\begin{aligned}\mathbf {J} ^{\mathsf {T}}\mathbf {J} &={\begin{pmatrix}m&\sum z&\sum z^{2}&\sum z^{3}\\\sum z&\sum z^{2}&\sum z^{3}&\sum z^{4}\\\sum z^{2}&\sum z^{3}&\sum z^{4}&\sum z^{5}\\\sum z^{3}&\sum z^{4}&\sum z^{5}&\sum z^{6}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}m&0&\sum z^{2}&0\\0&\sum z^{2}&0&\sum z^{4}\\\sum z^{2}&0&\sum z^{4}&0\\0&\sum z^{4}&0&\sum z^{6}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}5&0&10&0\\0&10&0&34\\10&0&34&0\\0&34&0&130\end{pmatrix}}\end{aligned}}}$$

Now, the normal equations can be factored into two separate sets of equations, by rearranging rows and columns, with

$${\displaystyle \mathbf {J} ^{\mathsf {T}}\mathbf {J} _{\text{even}}={\begin{pmatrix}5&10\\10&34\\\end{pmatrix}}\quad \mathrm {and} \quad \mathbf {J} ^{\mathsf {T}}\mathbf {J} _{\text{odd}}={\begin{pmatrix}10&34\\34&130\\\end{pmatrix}}}$$ Expressions for the inverse of each of these matrices can be obtained using Cramer's rule$${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)_{\text{even}}^{-1}={\frac {1}{70}}{\begin{pmatrix}34&-10\\-10&5\\\end{pmatrix}}\quad \mathrm {and} \quad \left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)_{\text{odd}}^{-1}={\frac {1}{144}}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}}$$ The normal equations become $${\displaystyle {\begin{pmatrix}a_{0}\\a_{2}\end{pmatrix}}_{j}={\frac {1}{70}}{\begin{pmatrix}34&-10\\-10&5\end{pmatrix}}{\begin{pmatrix}1&1&1&1&1\\4&1&0&1&4\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\end{pmatrix}}}$$ and $${\displaystyle {\begin{pmatrix}a_{1}\\a_{3}\end{pmatrix}}_{j}={\frac {1}{144}}{\begin{pmatrix}130&-34\\-34&10\end{pmatrix}}{\begin{pmatrix}-2&-1&0&1&2\\-8&-1&0&1&8\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\\\end{pmatrix}}}$$ Multiplying out and removing common factors, $${\displaystyle {\begin{aligned}a_{0,j}&={\frac {1}{35}}\left(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2}\right)\\[1ex]a_{1,j}&={\frac {1}{12}}\left(y_{j-2}-8y_{j-1}+8y_{j+1}-y_{j+2}\right)\\[1ex]a_{2,j}&={\frac {1}{14}}\left(2y_{j-2}-y_{j-1}-2y_{j}-y_{j+1}+2y_{j+2}\right)\\[1ex]a_{3,j}&={\frac {1}{12}}\left(-y_{j-2}+2y_{j-1}-2y_{j+1}+y_{j+2}\right)\end{aligned}}}$$

The coefficients of y in these expressions are known as convolution coefficients. They are elements of the matrix $${\displaystyle \mathbf {C} =\left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)^{-1}\mathbf {J} ^{\mathsf {T}}}$$ In general, $${\displaystyle (C\otimes y)_{j}\ =Y_{j}=\sum _{i=-s}^{s}C_{i}\,y_{j+i},\qquad s+1\leq j\leq n-s,\;\;s={\frac {m-1}{2}}}$$ In matrix notation this example is written as $${\displaystyle {\begin{pmatrix}Y_{3}\\Y_{4}\\Y_{5}\\\vdots \end{pmatrix}}={\frac {1}{35}}{\begin{pmatrix}-3&12&17&12&-3&0&0&\cdots \\0&-3&12&17&12&-3&0&\cdots \\0&0&-3&12&17&12&-3&\cdots \\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\ddots \end{pmatrix}}{\begin{pmatrix}y_{1}\\y_{2}\\y_{3}\\y_{4}\\y_{5}\\y_{6}\\y_{7}\\\vdots \end{pmatrix}}}$$

Tabellen mit Faltungskoeffizienten, die für m bis 25 auf dieselbe Weise berechnet wurden , wurden 1964 für den Savitzky-Golay-Glättungsfilter veröffentlicht ^{[ 3 ] [ 5 ].} _Der Wert des Mittelpunkts, z = 0, wird aus einem einzigen Koeffizientensatz ermittelt: 0 für die Glättung, 1 für die erste Ableitung usw. Die numerischen Ableitungen werden durch Differenzieren von Y erhalten _. Dies bedeutet, dass die Ableitungen für die geglättete Datenkurve berechnet werden . Für ein kubisches Polynom $${\displaystyle {\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}=a_{0}&{\text{ at }}z=0,x={\bar {x}}\\[1ex]{\frac {dY}{dx}}&={\frac {1}{h}}\left({a_{1}+2a_{2}z+3a_{3}z^{2}}\right)={\frac {1}{h}}a_{1}&{\text{ at }}z=0,x={\bar {x}}\\[1ex]{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}\left({2a_{2}+6a_{3}z}\right)={\frac {2}{h^{2}}}a_{2}&{\text{ at }}z=0,x={\bar {x}}\\[1ex]{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}}$$

Im Allgemeinen liefern Polynome vom Grad (0 und 1), ^{[ Anmerkung 3 ]} (2 und 3), (4 und 5) usw. dieselben Koeffizienten für Glättungs- und gerade Ableitungen. Polynome vom Grad (1 und 2), (3 und 4) usw. liefern dieselben Koeffizienten für ungerade Ableitungen.

Die Verwendung der Savitzky-Golay-Tabellen ist nicht immer erforderlich. Die Summen in der Matrix J T J lassen sich in geschlossener Form auswerten , sodass algebraische Formeln für die Faltungskoeffizienten hergeleitet werden können. [ 13 ] [ Anmerkung 4 ] Geeignete Funktionen für Kurven mit Wendepunkten sind: $${\displaystyle {\begin{aligned}\sum _{z=-s}^{s}z^{2}&={m(m^{2}-1) \over 12},&s={\frac {m-1}{2}}\\\sum _{z}z^{4}&={m(m^{2}-1)(3m^{2}-7) \over 240}\\\sum _{z}z^{6}&={m(m^{2}-1)(3m^{4}-18m^{2}+31) \over 1344}\end{aligned}}}$$

• Glättung, Polynomgrad 2,3 : (Der Wertebereich für i gilt auch für die folgenden Ausdrücke)${\displaystyle C_{0i}={\frac {\left({3m^{2}-7-20i^{2}}\right)/4}{m\left({m^{2}-4}\right)/3}};\quad {\frac {1-m}{2}}\leq i\leq {\frac {m-1}{2}}}$
• 1. Ableitung: Polynomgrad 3,4${\displaystyle C_{1i}={\frac {5\left({3m^{4}-18m^{2}+31}\right)i-28\left({3m^{2}-7}\right)i^{3}}{m\left({m^{2}-1}\right)\left({3m^{4}-39m^{2}+108}\right)/15}}}$
• 2. Ableitung: Polynomgrad 2,3${\displaystyle C_{2i}={\frac {12mi^{2}-m\left({m^{2}-1}\right)}{m^{2}\left({m^{2}-1}\right)\left({m^{2}-4}\right)/30}}}$
• 3. Ableitung: Polynomgrad 3,4${\displaystyle C_{3i}={\frac {-\left({3m^{2}-7}\right)i+20i^{3}}{m\left({m^{2}-1}\right)\left({3m^{4}-39m^{2}+108}\right)/420}}}$

Einfachere Ausdrücke, die für Kurven ohne Wendepunkt verwendet werden können, sind:

• Glättung, Polynomgrad 0,1 (gleitender Durchschnitt):${\displaystyle C_{0i}={\frac {1}{m}}}$
• 1. Ableitung, Polynomgrad 1,2:${\displaystyle C_{1i}={\frac {i}{m(m^{2}-1)/12}}}$

Höhere Ableitungen können ermittelt werden. Beispielsweise kann eine vierte Ableitung durch zweimalige Anwendung der zweiten Ableitungsfunktion erhalten werden. ^{[ 14 ]}

Eine Alternative zur Anpassung von m Datenpunkten durch ein einfaches Polynom in der Hilfsvariablen z ist die Verwendung orthogonaler Polynome . Dabei ist  P 0 , ...,  P k eine Menge paarig orthogonaler Polynome vom Grad 0, ..., k . Ausführliche Informationen zur Herleitung von Ausdrücken für die orthogonalen Polynome und zum Zusammenhang zwischen den Koeffizienten b und a finden sich bei Guest [ 2 ] . Ausdrücke für die Faltungskoeffizienten lassen sich leicht herleiten, da die Normalengleichungsmatrix J T J eine Diagonalmatrix ist , da das Produkt zweier orthogonaler Polynome aufgrund ihrer gegenseitigen Orthogonalität null ist. Daher ist jedes von null verschiedene Element ihrer Inversen einfach der Kehrwert des entsprechenden Elements in der Normalengleichungsmatrix. Die Berechnung wird weiter vereinfacht, indem orthogonale Gram-Polynome rekursiv konstruiert werden . Die gesamte Berechnung kann in wenigen Zeilen PASCAL- Code implementiert werden , einer Programmiersprache , die sich gut für rekursive Berechnungen eignet [ 15 ] .$${\displaystyle Y=b_{0}P^{0}(z)+b_{1}P^{1}(z)\cdots +b_{k}P^{k}(z).}$$

Savitzky-Golay-Filter werden am häufigsten verwendet, um den geglätteten Wert oder die Ableitung am Mittelpunkt z = 0 mithilfe eines einzigen Satzes von Faltungskoeffizienten zu erhalten. ( m  − 1)/2 Punkte am Anfang und Ende der Reihe können mit diesem Verfahren nicht berechnet werden. Es gibt verschiedene Strategien, um diese Einschränkung zu umgehen.

• Die Daten könnten künstlich erweitert werden, indem in umgekehrter Reihenfolge Kopien der ersten ( m  − 1)/2 Punkte am Anfang und Kopien der letzten ( m  − 1)/2 Punkte am Ende hinzugefügt werden. Beispielsweise werden bei m = 5 jeweils zwei  Punkte am Anfang und am Ende der Daten y₁ _, ..., yₙ hinzugefügt _.

y ₃ , y ₂ , y ₁ , ... , y _n , y _{n −1} , y _{n −2} .

• Wenn man sich das Anpassungspolynom noch einmal ansieht, wird deutlich, dass Daten für alle Werte von z berechnet werden können , indem man alle Sätze von Faltungskoeffizienten für ein einzelnes Polynom, a ₀ .. a _k , verwendet .
  Für ein kubisches Polynom$${\displaystyle {\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}\\{\frac {dY}{dx}}&={\frac {1}{h}}(a_{1}+2a_{2}z+3a_{3}z^{2})\\{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}(2a_{2}+6a_{3}z)\\{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}}$$
• Die Faltungskoeffizienten für die fehlenden ersten und letzten Punkte lassen sich ebenfalls leicht bestimmen. ^{[ 15 ]} Dies entspricht der Anpassung der ersten ( m  + 1)/2 Punkte mit demselben Polynom und analog für die letzten Punkte.

Die obige Vorgehensweise impliziert, dass alle Datenpunkte gleich gewichtet werden. Genauer gesagt, hat die im Verfahren der kleinsten Quadrate zu minimierende Zielfunktion Einheitsgewichte, w _{_i}  = 1. Sind die Gewichte nicht alle gleich, lauten die Normalgleichungen wie folgt: $${\displaystyle U=\sum _{i}w_{i}\left(Y_{i}-y_{i}\right)^{2}}$$

$${\displaystyle \mathbf {a} =\left(\mathbf {J} ^{\mathsf {T}}\mathbf {W} \mathbf {J} \right)^{-1}\mathbf {J} ^{\mathsf {T}}\mathbf {W} \mathbf {y} ,\qquad W_{i,i}\neq 1}$$

Verwendet man für alle Datenteilmengen denselben Satz von Diagonalgewichten, so lässt sich eine analytische Lösung der Normalgleichungen angeben. Zum Beispiel mit einem quadratischen Polynom: ${\displaystyle W={\text{diag}}(w_{1},w_{2},...,w_{m})}$

$${\displaystyle \mathbf {J} ^{\mathsf {T}}\mathbf {W} \mathbf {J} ={\begin{pmatrix}\sum w_{i}~~~&\sum w_{i}z_{i}&\sum w_{i}z_{i}^{2}\\\sum w_{i}z_{i}&\sum w_{i}z_{i}^{2}&\sum w_{i}z_{i}^{3}\\\sum w_{i}z_{i}^{2}&\sum w_{i}z_{i}^{3}&\sum w_{i}z_{i}^{4}\end{pmatrix}}}$$

Ein expliziter Ausdruck für die Inverse dieser Matrix lässt sich mithilfe der Cramerschen Regel gewinnen . Daraus kann dann ein Satz von Faltungskoeffizienten abgeleitet werden.

$${\displaystyle \mathbf {C} =\left(\mathbf {J} ^{\mathsf {T}}\mathbf {W} \mathbf {J} \right)^{-1}\mathbf {J} ^{\mathsf {T}}\mathbf {W} .}$$

Alternativ können die Koeffizienten C in einer Tabellenkalkulation berechnet werden , indem eine integrierte Matrixinversionsroutine verwendet wird, um die Inverse der Normalengleichungsmatrix zu erhalten. Dieser Koeffizientensatz kann nach der Berechnung und Speicherung für alle Berechnungen verwendet werden, die dasselbe Gewichtungsschema verwenden. Für jedes abweichende Gewichtungsschema wird ein anderer Koeffizientensatz benötigt.

Es wurde gezeigt, dass der Savitzky-Golay-Filter durch die Einführung von Gewichtungen verbessert werden kann, die an den Enden des Anpassungsintervalls abnehmen. ^{[ 16 ]}

Zweidimensionale Glättung und Differenzierung lassen sich auch auf Datentabellen anwenden, beispielsweise auf Intensitätswerte in einem fotografischen Bild, das aus einem rechteckigen Pixelraster besteht. ^{[ 17 ] [ 18 ]} Ein solches Raster wird als Kernel bezeichnet, die Datenpunkte des Kernels als Knoten. Der Trick besteht darin, den rechteckigen Kernel durch einfaches Ordnen der Knotenindizes in eine einzelne Zeile zu transformieren. Während die eindimensionalen Filterkoeffizienten durch Anpassen eines Polynoms in der Hilfsvariablen z an m Datenpunkte ermittelt werden, erhält man die zweidimensionalen Koeffizienten durch Anpassen eines Polynoms in den Hilfsvariablen v und w an die Werte der m × n Kernelknoten. Das folgende Beispiel für ein bivariates Polynom vom Grad 3 ( m = 7, n = 5) veranschaulicht den Prozess, der dem oben beschriebenen eindimensionalen Fall entspricht. ^{[ 19 ]} Der rechteckige Kern von 35 Datenwerten, d ₁  −  d ₃₅$${\displaystyle {\begin{aligned}v&={\frac {x-{\bar {x}}}{h(x)}};&w&={\frac {y-{\bar {y}}}{h(y)}}\end{aligned}}}$$$${\displaystyle Y=a_{00}+a_{10}v+a_{01}w+a_{20}v^{2}+a_{11}vw+a_{02}w^{2}+a_{30}v^{3}+a_{21}v^{2}w+a_{12}vw^{2}+a_{03}w^{3}}$$

wird zu einem Vektor, wenn die Zeilen nacheinander angeordnet werden.

Die Jacobi-Matrix besitzt 10 Spalten, eine für jeden der Parameter a _₀₀  −  a _₀₃ , und 35 Zeilen, eine für jedes Wertepaar von v und w . Jede Zeile hat die Form [ Formel einfügen] . Die Faltungskoeffizienten werden wie folgt berechnet : [Formel einfügen]. Die erste Zeile von C enthält 35 Faltungskoeffizienten, die jeweils mit den 35 Datenwerten multipliziert werden können, um den Polynomkoeffizienten zu erhalten. Dieser entspricht dem geglätteten Wert am zentralen Knoten des Kernels (d. h. am 18. Knoten der obigen Tabelle). Analog können die anderen Zeilen von C mit den 35 Werten multipliziert werden, um weitere Polynomkoeffizienten zu erhalten. Diese wiederum können verwendet werden, um geglättete Werte und verschiedene geglättete partielle Ableitungen an verschiedenen Knoten zu berechnen. $${\displaystyle J_{\text{row}}={\begin{pmatrix}1&v&w&v^{2}&vw&w^{2}&v^{3}&v^{2}w&vw^{2}&w^{3}\end{pmatrix}}}$$$${\displaystyle \mathbf {C} =\left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)^{-1}\mathbf {J} ^{\mathsf {T}}}$$${\displaystyle a_{00}}$

Nikitas und Pappa-Louisi zeigten, dass die Qualität der Glättung je nach Format des verwendeten Polynoms erheblich variieren kann. ^{[ 20 ]} Sie empfehlen die Verwendung eines Polynoms der Form

$${\displaystyle Y=\sum _{i=0}^{p}\sum _{j=0}^{q}a_{ij}v^{i}w^{j}}$$

Da solche Polynome sowohl im Zentrum als auch in den Randbereichen eines Kernels eine gute Glättung erzielen, können sie bedenkenlos zur Glättung sowohl der internen als auch der randnahen Datenpunkte eines abgetasteten Bereichs verwendet werden. Um eine schlechte Konditionierung bei der Lösung des Ausgleichsproblems zu vermeiden, gilt p < m und q < n . Software zur Berechnung der zweidimensionalen Koeffizienten sowie eine Datenbank solcher C -Werte finden Sie im Abschnitt über mehrdimensionale Faltungskoeffizienten weiter unten.

Die Idee der zweidimensionalen Faltungskoeffizienten lässt sich auf einfache Weise auch auf höhere räumliche Dimensionen erweitern ^{[ 17 ] [ 21 ]} , indem die mehrdimensionale Verteilung der Kernelknoten in einer einzigen Zeile angeordnet wird. In Anlehnung an die oben erwähnte Erkenntnis von Nikitas und Pappa-Louisi ^{[ 20 ]} für zweidimensionale Fälle wird für mehrdimensionale Fälle die Verwendung der folgenden Polynomform empfohlen:

$${\displaystyle Y=\sum _{i_{1}=0}^{p_{1}}\sum _{i_{2}=0}^{p_{2}}\cdots \sum _{i_{D}=0}^{p_{D}}(a_{i_{1}i_{2}\cdots i_{D}}\times u_{1}^{i_{1}}u_{2}^{i_{2}}\cdots u_{D}^{i_{D}})}$$

Dabei ist D die Dimension des Raumes, die Koeffizienten des Polynoms und die Koordinaten in den verschiedenen Raumrichtungen. Algebraische Ausdrücke für partielle Ableitungen beliebiger Ordnung, ob gemischt oder nicht, lassen sich leicht aus dem obigen Ausdruck ableiten. ^{[ 21 ]} Es ist zu beachten, dass C von der Anordnung der Kernknoten in einer Zeile und der Anordnung der Terme der entwickelten Form des obigen Polynoms bei der Erstellung der Jacobi-Matrix abhängt. ${\displaystyle a}$

Die genaue Berechnung von C in mehrdimensionalen Fällen wird schwierig, da die Genauigkeit gängiger Gleitkommazahlen in Programmiersprachen nicht mehr ausreicht. Die unzureichende Genauigkeit führt dazu, dass die Abschneidefehler der Gleitkommazahlen mit der Größenordnung einiger C- Elemente vergleichbar werden, was die Genauigkeit stark beeinträchtigt und die Berechnung unbrauchbar macht. Chandra Shekhar hat mit Advanced Convolution Coefficient Calculator (ACCC) und Precise Convolution Coefficient Calculator (PCCC) zwei Open-Source -Softwareprogramme entwickelt , die diese Genauigkeitsprobleme adäquat lösen. ACCC berechnet C iterativ mit Gleitkommazahlen. ^{[ 22 ]} Die Genauigkeit der Gleitkommazahlen wird in jeder Iteration mithilfe von GNU MPFR schrittweise erhöht . Sobald die in zwei aufeinanderfolgenden Iterationen erhaltenen C -Werte bis zu einem vordefinierten Abstand die gleiche Anzahl signifikanter Stellen aufweisen, gilt die Konvergenz als erreicht. Ist der Abstand ausreichend groß, liefert die Berechnung einen hochgenauen C- Wert . PCCC verwendet Berechnungen mit rationalen Zahlen unter Verwendung der GNU Multiple Precision Arithmetic Library und liefert ein vollständig genaues C im Format rationaler Zahlen . ^{[ 23 ]} Abschließend werden diese rationalen Zahlen in Gleitkommazahlen umgewandelt, bis eine vordefinierte Anzahl signifikanter Stellen erreicht ist.

Eine Datenbank mit C -Werten , die mittels ACCC für symmetrische Kerne sowie symmetrische und asymmetrische Polynome auf Kernknoten mit Einheitsabstand in ein-, zwei-, drei- und vierdimensionalen Räumen berechnet wurden, ist verfügbar. ^{[ 24 ]} Chandra Shekhar hat zudem ein mathematisches Rahmenwerk entwickelt, das die Verwendung von C-Werten , berechnet auf Kernknoten mit Einheitsabstand, zur Filterung und partiellen Differentiation (verschiedener Ordnungen) auf Kernknoten mit ungleichmäßigem Abstand beschreibt. ^{[ 21 ]} Dadurch können die in der oben genannten Datenbank enthaltenen C- Werte verwendet werden . Obwohl diese Methode nur Näherungswerte liefert, sind diese in den meisten technischen Anwendungen akzeptabel, sofern die Nicht-Uniformität der Kernknoten gering ist.

1. Die Summe der Faltungskoeffizienten zur Glättung ist gleich eins. Die Summe der Koeffizienten für ungerade Ableitungen ist null. ^{[ 25 ]}
2. Die Summe der quadrierten Faltungskoeffizienten zur Glättung entspricht dem Wert des zentralen Koeffizienten. ^{[ 26 ]}
3. Durch Glätten einer Funktion bleibt die Fläche unter der Funktion unverändert. ^{[ 25 ]}
4. Die Faltung einer symmetrischen Funktion mit Koeffizienten gerader Ableitungen erhält das Symmetriezentrum. ^{[ 25 ]}
5. Eigenschaften von Ableitungsfiltern. ^{[ 27 ]}

Es ist unvermeidlich, dass das Signal im Faltungsprozess verzerrt wird. Gemäß Eigenschaft 3 wird bei der Glättung von Daten mit einem Peak die Peakhöhe reduziert und die Halbwertsbreite erhöht. Sowohl das Ausmaß der Verzerrung als auch die Verbesserung des Signal-Rausch-Verhältnisses (S/N ) werden dadurch verbessert.

• sinkt mit zunehmendem Grad des Polynoms
• Sie nimmt zu, wenn die Breite m der Faltungsfunktion zunimmt.

Wenn beispielsweise das Rauschen in allen Datenpunkten unkorreliert ist und eine konstante Standardabweichung σ aufweist , wird die Standardabweichung des Rauschens durch Faltung mit einer m -Punkt-Glättungsfunktion auf ^{[ 26 ] [ Anmerkung 5 ] verringert.}

• Polynomgrad 0 oder 1: ( gleitender Durchschnitt )${\displaystyle {\sqrt {\frac {1}{m}}}\sigma }$
• Polynomgrad 2 oder 3: .${\displaystyle {\sqrt {\frac {3(3m^{2}-7)}{4m(m^{2}-4)}}}\sigma }$

Diese Funktionen sind in der Grafik rechts dargestellt. Beispielsweise werden mit einer linearen 9-Punkt-Funktion (gleitender Mittelwert) zwei Drittel des Rauschens entfernt, während mit einer quadratischen/kubischen 9-Punkt-Glättungsfunktion nur etwa die Hälfte des Rauschens entfernt wird. Der größte Teil des verbleibenden Rauschens ist niederfrequentes Rauschen (siehe Frequenzgang von Faltungsfiltern weiter unten).

Obwohl die gleitende Mittelwertfunktion eine bessere Rauschunterdrückung ermöglicht, ist sie ungeeignet zur Glättung von Daten mit Krümmung über m Punkte. Eine quadratische Filterfunktion ist ungeeignet, um die Ableitung einer Datenkurve mit einem Wendepunkt zu bestimmen , da ein quadratisches Polynom keinen Wendepunkt besitzt. Die optimale Wahl des Polynomgrades und der Anzahl der Faltungskoeffizienten stellt einen Kompromiss zwischen Rauschunterdrückung und Verzerrung dar. ^{[ 28 ]}

Eine Möglichkeit, Verzerrungen zu mindern und die Rauschunterdrückung zu verbessern, besteht darin, einen Filter mit geringerer Breite zu verwenden und diesen mehrfach zu falten. Zwei Durchläufe desselben Filters entsprechen einem Durchlauf eines Filters, der durch Faltung des ursprünglichen Filters mit sich selbst entsteht. ^{[ 29 ]} Beispielsweise sind zwei Durchläufe des Filters mit den Koeffizienten (1/3, 1/3, 1/3) äquivalent zu einem Durchlauf des Filters mit den Koeffizienten (1/9, 2/9, 3/9, 2/9, 1/9).

Der Nachteil der Mehrfachfilterung besteht darin, dass die äquivalente Filterbreite für Durchläufe einer n-Punkt-Funktion größer ist, wodurch die Mehrfachfilterung stärkeren Randeffekten unterliegt. Dennoch wurde die Mehrfachfilterung mit großem Erfolg eingesetzt. Beispielsweise lieferten etwa 40–80 Durchläufe von Daten mit einem Signal-Rausch-Verhältnis von nur 5 brauchbare Ergebnisse. ^{[ 30 ]} Die oben genannten Formeln zur Rauschunterdrückung sind nicht anwendbar, da die Korrelation zwischen den berechneten Datenpunkten mit jedem Durchlauf zunimmt. ${\displaystyle n}$${\displaystyle m}$${\displaystyle n(m-1)+1}$

Die Faltung entspricht der Multiplikation im Fourier -Koordinate . Die diskrete Fouriertransformierte eines Faltungsfilters ist eine reellwertige Funktion , die sich wie folgt darstellen lässt:

$${\displaystyle FT(\theta )=\sum _{j=-s}^{s}C_{j}\cos(j\theta ),\qquad s={\tfrac {m-1}{2}}}$$

θ läuft von 0 bis 180 Grad , danach wiederholt sich die Funktion. Der Verlauf einer quadratischen/kubischen Glättungsfunktion mit neun Punkten ist typisch. Bei sehr kleinen Winkeln ist der Verlauf nahezu flach, was bedeutet, dass niederfrequente Datenkomponenten durch die Glättung praktisch unverändert bleiben. Mit zunehmendem Winkel nimmt der Wert ab, sodass höherfrequente Komponenten immer stärker gedämpft werden. Dies zeigt, dass der Faltungsfilter als Tiefpassfilter beschrieben werden kann : Das entfernte Rauschen ist primär hochfrequentes Rauschen, während niederfrequentes Rauschen den Filter passiert. ^{[ 31 ]} Einige hochfrequente Rauschkomponenten werden stärker gedämpft als andere, was sich in Wellen in der Fouriertransformation bei großen Winkeln zeigt. Dies kann zu kleinen Oszillationen in den geglätteten Daten ^{[ 32 ]} und Phasenumkehr führen, d. h. hochfrequente Oszillationen in den Daten werden durch die Savitzky-Golay-Filterung invertiert. ^{[ 33 ]}

Die Faltung beeinflusst die Korrelation zwischen Fehlern in den Daten. Der Effekt der Faltung lässt sich als lineare Transformation ausdrücken.

$${\displaystyle Y_{j}\ =\sum _{i=-s}^{s}C_{i}\,y_{j+i},\qquad s={\tfrac {m-1}{2}}}$$

Gemäß dem Fehlerfortpflanzungsgesetz wird die Varianz-Kovarianz-Matrix der Daten, A , wie folgt in B transformiert: Um die praktische Anwendung zu verdeutlichen, betrachten wir den Effekt eines gleitenden 3-Punkte-Mittelwerts auf die ersten drei berechneten Punkte, Y₂  −  Y₄ , unter der Annahme _, dass die Datenpunkte die gleiche Varianz aufweisen und keine Korrelation zwischen ihnen besteht. A ist dann die Einheitsmatrix, _{multipliziert} mit einer Konstanten σ² ^, der Varianz an jedem Punkt. $${\displaystyle \mathbf {B} =\mathbf {C} \mathbf {A} \mathbf {C} ^{\mathsf {T}}}$$

$${\displaystyle {\begin{aligned}\mathbf {B} &={\frac {\sigma ^{2}}{9}}{\begin{pmatrix}1&1&1&0&0\\0&1&1&1&0\\0&0&1&1&1\end{pmatrix}}{\begin{pmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{pmatrix}}{\begin{pmatrix}1&0&0\\1&1&0\\1&1&1\\0&1&1\\0&0&1\end{pmatrix}}\\&={\frac {\sigma ^{2}}{9}}{\begin{pmatrix}3&2&1\\2&3&2\\1&2&3\end{pmatrix}}\end{aligned}}}$$

In diesem Fall ergeben sich die Korrelationskoeffizienten zwischen den berechneten Punkten i und j wie folgt : $${\displaystyle \rho _{ij}={\frac {B_{ij}}{\sqrt {B_{ii}B_{jj}}}},\ (i\neq j)}$$

$${\displaystyle \rho _{i,i+1}={2 \over 3}=0.66,\;\rho _{i,i+2}={1 \over 3}=0.33}$$

Im Allgemeinen sind die berechneten Werte korreliert, selbst wenn die beobachteten Werte nicht korreliert sind. Die Korrelation erstreckt sich jeweils über m  − 1 berechnete Punkte. ^{[ 34 ]}

Um den Einfluss mehrfacher Filterung auf das Rauschen und die Korrelation eines Datensatzes zu veranschaulichen, betrachten wir die Auswirkungen eines zweiten Durchlaufs eines gleitenden 3-Punkt-Mittelwertfilters. Für den zweiten Durchlauf ^{[ Anmerkung 6 ]}

$${\displaystyle {\begin{aligned}\mathbf {C} \mathbf {B} \mathbf {C} ^{\mathsf {T}}&={\frac {\sigma ^{2}}{81}}{\begin{pmatrix}1&1&1&0&0\\0&1&1&1&0\\0&0&1&1&1\end{pmatrix}}{\begin{pmatrix}3&2&1&0&0\\2&3&2&0&0\\1&2&3&0&0\\0&0&0&0&0\\0&0&0&0&0\end{pmatrix}}{\begin{pmatrix}1&0&0\\1&1&0\\1&1&1\\0&1&1\\0&0&1\end{pmatrix}}\\&={\frac {\sigma ^{2}}{81}}{\begin{pmatrix}19&16&10&4&1\\16&19&16&10&4\\10&16&19&16&10\\4&10&16&19&16\\1&4&10&16&19\end{pmatrix}}\end{aligned}}}$$ Nach zwei Durchläufen hat sich die Standardabweichung des Mittelwerts auf einen Wert von 0,58 σ verringert , verglichen mit 0,58 σ nach einem Durchlauf. Die Rauschreduzierung ist etwas geringer als bei einem Durchlauf eines gleitenden 5-Punkte-Mittelwerts, der unter denselben Bedingungen zu einer kleineren Standardabweichung der geglätteten Punkte von 0,45 σ führen würde . ${\displaystyle {\sqrt {\tfrac {19}{81}}}\sigma =0.48\sigma }$

Die Korrelation erstreckt sich nun über einen Bereich von 4 aufeinanderfolgenden Punkten mit Korrelationskoeffizienten.

$${\displaystyle {\begin{aligned}\rho _{i,i+1}&={\frac {16}{19}}=0.84,&\rho _{i,i+2}&={\frac {10}{19}}=0.53,\\[1ex]\rho _{i,i+3}&={\frac {4}{19}}=0.21,&\rho _{i,i+4}&={\frac {1}{19}}=0.05\end{aligned}}}$$

Der Vorteil, der sich durch die Durchführung von zwei Durchläufen mit der schmaleren Glättungsfunktion ergibt, besteht darin, dass dadurch weniger Verzerrungen in die berechneten Daten eingebracht werden.

Im Vergleich zu anderen Glättungsfiltern, z. B. der Faltung mit einer Gauß-Funktion oder der Mehrfachpass -Gleitmittelwertfilterung , weisen Savitzky-Golay-Filter im Frequenzbereich einen anfänglich flacheren Frequenzgang und eine schärfere Grenzfrequenz auf, insbesondere bei höheren Ordnungen des Anpassungspolynoms (siehe Frequenzcharakteristik ). Bei Daten mit begrenzter Signalbandbreite bedeutet dies, dass die Savitzky-Golay-Filterung ein besseres Signal-Rausch-Verhältnis als viele andere Filter liefern kann; so bleiben beispielsweise die Peak-Höhen der Spektren besser erhalten als bei anderen Filtern mit vergleichbarer Rauschunterdrückung. Nachteile der Savitzky-Golay-Filter sind die vergleichsweise schwache Unterdrückung einiger hoher Frequenzen (schlechte Sperrbereichsunterdrückung ) und Artefakte bei Verwendung von Polynomanpassungen für den ersten und letzten Messpunkt . ^{[ 16 ]}

Alternative Glättungsmethoden, die die Vorteile von Savitzky-Golay-Filtern nutzen und zumindest einige ihrer Nachteile ausgleichen, sind Savitzky-Golay-Filter mit geeignet gewählten alternativen Anpassungsgewichten , Whittaker-Henderson-Glättung und Hodrick-Prescott-Filter (äquivalente Methoden, die eng mit Glättungs- Splines verwandt sind ) sowie die Faltung mit einer Fenster -Sinc-Funktion . ^{[ 16 ]}

• sgolayfiltAus der Signal Processing Toolbox . Verfügbar seit Version vor R2006b.

• flattenaus dem Modul Lightkurve . Lightkurve ist die offizielle Bibliothek zur Analyse von Daten der Kepler- und TESS-Teleskope.
• scipy.signal.savgol_filteraus dem Modul SciPy . SciPy ist eine robuste Bibliothek, die in der akademischen Gemeinschaft für wissenschaftliches Rechnen weit verbreitet ist.

• Kernelglättung – Unterschiedliche Terminologie für viele der gleichen Prozesse, verwendet in der Statistik.
• Lokale Regression – die LOESS- und LOWESS-Methoden
• Numerische Differentiation – Anwendung auf die Differentiation von Funktionen
• Glättungs-Spline
• Stencil (Numerische Analysis) – Anwendung zur Lösung von Differentialgleichungen
• Hodrick-Prescott-Filter
• Kalman-Filter

Betrachten wir eine Menge von Datenpunkten .${\displaystyle (x_{j},y_{j})_{1\leq j\leq n}}$ Die Savitzky-Golay-Tabellen beziehen sich auf den Fall, dass die Schrittweite h konstant ist${\displaystyle x_{j}-x_{j-1}}$ . Beispiele für die Verwendung der sogenannten Faltungskoeffizienten mit einem kubischen Polynom und einer Fenstergröße m von 5 Punkten sind wie folgt.

Glättung

${\displaystyle Y_{j}={\frac {1}{35}}\left(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2}\right)}$ ;

1. Ableitung

${\displaystyle Y'_{j}={\frac {1}{12h}}\left(y_{j-2}-8y_{j-1}+0y_{j}+8y_{j+1}-y_{j+2}\right)}$ ;

2. Ableitung

${\displaystyle Y''_{j}={\frac {1}{7h^{2}}}\left(2y_{j-2}-y_{j-1}-2y_{j}-y_{j+1}+2y_{j+2}\right)}$Die

Ausgewählte Werte der Faltungskoeffizienten für Polynome vom Grad 1, 2, 3, 4 und 5 sind in den folgenden Tabellen angegeben ^{[ Anmerkung 7 ].} Die Werte wurden mit dem in Gorry bereitgestellten PASCAL-Code berechnet. ^{[ 15 ]}

• Gans, Peter (1992). Datenanpassung in den chemischen Wissenschaften: Mit der Methode der kleinsten Quadrate . Wiley. ISBN 9780471934127Die

• Erweiterter Faltungskoeffizientenrechner (ACCC) für mehrdimensionale Kleinste-Quadrate-Filter
• Savitzky-Golay-Filter in den Grundlagen der Statistik
• Ein breiteres Spektrum an Koeffizienten für verschiedene Datensatzgrößen, Anpassungsordnungen und Abweichungen vom Mittelpunkt