Filtr Savitzky'ego – Golaya vs. filtr liniowy IIR lub FIR

• Tradycyjny filtr IIR / FIR (dolnoprzepustowy w celu usunięcia oscylacji o wysokiej częstotliwości), np. Średnia ruchoma,

• lub filtr Savitzky-Golay

mogą być przydatne do wygładzenia sygnału, takiego jak sygnał obwiedni:

W jakim zastosowaniu filtr Savitzky-Golay byłby bardziej interesujący niż klasyczny dolnoprzepustowy?

Czym różni się od standardowego filtra i co dodaje w porównaniu ze standardowymi filtrami?

Czy dostosowuje się do danych wejściowych?

Czy lepiej jest zachować przejściowe zachowanie?

---

Czy byłeś kiedyś w sytuacji inżynierskiej, kiedy zdecydowałeś „Użyjmy filtra SG zamiast średniej ruchomej lub innego dolnoprzepustowego FIR! Lepiej, ponieważ to i to i to ...” ? To pytanie jest do Ciebie!

Ponieważ dyskusja w istniejących odpowiedziach i komentarzach koncentrowała się głównie na tym, czym faktycznie są filtry Savitzky-Golay (co było bardzo przydatne), postaram się dodać do istniejących odpowiedzi, podając informacje o tym, jak faktycznie wybrać filtr wygładzający, który jest, w moim rozumieniu, o co właściwie chodzi w tym pytaniu.

Przede wszystkim chciałbym powtórzyć to, co stało się jasne w dyskusji, powstałe z innych odpowiedzi: kategoryzację filtrów wygładzających w pytaniu na filtry liniowe i niezmienne czasowe (LTI) FIR / IIR z jednej strony, oraz Z drugiej strony filtry Savitzky-Golay są mylące. Filtr Savitkzy-Golay jest tylko standardowym filtrem FIR zaprojektowanym według określonego kryterium (lokalne przybliżenie wielomianowe). Więc wszystkie filtry wymienione w pytaniu są filtrami LTI.

Pozostaje pytanie, jak wybrać filtr wygładzający. Jeśli złożoność obliczeniowa i / lub pamięć stanowią problem, filtry IIR mogą być lepsze niż filtry FIR, ponieważ zazwyczaj osiągają porównywalne tłumienie szumów (tj. Tłumienie pasma zatrzymania) przy znacznie niższej kolejności filtrów niż filtry FIR. Należy jednak pamiętać, że jeśli konieczne jest przetwarzanie w czasie rzeczywistym, jedną z możliwych wad filtrów IIR jest to, że nie mogą one mieć dokładnie liniowej odpowiedzi fazowej. Zatem pożądany sygnał będzie podlegał pewnym zniekształceniom fazowym. W przypadku przetwarzania offline można uniknąć zniekształceń fazowych, nawet w przypadku filtrów IIR, stosując filtrowanie w fazie zerowej .

Oprócz rozważań omówionych w poprzednim akapicie, znaczenie ma przede wszystkim kryterium projektowe, nie tyle jeśli filtr jest FIR lub IIR, ponieważ każdy (stabilny) filtr IIR może być aproksymowany z dowolną dokładnością przez filtr FIR i dowolne Filtr FIR może być aproksymowany filtrem IIR, nawet jeśli ten ostatni może być znacznie trudniejszy. Odpowiednie kryterium projektowe oczywiście zależy od właściwości danych i hałasu. Jeśli chodzi o wygładzanie, zwykle zakładamy, że dane są wystarczająco nadpróbkowane (tj. Gładkie). Jeśli szum ma głównie składowe o wysokiej częstotliwości, tj. Jeśli dane w niewielkim stopniu pokrywają się z widmem, chcemy zmaksymalizować tłumienie pasma zatrzymania lub zminimalizować energię pasma zatrzymania, zachowując jednocześnie pożądany sygnał. W tym przypadku moglibyśmy wybrać filtr FIR z fazą liniową zaprojektowany zgodnie z kryterium minimax przy użyciu algorytmu Parks-McClellan. Możemy również zminimalizować energię pasma zatrzymania (tj. Zminimalizować moc szumu w paśmie zatrzymania), wybierając metodę najmniejszych kwadratów. Połączenie dwóch kryteriów (minimax i najmniejszych kwadratów) jest możliwe, wybierając aograniczona konstrukcja najmniejszych kwadratów , która minimalizuje energię pasma zatrzymania, jednocześnie ograniczając maksymalny błąd aproksymacji w paśmie przejścia.

Jeśli spektrum szumów w znacznym stopniu pokrywa się z widmem sygnału, konieczne jest bardziej ostrożne podejście, a tłumienie siły brutalnej nie będzie działać dobrze, ponieważ albo zostawisz zbyt dużo hałasu (wybierając zbyt wysoką częstotliwość odcięcia), albo zniekształcisz pożądane za dużo sygnału. W takim przypadku dobrym wyborem mogą być filtry Savitzky-Golay (SG). Cena do zapłaty to mierne tłumienie pasma, ale jedną zaletą jest to, że niektóre właściwości sygnału są zachowane bardzo dobrze. Ma to związek z faktem, że filtry SG mają płaską odpowiedź pasmową, tj.

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

gdzie jest rzędem aproksymującego wielomianu, a jest odpowiedzią częstotliwościową filtra. Właściwość gwarantuje, że pierwszy momenty sygnałem wejściowym zachowany na wyjściu, co oznacza, że szerokość i wysokość pików żądanego sygnału, są dobrze zakonserwowane.$r$$H(e^{j\omega})$$(1)$$r$

Oczywiście istnieje również kompromis między dwoma rodzajami filtrów wygładzających omówionych powyżej (wysokie tłumienie pasma i SG). Możemy zaprojektować filtr FIR o pewnym stopniu płaskości przy i wykorzystać pozostałe stopnie swobody, aby zmaksymalizować tłumienie pasma zatrzymania lub zminimalizować energię pasma zatrzymania. W przypadku filtrów FIR wynikający z tego problem projektowy jest wystarczająco prosty (i wypukły), a ogólne procedury optymalizacji dostępne w kilku pakietach oprogramowania można wykorzystać do uzyskania optymalnego filtra dla danej aplikacji.$\omega=0$

Dla osób zainteresowanych teorią filtrów SG najbardziej odpowiednie referencje, które mogę polecić, to:

• oryginalny artykuł Savitzky'ego i Golaya
• artykuł szkoleniowy Ronalda W. Schafera
• sekcja 8.3.5 książki DSP Orfanidisa

Jak zawsze, niektóre narzędzia są lepsze od innych.

Filtry średniej ruchomej (MA) mogą być używane do wygładzania danych i są FIR. Są one najprostszym filtrem, jaki możesz wymyślić, i sprawdzają się w wielu zadaniach, o ile nie próbujesz modelować nagłych skoków lub trendów wielomianowych. Pamiętaj jednak, że są to po prostu filtry dolnoprzepustowe, więc działają najlepiej, gdy dane, na których Ci zależy w sygnale, mają niską częstotliwość i dość wąskie pasmo.

Filtry Savitzky-Golay (SG) to specjalna grupa filtrów FIR, które zasadniczo dopasowują wielomian do twoich szeregów czasowych, gdy splot przesuwa się wzdłuż sygnału. Filtry SG są przydatne w przypadku sygnałów, w których ważne są niekoniecznie niskie częstotliwości i dość wąskie pasmo.

Myślę, że przekonasz się, że jeśli dokładnie przeczytasz stronę Wikipedii, do której odsyłałeś, to dość rygorystycznie wyjaśnia różnicę między filtrami SG i MA. Pamiętaj jednak: w końcu oba są tylko narzędziami do robienia podobnych rzeczy, od Ciebie zależy, kiedy wybierzesz odpowiednie narzędzie

EDYCJA :

Ponieważ wydaje się, że istnieje pewne zamieszanie, że filtry SG są w jakiś sposób „adaptacyjne” na poziomie podstawowym, zamieściłem prosty przykład MATLAB. Jak zauważył Dan, można je dostosować, ale ich podstawowe wdrożenie często nie jest. Sprawdzając kod, zobaczysz, że jest to po prostu wyszukiwanie matrycy ze specjalną obsługą. Nic w tym filtrze nie jest „adaptacyjne” w tradycyjnym znaczeniu, po prostu wybierasz dopasowanie wielomianowe i długość, na której ten wielomian będzie dopasowany do sygnału poprzez splot; SG jest niezbędnym FIR. Skrypt, który mam poniżej, tworzy tę fabułę:

Patrząc na tę figurę, można zauważyć, że MA i SG zasadniczo osiągają to samo, ale z kilkoma ważnymi różnicami:

1. MA świetnie sobie radzi z tłumieniem hałasu, ale słabo radzi sobie z przechwytywaniem przejściowych skoków sygnału. Możemy tłumić jeszcze więcej szumów, zwiększając długość filtra, ale wtedy będzie on działał jeszcze gorzej na stanach nieustalonych; efekt ten będzie postrzegany jako „rozmazywanie” w pobliżu wszelkich stanów przejściowych, które powinieneś zobaczyć na pokazanym rysunku.
2. SG świetnie sobie radzi z przechwytywaniem stanów nieustalonych sygnału, ale nie tak dobrze tłumi szumy (przynajmniej w porównaniu z MA tej samej wielkości). Możemy poprawić tłumienie szumów w pobliżu stanów nieustalonych przez zwiększenie długości ramki, ale wprowadzi to dzwonienie analogiczne do zjawiska Gibb'a w pobliżu jakichkolwiek stanów nieustalonych.

Aby lepiej zrozumieć działanie tych filtrów, zachęcam do wzięcia kodu tutaj i manipulowania nim oraz sprawdzenia, jak działają wszystkie poszczególne części filtra SG.

Kod dla przykładu MATLAB:

% 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');

UWAGA

moja poprzednia odpowiedź (przed tą edycją) oznaczająca filtr Savitzky-Golay (SG) jako nieliniową zależną od danych wejściowych zmienną w czasie była błędna z powodu przedwczesnej błędnej interpretacji sposobu, w jaki filtr Savitzky-Golay (SG) oblicza swoje dane wyjściowe zgodnie z podanym linkiem wiki. Więc teraz poprawiam to z korzyścią dla tych, którzy chcieliby zobaczyć, w jaki sposób filtry SG są implementowane przez filtrowanie FIR-LTI. Dzięki @MattL. za jego poprawkę, wielki link, który podał, i cierpliwość, której miał (czego nigdy nie mogłem wykazać) podczas mojego dochodzenia w tej sprawie. Chociaż szczerze wolę bardziej szczegółowy sprzeciw, który jednak nie jest konieczny. Należy również pamiętać, że poprawna odpowiedź to druga, ta służy jedynie dodatkowemu wyjaśnieniu właściwości LTI filtrów SG.

Nic dziwnego, że gdy ktoś (kto nigdy wcześniej nie korzystał z tych filtrów) stoi przed definicją filtra SG jako wielomianowego dopasowania LSE niskiego rzędu do danych , od razu doszedłby do wniosku, że są one zależne od danych, nieliniowe i zmienne w czasie (przesunięcie), filtry adaptacyjne.

Jednak procedura dopasowania wielomianu jest sprytnie interpretowana przez samych SG, tak że umożliwia całkowicie niezależne od danych, niezmienne w czasie, liniowe filtrowanie, dzięki czemu SG jest stałym filtrem LTI-FIR.

Poniżej znajduje się najkrótsze streszczenie z linku dostarczonego przez MattL. W celu uzyskania jakichkolwiek brakujących informacji zapoznaj się z oryginalnym dokumentem lub poproś o wyjaśnienie. Ale nie chciałbym tutaj reprodukować całego dokumentu.

Teraz rozważmy wartości danych wejściowych które są wyśrodkowane wokół i do którego chcemy dopasować wielomian rzędu , przy czym jest wskaźnikami czasu całkowitego:$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$$

Klasyczna procedura dopasowania wielomianu LSE oblicza te współczynniki aby znaleźć optymalny wielomian rzędu , który minimalizuje sumę kwadratów błędów$a_k$$N^{th}$$p[n]$

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

nad podanym wektorem danych .$x = \left[ x[-M],x[-M+1],...,x[0],x[1],...,x[M] \right]^{T}$

Te optymalne współczynniki wielomianowe są uzyskiwane przez ustawienie pochodnej na zero:$a_k$$\mathcal{E}$

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

Teraz dla tych, którzy są zaznajomieni z procedurą LIT Polyfit, po prostu napiszę wynikowe równanie macierzowe (z łącza), które definiuje optymalny zestaw współczynników:

$$\boxed{ a = (A^{T}A)^{-1} A^{T} x = H x } \tag{2}$$ gdzie to wektor danych wejściowych, macierz LSE polyfit i o macierzy jest macierzą Chwilę (moce momentach czasowych całkowitych ); tzn. należy zauważyć, że zarówno jak i są niezależne od wartości danych wejściowych, ponieważ jest podane przez:$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}$$

Teraz oprzyjmy się na chwilę i przedyskutujmy tutaj punkt.

Jak wyraźnie wskazuje równanie (2), chociaż i są niezależne od danych wejściowych i zależą tylko od wskaźników czasu , optymalne współczynniki wielomianowe LSE zależą od danych wejściowych. Ponadto, także zależy od wielkości okna i porządku wielomianu . Ponadto, ponieważ okno przesuwa się wzdłuż danych wejściowych , współczynnik powinien zostać ponownie obliczony (zaktualizowany), a zatem będzie również zależny od czasu. Dokładnie tak zdefiniowano filtr SG w linku na stronie poniżej:$A$$H$$n$$a_k$$M$$N$$x[n]$$a_k$$2^{nd}$

... To (polifunkcja LSE) można powtórzyć dla każdej próbki danych wejściowych, za każdym razem wytwarzając nowy wielomian i nową wartość sekwencji wyjściowej y [n] ...

Jak więc pokonać tę zagadkową niespodziankę? Poprzez interpretację i zdefiniowanie wyjściowego filtru SG jako:

Filtr SG rzędu , dla każdej próbki , przyjmuje zestaw wejściowy i wytwarza pojedynczą próbkę wyjściową zdefiniowaną jako wielomian obliczony przy ; to znaczy,$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 }$$

Oznacza to, że dla każdego zestawu wejściowego próbek (wyśrodkowanych wokół ) filtr SG wytwarza wynik oznaczony przez , który jest równoważny pojedynczemu współczynnikowi optymalnego LSE wielomian związany z tym konkretnym oknem próbek . Nawiasem mówiąc, gdy okno przesuwa się wzdłuż długości danych wejściowych, za każdym razem nowa próbka wyjściowa jest obliczana zgodnie z oknem próbek . Tutaj jest to filtr trzustkowy.$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]$

Nadszedł czas, aby pokazać, że współczynnik jest uzyskiwany jako liniowa kombinacja wartości sygnału wejściowego w oknie , a wyjście filtra jest zatem liniową kombinacją wartości wejściowych . I to jest właśnie definicja splotu (LTI) przez filtr FIR; wyjście w czasie jest liniową kombinacją jego wejścia i współczynników filtra . Ale jakie są współczynniki filtra dla tego filtra SG? Zobaczmy.$a_0$$x[n]$$y[n]$$x[n]$$n$$x[n]$$h[n]$

Ponownie rozważ obliczenia :$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}$$

z którego możemy zobaczyć, że pojedynczy współczynnik jest iloczynem iloczynu pierwszego wiersza z wektorem danych wejściowych ; to znaczy,$a_0$$H$$x$

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

gdzie w ostatniej równości zinterpretowaliśmy sumę iloczynu kropkowego jako sumę splotu, biorąc pod uwagę odpowiedź impulsową filtra SG na ,

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

dokładniej jest to odpowiedź impulsowa filtra SG rzędu o długości okna .$N$$2M+1$

I całkowite wyjście filtra SG N-tego rzędu o rozmiarze okna , dla wejścia o długości jest uzyskiwane przez pojedynczy splot LTI z odpowiedzią impulsową$y[n]$$2M+1$$x[n]$$L$$h_N[n]$

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

KOMENTARZ

Fakt, że współczynniki wielomianowe zależą od danych wejściowych, nie wyklucza, że ​​filtr może być FIR LTI. Ponieważ odpowiedź impulsowa może być zdefiniowana jako reprezentująca wynik do obliczenia z liniowych kombinacji próbek wejściowych. Liniowe kombinacje próbki sygnału są nieodłącznie dorozumiany przez iloczyn macierzy , który definiuje współczynniki optymalne z , a więc kombinacja liniowa również prowadzić w filtrze FIR LTI reprezentuje procedura dopasowania wielomianu LSE.$a_k$$h[n]$$y[n]$$x$$a = H x$$a_k$$p[n]$$a_k$$h[n]$

KOD MATLAB / OCTVE

Poniższego prostego MATLAB / OCTAVE można użyć do obliczenia odpowiedzi impulsowych filtra SG (Należy pamiętać, że jego wbudowany projektant SG może generować inny zestaw zgodnie z linkiem pdf)$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]');

Dane wyjściowe to:

Mam nadzieję, że to wyjaśnia problem.