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

Z Brain-wiki
(Utworzono nową stronę "==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 odpowi...")
 
 
(Nie pokazano 19 wersji utworzonych przez 2 użytkowników)
Linia 1: Linia 1:
==Filtry IIR==
+
=Notebook=
  
 +
[https://drive.google.com/open?id=0BzwQ_Lscn8yDdnphNl9JNjA3dHc notebook o filtrach]
 +
 +
Ustawiamy parametry wyświetlania i font.
 +
 +
<source lang="python"># encoding: utf-8
 +
% matplotlib inline
 +
import matplotlib
 +
 +
matplotlib.rcParams['figure.figsize'] = (7,7)#(10, 7)
 +
matplotlib.rcParams.update({'font.family': 'Arial'})
 +
matplotlib.rcParams.update({'font.size': 10})
 +
 +
from IPython.core.display import display, HTML
 +
display(HTML("<style>.container { width:100% !important; }</style>"))</source>
 +
Importujemy podstawowe moduły:
 +
 +
<source lang="python">import numpy as np
 +
import pylab as py</source>
 +
Importujemy funkcje specyficzne dla filtrowania w pythonie:
 +
 +
<source lang="python">from  scipy.signal import freqz, group_delay #funkcja obliczająca funkcję systemu
 +
from  scipy.signal import firwin, firwin2    # funkcje do projektowania filtrów FIR
 +
from  scipy.signal import butter, buttord    # funkcje do projektowania filtrów 
 +
from  scipy.signal import cheby1, cheb1ord    # funkcje do projektowania filtrów
 +
from  scipy.signal import cheby2, cheb2ord    # funkcje do projektowania filtrów
 +
from  scipy.signal import ellip, ellipord    # funkcje do projektowania filtrów eliptycznych
 +
from  scipy.signal import lfilter, filtfilt # funkcje do aplikowania filtrów</source>
 +
= Implementacja filtrowania: funkcja lfilter =
 +
 +
== Dla przypomnienia: ==
 +
 +
=== Działanie filtra w dziedzinie czasu ===
 +
 +
Najczęściej, wyjście filtra jest kombinacją liniową:
 +
 +
<math>y[n] = \sum_{i=0}^{n_b}b(i)x[n-i] - \sum_{i=1}^{n_a}a(i)y[n-i]</math>
 +
 +
gdzie:
 +
 +
* <math>n_b</math> liczba przeszłych próbek wejściowych <math>x</math>
 +
* <math>n_a</math> liczba przeszłych próbek wyjściowych <math>y</math>
 +
 +
użytych do obliczenia aktualnego wyjścia <math>y[n]</math>.
 +
 +
Większa z liczb <math>n_b</math> i <math>n_a</math> określa &quot;rząd&quot; filtra.
 +
 +
Zauważmy, że matematycznie operacje te odpowiadają splataniu próbek wejściowych z wektorem <math>b</math> i próbek wyjściowych z wektorem <math>a</math>.
 +
 +
=== Implementacja w pythonie ===
 +
 +
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:
 +
 +
<code>wy = scipy.signal.lfilter(b,a,we)</code>
 +
 +
gdzie b, a są to współczynniki z poprzedniego równania.
 +
 +
= Badanie własności filtra w dziedzinie czasu i częstości: =
 +
 +
Przy projektowaniu filtra musimy brać pod uwagę jego następujące własności: - w dziedzinie częstości: - moduł transmitancji (funkcji przenoszenia) - jest to zależność modułu funkcji systemu od częstości <math>\left|H\left(f\right)\right|</math> - opóźnienie grupowe - opóźnienie fazowe - w dziedzinie czasu: - funkcję odpowiedzi impulsowej - funkcję odpowiedzi schodkowej
 +
 +
== Zadanie: budujemy funkcję do ilustracji własności filtra ==
 +
 +
Nasza funkcja będzie przyjmowała na wejściu: * współczynniki filtra <code>a, b</code>, * wektor zawierający częstości <code>f</code>, dla których własności mają być policzone, * długość odcinka czasu, <code>T</code>, na którym mają być prezentowane własności czasowe filtra,<br />
 +
* oraz częstość próbkowania <code>Fs</code>.
 +
 +
Funkcja ta będzie rysowała wszystkie powyżej wspomniane charakterystyki filtra.
 +
 +
<source lang="python">def charkterystyki(a,b,f,T,Fs):
 +
    # przyda nam się oś czasu od -T do T sekund
 +
    t = np.arange(-T, T, 1/Fs)
 +
   
 +
    # oś częstości przeliczamy na radiany
 +
    w = 2*np.pi* f/Fs
 +
   
 +
    # obliczamy transmitancję
 +
    w, h = freqz(...)
 +
   
 +
    # obliczamy moduł transmitancji
 +
    m = np.abs(h)
 +
       
 +
    # obliczamy fazę i "rozwijamy" ją 
 +
    faza = np.unwrap(np.angle(h))
 +
 +
    # obliczamy opóźnienie fazowe
 +
    opoznienieFazowe = - faza/w
 +
   
 +
    # obliczamy opóźnienie grupowe
 +
    df = np.diff(faza)
 +
    idx, = np.where(np.abs(df-np.pi)<0.05) #to zabezpieczenie na błędy przy "rozwijaniu" fazy
 +
    df[idx] = (df[idx+1]+df[idx-1])/2
 +
    grupowe = - df/np.diff(w)
 +
   
 +
    # obliczamy odpowiedź impulsową
 +
    x = np.zeros(len(t))
 +
    x[len(t)//2] = 1 # impuls
 +
    y = lfilter(b,a,x)
 +
   
 +
    # obliczamy odpowiedź schodkową
 +
    s = np.zeros(len(t))
 +
    s[len(t)//2:] = 1 # schodek
 +
    ys = lfilter(b,a,s) # przepuszczamy impuls przez filtr i obserwujemy odpowiedź impulsową
 +
   
 +
    # rysujemy
 +
    fig = py.figure()
 +
    py.subplot(3,2,1)
 +
    py.title('moduł transmitancji')
 +
    py.plot(f,20*np.log10(m))
 +
    py.ylabel('[dB]')
 +
    py.grid('on')
 +
   
 +
    py.subplot(3,2,3)
 +
    py.title('opóźnienie fazowe')
 +
    py.plot(f, opoznienieFazowe)
 +
    py.ylabel('próbki')
 +
    py.grid('on')
 +
   
 +
    py.subplot(3,2,5)
 +
    py.title('opóźnienie grupowe')
 +
    py.plot(f[:-1],grupowe)
 +
    py.ylabel('próbki')
 +
    py.xlabel('Częstość [Hz]')
 +
    py.grid('on')
 +
    py.ylim([0, np.max(grupowe)+1])
 +
   
 +
    py.subplot(3,2,2)
 +
    py.title('odpowiedź impulsowa')
 +
    py.plot(t, x)
 +
    py.plot(t, y)
 +
    py.xlim([-T/2,T])
 +
    py.grid('on')
 +
   
 +
    py.subplot(3,2,4)
 +
    py.title('odpowiedź schodkowa')
 +
    py.plot(t, s)
 +
    py.plot(t, ys)
 +
    py.xlim([-T/2,T])
 +
    py.xlabel('Czas [s]')
 +
    py.grid('on')
 +
   
 +
    fig.subplots_adjust(hspace=.5)
 +
    py.show()
 +
    </source>
 +
Testujemy:
 +
 +
<source lang="python">b = firwin(54,0.5)# licznik
 +
a = np.array([1.0]) # mianownik
 +
Fs = 100
 +
T = 1
 +
f = np.arange(0.01,Fs/2,0.01)
 +
charkterystyki(a,b,f,T,Fs)</source>
 +
= Funkcje do projektowania filtrów =
 +
 +
W bibliotece <code>scipy.signal</code> jest kilka funkcji do projektowania filtrów o zadanych parametrach. My skupimy się na dwóch zasadniczych grupach: * FIR - filtry o skończonej odpowiedzi impulsowej * klasyczne IIR - filtry o nieskończonej odpowiedzi impulsowej
 +
 +
== FIR ==
 +
 +
Filtry typu FIR zwykle wymagają znacznie wyższych rzędów aby osiągnąć transmitancję o porządanej formie. Mają jednak dwie podstawowe zalety: * ich funkcja odpowiedzi jest skończona opisana wektorem <code>b</code> - efekty brzegowe sięgają z obu końców filtrowanego sygnału na dokładnie połowę długości wektora <code>b</code> * mają liniową zależnaość fazy od częstości. Z tego powodu opóżnienie grupowe dla wszystkich częstości jest takie samo.
 +
 +
W module <tt>scipy.signal</tt> mamy kilka funkcji do projektowania filtrów typu FIR: <code>firwin</code> i <code>firwin2</code>.
 +
 +
=== firwin ===
 +
 +
Najprostszą koncepcyjnie metodą projektowania filtrów FIR jest metoda okienkowa.<br />
 +
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
 +
 +
* <code>numtaps</code>: 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łady: Zbadaj włsności przykładowych projektów ===
 +
 +
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:
 +
* <code>firwin(21, 40, nyq= 128)</code>
 +
 +
<source lang="python">Fs = 128
 +
T = 0.2
 +
f = np.arange(0.01,Fs/2,0.01)
 +
b = firwin(21, 40, nyq= 128)
 +
charkterystyki(1,b,f,T ,Fs)</source>
 +
* filtr górnoprzepustowy rzędu 15 z częstością odcięcia 5 Hz:
 +
* <code>firwin(17, 15, nyq= 128, pass_zero=False)</code>
 +
 +
<source lang="python">Fs = 256
 +
T = 1
 +
f = np.arange(0.01,Fs/2,0.01)
 +
 +
b = firwin(17, 15, nyq= 128, pass_zero=False)
 +
charkterystyki(1,b,f,T = 0.2,Fs=Fs)</source>
 +
* pasmowo przepustowy 51 rzędu przenoszący częstości pomiędzy 8 a 14 Hz:
 +
* <code>firwin(51, [8, 14], nyq= 128, pass_zero=False)</code>
 +
 +
<source lang="python">Fs = 256
 +
f = np.arange(0.01,Fs/2,0.01)
 +
 +
b=firwin(51, [8, 14], nyq= Fs/2, pass_zero=False)
 +
charkterystyki(1,b,f,T = 0.2,Fs=Fs)</source>
 +
=== Zadanie: Zaprojektuj i zbadaj własności filtra: ===
 +
 +
FIR dolno z pasmem przenoszenia od 30 Hz dla sygnału próbkowanego 256 Hz
 +
 +
<source lang="python"></source>
 +
=== Zadanie: Znajdź rząd filtra FIR: ===
 +
 +
dolnoprzepustowego z pasmem przenoszenia do 40 Hz dla sygnału próbkowanego 256 Hz, tak aby dla częstości powyżej 45 Hz jego tłumienie było nie mniejsze niż 20dB.
 +
 +
<source lang="python"></source>
 +
=== firwin2 ===
 +
 +
Funkcja
 +
 +
<code>scipy.signal.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0)</code>
 +
 +
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: * <code>freq</code> i Wektor freq definiuje punkty w częstości (jednostki takie same jak <code>nyq</code>, muszą zawierać 0 i <code>nyq</code>) dla których znana jest wartość pożądanego przenoszenia. Wartości <code>freq</code> 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. * <code>gain</code> Pożądane wartości przenoszenia odpowiadające kolejnym częstościom definiowane są w <code>gain</code>. * Znaczenie pozostałych parametrów jest takie samo jak dla ``firwin.
 +
 +
=== Zadanie: filtr wielopasmowy ===
 +
 +
Zaprojektuj filtr przenoszący częstości w pasmach pomiędzy : 10-11, 20-21 i 30-31 Hz, który w paśmie zaporowym ma co najmniej 60 dB tłumienia.
 +
 +
<source lang="python">Fs =100
 +
T = 2
 +
f = np.arange(0.01,Fs/2,0.01)
 +
freq = np.array([0, 10, 10, 11, 11, 20, 20, 21, 21, 30, 30, 31, 31, 50])
 +
gain = np.array([0,  0,  1,  1,  0,  0,  1,  1,  0,  0,  1,  1,  0,  0])
 +
b = firwin2(100, freq, gain, nyq= 50)
 +
charkterystyki(1,b,f,T,Fs)</source>
 +
= 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.
  
W module <tt>scipy.signal</tt> mamy zaimplementowane kilka funkcji do projektowania &bdquo;optymalnych&rdquo; pod różnymi względami filtrów w klasycznych konfiguracjach:
+
W module <tt>scipy.signal</tt> 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.
dolno- albo górnoprzepustowe i pasmowo-przepustowe albo pasmowo-zaporowe.
+
 
 +
== Funkcje do projektowania filtrów IIR dostępne w module <tt>scipy.signal</tt> ==
 +
 
 +
W module <tt>scipy.signal</tt> 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: * <tt>wp</tt>, <tt>ws</tt> — 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: <tt>wp = 0.2</tt>, <tt>ws = 0.3</tt> * * górno-przepustowy: <tt>wp = 0.3</tt>, <tt>ws = 0.2</tt> * * pasmowo-przepustowy: <tt>wp = [0.2, 0.5]</tt>, <tt>ws = [0.1, 0.6]</tt> * * pasmowo-zaporowy: <tt>wp = [0.1, 0.6]</tt>, <tt>ws = [0.2, 0.5]</tt> * <tt>gpass</tt> — maksymalna dopuszczalna strata w pasmie przenoszenia (w funkcjach projektujących filtry jest to <tt>rp</tt>) (dB). * <tt>gstop</tt> — minimalne wymagane tłumienie w pasmie tłumienia (w funkcjach projektujących filtry jest to <tt>rs</tt>) (dB). * <tt>btype</tt> — typ filtra (<tt>'lowpass'</tt>, <tt>'highpass'</tt>, <tt>'bandpass'</tt>, <tt>'bandstop'</tt>).
  
==Funkcje do projektowania filtrów IIR dostępne w module <tt>scipy.signal</tt>==
+
Funkcje do projektowania filtrów są zaimplementowane parami: * jedna pomaga dobierać rząd filtru do wymagań projektowych, * druga oblicza współczynniki filtru.
W module <tt>scipy.signal</tt> dostępne 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:
 
* <tt>wp</tt>, <tt>ws</tt> &mdash; 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:  <tt>wp = 0.2</tt>,          <tt>ws = 0.3</tt>
 
** górno-przepustowy:  <tt>wp = 0.3</tt>,          <tt>ws = 0.2</tt>
 
** pasmowo-przepustowy:  <tt>wp = [0.2, 0.5]</tt>,  <tt>ws = [0.1, 0.6]</tt>
 
** pasmowo-zaporowy:  <tt>wp = [0.1, 0.6]</tt>,  <tt>ws = [0.2, 0.5]</tt>
 
* <tt>gpass</tt> &mdash; maksymalna dopuszczalna strata w pasmie przenoszenia (w funkcjach projektujących filtry jest to <tt>rp</tt>) (dB).
 
* <tt>gstop</tt> &mdash; minimalne wymagane tłumienie w pasmie tłumienia (w funkcjach projektujących filtry jest to <tt>rs</tt>) (dB).
 
* <tt>btype</tt> &mdash; typ filtra (<tt>'lowpass'</tt>, <tt>'highpass'</tt>, <tt>'bandpass'</tt>, <tt>'bandstop'</tt>).
 
  
Funkcje do projektowania filtrów są zaimplementowane parami: jedna pomaga dobierać rząd filtru do wymagań projektowych, a druga oblicza współczynniki filtru.
+
W poniższych przykładach przyjmiemy:
  
 +
<source lang="python">T = 0.3
 +
Fs = 100 # Hz
 +
f = np.arange(0.01,Fs/2,0.01)
  
'''Filtr Butterwortha'''
+
wp = 10/(Fs/2)
daje gładką (bez tętnień) funkcję przenoszenia w całym zakresie częstości
+
ws = 30/(Fs/2)
<source lang=python>>
+
gpass = 1
[n,Wn]=buttord(wp, ws, gpass, gstop, analog=0)
+
gstop = 25
 +
analog=0
 +
rp = gpass
 +
rs = gstop</source>
 +
=== Filtr Butterwortha ===
 +
 
 +
daje gładką (bez tętnień) funkcję przenoszenia w całym zakresie częstości.
 +
 
 +
<source lang="python">[n,Wn]=buttord(wp, ws, gpass, gstop)
 
[b,a]=butter(n,Wn)
 
[b,a]=butter(n,Wn)
</source>
+
charkterystyki(a,b,f,T,Fs)
<tt>n</tt> &mdash; rząd filtra, <tt>Wn</tt> &mdash; krawędzie pasm
+
print('Filtr Butterwortha, rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))</source>
 +
=== Filtr Czebyszewa I rodzaju ===
 +
 
 +
gładka funkcja przenoszenia w paśmie tłumienia, minimalizowane są tętnienia w paśmie przenoszenia
 +
 
 +
<source lang="python">[n,Wn]=cheb1ord(wp, ws, gpass, gstop, analog)
 +
[b,a]=cheby1(n, rp, Wn, btype='low',  output='ba')
 +
charkterystyki(a,b,f,T,Fs)
 +
print('Czebyszewa I Typu: rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))</source>
 +
=== Filtr Czebyszewa II rodzaju ===
 +
 
 +
gładka funkcja przenoszenia w paśmie przenoszenia, minimalizowane tętnienia w paśmie tłumienia
  
'''Filtr Czebyszewa I rodzaju''' &mdash; gładka funkcja przenoszenia w paśmie tłumienia, minimalizowane są tętnienia  w paśmie przenoszenia
+
<source lang="python">[n,Wn]=cheb2ord(wp, ws, gpass, gstop, analog=0);
<source lang=python>>
 
[n,Wn]=cheb1ord(wp, ws, gpass, gstop, analog=0);
 
[b,a]=cheby1(n, rp, Wn, btype='low', analog=0, output='ba')
 
</source>
 
'''Filtr Czebyszewa II''' rodzaju &mdash; gładka funkcja przenoszenia w paśmie przenoszenia, minimalizowane tętnienia w paśmie tłumienia
 
<source lang=python>
 
[n,Wn]=cheb2ord(wp, ws, gpass, gstop, analog=0);
 
 
[b,a]=cheby2(n, rs, Wn, btype='low', analog=0, output='ba')
 
[b,a]=cheby2(n, rs, Wn, btype='low', analog=0, output='ba')
</source>
+
charkterystyki(a,b,f,T,Fs)
'''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
+
print('Czebyszewa II Typu: rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))</source>
<source lang=python>
+
=== Filtr eliptyczny ===
[n,Wn]=ellipord(Wp, Ws, Rp,Rs);
+
 
 +
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
 +
 
 +
<source lang="python">[n,Wn]=ellipord(wp, ws, rp,rs);
 
[b,a]=ellip(n, rp, rs, Wn, btype='low', analog=0, output='ba')
 
[b,a]=ellip(n, rp, rs, Wn, btype='low', analog=0, output='ba')
</source>
+
charkterystyki(a,b,f,T,Fs)
 
+
print('Eliptyczny: rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))</source>
==Filtrowanie z zerowym przesunięciem fazowym==
+
== 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>.
+
=== Zadanie ===
* 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.
+
Filtrowanie sygnałów off-line można zrealizować tak, aby sygnał wyjściowy nie miał przesunięcia fazowego. 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:
=== Przykład ===
 
Procedura powyższa zaimplementowana jest w funkcji: <tt>scipy.signal.filtfilt</tt>.
 
Jej działanie i porównanie z efektami funkcji <tt>lfilter</tt> przedstawia poniższy przykład:
 
<source lang = 'py'>
 
from scipy.signal import filtfilt, butter, freqz,lfilter
 
import numpy as np
 
import pylab as py
 
  
# częstość próbkowania
+
<source lang="python"># częstość próbkowania
 
Fs = 100.0
 
Fs = 100.0
 
# projekt dolnoprzepustowego filtra Butterwortha 5 rzędu
 
# projekt dolnoprzepustowego filtra Butterwortha 5 rzędu
 
# o częstości odcięcia 10 Hz   
 
# o częstości odcięcia 10 Hz   
[b,a] = butter(5,10.0/(Fs/2.0))
+
[b,a] = butter(...)
 
+
 
# obliczamy funkcję przenoszenia
 
# obliczamy funkcję przenoszenia
w,h = freqz(b,a,500)
+
w,h = freqz(...)
transmitancja = np.abs(h)
+
transmitancja = ...
 
+
 
#opóźnienie grupowe
 
#opóźnienie grupowe
 
grupowe = -np.diff(np.unwrap(np.angle(h)))/np.diff(w)   
 
grupowe = -np.diff(np.unwrap(np.angle(h)))/np.diff(w)   
 
+
# przeliczamy skalę częstości na Hz  
+
# przeliczamy skalę w (radiany) na częstości w Hz  
f = w/(np.pi)*Fs/2.0
+
f = ...
 
+
 
# generujemy sygnał
 
# generujemy sygnał
 
t = np.arange(0,1,1/Fs)
 
t = np.arange(0,1,1/Fs)
 
s = np.sin(2*np.pi*5*t)*np.hanning(len(t))
 
s = np.sin(2*np.pi*5*t)*np.hanning(len(t))
 
+
 
# Filtrowanie z zerowym opoznieniem fazowym
 
# Filtrowanie z zerowym opoznieniem fazowym
y = filtfilt(b,a,s)
+
y = filtfilt(...)
 
+
 
# Filtrowanie standardowe
 
# Filtrowanie standardowe
 
y1 = lfilter(b,a,s)
 
y1 = lfilter(b,a,s)
 
+
 
 
 
# WYKRESY
 
# WYKRESY
 
py.subplot(2,2,1)
 
py.subplot(2,2,1)
 
py.plot(f,20*np.log10(transmitancja)) # przeliczenie modułu transmitancji na dB
 
py.plot(f,20*np.log10(transmitancja)) # przeliczenie modułu transmitancji na dB
py.title('transmitancja')
+
py.title('moduł transmitancji')
py.xlabel('[Hz]')
 
 
py.ylabel('[dB]')
 
py.ylabel('[dB]')
 
+
 
py.subplot(2,2,3)
 
py.subplot(2,2,3)
 
py.plot(f[:-1], grupowe )
 
py.plot(f[:-1], grupowe )
Linia 116: Linia 339:
 
py.xlabel('[Hz]')
 
py.xlabel('[Hz]')
 
py.ylabel('punkty')
 
py.ylabel('punkty')
 
+
 
 
 
py.subplot(1,2,2)
 
py.subplot(1,2,2)
 
py.plot(t,s)
 
py.plot(t,s)
Linia 125: Linia 347:
 
py.xlabel('[s]')
 
py.xlabel('[s]')
 
py.title('sygnaly')
 
py.title('sygnaly')
py.show()
+
py.show()</source>
 +
=== Zadanie ===
 +
 
 +
Skonstruować filtry dolnoprzepustowe rzędu n=5, 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.
 +
 
 +
<source lang="python"></source>
 +
=== Zadanie ===
 +
 
 +
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.
 +
 
 +
<source lang="python"></source>
 +
=== Zadanie: filtr pasmowy do wyszukiwania wrzecion snu ===
 +
 
 +
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. Działanie filtra przetestować na sygnale:
 +
 
 +
http://www.fuw.edu.pl/~jarekz/SYGNALY/TF/c4spin.txt
 +
 
 +
Sygnał ten to fragment zapisu EEG z II stadium snu elektroda C4 próbkownie 128Hz.
 +
 
 +
<source lang="python">s = np.loadtxt('c4spin.txt') # wczytujemy sygnał z pliku tekstowego
 +
Fs = 128
 +
t = np.arange(0,len(s))/Fs # przygotowujemy oś czasu
 +
 
 +
py.plot(t, s)
 +
[b,a] = butter(...)
 +
sf = filtfilt(...)
 +
py.plot(t,sf)
 +
py.show()</source>
 +
=== Zadanie 5: uwaga na odpowiedź impulsową ===
 +
 
 +
Przydadzą nam się pliki: * https://drive.google.com/open?id=0BzwQ_Lscn8yDNGc0aU5jSDFFMmc Plik z sygnałem EKG * https://drive.google.com/open?id=0BzwQ_Lscn8yDOF9jX0pjcG9LSGc Plik z metadanymi do tego sygnału
 +
 
 +
Proszę: * wykreślić pierwsze 10 sekund sygnału * zastosować filtr górnoprzepustowy Butterwartha o częstościach odcięcia: 0.01, 0.1, 0.5 -&gt; zaobserwuj jak długo stabilizuje się sygnał * Zastosuj filtr pasmowoprzepustowy (11 - 14 Hz) i wykreśl wynik jego zastosowania na tle poprzedniej wersji sygnału. Zinterpretuj wynik w kontekście funkcji odpowiedzi impulsowej tego filtra.
  
</source>
+
<source lang="python">s = np.fromfile('EKG.bin', dtype='<f') # tworzymy tablicę sig o typie określonym przez ''dtype''                               
 +
# ustawiamy częstość próbkowania
 +
Fs = 128
 +
# tworzymy wektor czasu
 +
t = np.arange(0,len(s))/Fs
  
==Zadania==
+
# ustalamy zakres indeksów sygnału i czasu do wyświetlania
 +
zakres = np.logical_and(0<t, t<10)
  
 +
py.plot(t[zakres],s[zakres])
 +
# filtr górnoprzepustowy (0.01, 0.1, 0.5)
 +
[b,a] = butter(... )
 +
sf = filtfilt(...)
 +
py.plot(t[zakres],sf[zakres])
  
===Zadania===
+
# filtr pasmowy
 +
[bl,al] = butter(... )
 +
sf_l = filtfilt(bl,al,sf)
 +
py.plot(t[zakres],sf_l[zakres])
  
* 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.
+
py.show()</source>
* Dobrać rząd i zaprojektować, a następnie zbadać własności otrzymanego filtru Butterwortha spełniającego poniższe kryteria:
+
=== Zadanie: Uwaga na odpowiedź schodkową ===
**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. [http://brain.fuw.edu.pl/~jarek/SYGNALY/TF/c4spin.txt]
 
  
==Przepróbkowywanie==
+
Wykorzystajmy fragment sygnału EKG z poprzedniego zadania (pomiędzy 12 a 40 -tą sekundą). * wykreśl ten fragment * zaprojektuj filtr górnoprzepustowy Butterwortha o częstości odcięcia 0.1 (potem 0.5), rzedu 1 (potem 5), * przefiltruj sygnał z tymi współczynnikami za pomocą funkcji <code>filtfilt</code> i <code>lfilter</code> , * dodaj do sygnału z lewej strony jego kopię odwróconą w czasie, * ten sygnał przefiltruj funkcją <code>lfilter</code> i wykreśl jego drugą połowę, * zinterpretuj uzyskane wyniki w kontekście funkcji odpowiedzi impulsowej.
  
 +
<source lang="python">zakres =  np.logical_and(10<t, t<40)
 +
t_z = t[zakres]
 +
s_z = s[zakres]
 +
py.plot(t_z,s_z, label = 'oryginalny')
 +
py.grid('on')
  
===Do góry: Zwiększamy częstość prókowania całkowitą ilość razy <math>P</math>. ===
+
# filtr górnoprzepustowy ( 0.1, 0.5)
 +
[b,a] = butter(...)
 +
sff = filtfilt(...)
 +
py.plot(t_z, sff, label = 'filtfilt')
  
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 &mdash; rozciąganie sygnału nie dokłada do niego nowej informacji więc i tak nic nie tracimy.
+
# lfilter
 +
sfl = lfilter(...)
 +
py.plot(t_z,...,label = 'lfilter')
  
 +
# lfilter z przedłużeniem
 +
x=np.concatenate((s_z[::-1],s_z))
 +
sfl_p = lfilter(...)
 +
py.plot(t_z,sfl_p[len(t_z):],label = 'lfilter z przedłużaniem')
 +
py.legend()
 +
py.show()</source>
 +
== Przepróbkowywanie ==
  
 +
=== Przepróbkowywanie do góry: ===
  
===Przykład przepróbkowania do góry:===
+
Zwiększamy częstość prókowania całkowitą ilość razy P
<source lang=python>>
 
from scipy.signal import filtfilt, butter
 
import numpy as np
 
import pylab as py
 
  
t = np.arange(0,0.05,0.001)# czas
+
Najpowszechniej stosowana metoda polega na dodaniu P zer pomiędzy istniejące próbki sygnału tak aby osiągnął on P-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:
# sygnal
 
x = np.sin(2*np.pi*30*t) + np.sin(2*np.pi*60*t)
 
  
 +
<source lang="python">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.subplot(3,1,1)
 
py.plot(x,'.')
 
py.plot(x,'.')
 +
py.title('Sygnał oryginalny')
 +
 
py.subplot(3,1,2)
 
py.subplot(3,1,2)
 
X = np.zeros(4*len(x))
 
X = np.zeros(4*len(x))
 
X[::4] = x
 
X[::4] = x
 
py.plot(X,'.')
 
py.plot(X,'.')
[b,a]=butter(8,1.0/4)
+
py.title('Sygnał poprzetykany zerami')
y=filtfilt(b,a,X);
+
 
 +
[b,a] = butter(8,...) # filtr powinien przepuszczać tylko te częstości,
 +
                        # które były w oryginalnym sygnale tzn. poniżej pierwotnego Nyqista
 +
y = filtfilt(b,a,X); # filtrujemy dolnoprzepustowo 
 
py.subplot(3,1,3)
 
py.subplot(3,1,3)
 
py.plot(y,'.')
 
py.plot(y,'.')
+
py.show()</source>
 +
=== Przepróbkowanie do dołu: ===
  
py.show()
+
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.
</source>
 
 
 
===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.
 
<source lang=python>>
 
from scipy.signal import filtfilt, butter
 
from numpy import pi, arange, sin
 
from pylab import plot, subplot, show, title
 
 
 
Fs1=128.0 #perwitna częstość próbkowaniaHz
 
FN1=Fs1/2 # pierwotna częstość Nyquista
 
  
 +
<source lang="python">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
 
t = arange(0,1,1.0/Fs1) # czas probkowany 1/Fs1
 
f1 = 6 # Hz
 
f1 = 6 # Hz
Linia 195: Linia 469:
 
plot(t,s,'.-')
 
plot(t,s,'.-')
 
title(u'sygnał pierwotny')
 
title(u'sygnał pierwotny')
#obnizamy czestosc probkowania k razy
+
# obnizamy czestosc probkowania k razy
 
k = 2
 
k = 2
 
Fs2 = Fs1/k # nowa czestosc probkowania jest k razy niższa
 
Fs2 = Fs1/k # nowa czestosc probkowania jest k razy niższa
Linia 202: Linia 476:
 
                           # tak aby nic nie zostało powyzej
 
                           # tak aby nic nie zostało powyzej
 
                           # docelowej częstości Nyquista  
 
                           # docelowej częstości Nyquista  
ss=filtfilt(b,a,s);
+
ss = filtfilt(b,a,s);
t2=arange(0,1,1.0/Fs2)  
+
t2 = arange(0,1,1.0/Fs2)  
 
subplot(4,1,2)
 
subplot(4,1,2)
 
plot(t,ss,'.-')
 
plot(t,ss,'.-')
 
title(u'sygnał przefiltrowany dolnoprzepustowy')
 
title(u'sygnał przefiltrowany dolnoprzepustowy')
 
+
 
subplot(4,1,3)
 
subplot(4,1,3)
ss2=ss[::k]
+
ss2 = ss[::k]
 
plot(t2,ss2,'.-')
 
plot(t2,ss2,'.-')
 
title(u'sygnał przepróbkowany prawidłowo')
 
title(u'sygnał przepróbkowany prawidłowo')
 
+
 
subplot(4,1,4)
 
subplot(4,1,4)
ss3=s[::k]
+
ss3 = s[::k]
 
plot(t2,ss3,'.-')
 
plot(t2,ss3,'.-')
 
title(u'sygnał przepróbkowany nieprawidłowo, bez filtrowania dolnoprzepustowego')
 
title(u'sygnał przepróbkowany nieprawidłowo, bez filtrowania dolnoprzepustowego')
show()
+
show()</source>
</source>
+
=== Zmiana częstości o wymierną ilość razy: ===
 +
 
 +
Zmieniamy częstość próbkowania o wymierną <math>\frac{P}{Q}</math> liczbę razy — uzyskujemy składając powyższe kroki tzn. najpierw zwiększamy częstość P-krotnie, a następnie zmniejszamy Q-krotnie.
 +
 
 +
=== Zadanie 6 ===
 +
 
 +
Proszę napisać funkcję, która będzie przepróbkowywać sygnał o wymierną liczbę razy. Funkcja powinna przyjmować sygnał, częstość próbkowania, parametry P i Q i zwracać przepróbkowany sygnał i nową częstość próbkowania
  
===Zmiana częstości o wymierną ilość razy:===
+
<source lang="python">def resample(s,fs,P=1,Q=1):
 +
    if P>1 and isinstance(P,int):
 +
        sP = np.zeros(P*len(s))
 +
        sP[...] = s
 +
        fs = ...
 +
        [b,a] = butter(...)
 +
        s = filtfilt(...)
 +
    if Q>1 and isinstance(Q,int):
 +
        fs = fs/Q
 +
        [b,a] = butter(...)
 +
        s = filtfilt(...)
 +
        s = s[...]
 +
    return s,fs
  
Zmieniamy częstość próbkowania o wymierną <math>\frac{P}{Q}</math> ilość razy &mdash; 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.
+
fs = 1000
 +
t = np.arange(0,0.05,0.001) # czas
 +
s1 = np.sin(2*np.pi*30*t) + np.sin(2*np.pi*60*t) # sygnał
 +
 +
py.subplot(3,1,1)
 +
py.plot(s1,'.')
 +
y,fs2 = resample(s1,fs,P=10,Q=1)
 +
py.subplot(3,1,2)
 +
py.plot(y,'.')
 +
py.show()
 +
 
 +
fs = 128.0
 +
t = np.arange(0,1,1.0/fs)
 +
f1 = 6 # Hz
 +
f2 = 40
 +
fi = np.pi/2
 +
s2 = np.sin(2*np.pi*t*f1+fi) + np.sin(2*np.pi*t*f2+fi)
 +
py.subplot(3,1,1)
 +
py.plot(s2,'.-')
 +
y,fs2 = resample(s2,fs,P=1,Q=2)
 +
py.subplot(3,1,2)
 +
py.plot(y,'.-')
 +
py.show()</source>

Aktualna wersja na dzień 12:54, 12 gru 2016

Notebook

notebook o filtrach

Ustawiamy parametry wyświetlania i font.

# encoding: utf-8
% matplotlib inline
import matplotlib

matplotlib.rcParams['figure.figsize'] = (7,7)#(10, 7)
matplotlib.rcParams.update({'font.family': 'Arial'})
matplotlib.rcParams.update({'font.size': 10})

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

Importujemy podstawowe moduły:

import numpy as np
import pylab as py

Importujemy funkcje specyficzne dla filtrowania w pythonie:

from  scipy.signal import freqz, group_delay #funkcja obliczająca funkcję systemu
from  scipy.signal import firwin, firwin2     # funkcje do projektowania filtrów FIR
from  scipy.signal import butter, buttord     # funkcje do projektowania filtrów  
from  scipy.signal import cheby1, cheb1ord    # funkcje do projektowania filtrów 
from  scipy.signal import cheby2, cheb2ord    # funkcje do projektowania filtrów 
from  scipy.signal import ellip, ellipord     # funkcje do projektowania filtrów eliptycznych
from  scipy.signal import lfilter, filtfilt # funkcje do aplikowania filtrów

Implementacja filtrowania: funkcja lfilter

Dla przypomnienia:

Działanie filtra w dziedzinie czasu

Najczęściej, wyjście filtra jest kombinacją liniową:

[math]y[n] = \sum_{i=0}^{n_b}b(i)x[n-i] - \sum_{i=1}^{n_a}a(i)y[n-i][/math]

gdzie:

  • [math]n_b[/math] liczba przeszłych próbek wejściowych [math]x[/math]
  • [math]n_a[/math] liczba przeszłych próbek wyjściowych [math]y[/math]

użytych do obliczenia aktualnego wyjścia [math]y[n][/math].

Większa z liczb [math]n_b[/math] i [math]n_a[/math] określa "rząd" filtra.

Zauważmy, że matematycznie operacje te odpowiadają splataniu próbek wejściowych z wektorem [math]b[/math] i próbek wyjściowych z wektorem [math]a[/math].

Implementacja w pythonie

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.

Badanie własności filtra w dziedzinie czasu i częstości:

Przy projektowaniu filtra musimy brać pod uwagę jego następujące własności: - w dziedzinie częstości: - moduł transmitancji (funkcji przenoszenia) - jest to zależność modułu funkcji systemu od częstości [math]\left|H\left(f\right)\right|[/math] - opóźnienie grupowe - opóźnienie fazowe - w dziedzinie czasu: - funkcję odpowiedzi impulsowej - funkcję odpowiedzi schodkowej

Zadanie: budujemy funkcję do ilustracji własności filtra

Nasza funkcja będzie przyjmowała na wejściu: * współczynniki filtra a, b, * wektor zawierający częstości f, dla których własności mają być policzone, * długość odcinka czasu, T, na którym mają być prezentowane własności czasowe filtra,

  • oraz częstość próbkowania Fs.

Funkcja ta będzie rysowała wszystkie powyżej wspomniane charakterystyki filtra.

def charkterystyki(a,b,f,T,Fs):
    # przyda nam się oś czasu od -T do T sekund
    t = np.arange(-T, T, 1/Fs)
    
    # oś częstości przeliczamy na radiany
    w = 2*np.pi* f/Fs
    
    # obliczamy transmitancję
    w, h = freqz(...) 
    
    # obliczamy moduł transmitancji
    m = np.abs(h)
        
    # obliczamy fazę i "rozwijamy" ją   
    faza = np.unwrap(np.angle(h))

    # obliczamy opóźnienie fazowe
    opoznienieFazowe = - faza/w
    
    # obliczamy opóźnienie grupowe
    df = np.diff(faza)
    idx, = np.where(np.abs(df-np.pi)<0.05) #to zabezpieczenie na błędy przy "rozwijaniu" fazy
    df[idx] = (df[idx+1]+df[idx-1])/2
    grupowe = - df/np.diff(w)
    
    # obliczamy odpowiedź impulsową
    x = np.zeros(len(t))
    x[len(t)//2] = 1 # impuls
    y = lfilter(b,a,x)
    
    # obliczamy odpowiedź schodkową
    s = np.zeros(len(t))
    s[len(t)//2:] = 1 # schodek
    ys = lfilter(b,a,s) # przepuszczamy impuls przez filtr i obserwujemy odpowiedź impulsową
    
    # rysujemy
    fig = py.figure()
    py.subplot(3,2,1)
    py.title('moduł transmitancji')
    py.plot(f,20*np.log10(m))
    py.ylabel('[dB]')
    py.grid('on')
    
    py.subplot(3,2,3)
    py.title('opóźnienie fazowe')
    py.plot(f, opoznienieFazowe)
    py.ylabel('próbki')
    py.grid('on')
    
    py.subplot(3,2,5)
    py.title('opóźnienie grupowe')
    py.plot(f[:-1],grupowe)
    py.ylabel('próbki')
    py.xlabel('Częstość [Hz]')
    py.grid('on')
    py.ylim([0, np.max(grupowe)+1])
    
    py.subplot(3,2,2)
    py.title('odpowiedź impulsowa')
    py.plot(t, x)
    py.plot(t, y)
    py.xlim([-T/2,T])
    py.grid('on')
    
    py.subplot(3,2,4)
    py.title('odpowiedź schodkowa')
    py.plot(t, s)
    py.plot(t, ys)
    py.xlim([-T/2,T])
    py.xlabel('Czas [s]')
    py.grid('on')
    
    fig.subplots_adjust(hspace=.5)
    py.show()

Testujemy:

b = firwin(54,0.5)# licznik
a = np.array([1.0]) # mianownik
Fs = 100
T = 1
f = np.arange(0.01,Fs/2,0.01)
charkterystyki(a,b,f,T,Fs)

Funkcje do projektowania filtrów

W bibliotece scipy.signal jest kilka funkcji do projektowania filtrów o zadanych parametrach. My skupimy się na dwóch zasadniczych grupach: * FIR - filtry o skończonej odpowiedzi impulsowej * klasyczne IIR - filtry o nieskończonej odpowiedzi impulsowej

FIR

Filtry typu FIR zwykle wymagają znacznie wyższych rzędów aby osiągnąć transmitancję o porządanej formie. Mają jednak dwie podstawowe zalety: * ich funkcja odpowiedzi jest skończona opisana wektorem b - efekty brzegowe sięgają z obu końców filtrowanego sygnału na dokładnie połowę długości wektora b * mają liniową zależnaość fazy od częstości. Z tego powodu opóżnienie grupowe dla wszystkich częstości jest takie samo.

W module scipy.signal mamy kilka funkcji do projektowania filtrów typu FIR: firwin i firwin2.

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) http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html#scipy.signal.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łady: Zbadaj włsności przykładowych projektów

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, nyq= 128)
Fs = 128
T = 0.2
f = np.arange(0.01,Fs/2,0.01)
b = firwin(21, 40, nyq= 128)
charkterystyki(1,b,f,T ,Fs)
  • filtr górnoprzepustowy rzędu 15 z częstością odcięcia 5 Hz:
  • firwin(17, 15, nyq= 128, pass_zero=False)
Fs = 256
T = 1
f = np.arange(0.01,Fs/2,0.01)

b = firwin(17, 15, nyq= 128, pass_zero=False)
charkterystyki(1,b,f,T = 0.2,Fs=Fs)
  • pasmowo przepustowy 51 rzędu przenoszący częstości pomiędzy 8 a 14 Hz:
  • firwin(51, [8, 14], nyq= 128, pass_zero=False)
Fs = 256
f = np.arange(0.01,Fs/2,0.01) 

b=firwin(51, [8, 14], nyq= Fs/2, pass_zero=False)
charkterystyki(1,b,f,T = 0.2,Fs=Fs)

Zadanie: Zaprojektuj i zbadaj własności filtra:

FIR dolno z pasmem przenoszenia od 30 Hz dla sygnału próbkowanego 256 Hz

Zadanie: Znajdź rząd filtra FIR:

dolnoprzepustowego z pasmem przenoszenia do 40 Hz dla sygnału próbkowanego 256 Hz, tak aby dla częstości powyżej 45 Hz jego tłumienie było nie mniejsze niż 20dB.

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 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. 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. * gain Pożądane wartości przenoszenia odpowiadające kolejnym częstościom definiowane są w gain. * Znaczenie pozostałych parametrów jest takie samo jak dla ``firwin.

Zadanie: filtr wielopasmowy

Zaprojektuj filtr przenoszący częstości w pasmach pomiędzy : 10-11, 20-21 i 30-31 Hz, który w paśmie zaporowym ma co najmniej 60 dB tłumienia.

Fs =100
T = 2
f = np.arange(0.01,Fs/2,0.01) 
freq = np.array([0, 10, 10, 11, 11, 20, 20, 21, 21, 30, 30, 31, 31, 50])
gain = np.array([0,  0,  1,  1,  0,  0,  1,  1,  0,  0,  1,  1,  0,  0])
b = firwin2(100, freq, gain, nyq= 50)
charkterystyki(1,b,f,T,Fs)

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, * druga oblicza współczynniki filtru.

W poniższych przykładach przyjmiemy:

T = 0.3
Fs = 100 # Hz
f = np.arange(0.01,Fs/2,0.01) 

wp = 10/(Fs/2)
ws = 30/(Fs/2)
gpass = 1 
gstop = 25
analog=0
rp = gpass
rs = gstop

Filtr Butterwortha

daje gładką (bez tętnień) funkcję przenoszenia w całym zakresie częstości.

[n,Wn]=buttord(wp, ws, gpass, gstop)
[b,a]=butter(n,Wn)
charkterystyki(a,b,f,T,Fs)
print('Filtr Butterwortha, rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))

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)
[b,a]=cheby1(n, rp, Wn, btype='low',  output='ba')
charkterystyki(a,b,f,T,Fs)
print('Czebyszewa I Typu: rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))

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')
charkterystyki(a,b,f,T,Fs)
print('Czebyszewa II Typu: rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))

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')
charkterystyki(a,b,f,T,Fs)
print('Eliptyczny: rząd: {}, częstość odcięcia {:.3f}'.format(n,Wn*Fs/2))

Filtrowanie z zerowym przesunięciem fazowym

Zadanie

Filtrowanie sygnałów off-line można zrealizować tak, aby sygnał wyjściowy nie miał przesunięcia fazowego. 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:

# 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(...)
 
# obliczamy funkcję przenoszenia
w,h = freqz(...)
transmitancja = ...
 
#opóźnienie grupowe
grupowe = -np.diff(np.unwrap(np.angle(h)))/np.diff(w)  
 
# przeliczamy skalę w (radiany) na częstości w Hz 
f = ...
 
# 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(...)
 
# 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('moduł transmitancji')
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()

Zadanie

Skonstruować filtry dolnoprzepustowe rzędu n=5, 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.

Zadanie

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.

Zadanie: filtr pasmowy do wyszukiwania wrzecion snu

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. Działanie filtra przetestować na sygnale:

http://www.fuw.edu.pl/~jarekz/SYGNALY/TF/c4spin.txt

Sygnał ten to fragment zapisu EEG z II stadium snu elektroda C4 próbkownie 128Hz.

s = np.loadtxt('c4spin.txt') # wczytujemy sygnał z pliku tekstowego
Fs = 128
t = np.arange(0,len(s))/Fs # przygotowujemy oś czasu

py.plot(t, s)
[b,a] = butter(...)
sf = filtfilt(...)
py.plot(t,sf)
py.show()

Zadanie 5: uwaga na odpowiedź impulsową

Przydadzą nam się pliki: * https://drive.google.com/open?id=0BzwQ_Lscn8yDNGc0aU5jSDFFMmc Plik z sygnałem EKG * https://drive.google.com/open?id=0BzwQ_Lscn8yDOF9jX0pjcG9LSGc Plik z metadanymi do tego sygnału

Proszę: * wykreślić pierwsze 10 sekund sygnału * zastosować filtr górnoprzepustowy Butterwartha o częstościach odcięcia: 0.01, 0.1, 0.5 -> zaobserwuj jak długo stabilizuje się sygnał * Zastosuj filtr pasmowoprzepustowy (11 - 14 Hz) i wykreśl wynik jego zastosowania na tle poprzedniej wersji sygnału. Zinterpretuj wynik w kontekście funkcji odpowiedzi impulsowej tego filtra.

s = np.fromfile('EKG.bin', dtype='<f') # tworzymy tablicę sig o typie określonym przez ''dtype''                                
# ustawiamy częstość próbkowania
Fs = 128
# tworzymy wektor czasu
t = np.arange(0,len(s))/Fs

# ustalamy zakres indeksów sygnału i czasu do wyświetlania
zakres =  np.logical_and(0<t, t<10)

py.plot(t[zakres],s[zakres])
# filtr górnoprzepustowy (0.01, 0.1, 0.5)
[b,a] = butter(... )
sf = filtfilt(...)
py.plot(t[zakres],sf[zakres])

# filtr pasmowy 
[bl,al] = butter(... )
sf_l = filtfilt(bl,al,sf)
py.plot(t[zakres],sf_l[zakres])

py.show()

Zadanie: Uwaga na odpowiedź schodkową

Wykorzystajmy fragment sygnału EKG z poprzedniego zadania (pomiędzy 12 a 40 -tą sekundą). * wykreśl ten fragment * zaprojektuj filtr górnoprzepustowy Butterwortha o częstości odcięcia 0.1 (potem 0.5), rzedu 1 (potem 5), * przefiltruj sygnał z tymi współczynnikami za pomocą funkcji filtfilt i lfilter , * dodaj do sygnału z lewej strony jego kopię odwróconą w czasie, * ten sygnał przefiltruj funkcją lfilter i wykreśl jego drugą połowę, * zinterpretuj uzyskane wyniki w kontekście funkcji odpowiedzi impulsowej.

zakres =  np.logical_and(10<t, t<40)
t_z = t[zakres]
s_z = s[zakres]
py.plot(t_z,s_z, label = 'oryginalny')
py.grid('on')

# filtr górnoprzepustowy ( 0.1, 0.5)
[b,a] = butter(...)
sff = filtfilt(...)
py.plot(t_z, sff, label = 'filtfilt')

# lfilter
sfl = lfilter(...)
py.plot(t_z,...,label = 'lfilter')

# lfilter z przedłużeniem 
x=np.concatenate((s_z[::-1],s_z))
sfl_p = lfilter(...)
py.plot(t_z,sfl_p[len(t_z):],label = 'lfilter z przedłużaniem')
py.legend()
py.show()

Przepróbkowywanie

Przepróbkowywanie do góry:

Zwiększamy częstość prókowania całkowitą ilość razy P

Najpowszechniej stosowana metoda polega na dodaniu P zer pomiędzy istniejące próbki sygnału tak aby osiągnął on P-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:

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.title('Sygnał oryginalny')

py.subplot(3,1,2)
X = np.zeros(4*len(x))
X[::4] = x
py.plot(X,'.')
py.title('Sygnał poprzetykany zerami')

[b,a] = butter(8,...) # filtr powinien przepuszczać tylko te częstości, 
                        # które były w oryginalnym sygnale tzn. poniżej pierwotnego Nyqista
y = filtfilt(b,a,X); # filtrujemy dolnoprzepustowo  
py.subplot(3,1,3)
py.plot(y,'.')
py.show()

Przepróbkowanie 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.

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] liczbę razy — uzyskujemy składając powyższe kroki tzn. najpierw zwiększamy częstość P-krotnie, a następnie zmniejszamy Q-krotnie.

Zadanie 6

Proszę napisać funkcję, która będzie przepróbkowywać sygnał o wymierną liczbę razy. Funkcja powinna przyjmować sygnał, częstość próbkowania, parametry P i Q i zwracać przepróbkowany sygnał i nową częstość próbkowania

def resample(s,fs,P=1,Q=1):
    if P>1 and isinstance(P,int):
        sP = np.zeros(P*len(s))
        sP[...] = s
        fs = ...
        [b,a] = butter(...)
        s = filtfilt(...)
    if Q>1 and isinstance(Q,int):
        fs = fs/Q 
        [b,a] = butter(...) 
        s = filtfilt(...)
        s = s[...]
    return s,fs

fs = 1000
t = np.arange(0,0.05,0.001) # czas
s1 = np.sin(2*np.pi*30*t) + np.sin(2*np.pi*60*t) # sygnał
 
py.subplot(3,1,1)
py.plot(s1,'.')
y,fs2 = resample(s1,fs,P=10,Q=1)
py.subplot(3,1,2)
py.plot(y,'.')
py.show()

fs = 128.0
t = np.arange(0,1,1.0/fs)
f1 = 6 # Hz
f2 = 40
fi = np.pi/2
s2 = np.sin(2*np.pi*t*f1+fi) + np.sin(2*np.pi*t*f2+fi)
py.subplot(3,1,1)
py.plot(s2,'.-')
y,fs2 = resample(s2,fs,P=1,Q=2)
py.subplot(3,1,2)
py.plot(y,'.-')
py.show()