Czy implementacja filtrowania jako zwielokrotnienia lub splotu jest bardziej stabilna numerycznie?

Piszę program do filtrowania w trybie offline sygnału 20 000 próbek przy użyciu filtra Butterwortha piątego rzędu. Mam dostęp do implementacji FFT. Wydaje się, że istnieją dwie możliwości wdrożenia filtrowania:

• splot sygnału z odpowiedzią impulsową w dziedzinie czasu, lub
• zwielokrotnienie sygnału przez odpowiedź impulsową w dziedzinie częstotliwości i odwrotne przekształcenie wyniku

Metody te byłyby identyczne w teoretycznym przypadku FT. Robię to w prawdziwym życiu z DFT, ale chyba wszystko jest inne. Czy jedna z metod jest bardziej stabilna numerycznie? Czy są jeszcze jakieś problemy, o których powinienem wiedzieć? Liczba obliczeń nie jest ważna.

Dzięki splotowi nie napotkasz żadnych problemów ze stabilnością, ponieważ nie ma filtrowania rekurencyjnego, więc nie będziesz kumulować żadnych błędów. Innymi słowy, w systemie są wszystkie zera, bez biegunów. Słyszałem anegdotycznie, ale nie zweryfikowałem dla siebie, że splot oparty na FFT ma mniejszy błąd niż splot w dziedzinie czasu, po prostu dlatego, że ma operacje arytmetyczne O (n log n) zamiast O (n ^ 2).

Zazwyczaj, o ile mi wiadomo, filtry Butterworth są implementowane jako filtry rekurencyjne (IIR), więc jest to inny temat. Filtry IIR mają bieguny i zera, więc w praktyce mogą wystąpić problemy ze stabilnością. Ponadto, w przypadku filtrów IIR, metody oparte na FFT nie są opcją, ale z drugiej strony, filtry IIR wydają się być bardzo niskiego rzędu.

Jeśli chodzi o problemy ze stabilnością z filtrami IIR, mają one zwykle problemy przy wyższych rzędach - po prostu wyrzucę liczbę i powiem, że popycha ją około 6. rzędu. Zamiast tego są zwykle implementowane jako kaskadowe biquady (sekcje filtrów drugiego rzędu). Dla filtra piątego rzędu napisz go jako funkcję transferu domeny z (będzie to funkcja racjonalna 5. stopnia), a następnie podziel go na 5 biegunów i 5 zer. Zbierz złożone koniugaty, a otrzymasz dwie bikwady i jeden filtr pierwszego rzędu. Ogólnie problemy ze stabilnością pojawiają się, gdy bieguny zbliżają się do okręgu jednostki.

Mogą występować również problemy z szumem i cyklem granicznym w filtrach IIR, więc istnieją różne topologie filtrów (tj. Bezpośrednia forma I, bezpośrednia forma II), które mają różne właściwości liczbowe, ale nie zastanowiłbym się nad tym - po prostu użyj podwójnego- precyzja i prawie na pewno będzie wystarczająco dobra.

W prawie wszystkich przypadkach najlepszym wyborem nie jest ani splot, ani FFT, ale po prostu zastosowanie filtra IIR bezpośrednio (za pomocą np. Funkcji sosfilt ()). Będzie to znacznie bardziej wydajne pod względem zużycia procesora i pamięci.

To, czy ma to znaczenie liczbowe, zależy od konkretnego filtra. Jedynym przypadkiem, w którym może się wkraść pewna różnica, jest to, że bieguny są bardzo, bardzo blisko koła jednostki. Nawet tam jest kilka sztuczek, które mogą pomóc. NIE UŻYWAJ reprezentacji funkcji przesyłania i filter (), ale używaj biegunów i zer z sosfilt (). Oto przykład różnicy.

n = 2^16;  % filter length
fs = 44100; % sample rate
x = zeros(n,1); x(1) = 1;
f0 = 15; % cutoff frequency in Hz
% design with poles and zeroes
[z,p,k] = butter(5,f0*2/fs);
clf
plot(sosfilt(zp2sos(z,p,k),x));
% design with transfer function
[b,a] = butter(5,f0*2/fs);
hold on
plot(filter(b,a,x),'k');

filter () psuje się przy odcięciu około 15 Hz przy 44,1 kHz. W przypadku sosfilt () odcięcie może być znacznie poniżej 1/100 Hz przy 44,1 kHz bez żadnych problemów.

JEŚLI masz problemy ze stabilnością, FFT też niewiele pomaga. Ponieważ twój filtr jest filtrem IIR, odpowiedź impulsowa jest nieskończona i musiałaby zostać najpierw obcięta. Przy tak bardzo niskiej częstotliwości odpowiedź impulsowa staje się tak długa, że ​​FFT również staje się niepraktyczne.

Na przykład, jeśli chcesz mieć wartość graniczną 1/100 Hz przy 44,1 kHz i chcesz zakres dynamiczny w odpowiedzi impulsowej 100 dB, potrzebujesz około 25 milionów próbek !!! To prawie 10 minut przy 44,1 kHz i wiele, wiele razy dłużej niż oryginalny sygnał

Jak myślisz, dlaczego będzie inaczej? Koncepcje teoretyczne powinny przełożyć się na praktyczne zastosowania, z tą różnicą, że problemami zmiennoprzecinkowymi są problemy, których nie możemy uniknąć. Możesz łatwo zweryfikować za pomocą prostego przykładu w MATLAB:

x=randn(5,1);
y=randn(5,1);
X=fft(x,length(x)+length(y)-1);
Y=fft(y,length(x)+length(y)-1);

z1=conv(x,y);z2=ifft(X.*Y);
z1-z2

ans =

   1.0e-15 *

   -0.4441
   -0.6661
         0
   -0.2220
    0.8882
   -0.2220
         0
   -0.4441
    0.8882

Jak widać, błędy są rzędu dokładności maszyny. Nie powinno być żadnego powodu, aby nie używać metody FFT. Jak wspomniano Endolith, znacznie częściej stosuje się metodę FFT do filtrowania / zwijania / korelacji krzyżowej itp. I znacznie szybciej, z wyjątkiem bardzo małych próbek (jak w tym przykładzie). Nie chodzi o to, że przetwarzanie w dziedzinie czasu nigdy nie jest wykonywane ... wszystko sprowadza się do aplikacji, potrzeb i ograniczeń.