Savitzky-Golay-Filter vs. IIR- oder FIR-Linearfilter

• Ein traditionelles IIR / FIR-Filter (Tiefpass zum Entfernen der Hochfrequenzschwingungen), z. B. gleitender Durchschnitt,

• oder ein Savitzky-Golay- Filter

kann alle nützlich sein, um ein Signal zu glätten, z. B. ein Hüllkurvensignal:

Für welche Anwendung wäre ein Savitzky-Golay-Filter interessanter als ein klassischer Tiefpass?

Was unterscheidet es von einem Standardfilter und was fügt es im Vergleich zu Standardfiltern hinzu?

Passt es sich an die Eingabedaten an?

Ist es besser für die vorübergehende Konservierung?

---

Waren Sie eines Tages jemals in einer technischen Situation, als Sie beschlossen haben, "einen SG-Filter anstelle eines gleitenden Durchschnitts oder eines anderen FIR-Tiefpasses zu verwenden! Es ist besser, weil dies und das und das ..." ? Dann ist diese Frage für Sie!

Da sich die Diskussion in den vorhandenen Antworten und Kommentaren hauptsächlich auf die tatsächlichen Savitzky-Golay-Filter konzentriert hat (was sehr nützlich war), werde ich versuchen, die vorhandenen Antworten zu ergänzen, indem ich einige Informationen zur tatsächlichen Auswahl eines Glättungsfilters bereitstelle ist nach meinem Verständnis, worum es bei der Frage eigentlich geht.

Zunächst möchte ich wiederholen, was in der Diskussion deutlich geworden ist, die sich aus den anderen Antworten ergibt: die Kategorisierung der Glättungsfilter in der Frage in lineare und zeitinvariante (LTI) FIR / IIR-Filter einerseits und Savitzky-Golay-Filter sind dagegen irreführend. Ein Savitkzy-Golay-Filter ist nur ein Standard-FIR-Filter, der nach einem bestimmten Kriterium (lokale Polynomnäherung) entwickelt wurde. Alle in der Frage genannten Filter sind also LTI-Filter.

Die verbleibende Frage ist, wie ein Glättungsfilter ausgewählt wird. Wenn Rechenkomplexität und / oder Speicher ein Problem darstellen, können IIR-Filter gegenüber FIR-Filtern vorzuziehen sein, da sie typischerweise eine vergleichbare Rauschunterdrückung (dh Stoppbanddämpfung) mit einer viel niedrigeren Filterreihenfolge als FIR-Filter erzielen. Beachten Sie jedoch, dass ein möglicher Nachteil von IIR-Filtern darin besteht, dass sie keine genau lineare Phasenantwort haben können, wenn eine Echtzeitverarbeitung erforderlich ist. Das gewünschte Signal leidet also unter einigen Phasenverzerrungen. Bei der Offline-Verarbeitung können die Phasenverzerrungen auch bei IIR-Filtern durch Anwendung der Nullphasenfilterung vermieden werden .

Abgesehen von den im vorherigen Absatz erörterten Überlegungen ist hauptsächlich das Entwurfskriterium von Bedeutung, nicht so sehr, wenn es sich bei dem Filter um FIR oder IIR handelt, da jedes (stabile) IIR-Filter durch ein FIR-Filter und jedes beliebige mit beliebiger Genauigkeit angenähert werden kann Das FIR-Filter kann durch ein IIR-Filter angenähert werden, obwohl letzteres viel schwieriger sein kann. Das geeignete Entwurfskriterium hängt offensichtlich von den Eigenschaften der Daten und dem Rauschen ab. Wenn es um die Glättung geht, gehen wir normalerweise von ausreichend überabgetasteten (dh glättenden) Daten aus. Wenn das Rauschen hauptsächlich Hochfrequenzkomponenten aufweist, dh wenn es nur eine geringe spektrale Überlappung zwischen den Daten und dem Rauschen gibt, möchten wir die Stoppbanddämpfung maximieren oder die Stoppbandenergie minimieren, während das gewünschte Signal so gut wie möglich erhalten bleibt. In diesem Fall könnten wir ein linearphasiges FIR-Filter wählen, das nach einem Minimax-Kriterium unter Verwendung des Parks-McClellan-Algorithmus entworfen wurde. Wir könnten auch die Stoppbandenergie minimieren (dh die Rauschleistung im Stoppband minimieren), indem wir eine Methode der kleinsten Quadrate wählen. Eine Mischung zwischen den beiden Kriterien (Minimax und kleinste Quadrate) ist durch Auswahl von a möglicheingeschränktes Design der kleinsten Quadrate , das die Stoppbandenergie minimiert und gleichzeitig den maximalen Approximationsfehler im Durchlassband einschränkt.

Wenn sich das Rauschspektrum erheblich mit dem Signalspektrum überschneidet, ist ein vorsichtigerer Ansatz erforderlich, und die Brute-Force-Dämpfung funktioniert nicht gut, da Sie entweder zu viel Rauschen hinterlassen (indem Sie die Grenzfrequenz zu hoch wählen) oder das gewünschte verzerren Signal zu viel. In diesem Fall können Savitzky-Golay (SG) -Filter eine gute Wahl sein. Der zu zahlende Preis ist eine mittelmäßige Sperrdämpfung, aber ein Vorteil ist, dass einige Signaleigenschaften sehr gut erhalten bleiben. Dies hat mit der Tatsache zu tun, dass SG-Filter eine Flat-Pass-Band-Antwort haben, dh

$$\frac{d^kH(e^{j\omega})}{d\omega^k}{\huge |}_{\omega=0}=0\quad k=1,2,\ldots,r\tag{1}$$

Dabei ist die Reihenfolge des approximierenden Polynoms und der Frequenzgang des Filters. Die Eigenschaft garantiert, dass die ersten Momente des Eingangssignals im Ausgang erhalten bleiben, was bedeutet, dass die Breite und Höhe der Spitzen im gewünschten Signal gut erhalten bleiben.$r$$H(e^{j\omega})$$(1)$$r$

Natürlich gibt es auch einen Kompromiss zwischen den beiden oben diskutierten Arten von Glättungsfiltern (hohe Sperrbanddämpfung und SG). Wir könnten ein FIR-Filter mit einem bestimmten Grad an Ebenheit bei entwerfen und die verbleibenden Freiheitsgrade verwenden, um die Stoppbanddämpfung zu maximieren oder die Stoppbandenergie zu minimieren. Im Fall von FIR-Filtern ist das resultierende Entwurfsproblem ausreichend einfach (und konvex), und allgemeine Optimierungsroutinen, die in mehreren Softwarepaketen verfügbar sind, können verwendet werden, um den optimalen Filter für die gegebene Anwendung zu erhalten.$\omega=0$

Für diejenigen, die sich für die Theorie der SG-Filter interessieren, sind die wichtigsten Referenzen, die ich empfehlen kann, die folgenden:

• das Originalpapier von Savitzky und Golay
• ein Tutorial von Ronald W. Schafer
• Abschnitt 8.3.5 des DSP-Buches von Orfanidis

Wie bei allem sind manchmal bestimmte Werkzeuge besser als andere.

Filter mit gleitendem Durchschnitt (MA) können zum Glätten von Daten verwendet werden und sind FIR. Sie sind so ziemlich der einfachste Filter, den Sie sich vorstellen können, und sie eignen sich gut für viele Aufgaben, solange Sie nicht versuchen, plötzliche Sprünge oder Polynomtrends zu modellieren. Beachten Sie jedoch, dass es sich im Wesentlichen nur um Tiefpassfilter handelt. Sie funktionieren daher am besten, wenn die Daten, die Sie im Signal interessieren, niederfrequent und ziemlich schmalbandig sind.

Savitzky-Golay (SG) -Filter sind eine spezielle Gruppe von FIR-Filtern, die im Wesentlichen ein Polynom zu Ihrer Zeitreihe passen, wenn die Faltung entlang des Signals gleitet. SG-Filter sind nützlich für Signale, bei denen die Dinge, die Sie interessieren, nicht unbedingt niederfrequent und ziemlich schmalbandig sind.

Ich denke, Sie werden feststellen, dass wenn Sie die Wikipedia-Seite, die Sie verlinkt haben, ziemlich gründlich lesen, dies den Unterschied zwischen SG- und MA-Filtern sehr genau erklärt. Denken Sie jedoch daran: Letztendlich sind beide nur Werkzeuge, um ähnliche Dinge zu tun. Es liegt an Ihnen, herauszufinden, wann Sie das richtige Werkzeug verwenden müssen

EDIT :

Da es einige Verwirrung zu geben scheint, dass SG-Filter auf einer grundlegenden Ebene in irgendeiner Weise "adaptiv" sind, habe ich ein einfaches MATLAB-Beispiel beigefügt. Wie Dan betonte, können diese anpassungsfähig gemacht werden, ihre grundlegende Implementierung jedoch häufig nicht. Wenn Sie den Code überprüfen, werden Sie feststellen, dass dies einfach nur eine Matrix-Suche mit einer speziellen Behandlung ist. Nichts an diesem Filter ist im herkömmlichen Sinne "adaptiv". Sie wählen einfach eine Polynomanpassung und die Länge, über die dieses Polynom durch Faltung auf das Signal angepasst wird. SG ist wesentlich FIR. Das Skript, das ich unten habe, erzeugt diese Handlung:

Wenn Sie sich diese Abbildung ansehen, können Sie sehen, dass MA und SG im Wesentlichen dasselbe erreichen, jedoch mit einigen wichtigen Unterschieden:

1. Der MA kann das Rauschen hervorragend unterdrücken, aber die vorübergehenden Sprünge im Signal werden schlecht erfasst. Wir können noch mehr Rauschen unterdrücken, indem wir die Länge des Filters erhöhen, aber dann wird es bei den Transienten noch schlechter abschneiden. Dieser Effekt wird als "Verschmieren" in der Nähe von Transienten angesehen, die Sie in der gezeigten Abbildung sehen sollten.
2. Der SG erfasst die Transienten des Signals hervorragend, unterdrückt jedoch nicht so gut das Rauschen (zumindest im Vergleich zu einem MA gleicher Größe). Wir können die Rauschunterdrückung in der Nähe von Nicht-Transienten verbessern, indem wir die Rahmenlänge erhöhen. Dies führt jedoch zu einem Klingeln analog zum Gibb-Phänomen in der Nähe von Transienten.

Um ein besseres Verständnis der Funktionsweise dieser Filter zu erhalten, empfehlen wir Ihnen, den Code hier zu verwenden und zu bearbeiten und zu sehen, wie alle einzelnen Teile des SG-Filters funktionieren.

Code für das MATLAB-Beispiel:

% Generate a signal "s" that has square waves, and scale it with a
% polynomial of order 5
up = 1*ones(1,100);
down = zeros(1,150);
s = [down down up up down up down up down up up up down down down down down];
n = numel(s);
nn = linspace(0,4,numel(s));
s = s .* (nn .^5);
sn = (s + 4*randn(size(s))).';

% smooth it with SMA of length 15
sz = 15;
h = 1/sz * ones(1,sz);
sn_sma = conv(sn,h,'same');

% smooth it with sgolay of frame length 15
m = (sz-1)/2;
% look up the SG matrix for this order and size
B = sgolay(5, sz);

% compute the steady state response for the signal, i.e. everywhere that
% isnt the first or last "frame" for the SG filter
steady = conv(sn, B(m+1,:), 'same');
% handle the transiet portion at the start of the signal
y_st   = B(1:m,:)*sn(1:sz);
% handle the transient portion at the end of the signal
y_en   = B(sz -m+1 : sz, :) * sn(n - sz+1:n);

% combine our results
sn_sg  = steady;
sn_sg(1:m) = y_st;
sn_sg(n-m+1:n) = y_en;

% plots
figure(1);
hold off;
plot(sn,'Color',[0.75 0.75 0.75]);
hold on;
plot(sn_sma,'b');
plot(sn_sg,'r');

legend('Noisy Signal','MA Smoothing','SG Smoothing, order 5','Location','NorthWest');

HINWEIS

Meine vorherige Antwort (vor dieser Bearbeitung), in der der Savitzky-Golay (SG) -Filter als nichtlineares, zeitvariables Eingabedaten abhängig bezeichnet wurde, war falsch, da vorzeitig falsch interpretiert wurde, wie ein Savitzky-Golay (SG) -Filter seine Ausgabe berechnet gemäß dem bereitgestellten Wiki-Link. Jetzt korrigiere ich es zum Nutzen derer, die auch sehen würden, wie SG-Filter durch FIR-LTI-Filterung implementiert werden können. Danke an @MattL. für seine Korrektur, die großartige Verbindung, die er bereitgestellt hat, und die Geduld, die er (die ich niemals hätte zeigen können) während meiner Untersuchung des Problems hatte. Obwohl ich ehrlich gesagt einen ausführlicheren Einwand vorziehen würde, der dennoch eindeutig nicht notwendig ist. Bitte beachten Sie auch, dass die richtige Antwort die andere ist. Diese dient nur zur zusätzlichen Verdeutlichung der LTI-Eigenschaft von SG-Filtern.

Nun ist es nicht verwunderlich, dass jemand, der (der diese Filter noch nie zuvor verwendet hat) der Definition des SG-Filters als LSE-Polynom niedriger Ordnung für bestimmte Daten gegenübersteht, sofort zu dem Schluss kommt, dass diese datenabhängig, nichtlinear und sind zeitvariable, adaptive Filter.

Das Polynomanpassungsverfahren wird jedoch von SG selbst geschickt interpretiert, so dass es eine vollständig datenunabhängige, zeitinvariante, lineare Filterung ermöglicht, wodurch SG als festes LTI-FIR-Filter möglich wird.

Das ist eine kürzeste Zusammenfassung von der Verbindung von mattl zur Verfügung gestellt. Für Details, die zu fehlen scheinen, konsultieren Sie bitte das Originaldokument oder bitten Sie um Klärung. Aber ich möchte hier nicht das gesamte Dokument nachproduzieren.

Betrachten Sie nun Eingabedatenwerte die um zentriert sind und in die wir ein Polynom der Ordnung wollen , wobei die ganzzahligen :$2M+1$$x[-M],x[-M+1],...,x[0],x[1],...,x[M]$$n=0$$p[n]$$N$$n = -M,-M+1,...,-1,0,1,...M$

$$p[n] =\sum_{k=0}^{N} a_k n^k = a_0 + a_1 n + a_2 n^2 + ... + a_N n^N$$

Das klassische LSE-Polynomanpassungsverfahren berechnet diese Koeffizienten , um das optimale Polynom Ordnung herauszufinden, das die Summe der Fehlerquadrate minimiert$a_k$$N^{th}$$p[n]$

$$\mathcal{E} = \sum_{-M}^{M} (p[n]-x[n])^2$$

über den gegebenen Datenvektor .$x = \left[ x[-M],x[-M+1],...,x[0],x[1],...,x[M] \right]^{T}$

Diese optimalen Polynomkoeffizienten werden erhalten, indem die Ableitung von auf Null gesetzt wird:$a_k$$\mathcal{E}$

$$\frac{ \partial \mathcal{E} }{\partial a_i} = 0 ~~~,~~~ \text{for }~~~ i=0,1,..,N \tag{1}$$

Für diejenigen, die mit dem LSE-Polyfit-Verfahren vertraut sind, schreibe ich einfach die resultierende Matrixgleichung (aus dem Link), die den optimalen Koeffizientensatz definiert:

$$\boxed{ a = (A^{T}A)^{-1} A^{T} x = H x } \tag{2}$$ wobei der Eingangsdatenvektor ist die LSE-Polyfit-Matrix und die durch Matrix ist die Zeit-Momentan-Matrix (Potenzen von ganzzahligen Zeitpunkten ); das heißt, zu beachten , daß sowohl und folglich sind von Eingangsdaten indpendent Werte als gegeben ist durch :$x$$(2M+1) \times 1$$H$$2M+1$$N$$A$$n$$A$$H$$A$

$$A = \begin{bmatrix} \alpha_{n,i} \end{bmatrix} = \begin{bmatrix} (-M)^0 & (-M)^1 & ... & (-M)^N \\ (-M+1)^0 & (-M+1)^1 & ... & (-M+1)^N \\ &... \\ (0)^0 & (0)^1 & ... & (0)^N \\ &...\\ (M)^0 & (M)^1 & ... & (M)^N \\ \end{bmatrix}$$

Jetzt lehnen wir uns für einen Moment zurück und diskutieren hier einen Punkt.

Wie aus Gleichung (2) klar hervorgeht, sind die optimalen LSE-Polynomkoeffizienten von den Eingabedaten abhängig , obwohl und unabhängig von den Eingabedaten sind und nur von den Zeitindizes abhängen. Darüber hinaus hängen sie auch von der Größe des Fensters und der Reihenfolge des Polynoms . Wenn das Fenster entlang der Eingabedaten gleitet , sollte der Koeffizient neu berechnet (aktualisiert) werden und ist daher auch zeitabhängig. Genau so wird der SG-Filter im Link auf der Seite definiert, die unten aufgeführt ist:$A$$H$$n$$a_k$$M$$N$$x[n]$$a_k$$2^{nd}$

... Dies (das LSE-Polyfit) kann bei jedem Abtastwert der Eingabe wiederholt werden, wobei jedes Mal ein neues Polynom und ein neuer Wert der Ausgabesequenz y [n] erzeugt werden ...

Wie überwinden wir diese rätselhafte Überraschung? Durch Interpretieren und Definieren der SG-Filterausgabe wie folgt:

Das SG-Filter der Ordnung akzeptiert für jede Zeitabtastung einen Eingabesatz und erzeugt eine einzelne Ausgangsabtastung , die als das bei bewertete Polynom definiert ist ; dh$N$$n$$x[n]$$y[n]$$p[n]$$n=0$

$$\boxed{ y[n] = y[0] = \sum_{m=0}^{N} a_m n^m = a_0 }$$

Das heißt, für jeden Eingangssatz von Abtastwerten von (zentriert um ) erzeugt das SG-Filter die mit bezeichnete Ausgabe, die dem einzelnen Koeffizienten der optimalen LSE entspricht Polynom , das diesem bestimmten Fenster der Abtastwerte . Während das Fenster entlang der Eingangsdatenlänge gleitet, wird übrigens zu jedem Zeitpunkt eine neue Ausgangsabtastung gemäß dem Fenster der Abtastwerte berechnet. . Hier ist dies ein Nancusal-Filter.$2M+1$$x[n]$$n=d$$y[n]$$a_0$$p[n]$$x[n]$$n=d$$y[d]$$x[d-M],x[d-M+1],...,x[d-1],x[d],x[d+1],...x[d+M]$

Es ist nun an der Zeit zu zeigen, dass der Koeffizient als lineare Kombination von Eingangssignalwerten im Fenster erhalten wird und der Filterausgang daher eine lineare Kombination von Eingangswerten . Und dies ist genau die Definition der (LTI) -Faltung über ein FIR-Filter; Die Ausgabe zum Zeitpunkt ist die lineare Kombination ihrer Eingabe und der Filterkoeffizienten . Aber was sind dann die Filterkoeffizienten für diesen SG-Filter? Mal sehn.$a_0$$x[n]$$y[n]$$x[n]$$n$$x[n]$$h[n]$

Betrachten Sie noch einmal die Berechnung von :$a_k$

$$a = H x$$

$$\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_N \\ \end{bmatrix} = \begin{bmatrix} h(0,0) & h(0,1) & ... & h(0,2M) \\ h(1,0) & h(1,1) & ... & h(1,2M) \\ &... \\ h(N,0) & h(0,1) & ... & h(0,2M) \\ \end{bmatrix} \cdot \begin{bmatrix} x[-M] \\ x[-M+1] \\ ... \\ x[M] \\ \end{bmatrix}$$

woraus wir sehen können, dass der einzelne Koeffizient durch das Punktprodukt der ersten Zeile von mit dem Eingabedatenvektor ; dh$a_0$$H$$x$

$$a_0 = H(0,n) \cdot x = \sum H(0,k) x[k] = H(0,-n) \star x[n]$$

wo in der letzten Gleichheit wir die Punktproduktsumme als die Faltungssumme interpretiert haben, indem wir die Impulsantwort des SG-Filters als ,

$$\boxed{ h[n] = H(0,-n) }$$

Insbesondere ist es die Impulsantwort des SG-Filters der Ordnung mit einer Fensterlänge von .$N$$2M+1$

Und die vollständige Ausgabe des SG-Filters N-ter Ordnung mit einer Fenstergröße von für eine Eingabe der Länge wird durch eine einzelne LTI-Faltung mit der Impulsantwort$y[n]$$2M+1$$x[n]$$L$$h_N[n]$

$$\boxed{ y[n] = x[n] \star h_N[n] }$$

KOMMENTAR

Die Tatsache, dass die Polynomkoeffizienten von den Eingabedaten abhängen, verhindert nicht, dass der Filter eine LTI-FIR ist. Weil eine Impulsantwort definiert werden kann, um den Ausgang darzustellen , der aus einer linearen Kombination von Eingangsabtastwerten berechnet werden soll. Die linearen Kombinationen von Eingangsabtastwerten werden von Natur aus durch das Matrixprodukt impliziert , das die optimalen Koeffizienten von . Daher würde jede lineare Kombination von auch dazu führen, dass ein FIR-LTI-Filter dargestellt wird das LSE-Polynomanpassungsverfahren.$a_k$$h[n]$$y[n]$$x$$a = H x$$a_k$$p[n]$$a_k$$h[n]$

MATLAB / OKTVE-CODE

Das folgende einfache MATLAB / OCTAVE kann verwendet werden, um diese SG-Filter-Impulsantworten zu berechnen (Beachten Sie, dass der eingebaute SG-Designer möglicherweise einen anderen Satz von wie im verknüpften PDF beschrieben).$h[n]$$h[n]$

% Savitzky-Golay Filter
% 
clc; clear all; close all;

N = 3;                      % a0,a1,a2,a3 : 3rd order polynomial
M = 4;                      % x[-M],..x[M] . 2M + 1 data

A = zeros(2*M+1,N+1);
for n = -M:M
    A(n+M+1,:) = n.^[0:N];
end

H = (A'*A)^(-1)* A';        % LSE fit matrix

h = H(1,:);                 % S-G filter impulse response (nancausal symmetric FIR)

figure,subplot(2,1,1)
stem([-M:M],h);
title(['Impulse response h[n] of Savitzky-Golay filter of order N = ' num2str(N), ' and window size 2M+1 =  ' , num2str(2*M+1)]);

subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('Frequency response magnitude of h[n]');

Die Ausgabe ist:

Hoffe, das klärt das Problem.