Ćwiczenia 7

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

1. Zaprojektuj i zbadaj własności filtru:
   • FIR dolno z pasmem przenoszenia od 30 Hz dla sygnału próbkowanego 256 Hz
2. 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.