%matplotlib inline
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
from util import set_style
set_style()
# częstotliwość próbkowania
fs = 48000
Filtry o skończonej odpowiedzi impulsowej (FIR)¶
Filtry są to układy lub algorytmy, które modyfikują sygnał. Filtr cyfrowy to algorytm, który na podstawie pewnej deterministycznej funkcji przekształca wejściowy ciąg próbek w zmodyfikowany ciąg wyjściowy. Najczęściej stosowane są filtry przepustowo-zaporowe, których zadaniem jest przepuszczenie pewnych zakresów częstotliwości sygnału i wytłumienie pozostałych. Istnieją jednak również inne typy filtrów, np. wszechprzepustowe, Hilberta, itp.
Filtry o skończonej odpowiedzi impulsowej (SOI), określane angielskim skrótem FIR (Finite Impulse Response) są cyfrową realizacją filtrów, nie mają one odpowiednika w technice analogowej. Są to bardzo proste algorytmy. Zasadę działania filtrów FIR można opisać jednym zdaniem: $L$ ostatnich próbek sygnału jest mnożonych przez współczynniki filtru (filter coefficients lub taps), a wyniki mnożenia są sumowane. Równanie różnicowe filtru FIR jest następujące:
$$y(n) = b_0 x(n) + b_1 x(n-1) + b_2 x(n-2) + ... + b_N x(n-N)$$
Filtr FIR oblicza więc ważoną sumę pewnej liczby ostatnich próbek sygnału, ze współczynnikami $b_i$ jako wagami sumowania. Odpowiedź impulsowa filtru FIR (impulse response) jest po prostu zbiorem współczynników $b_i$. Transformata Fouriera odpowiedzi impulsowej jest charakterystyką częstotliwościową filtru (frequency response). Transmitancja filtru FIR w dziedzinie $z$ wyraża się wzorem:
$$H(z) =\sum_{i=0}^{N} h(n) z^{-i} = \sum_{i=0}^{N} b_i z^{-i}$$
Wartość $N$ nazywa się rzędem filtru (filter order). Rząd filtru FIR jest zawsze o jeden mniejszy niż liczba współczynników.
Celem projektowania filtru FIR (filter design) jest obliczenie współczynników $b_i$ takich, aby uzyskać pożądany kształt charakterystyki częstotliwościowej.
Filtr dolnoprzepustowy i metoda okienkowania¶
Jako pierwszy rozpatrzymy filtr dolnoprzepustowy (DP) (high pass, HP). Filtr ten przepuszcza składowe częstotliwościowe do ustalonej częstotliwości granicznej ($f_c$), nazywanej też częstotliwością odcięcia (cutoff frequency). Do częstotliwości granicznej mamy pasmo przepustowe (pass band), powyżej - pasmo zaporowe (stop band). Idealna charakterystyka filtru dolnoprzepustowego o częstotliwości granicznej 1 kHz wygląda tak jak na poniższym wykresie.
char = np.zeros(5001)
f = np.arange(5001)
char[f<=1000] = 1
plt.plot(f, char)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('wzmocnienie')
plt.title('Charakterystyka idealnego filtru dolnoprzepustowego')
plt.show()
Jeżeli obliczymy odwrotne przekształcenie Fouriera powyższej charakterystyki, otrzymamy odpowiedź impulsową filtru. Nie da się jednak wykonać takiego filtru (formalne określenie: filtr nierealizowalny), z dwóch powodów:
- liczba współczynników filtru jest nieskończona,
- odpowiedź impulsowa filtru jest niezerowa dla ujemnych wartości czasu, czyli zanim wprowadzone zostaje pobudzenie; taki filtr jest nieprzyczynowy (non-casual).
Pierwszy problem można rozwiązać wycinając wybraną liczbę współczynników wokół zerowej wartości czasu. Np. biorąc współczynniki od $n=-15$ do $n=15$ dostaniemy 31 współczynników filtru (jeżeli w punkcie zerowym też jest współczynnik). Obcinamy więc nieskończoną odpowiedź filtru do skończonej liczby współczynników. Mówiąc inaczej, nakładamy okno czasowe na odpowiedź impulsową filtru, centralnie wokół zerowego punktu na osi czasu. Stąd nazwa tej metody projektowania filtru: metoda okienkowania (windowing). Podobnie jak w przypadku analizy częstotliwościowej, należy używać specjalnych funkcji okna, aby ograniczyć zafalowania charakterystyki (listki boczne). Możemy zastosować te same funkcje, co w przypadku FFT, a więc okno Hamminga lub Blackmana. Istnieją również specyficzne typy okien do projektowania filtrów, np. okno Kaisera, które posiada parametr umożliwiający regulowanie stopnia zafalowania charakterystyki w wybranym paśmie częstotliwości.
Problem nieprzyczynowości filtru rozwiązujemy przesuwając okno o połowę jego długości, tak że pierwsza wartość odpowiedzi impulsowej znajduje się w punkcie zerowym. To przesunięcie ma praktyczne znaczenie dla działania filtru: wyznacza ono opóźnienie filtru.
W celu zaprojektowania filtru metodą okienkową w języku Python musimy sięgnąć po moduł scipy.signal, który w poniższych przykładach został zaimportowany pod skrótem sig. Do projektowania filtru użyjemy funkcji firwin. Załóżmy długość filtru $L=30$ i częstotliwość graniczną 1 kHz.
dp = sig.firwin(30, 1000, window='hamming', pass_zero=True, fs=fs)
plt.stem(dp)
plt.xlabel('współczynnik filtru')
plt.title('Odpowiedź impulsowa filtru DP, L=30')
plt.show()
Wyjaśnienia wymaga parametr pass_zero. Należy go ustawić na True jeżeli zerowa częstotliwość ma znajdować się w paśmie przepustowym, False jeżeli w paśmie zaporowym. W przypadku filtru dolnoprzepustowego podajemy True (jest to wartość domyślna).
Na wykresie odpowiedzi impulsowej widzimy że, jest ona symetryczna. Ponieważ wybraliśmy parzystą liczbę współczynników (30), środek symetrii leży pomiędzy środkowymi próbkami. Gdybyśmy wzięli nieparzystą długość filtru, środkowa próbka wyznaczałaby środek symetrii.
Aby wykreślić charakterystykę częstotliwościową filtru, możemy obliczyć ją za pomocą funkcji freqz. Domyślnie oblicza ona 512 punktów charakterystyki, liczbę punktów można zmienić dodając parametr worN. Funkcja zwraca dwie wartości: tablicę pulsacji punktów (od 0 do $\pi$) oraz tablicę zespolonych amplitud charakterystyki częstotliwościowej w poszczególnych punktach.
w, h = sig.freqz(dp)
f = w * fs / (2 * np.pi)
idealna = np.zeros(len(f))
idealna[f > 1000] = -100
plt.plot(f, idealna, '#808080')
plt.plot(f, 20 * np.log10(np.abs(h)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakerystyka częstotliwościowa filtru DP, L=30')
plt.show()
Wpływ długości filtru na jego charakterystykę¶
Na wykresie charakterystyki filtru o długości 30 widać że charakterystyka ta bardzo wyraźnie odbiega od idealnej: przejście do pasma zaporowego jest bardzo szerokie, a w paśmie zaporowym są widoczne zafalowania charakterystyki. Jak zmieni się kształt charakterystyki dla większej długości filtru? Większe okno obejmie dłuższy fragment idealnej odpowiedzi impulsowej filtru, a więc można się spodziewać że charakterystyka filtru będzie "dokładniejsza". Zilustrujmy to na wykresie charakterystyk filtrów o różnej długości i stałej częstotliwości granicznej 1 kHz.
dp_30 = sig.firwin(30, 1000, fs=fs)
_, char_30 = sig.freqz(dp_30)
dp_50 = sig.firwin(50, 1000, fs=fs)
_, char_50 = sig.freqz(dp_50)
dp_75 = sig.firwin(75, 1000, fs=fs)
_, char_75 = sig.freqz(dp_75)
dp_100 = sig.firwin(100, 1000, fs=fs)
_, char_100 = sig.freqz(dp_100)
plt.plot(f, idealna, '#808080', label='idealna')
plt.plot(f, 20 * np.log10(np.abs(char_30)), label='L=30')
plt.plot(f, 20 * np.log10(np.abs(char_50)), label='L=50')
plt.plot(f, 20 * np.log10(np.abs(char_75)), label='L=75')
plt.plot(f, 20 * np.log10(np.abs(char_100)), label='L=100')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.ylim(-100, 0)
plt.legend()
plt.title('Porównanie charakterystyk filtrów DP o różnej długości')
plt.show()
Ten sam wykres, ale z powiększeniem skali częstotliwości:
plt.plot(f, idealna, '#808080', label='idealna')
plt.plot(f, 20 * np.log10(np.abs(char_30)), label='L=30')
plt.plot(f, 20 * np.log10(np.abs(char_50)), label='L=50')
plt.plot(f, 20 * np.log10(np.abs(char_75)), label='L=75')
plt.plot(f, 20 * np.log10(np.abs(char_100)), label='L=100')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.legend()
plt.xlim(0, 5000)
plt.ylim(-100, 0)
plt.title('Porównanie charakterystyk filtrów DP o różnej długości')
plt.show()
Z wykresów wynika, że w "nieidealnym" filtrze FIR zawsze występuje pasmo przejściowe (transition band) w pobliżu częstotliwości granicznej. Zgodnie z przewidywaniem, większa długość filtru powoduje zmniejszenie szerokości pasma przejściowego. W paśmie zaporowym występują zafalowania, ich kształt również zależy od długości filtru.
Zwiększenie liczby współczynników filtru poprawia "dokładność" filtracji. Trzeba jednak zdawać sobie sprawę z tego, że większa liczba współczynników oznacza również:
- większą liczbę operacji do wykonania (każdy współczynnik to dodatkowe mnożenie i dodawanie),
- większą pamięć potrzebną do przechowania wartości współczynników i próbek sygnału,
- większe opóźnienie wprowadzane przez filtr.
Powyższe aspekty mają znaczenie przy implementacji filtrów FIR w systemach o ograniczonych zasobach, np. na procesorach sygnałowych. Opóźnienie wprowadzane przez filtr ma też znaczenie w systemach, które mają działać w czasie rzeczywistym.
Praktycznie stosowane filtry FIR mają często długość ponad 100. Należy zawsze znaleźć kompromis pomiędzy dokładnością filtracji oraz wykorzystaniem procesora i pamięci oraz czasem wykonywania obliczeń.
Filtracja sygnału za pomocą filtru FIR¶
Mając zaprojektowany filtr FIR, możemy za jego pomocą przetworzyć sygnał. Szum biały to sygnał losowy o płaskim widmie częstotliwościowym, więc powinniśmy wyraźnie zauważyć efekt filtracji. Aby wygenerować szum biały, wykorzystamy generator liczb losowych o rozkładzie równomiernym.
Do przeprowadzenia filtracji użyjemy funkcji lfilter z modułu scipy.signal. Pierwszym parameterem są obliczone współczynniki filtru, drugi parametr ma dla filtru FIR zawsze wartość 1, trzecim parametrem jest sygnał poddawany filtracji.
szum = np.random.random(2048) * 2 - 1
przefiltrowany = sig.lfilter(dp_100, 1, szum)
widmo_szum = 20 * np.log10(np.abs(np.fft.rfft(szum * np.hamming(2048))) / 1024)
widmo_przefiltrowany = 20 * np.log10(np.abs(np.fft.rfft(przefiltrowany * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f, widmo_szum, '#c0c0c0', label='oryginalny')
plt.plot(f, widmo_przefiltrowany, label='po filtracji')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.legend(loc='upper right')
plt.title('Wynik filtracji szumu białego przez filtr DP 1 kHz')
plt.show()
Widzimy efekt filtracji dolnoprzepustowej.
W drugim przykładzie "przepuścimy" przez filtr sygnał złożony z trzech sinusów o częstotliwościach 1, 3 i 5 kHz. Chcemy odfiltrować trzeci prążek. Zastosujemy filtr o częstotliwości granicznej 4 kHz i długości 100.
n = np.arange(2048)
trzysin = (np.sin(2 * np.pi * n * 1000 / fs) +
np.sin(2 * np.pi * n * 3000 / fs) +
np.sin(2 * np.pi * n * 5000 / fs)) / 3
widmo_trzysin = 20 * np.log10(np.abs(np.fft.rfft(trzysin * np.hamming(2048))) / 1024)
dp_4 = sig.firwin(100, 4000, fs=fs)
trzysin_filtr = sig.lfilter(dp_4, 1, trzysin)
widmo_filtr = 20 * np.log10(np.abs(np.fft.rfft(trzysin_filtr * np.hamming(2048))) / 1024)
_, ax = plt.subplots(2)
ax[0].plot(f, widmo_trzysin, label='oryginalny')
ax[1].plot(f, widmo_filtr, label='po filtracji')
for a in ax:
a.set_xlabel('częstotliwość [Hz]')
a.set_ylabel('amplituda [dB]')
a.legend()
a.set_xlim(0, 10000)
ax[0].set_title('Filtracja DP sumy trzech sinusów')
plt.show()
Prążek na częstotliwości 5 kHz został przez filt wytłumiony, amplituda pozostałych prążków nie uległa zmianie.
Filtr górnoprzepustowy¶
Filtr górnoprzepustowy (GP), high pass (HP), działa "odwrotnie" do filtru dolnoprzepustowego: tłumi czestotliwości poniżej częstotliwości granicznej, przepuszcza pozostałą część pasma. Zasada projektowania jest taka jak dla filtru dolnoprzepustowego, z dwoma różnicami. Po pierwsze, parametr pass_zero w funkcji firwin musi być ustawiony na False, ponieważ zerowa częstotliwość leży w paśmie zaporowym. Po drugie, jeżeli spróbujemy obliczyć współczynniki filtru GP o parzystej długości za pomocą funkcji firwin, otrzymamy komunikat o błędzie:
ValueError: A filter with an even number of coefficients must have zero response at the Nyquist frequency.O co chodzi? Projektując filtr FIR metodą okienkowania otrzymujemy symetryczną odpowiedź impulsową. Jeżeli wybierzemy nieparzystą liczbę współczynników, środkowa wartość odpowiedzi impulsowej będzie leżała w środku symetrii. Tak zaprojektowany filtr FIR nazywamy filtrem typu I. Typ I nadaje się do zaprojektowania filtru przepustowo-zaporowego o dowolnej charakterystyce. Natomiast jeżeli liczba współczynników będzie parzysta, środek symetrii wypadnie pomiędzy dwoma środkowymi wartościami odpowiedzi impulsowej. Jest to filtr typu II. Parzyste filtry są przydatne w implementacji np. na procesorach sygnałowych, ponieważ możemy wykorzystać symetrię odpowiedzi impulsowej do zredukowania liczby mnożeń. Ale typ II ma istotne ograniczenie: częstotliwość Nyquista ($f_s/2$) nie może leżeć w paśmie przepustowym. Zatem nie zaprojektujemy "parzystego" górnoprzepustowego filtru FIR, musimy wziąć nieparzystą liczbę współczynników.
Na marginesie, można zaprojektować filtr FIR tak, że odpowiedź impulsowa jest antysymetryczna. W ten sposób powstają filtry typu III (nieparzyste) i IV (parzyste). Te filtry stosuje się do specyficznych zastosowań (np. filtry Hilberta).
Charakterystyka zaprojektowanego filtru górnoprzepustowego o 101 współczynnikach, częstotliwość graniczna 4 kHz:
gp = sig.firwin(101, 4000, pass_zero=False, fs=fs)
w, h_gp = sig.freqz(gp)
f = w * fs / (2 * np.pi)
idealna = np.zeros(len(f))
idealna[f < 4000] = -100
plt.plot(f, idealna, '#808080')
plt.plot(f, 20 * np.log10(np.abs(h_gp)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakterystyka filtru GP 4 kHz, L=101')
plt.show()
Wynik filtrowania szumu białego oraz sygnału "trzy sinusy" - dwa pierwsze prążki zostały wytłumione przez filtr (drugi prążek przetrwał w "szczątkowej" formie, ponieważ był objęty krańcem pasma przejściowego).
szum = np.random.random(2048) * 2 - 1
przefiltrowany = sig.lfilter(gp, 1, szum)
widmo_szum = 20 * np.log10(np.abs(np.fft.rfft(szum * np.hamming(2048))) / 1024)
widmo_przefiltrowany = 20 * np.log10(np.abs(np.fft.rfft(przefiltrowany * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f, widmo_szum, '#c0c0c0', label='oryginalny')
plt.plot(f, widmo_przefiltrowany, label='po filtracji')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.legend()
plt.title('Wynik filtracji szumu białego przez filtr GP 4 kHz')
plt.show()
trzysin_filtr_gp = sig.lfilter(gp, 1, trzysin)
widmo_filtr_gp = 20 * np.log10(np.abs(np.fft.rfft(trzysin_filtr_gp * np.hamming(2048))) / 1024)
_, ax = plt.subplots(2)
ax[0].plot(f, widmo_trzysin, label='oryginalny')
ax[1].plot(f, widmo_filtr_gp, label='po filtracji')
for a in ax:
a.set_xlabel('częstotliwość [Hz]')
a.set_ylabel('amplituda [dB]')
a.legend()
a.set_xlim(0, 10000)
ax[0].set_title('Filtracja GP sumy trzech sinusów')
plt.show()
Filtr środkowoprzepustowy i środkowozaporowy¶
Filtr środkowoprzepustowy (SP), band pass (BP), nazywany też pasmowo-przepustowym, przepuszcza częstotliwości w zakresie leżącym pośrodku pełnego pasma. Ma on dwa pasma zaporowe i jedno przepustowe. Jest on więc opisany dwoma częstotliwościami granicznymi - górną i dolną. Filtr ma zerowe wzmocnienie zarówno dla zerowej częstotliwości (pass_zero=False), jak i dla częstotliwości Nyquista, zatem możemy stosować typ I albo typ II. Poniżej pokazano charakterystykę filtru o 100 współczynnikach i paśmie przepustowym od 1 kHz do 5 kHz. Przy projektowaniu filtru za pomocą funkcji firwin należy podać obie częstotliwości graniczne w nawiasach () albo [].
pp = sig.firwin(100, (1000, 5000), pass_zero=False, fs=fs)
w, h_pp = sig.freqz(pp)
f = w * fs / (2 * np.pi)
idealna = np.zeros(len(f))
idealna[(f < 1000) | (f > 5000)] = -100
plt.plot(f, idealna, '#808080')
plt.plot(f, 20 * np.log10(np.abs(h_pp)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakterystyka czestotliwościowa filtru SP 1-5 kHz, L=100')
plt.show()
Poniżej wynik filtracji szumu białego za pomocą tego samego filtru oraz wynik filtracji sygnału "trzy sinusy" za pomocą filtru o paśmie od 2 do 4 KHz, którego zadaniem jest przepuszczenie tylko środkowego prążka.
przefiltrowany = sig.lfilter(pp, 1, szum)
widmo_przefiltrowany = 20 * np.log10(np.abs(np.fft.rfft(przefiltrowany * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f, widmo_szum, '#c0c0c0', label='oryginalny')
plt.plot(f, widmo_przefiltrowany, label='po filtracji')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.legend()
plt.title('Wynik filtracji szumu białego przez filtr SP 1-5 kHz')
plt.show()
pp_24 = sig.firwin(100, (2000, 4000), pass_zero=False, fs=fs)
trzysin_filtr_pp = sig.lfilter(pp_24, 1, trzysin)
widmo_filtr_pp = 20 * np.log10(np.abs(np.fft.rfft(trzysin_filtr_pp * np.hamming(2048))) / 1024)
_, ax = plt.subplots(2)
ax[0].plot(f, widmo_trzysin, label='oryginalny')
ax[1].plot(f, widmo_filtr_pp, label='po filtracji')
for a in ax:
a.set_xlabel('częstotliwość [Hz]')
a.set_ylabel('amplituda [dB]')
a.legend()
a.set_xlim(0, 10000)
ax[0].set_title('Filtracja SP sumy trzech sinusów')
plt.show()
Działania filtru środkowozaporowego (SZ), band stop (BS), nazywanego też pasmowo-zaporowym, łatwo się domyślić: jest ono przeciwne do filtru środkowoprzepustowego. Posiada pasmo zaporowe "w środku" i dwa pasma przepustowe. Częstotliwość zerowa leży w paśmie przepustowym (pass_zero=True), czestotliwość Nyquista - również, a więc nie możemy wykorzystać typu II o parzystej liczbie współczynników. Charakterystyka częstotliwościowa filtru SZ o paśmie zaporowym 1 kHz do 5 kHz i długości filtru 101 wygląda następująco.
pz = sig.firwin(101, (1000, 5000), fs=fs)
w, h_pz = sig.freqz(pz)
f = w * fs / (2 * np.pi)
idealna = np.zeros(len(f))
idealna[(f > 1000) & (f < 5000)] = -100
plt.plot(f, idealna, '#808080')
plt.plot(f, 20 * np.log10(np.abs(h_pz)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakterystyka częstotliwościowa filtru SZ 1-5 kHz, L=101')
plt.show()
Wynik filtracji szumu białego za pomocą tego filtru:
przefiltrowany = sig.lfilter(pz, 1, szum)
widmo_przefiltrowany = 20 * np.log10(np.abs(np.fft.rfft(przefiltrowany * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / fs)
plt.plot(f, widmo_szum, '#c0c0c0', label='oryginalny')
plt.plot(f, widmo_przefiltrowany, label='po filtracji')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.legend()
plt.title('Wynik filtracji szumu białego przez filtr SZ 1-5 kHz')
plt.show()
Filtry środkowozaporowe o bardzo wąskim paśmie zaporowym nazywa się filtrami notch. Stosuje się je wtedy, gdy istnieje potrzeba usunięcia z widma sygnału konkretnej składowej. Muszą one mieć dużą liczbę współczynników aby zapewnić odpowiednie tłumienie w paśmie zaporowym, ponieważ częstotliwości graniczne są położone blisko siebie. Poniżej pokazano charakterystykę filtru o długości 171 i paśmie zaporowym od 2,5 do 3,5 kHz oraz wynik filtracji sygnału "trzy sinusy" w celu usunięcia drugiego prążka (o częstotliwości 3 kHz). Filtr o długości 101 nie był wystarczający do wytłumienia prążka. Wartość 171 dobrano eksperymentalnie, dla większych długości zaczynają się pojawiać nadmierne zafalowania w paśmie zaporowym.
pz_notch = sig.firwin(171, (2500, 3500), fs=fs)
w, h_pz_notch = sig.freqz(pz_notch)
f = w * fs / (2 * np.pi)
plt.plot(f, 20 * np.log10(np.abs(h_pz_notch)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakterystyka filtru SZ 2,5-3,5 kHz, L=171')
plt.xlim(0, 5000)
plt.show()
trzysin_filtr_notch = sig.lfilter(pz_notch, 1, trzysin)
widmo_filtr_notch = 20 * np.log10(np.abs(np.fft.rfft(trzysin_filtr_notch * np.hamming(2048))) / 1024)
f = np.fft.rfftfreq(2048, 1 / fs)
_, ax = plt.subplots(2)
ax[0].plot(f, widmo_trzysin, label='oryginalny')
ax[1].plot(f, widmo_filtr_notch, label='po filtracji')
for a in ax:
a.set_xlabel('częstotliwość [Hz]')
a.set_ylabel('amplituda [dB]')
a.legend()
a.set_xlim(0, 10000)
ax[0].set_title('Filtracja SZ notch sumy trzech sinusów')
plt.show()
Liniowofazowość filtrów FIR, opóźnienie grupowe¶
Do tej pory rozpatrywaliśmy charakterystyki amplitudowe filrów FIR, teraz weźmy pod uwagę charakterystykę fazową. Dla filtru dolnoprzepustowego o 31 współczynnikach i częstotliwości granicznej 1 kHz wygląda ona następująco.
dp = sig.firwin(31, 1000, fs=fs)
w, h = sig.freqz(dp)
f = w * fs / (2 * np.pi)
plt.plot(f, np.degrees(np.angle(h)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('faza [st]')
plt.title('Charakterystyka fazowa filtru DP 1 kHz, L=31')
plt.show()
Charakterystyka fazowa jest odcinkami liniowa, z przeskokami fazy. Sama faza filtru niewiele nam tu powie. Bardziej istotna jest charakterystyka opóźnienia grupowego (group delay, oznaczane symbolem $\tau$), które jest pochodną charakterystyki fazowej względem pulsacji, z ujemnym znakiem. Do obliczenia opóźnienia grupowego użyjemy funkcji scipy.signal.group_delay. Zwraca ona dwie wartości: tablicę pulsacji punktów charakterystyki (od 0 do $\pi$) oraz wartości opóźnienia grupowego dla tych punktów. Podwójne nawiasy w wywołaniu funkcji są konieczne.
w, tau = sig.group_delay((dp, 1))
plt.plot(w * fs / (2 * np.pi), tau)
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('opóźnienie grupowe [próbki]')
plt.title('Opóźnienie grupowe filtru DP 1 kHz, L=31')
plt.ylim(0, 30)
plt.show()
Opóźnienie grupowe dla tego filtru jest stałe i wynosi 15 próbek. Oznacza to, że każda składowa częstotliwościowa sygnału jest opóźniana o tę samą liczbę próbek. Na przykład, jeżeli w podawanym wcześniej przykładzie mamy sygnał złożony z trzech sinusów, to "czas przejścia" tych składowych przez filtr jest taki sam.
Filtry, które mają (odcinkami) liniową charakterystykę fazową, a w konsekwencji stałą charakterystykę opóźnienia grupowego, nazywamy filtrami liniowofazowymi (linear phase filters). Zaletą tych filtrów jest to, że nie zaburzają zależności fazowych pomiędzy składowymi częstotliwościami sygnału. Ma to znaczenie w sytuacjach, w których zależności fazowe są istotną cechą sygnału, przykładem jest stereofoniczny sygnał dźwiękowy.
Filtry FIR typu I/II/III/IV, czyli o symetrycznej/antysymetrycznej odpowiedzi impulsowej i o nieparzystej/parzystej liczbie współczynników, są zawsze liniowofazowe. Z tego powodu te typy filtrów są stosowane najczęściej, jednak projektuje się również filtry FIR nie mające liniowej charakterystyki fazowej (np. filtry minimalnofazowe), jeżeli wymagania dotyczące filtru nie obejmują liniowości fazy, a potrzebne jest inne działanie filtru.
W powyższym przykładzie, opóźnienie grupowe filtru wynosi 15 próbek, czyli połowę wartości rzędu filtru. Taka zależność zachodzi dla każdego liniowofazowego filtru FIR. Dla filtrów o nieparzystej liczbie współczynników (typy I i III), opóźnienie grupowe jest zawsze liczbą całkowitą, natomiast dla typów II i IV (o parzystej liczbie współczynników), opóźnienie grupowe jest niecałkowite ("połówkowe").
W praktyce zdarza się, że potrzebujemy do przetwarzania zarówno sygnału przefiltrowanego, jak i oryginalnego. Aby zilustrować wpływ opóźnienia grupowego na sygnał, pokażmy taki przypadek dla sygnału świergotowego (chirp) i filtru dolnoprzepustowego pokazanego powyżej.
chirp = sig.chirp(np.linspace(0, 1, 2001), 20, 1, 1000)
chirp_filtr = sig.lfilter(dp, 1, chirp)
plt.plot(chirp[:200], label='oryginalny')
plt.plot(chirp_filtr[:200], label='po filtracji')
plt.xlabel('nr próbki')
plt.ylabel('amplituda [dB]')
plt.legend()
plt.title('Sygnał oryginalny i przefiltrowany')
plt.show()
Jak widać, sygnał przefiltrowany "spóźnia się" względem oryginalnego. Aby zsynchronizować oba sygnały, musimy skompensować opóźnienie wprowadzone przez filtr, opóźniając sygnał oryginalny. O ile musimy go opóźnić? Oczywiście o wartość opóźnienia grupowego filtru, czyli o 15 próbek.
chirp_op = np.insert(chirp, 0, np.zeros(15))
plt.plot(chirp_op[:200], label='oryginalny')
plt.plot(chirp_filtr[:200], label='po filtracji')
plt.xlabel('nr próbki')
plt.ylabel('amplituda [dB]')
plt.legend(loc='upper right')
plt.title('Wyrównane czasowo sygnały: przed i po filtracji')
plt.show()
Teraz oba sygnały są wyrównane w czasie. Zauważmy, że jeżeli stosujemy filtr typu II, opóźnienie grupowe nie będzie liczbą całkowitą. W takiej sytuacji albo musimy pogodzić się z opóźnieniem o 0,5 próbki, albo musimy zastosować ułamkowe opóźnienie (co wymaga interpolacji sygnału).
Filtracja blokowa, pamięć filtru¶
W poprzednich przykładach przetwarzaliśmy sygnał w całości. W praktycznych zastosowaniach filtrowany sygnał jest zazwyczaj ciągły (stale napływa na wejście układu) i musimy go przetwarzać w blokach próbek. Nie możemy po prostu przefiltrować każdego bloku osobno, ponieważ w takim przypadku dla każdego bloku przyjmowalibyśmy zerowe warunki początkowe (filtr zaczynałby pracę zawsze od początku). Jeżeli np. rozpoczynamy filtrację drugiego bloku, potrzebujemy próbek z pierwszego bloku. Dlatego dane potrzebne do filtracji, w przypadku filtrów FIR: $N$ poprzednich próbek, muszą być przechowywane w pamięci filtru, utrzymywanej pomiędzy blokami filtrowanych próbek. Dane z pamieci filtru ustalają warunki początkowe (initial conditions) dla każdego wywołania filtru.
Funkcja lfilter umożliwia podanie zapamiętanego stanu filtru jako parametr zi. Zmienia się wtedy sposób działania funkcji: zwraca ona dwa wyniki. Pierwszym są wyniki filtracji, drugim — końcowy stan pamięci filtru, który powinniśmy podać przy kolejnym wywołaniu funkcji. Dla pierwszego bloku musimy obliczyć początkowy stan filtru, przyjmując zerowe wartości na wejściu i na wyjściu filtru. Służy do tego funkcja scipy.signal.lfiltic. Podajemy w niej współczynniki b filtru, wartość 1 dla filtru FIR oraz warunki początkowe — dla zerowych warunków wystarczy podać [0].
W poniższym przykładzie najpierw obliczamy zerowy stan początkowy. Następnie symulujemy podział sygnału na bloki po 128 próbek, każdy blok przetwarzamy za pomocą filtru, podając poprzedni stan i otrzymując przefiltrowane próbki bloku oraz nowy stan filtru.
pamiec = sig.lfiltic(dp, 1, [0])
trzysin_bloki = trzysin.reshape(-1, 128)
for blok in trzysin_bloki:
wynik, pamiec = sig.lfilter(dp, 1, blok, zi=pamiec)
Okno Kaisera¶
Do tej pory projektowaliśmy filtry FIR podając częstotliwości graniczne i liczbę współczynników. Często jednak chcemy, aby filtr spełniał określone wymagania, np. aby tłumienie w paśmie zaporowym wynosiło co najmniej 70 dB. Możemy próbować różnych długości filtru aż otrzymamy pożądane tłumienie, ale są lepsze metody uzyskania tego efektu. Tutaj opiszemy wykorzystanie okna Kaisera.
W przeciwieństwie do typowych okien czasowych, okno Kaisera jest parametryczne. Posiada parametr $\beta$, który steruje kształtem okna - większe wartości $\beta$ powodują, że okno staje się coraz węższe, przez co po zastosowaniu okna, zwiększa się tłumienie listków bocznych i poszerza się listek główny. Możemy obliczyć wartość $\beta$ dającą pożądane tłumienie. Możliwe jest także otrzymanie charakterystyki o ustalonej szerokości pasma przejściowego. Jak pamietamy, szerokość pasma przejściowego zależy od rzędu filtru. Funkcja scipy.signal.kaiserord oblicza liczbę współczynników filtru (czyli rozmiar okna Kaisera) oraz parametr $\beta$ dla podanej wartości tłumienia oraz szerokość pasma przejściowego - tę drugą wartość podajemy nie w Hz, ale jako ułamek częstotliwości Nyquista. Np. dla tłumienia 70 dB i szerokości pasma przejściowego 500 Hz mamy:
L, beta = sig.kaiserord(70, 500 / (fs/2))
print('Długość filtru = {}, parametr beta = {}'.format(L, beta))
Obliczamy współczynniki filtru tak jak wcześniej, ale oprócz typu okna musimy podać też wartość $\beta$.
dp_kaiser = sig.firwin(L, 1000, window=('kaiser', beta), fs=fs)
w, h = sig.freqz(dp_kaiser)
plt.plot(w * fs / (2 *np.pi), 20 * np.log10(np.abs(h)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Filtr DP zaprojektowany za pomocą okna Kaisera')
plt.xlim(0, 5000)
plt.show()
Alternatywna metoda projektowania FIR metodą okienkowania¶
W module scipy.signal jest również dostępna funkcja firwin2, która umożliwia zaprojektowanie filtru FIR metodą okienkowania w nieco inny sposób. Nie podaje się częstotliwości granicznych, zamiast tego podaje się punkty na osi częstotliwości oraz wartości wzmocnienia w tych punktach. Na przykład, chcemy uzyskać filtr dolnoprzepustowy o paśmie przepustowym do 1000 Hz, paśmie przejściowym 1000-1500 Hz i paśmie zaporowym powyżej 1500 Hz. Wzmocnienie filtru ma wynosić 1 w paśmie przepustowym i 0 w paśmie zaporowym. Zakładając że liczba współczynników ma wynosić 101, filtr projektujemy następująco. Lista częstotliwości (drugi parametr) musi zaczynać się od 0 i kończyć częstotliwością Nyqusita.
dp2 = sig.firwin2(101, [0, 1000, 1500, 24000], [1, 1, 0, 0], fs=fs)
w, h = sig.freqz(dp2)
plt.plot(w * fs / (2 *np.pi), 20 * np.log10(np.abs(h)))
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Charakterystyka filtru DP zaprojektowana za pomocą firwin2')
plt.xlim(0, 5000)
plt.ylim(-90, 10)
plt.show()
Zaletą tej funkcji jest nieco bardziej intuicyjne podawanie specyfikacji (nie trzeba używać parametru pass_zero), możliwość ustalania szerokości pasma przejściowego, możliwość elastycznego ustalania punktów charakterystyki, a także możliwość projektowania filtrów typu III i IV przez podanie parametru antisymmetric=True. Wybór funkcji firwin lub firwin2 zależy od użytkownika.
Metody optymalnego projektowania filtrów FIR - najmniejszych kwadratów i Parksa-McClellana¶
Metoda okienkowania jest prostą metodą projektowania filtrów FIR, natomiast nie jest metodą optymalną, tzn. nie daje najlepszego możliwego dopasowania charakterystyki zaprojektowanego filtru do idealnej charakterystyki. Pierwszą z metod optymalnego projektowania filtrów FIR jest metoda najmniejszych kwadratów. Celem metody jest uzyskanie jak najmniejszej ważonej różnicy pomiędzy zaprojektowaną a idealną charakterystyką w paśmie przepustowym oraz zaporowym. Jest ona realizowana za pomocą funkcji scipy.signal.firls. Specyfikacja charakterystyki filtru jest taka jak dla funkcji firwin2.
dp_ls = sig.firls(101, [0, 1000, 1500, 24000], [1, 1, 0, 0], fs=fs)
Komputery bardzo dobrze radzą sobie z rozwiązywaniem problemów optymalizacyjnych. Dajemy komputerowi pewien problem do rozwiązania, np. znaleźć współczynniki filtru o podanym rzędzie, tak aby różnica w stosunku do idealnej charakterystyki była jak najmniejsza. Komputer wykonuje obliczenia w sposób iteracyjny, tak aby z każdą kolejną iteracją przybliżać się do optymalnego rozwiązania. Metoda Parksa-McClellana jest najczęściej stosowaną iteracyjną metodą projektowania filtrów FIR. Podobnie jak metoda najmniejszych kwadratów, minimalizuje ona różnicę między zaprojektowaną a idealną charakterystyką filtru w paśmie przepustowym i zaporowym. Metoda ta opiera się na algorytmie Remeza, dlatego funkcjonuje ona również pod tą nazwą (ponieważ jest prostsza do napisania).
W module scipy.signal mamy funkcję właśnie o nazwie remez. Specyfikację filtru podaje się jako listę punktów na osi częstotliwości oraz - inaczej niż w poprzednich metodach - listę pożądanych wartości wzmocnienia w kolejnych pasmach przepustowych i zaporowych. W przypadku filtru dolnoprzepustowego podajemy wzmocnienie 1 w pierwszym paśmie i 0 w drugim.
dp_rem = sig.remez(101, [0, 1000, 1500, 24000], [1, 0], fs=fs)
Wykreślmy teraz charakterystyki trzech filtrów o tych samych parametrach, zaprojektowanych trzema metodami: okienkową, najmniejszych kwadratów i Remeza.
w, h_ls = sig.freqz(dp_ls)
w, h_rem = sig.freqz(dp_rem)
f = w * fs / (2 *np.pi)
plt.plot(f, 20 * np.log10(np.abs(h)), label='firwin2')
plt.plot(f, 20 * np.log10(np.abs(h_ls)), label='firls')
plt.plot(f, 20 * np.log10(np.abs(h_rem)), label='remez')
plt.xlabel('częstotliwość [Hz]')
plt.ylabel('amplituda [dB]')
plt.title('Porównanie charakterystyk filtrów DP zaprojektowanych różnymi metodami')
plt.xlim(0, 5000)
plt.ylim(-90, 10)
plt.legend()
plt.show()
Z powyższego porównania widać, że co prawda filtr zaprojektowany metodą okienkową ma największe tłumienie w paśmie zaporowym (60 dB), ale szerokość pasma przejściowego jest znacznie większa niż zakładano - pierwsze minimum wypada na częstotliwości ok. 2100 Hz zamiast 1500 Hz. Filtry zaprojektowane metodami optymalnymi mają słabsze tłumienie w paśmie zaporowym - ok. 25 dB, przy czym dla metody najmniejszych kwadratów tłumienie wzrasta z częstotliwością, a dla metody Remeza jest stałe. Natomiast szerokość pasma przejściowego jest znacznie mniejsza niż dla metody okienkowej i mniej więcej zgodna ze specyfikacją - pierwsze minimum w paśmie zaporowym znajduje się na ok. 1500 Hz.
Należy pamiętać, że obie metody optymalnego projektowania filtrów dopasowują charakterystykę do idealnej tylko w pasmach zaporowych i przepustowych. Kształt charakterystyki filtru w pasmach przejściowych nie jest zdefiniowany. W niektórych przypadkach może zdarzyć się że charakterystyka w paśmie przejściowym nie jest gładka, a pojawiają się w niej "przerzuty". Ponadto obie metody umożliwiają nadanie względnych wag poszczególnym pasmom, co przydaje się jeżeli zależy nam na dokładniejszym dopasowaniu charakterystyki w określonym paśmie częstotliwości.
Podsumowanie¶
- Filtry FIR to cyfrowe realizacje układów modyfikujących sygnały.
- Najczęściej używane filtry: dolnoprzepustowy, górnoprzepustowy, środkowoprzepustowy, środkowozaporowy.
- Długość filtru to liczba współczynników filtru, rząd filtru jest o jeden mniejszy niż jego długość.
- Metoda okienkowania jest podstawową metodą projektowania filtrów FIR. Specyfikacja filtru objemuje częstotliwości graniczne oraz liczbę współczynników.
- Filtry FIR projektowane metodą okienkową są zawsze liniowofazowe - mają stałe opóźnienie równe połowie rzędu filtru, dla wszystkich częstotliwości.
- Filtry typu I wymagają nieparzystej liczby współczynników, nadają się do projektowania każdej charakterystyki filtru.
- Filtry typu II wymagają parzystej liczby współczynników, nie nadają się do charakterystyk mających niezerowe wzmocnienie na częstotliwości Nyquista (górnoprzepustowych i pasmowo-zaporowych).
- Użycie okna Kaisera umożliwia uzyskanie pożądanej szerokości pasma przejściowego i tłumienia w paśmie zaporowym.
- Metody optymalnego projektowania filtrów FIR - metoda najmniejszych kwadratów oraz metoda Parksa-McClellana (Remeza) - umożliwiają lepsze dopasowanie charakterystyki filtru w paśmie przepustowym i zaporowym do idealnej charakterystyki. W porównaniu do metody okienkowania: zmniejszają szerokość pasma przejściowego, ale mogą zmniejszyć tłumienie w paśmie zaporowym.