Ćwiczenia 7: Różnice pomiędzy wersjami

Z Brain-wiki
Linia 1: Linia 1:
==Filtry IIR==
+
=Implementacja filtrowania: funkcja lfilter=
 +
Filtrowanie zgodne z powyższymi równaniami zaimplementowane jest w pythonie w module <tt>scipy.signal</tt> jako funkcja [http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html lfilter].
 +
 
 +
Podstawowe wywołanie tej funkcji dla sygnału <tt>we</tt> wygląda następująco:
 +
<tt>wy = scipy.signal.lfilter(b,a,we)</tt> gdzie <tt>b, a</tt> są to współczynniki z poprzedniego równania.
 +
np:
 +
<source lang = python>
 +
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()
 +
</source>
 +
 
 +
==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 <tt>scipy.signal</tt> 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
 +
<tt>scipy.signal.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=1.0)</tt>.
 +
Pozwala ona projektować filtry dolno- i górno- przepustowe oraz pasmowo przepustowe i pasmowo zaporowe metodą okienkową.
 +
 
 +
Najważniejsze parametry (kompletny opis w dokumentacji [http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html#scipy.signal.firwin firwin]):
 +
* <tt>numtaps</tt>: int, ilość współczynników filtru (rząd filtru+1).  Liczba ta musi być parzysta jeśli pasmo przenoszenia ma zawierać częstość Nyquista.
 +
* <tt>cutoff</tt>: 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 <tt>nyq</tt> i musi zawierać się pomiędzy 0 a <tt>nyq</tt>.
 +
* <tt>window</tt>: napis lub krotka: określa jakiego okna użyć do projektu filtru. Może to być dowolne okno spośród opisanych w <tt>scipy.signal.get_window</tt>
 +
* <tt>pass_zero</tt>: 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 <tt>cutoff</tt> - czy ma to być pasmo przenoszone czy tłumione.
 +
* <tt>nyq</tt>: float. Częstość Nyquista.
 +
* Zwraca: współczynniki <tt>b</tt>
 +
 
 +
===== 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:
 +
<tt>firwin(21, 40/128.0, nyq= 1)</tt>
 +
<tt>firwin(21, 40, nyq= 128)</tt>
 +
* filtr górnoprzepustowy rzędu 15 z częstością odcięcia 5 Hz:
 +
<tt>firwin(16, 5, nyq= 128, pass_zero=False)</tt>
 +
* pasmowo przepustowy 5 rzędu przenoszący częstości pomiędzy 8 a 14 Hz:
 +
<tt>firwin(6, [8 14], nyq= 128, pass_zero=False)</tt>
 +
* pasmowo zaporowy 5 rzędu tłumiący częstości pomiędzy 8 a 14 Hz:
 +
<tt>firwin(6, [8 14], nyq= 128, pass_zero=True)</tt>
 +
 
 +
Demo własności "audio" filtrów tego typu: [http://www.falstad.com/dfilter/ applet ]
 +
 
 +
==== firwin2 ====
 +
Funkcja [http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin2.html#scipy.signal.firwin2 <tt>scipy.signal.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0)</tt>] 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:
 +
<tt>freq</tt> i <tt>gain</tt>. Wektor <tt>freq</tt> definiuje punkty w częstości (jednostki takie same jak <tt>nyq</tt>, muszą zawierać 0 i <tt>nyq</tt>) dla których znana jest wartość pożądanego przenoszenia. Pożądane wartości przenoszenia odpowiadające kolejnym częstościom definiowane są w <tt>gain</tt>. Wartości <tt>freq</tt> 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 <tt>gain</tt> aby zdefiniować nieciągłość funkcji przenoszenia. Znaczenie pozostałych parametrów jest takie samo jak dla <tt>firwin</tt>.
 +
 
 +
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:
 +
<source lang = python>
 +
f = np.array([0, 20, 40, 128])
 +
g = np.array([0, 0,  1,  1])
 +
b = firwin2(50,f,g,nfreqs, nyq=128)# licznik
 +
</source>
 +
 
 +
===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 <tt>scipy.signal</tt> zaimplementowana jest funkcja <tt>freqz</tt>, która wylicza wartości funkcję przenoszenia filtru zadanego współczynnikami <tt>b, a</tt> w zadanej ilości dyskretnych częstości. Jej użycie ilustruje poniższy przykład:
 +
<!--
 +
<source lang= python>
 +
import numpy as np
 +
import pylab as py
 +
from  scipy.signal import freqz
 +
 
 +
 
 +
b = np.array([0.7,  0.5]) # licznik
 +
a = np.array([1.0, -0.9]) # mianownik
 +
n = 100 # 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(2,1,1)
 +
py.semilogy(w,m)
 +
py.ylabel('[dB]')
 +
py.title('moduł transmitancji')
 +
py.subplot(2,1,2)
 +
py.plot(w,phi)
 +
py.title('faza transmitancji')
 +
py.ylabel('rad')
 +
py.xlabel('częstość [rad]')
 +
py.show()
 +
 
 +
</source>
 +
-->
 +
====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 =====
 +
<source lang = python>
 +
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()
 +
</source>
 +
 
 +
=====W dziedzinie czasu =====
 +
Porównajmy działanie tego filtra w czasie z funkcjami przesunięć fazowych i grupowych:
 +
 
 +
<source lang = python>
 +
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()
 +
</source>
 +
 
 +
==Zadania==
 +
<!--# Zbadaj własności filtrów podanych jako [[Ćwiczenia_6#Przyk.C5.82adowe_projekty|przykładowe projekty]] -->
 +
# Zaprojektuj i zbadaj własności filtru:
 +
<!--#* FIR <math>48</math> rzędu z pasmem przenoszenia 20 -40 Hz dla sygnału próbkowanego 250 Hz
 +
-->
 +
#* 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.
 
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.

Wersja z 16:29, 30 lis 2016

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.