Ćwiczenia 7
Spis treści
Notebook
Implementacja filtrowania: funkcja lfilter
Filtrowanie zgodne z powyższymi równaniami zaimplementowane jest w pythonie w module scipy.signal jako funkcja lfilter.
Podstawowe wywołanie tej funkcji dla sygnału we wygląda następująco: wy = scipy.signal.lfilter(b,a,we) gdzie b, a są to współczynniki z poprzedniego równania. np:
import numpy as np
import pylab as py
from scipy.signal import lfilter
b = np.array([0.7, 0.5]) # licznik
a = np.array([1.0, -0.9]) # mianownik
we = np.random.randn(100)
wy=lfilter(b,a,we);
py.subplot(2,1,1)
py.plot(we)
py.subplot(2,1,2)
py.plot(wy)
py.show()
Projektowanie filtru
Specyfikacja własności filtru
Filtr specyfikujemy opisując moduł jego pożądanej funkcji przenoszenia:
- ogólne określenie granic pasma przenoszenia np: dla sygnału próbkowanego 128 Hz zaprojektować filtr dolnoprzepustowy 30 Hz, oznacza, że w idealnej sytuacji filtr będzie przenosił bez zmian częstości od 0 do 30 Hz a od 30Hz do 64Hz będzie całkowicie tłumił.
- w bardziej rygorystycznym opisie możemy dodatkowo podać:
- dopuszczalną amplitudę oscylacji Rp w paśmie przenoszenia (pass band),
- minimalne tłumienie Rs pasma tłumienia (stop band),
- szerokość pasma przejściowego
Oznaczenia jak na poniższym rysunku.
Funkcje do projektownaia filtrów FIR
W module scipy.signal mamy kilka funkcji do projektowania filtrów typu FIR.
firwin
Najprostszą koncepcyjnie metodą projektowania filtrów FIR jest metoda okienkowa. Metoda składa się z następujących kroków: w dziedzinie częstości projektowana jest idealna funkcja przenoszenia, obliczana jest od niej odwrotna transformata Fouriera, następnie otrzymana sekwencja czasowa (odpowiedź impulsowa) jest przemnażana przez wybrane okno.
Metoda ta zaimplementowana jest w funkcji scipy.signal.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=1.0). Pozwala ona projektować filtry dolno- i górno- przepustowe oraz pasmowo przepustowe i pasmowo zaporowe metodą okienkową.
Najważniejsze parametry (kompletny opis w dokumentacji firwin):
- numtaps: int, ilość współczynników filtru (rząd filtru+1). Liczba ta musi być parzysta jeśli pasmo przenoszenia ma zawierać częstość Nyquista.
- cutoff: częstość odcięcia filtru. Może być jedną liczbą zmiennoprzecinkową dla filtru dolno- lub górno- przepustowego lub tablicą dla filtrów pasmowych. Wyrażamy ją w tych samych jednostkach co nyq i musi zawierać się pomiędzy 0 a nyq.
- window: napis lub krotka: określa jakiego okna użyć do projektu filtru. Może to być dowolne okno spośród opisanych w scipy.signal.get_window
- pass_zero: bool, Jeśli True to zero jest przenoszone, jeśli False to nie jest. Ten parametr decyduje jak jest interpretowane pierwsze pasmo od 0 do cutoff - czy ma to być pasmo przenoszone czy tłumione.
- nyq: float. Częstość Nyquista.
- Zwraca: współczynniki b
Przykładowe projekty
We wszystkich poniższych przykładach zakładamy, że częstość próbkowania wynosi 256Hz:
- filtr dolnoprzepustowy rzędu 20 z częstością odcięcia 40Hz:
firwin(21, 40/128.0, nyq= 1) firwin(21, 40, nyq= 128)
- filtr górnoprzepustowy rzędu 15 z częstością odcięcia 5 Hz:
firwin(16, 5, nyq= 128, pass_zero=False)
- pasmowo przepustowy 5 rzędu przenoszący częstości pomiędzy 8 a 14 Hz:
firwin(6, [8 14], nyq= 128, pass_zero=False)
- pasmowo zaporowy 5 rzędu tłumiący częstości pomiędzy 8 a 14 Hz:
firwin(6, [8 14], nyq= 128, pass_zero=True)
Demo własności "audio" filtrów tego typu: applet
firwin2
Funkcja scipy.signal.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0) również implementuje okienkową metodę projektowania filtrów FIR. Daje ona nieco większą swobodę w kształtowaniu idealnej funkcji przenoszenia. Zadaje się ją przez podanie dwóch wektorów: freq i gain. Wektor freq definiuje punkty w częstości (jednostki takie same jak nyq, muszą zawierać 0 i nyq) dla których znana jest wartość pożądanego przenoszenia. Pożądane wartości przenoszenia odpowiadające kolejnym częstościom definiowane są w gain. Wartości freq muszą być ułożone w kolejności rosnącej, dopuszczalne jest powtórzenie tej samej wartości częstości i odpowiadających im różnych wartości gain aby zdefiniować nieciągłość funkcji przenoszenia. Znaczenie pozostałych parametrów jest takie samo jak dla firwin.
funkcja dostepna jest począwszy od wersji scipy 0.9.0
Przykład: filtr górno przepustowy. Liniowo narastające przenoszenie pomiędzy 20 a 40 Hz:
f = np.array([0, 20, 40, 128])
g = np.array([0, 0, 1, 1])
b = firwin2(50,f,g,nfreqs, nyq=128)# licznik
Badanie własności filtru w dziedzinie częstości
Podstawowe własności filtru możemy łatwo odczytać z wykresu przedstawiającego moduł funkcji przenoszenia i jej fazę w zależności od częstości. Moduł funkcji przenoszenia najczęściej wykreśla się w skali pół logarytmicznej.
W module scipy.signal zaimplementowana jest funkcja freqz, która wylicza wartości funkcję przenoszenia filtru zadanego współczynnikami b, a w zadanej ilości dyskretnych częstości. Jej użycie ilustruje poniższy przykład:
Przykład
Poniższy przykład przedstawia własności filtra dolnoprzepustowego rzędu 50 o częstości odcięcia równej połowie częstości Nyquista.
W dziedzinie częstości
import numpy as np
import pylab as py
from scipy.signal import freqz, firwin
b = firwin(50,0.5)# licznik
a = np.array([1.0]) # mianownik
n = 1000 # n ilość punktów na których będzie obliczona funkcja h
w, h = freqz(b,a,n)
m = np.abs(h) # moduł transmitancji
phi = np.unwrap(np.angle(h)) # faza
py.subplot(4,1,1)
py.plot(w,20*np.log10(m))
py.ylabel('[dB]')
py.title('moduł transmitancji')
py.subplot(4,1,2)
py.plot(w,phi)
py.title('faza transmitancji')
py.ylabel('rad')
py.xlabel('rad/próbki')
fazowe = - phi/w
py.subplot(4,1,3)
py.plot(w,fazowe)
py.ylim([20,30])
py.title('opóźnienie fazowe')
py.ylabel('próbki')
py.xlabel('rad')
grupowe = - np.diff(phi)/np.diff(w)
py.subplot(4,1,4)
py.plot(w[:-1],grupowe)
py.ylim([20,30])
py.title('opóźnienie grupowe')
py.ylabel('próbki')
py.xlabel('rad')
py.show()
W dziedzinie czasu
Porównajmy działanie tego filtra w czasie z funkcjami przesunięć fazowych i grupowych:
import numpy as np
import pylab as py
from scipy.signal import lfilter, firwin
b = firwin(50,0.5)# licznik
py.subplot(3,1,1)
py.plot(b,'o')
s1 = np.zeros(100)
t = np.arange(0,1,0.01)
s2 = np.sin(2*np.pi*10*t)
s = np.concatenate((s1,s2))
py.subplot(3,1,2)
py.plot(s)
py.subplot(3,1,3)
y = lfilter(b,1,s)
py.plot(y)
py.show()
Zadania
- Zaprojektuj i zbadaj własności filtru:
- FIR dolno z pasmem przenoszenia od 30 Hz dla sygnału próbkowanego 256 Hz
- Znajdź rząd filtru FIR:
- dolnoprzepustowego z pasmem przenoszenia do 40 Hz dla sygnału próbkowanego [math]256[/math] Hz, tak aby dla częstości powyżej 45 Hz jego tłumienie było nie mniejsze niż 20dB.
Filtry IIR
Filtry o nieskończonej odpowiedzi impulsowej (infinite impulse response, IIR) mają zazwyczaj dużo niższe rzędy niż filtry o skończonej odpowiedzi impulsowej (finite impulse response, FIR) z analogicznym poziomem tłumienia i szerokością pasma przejściowego.
W module scipy.signal mamy zaimplementowane kilka funkcji do projektowania „optymalnych” pod różnymi względami filtrów w klasycznych konfiguracjach: dolno- albo górnoprzepustowe i pasmowo-przepustowe albo pasmowo-zaporowe.
Funkcje do projektowania filtrów IIR dostępne w module scipy.signal
W module scipy.signal dostępne są funkcje do projektowania czterech typów filtrów: Butterwortha, Czebyszewa typu I i II, oraz eliptyczny. Do opisu wymagań projektowych funkcje te wykorzystują następujące pojęcia:
- wp, ws — krawędzie pasma przenoszenia i tłumienia. Częstości są znormalizowane do zakresu od 0 do 1 (1 odpowiada częstości Nyquista) przykładowo:
- dolno-przepustowy: wp = 0.2, ws = 0.3
- górno-przepustowy: wp = 0.3, ws = 0.2
- pasmowo-przepustowy: wp = [0.2, 0.5], ws = [0.1, 0.6]
- pasmowo-zaporowy: wp = [0.1, 0.6], ws = [0.2, 0.5]
- gpass — maksymalna dopuszczalna strata w pasmie przenoszenia (w funkcjach projektujących filtry jest to rp) (dB).
- gstop — minimalne wymagane tłumienie w pasmie tłumienia (w funkcjach projektujących filtry jest to rs) (dB).
- btype — typ filtra ('lowpass', 'highpass', 'bandpass', 'bandstop').
Funkcje do projektowania filtrów są zaimplementowane parami: jedna pomaga dobierać rząd filtru do wymagań projektowych, a druga oblicza współczynniki filtru.
Filtr Butterwortha daje gładką (bez tętnień) funkcję przenoszenia w całym zakresie częstości
[n,Wn]=buttord(wp, ws, gpass, gstop, analog=0)
[b,a]=butter(n,Wn)
n — rząd filtra, Wn — krawędzie pasm
Filtr Czebyszewa I rodzaju — gładka funkcja przenoszenia w paśmie tłumienia, minimalizowane są tętnienia w paśmie przenoszenia
[n,Wn]=cheb1ord(wp, ws, gpass, gstop, analog=0);
[b,a]=cheby1(n, rp, Wn, btype='low', analog=0, output='ba')
Filtr Czebyszewa II rodzaju — gładka funkcja przenoszenia w paśmie przenoszenia, minimalizowane tętnienia w paśmie tłumienia
[n,Wn]=cheb2ord(wp, ws, gpass, gstop, analog=0);
[b,a]=cheby2(n, rs, Wn, btype='low', analog=0, output='ba')
Filtr eliptyczny daje najostrzejsze przejście pomiędzy pasmem tłumienia i przenoszenia przy tym samym rzędzie w porównaniu z wyżej wymienionymi filtrami, tętnienia są obecne zarówno w paśmie przenoszenia jak i w paśmie tłumienia
[n,Wn]=ellipord(Wp, Ws, Rp,Rs);
[b,a]=ellip(n, rp, rs, Wn, btype='low', analog=0, output='ba')
Filtrowanie z zerowym przesunięciem fazowym
Filtrowanie sygnałów off-line można zrealizować tak, aby sygnał wyjściowy nie miał przesunięcia fazowego. Metodę tę zaprezentujemy na poniższym przykładzie. Rozważymy co dzieje się z pojedynczą składową fourierowską sygnału o częstości [math]\omega[/math].
- Niech ta składowa nazywa się
- [math]s_{\omega} (t) = A_{\omega} e^{i \omega t} [/math],
- natomiast transmitancja filtru dla częstości
- [math]H(\omega) = |H_{\omega}| e^{i \phi_{\omega}}[/math].
- Podając tę składową na wejście filtru, na jego wyjściu otrzymujemy [math]s^'_{\omega}(t)[/math]
- [math]s^'_{\omega}(t) = A_{\omega}|H_{\omega}| e^{i \omega t}e^{i \phi_{\omega}} [/math]
- korzystając z definicji przesunięcia fazowego możemy to zapisać:
- [math]s^'_{\omega}(t) = A_{\omega}|H_{\omega}| e^{i \omega( t- \tau_{\omega})} [/math]
- w sygnale tym odwracamy sekwencję czasową (kierunek czasu):
- [math]s^{''}_{\omega}(t) =s^{'}_{\omega}(-t) = A_{\omega}|H_{\omega}| e^{-i \omega( t- \tau_{\omega})} [/math]
- otrzymany sygnał podajemy na wejście filtru, wówczas na wyjściu pojawi się następujący sygnał:
- [math]s^{'''}_{\omega}(t) = A_{\omega}|H_{\omega}||H_{\omega}| e^{-i \omega( t- \tau_{\omega})}e^{i \phi_{\omega}} [/math]
- ponownie korzystając z definicji przesunięcia fazowego możemy napisać:
- [math]s^{'''}_{\omega}(t) = A_{\omega}|H_{\omega}|^2 e^{-i \omega t + \omega\tau_{\omega}- \omega\tau_{\omega}}= A_{\omega}|H_{\omega}|^2 e^{-i \omega t } [/math]
- w tym sygnale ponownie odwracamy sekwencję czasową i otrzymujemy:
- [math]y(t) = s^{'''}_{\omega}(-t) =A_{\omega}|H_{\omega}|^2 e^{i \omega t } [/math]
- Jest to sygnał wyjściowy.
- Zauważmy, że względem sygnału wejściowego ma on amplitudę zmienioną [math]|H_{\omega}|^2[/math] razy i nie zmienioną fazę.
Powyższe rozważania są prawdziwe dla dowolnego sygnału, gdyż z poprzednich zajęć wiemy, że działanie filtra (systemu LTI) można analizować niezależnie dla każdej składowej fourierowskiej.
Przykład
Procedura powyższa zaimplementowana jest w funkcji: scipy.signal.filtfilt. Jej działanie i porównanie z efektami funkcji lfilter przedstawia poniższy przykład:
from scipy.signal import filtfilt, butter, freqz, lfilter
import numpy as np
import pylab as py
# częstość próbkowania
Fs = 100.0
# projekt dolnoprzepustowego filtra Butterwortha 5 rzędu
# o częstości odcięcia 10 Hz
[b,a] = butter(5,10.0/(Fs/2.0))
# obliczamy funkcję przenoszenia
w,h = freqz(b,a,500)
transmitancja = np.abs(h)
#opóźnienie grupowe
grupowe = -np.diff(np.unwrap(np.angle(h)))/np.diff(w)
# przeliczamy skalę częstości na Hz
f = w/(np.pi)*Fs/2.0
# generujemy sygnał
t = np.arange(0,1,1/Fs)
s = np.sin(2*np.pi*5*t)*np.hanning(len(t))
# Filtrowanie z zerowym opoznieniem fazowym
y = filtfilt(b,a,s)
# Filtrowanie standardowe
y1 = lfilter(b,a,s)
# WYKRESY
py.subplot(2,2,1)
py.plot(f,20*np.log10(transmitancja)) # przeliczenie modułu transmitancji na dB
py.title('transmitancja')
py.xlabel('[Hz]')
py.ylabel('[dB]')
py.subplot(2,2,3)
py.plot(f[:-1], grupowe )
py.title('opoznienie grupowe')
py.xlabel('[Hz]')
py.ylabel('punkty')
py.subplot(1,2,2)
py.plot(t,s)
py.plot(t,y,'g.')
py.plot(t,y1,'r')
py.legend(('s','filtfilt','lfilter'))
py.xlabel('[s]')
py.title('sygnaly')
py.show()
Zadania
Cel: zaobserwowanie konieczności kompromisu pomiędzy ostrym cięciem w dziedzinie częstości a "dobrymi" parametrami w dziedzinie czasu - funkcja odpowiedzi impulsowej i schodkowej.
Zadania
- Skonstruować filtry dolnoprzepustowe rzędu [math]n=5[/math], o częstości odcięcia 30 Hz przy częstości próbkowania sygnału 128 Hz, rp = 0,5 dB, rs = 20 dB, przy pomocy wszystkich podanych powyżej funkcji i porównać ich własności.
- Dobrać rząd i zaprojektować, a następnie zbadać własności otrzymanego filtru Butterwortha spełniającego poniższe kryteria:
- pasmo przenoszenia 1000-2000 Hz, pasmo tłumienia zaczyna się 500 Hz od każdego z brzegów pasma przenoszenia,
- próbkowanie 10 kHz,
- najwyżej 1 dB tętnienia w paśmie przenoszenia,
- co najmniej 60 dB tłumienia w paśmie tłumienia.
- Zaprojektować filtr do wyławiania wrzecion snu z sygnału. Wrzeciona snu to struktury w sygnale EEG rejestrowanym w czasie snu zawierające się w paśmie 11-15 Hz. [1]. Filtr zaprojektuj w programie Svarog.
Przepróbkowywanie
Do góry: Zwiększamy częstość prókowania całkowitą ilość razy [math]P[/math].
Najpowszechniej stosowana metoda polega na dodaniu [math]P[/math] zer pomiędzy istniejące próbki sygnału tak aby osiągnął on [math]P[/math]-krotnie większą długość. Następnie taki rozciągnięty sygnał filtrujemy filtrem dolnoprzepustowym o częstości odcięcia nie większej niż częstość Nyquista oryginalnego sygnału — rozciąganie sygnału nie dokłada do niego nowej informacji więc i tak nic nie tracimy.
Przykład przepróbkowania do góry:
from scipy.signal import filtfilt, butter
import numpy as np
import pylab as py
t = np.arange(0,0.05,0.001) # czas
x = np.sin(2*np.pi*30*t) + np.sin(2*np.pi*60*t) # sygnał
py.subplot(3,1,1)
py.plot(x,'.')
py.subplot(3,1,2)
X = np.zeros(4*len(x))
X[::4] = x
py.plot(X,'.')
[b,a] = butter(8,1.0/4)
y = filtfilt(b,a,X);
py.subplot(3,1,3)
py.plot(y,'.')
py.show()
Do dołu: Zmniejszamy częstość próbkowania całkowitą ilość razy.
Musimy pamiętać o tym, żeby wyfiltrować to, co było w oryginalnym sygnale powyżej docelowego Nyquista, żeby uniknąć aliasingu w wynikowym sygnale.
from scipy.signal import filtfilt, butter
from numpy import pi, arange, sin
from pylab import plot, subplot, show, title
Fs1 = 128.0 # pierwotna częstość próbkowania [Hz]
FN1 = Fs1/2 # pierwotna częstość Nyquista
t = arange(0,1,1.0/Fs1) # czas probkowany 1/Fs1
f1 = 6 # Hz
f2 = 60
fi = pi/2
s = sin(2*pi*t*f1+fi) + sin(2*pi*t*f2+fi)
subplot(4,1,1)
plot(t,s,'.-')
title(u'sygnał pierwotny')
# obnizamy czestosc probkowania k razy
k = 2
Fs2 = Fs1/k # nowa czestosc probkowania jest k razy niższa
FN2 = Fs2/2 # nowa częstość Nyquista
[b,a] = butter(8,FN2/FN1) # przefiltrujemy filtrem dolnoprzepustowym
# tak aby nic nie zostało powyzej
# docelowej częstości Nyquista
ss = filtfilt(b,a,s);
t2 = arange(0,1,1.0/Fs2)
subplot(4,1,2)
plot(t,ss,'.-')
title(u'sygnał przefiltrowany dolnoprzepustowy')
subplot(4,1,3)
ss2 = ss[::k]
plot(t2,ss2,'.-')
title(u'sygnał przepróbkowany prawidłowo')
subplot(4,1,4)
ss3 = s[::k]
plot(t2,ss3,'.-')
title(u'sygnał przepróbkowany nieprawidłowo, bez filtrowania dolnoprzepustowego')
show()
Zmiana częstości o wymierną ilość razy:
Zmieniamy częstość próbkowania o wymierną [math]\frac{P}{Q}[/math] ilość razy — uzyskujemy składając powyższe kroki tzn. najpierw zwiększamy częstość [math]P[/math]-krotnie, a następnie zmniejszamy [math]Q[/math]-krotnie.
