Ćwiczenia 7

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 dziedzinei czasu - funkcja odpowiedzi impulsjowej 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]

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
# sygnal
x = np.sin(2*np.pi*30*t) + np.sin(2*np.pi*60*t)

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.