Ćwiczenia 7

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()