This is an old revision of the document!

Laboratorul 11.

Filtre FIR trece bandă și trece-sus, filtre IIR

În acest laborator ne vom înbunătății cunoștințele despre filtrele FIR, experimentând cu filtre trece bandă și trece-sus și vom învăta despre implementarea filtrelor IIR( filtre cu răspuns infinit la impuls), care pot mult mai eficiente decât filtrele FIR (ca și număr de operații de înmulțire pentru performanța similară).

Mareriale utile:

Vedeți slide-urile din R. Lyons aici

Exercițiul 1 -- Filtre FIR trece bandă și trece-sus [4p]

Similar cu al doilea exercițiu din laboratorul precedent, vom utiliza metoda proiectării cu fereastră pentru a crea filtre FIR trece-banda și trece-sus. Precum am văzut la curs, putem folosi același principiu pentru a crea filtre trece-jos, trece-banda sau trece-sus. Tot ce trebuie sa facem este să înmulțim coeficienții filtrelor (adică secvența h(k)) cu volorile unei sinusoide cu anume frecvență ( centrul frecvențelor pentru filtrul trece-bandă).

Pentru a crea un filtru trece-bandă cu frecvența centrală $f_B$ ar trebui procedați în felul următor:

Generați filtrul în timp h(k) precum în laboratorul precedent (creați filtrul în frecvență, treceți în domeniul timp, înmulțiti cu o fereastră precum Blackman sau alta)
Înmulțiti h(k) element cu element cu secvența $\cos(\frac{2\pi f_B n}{f_s})$, unde $f_s$ este frecvența de eșantionare.

Câteva cazuri particulare:

$f_B = \frac{f_s}{4}$, în acest caz secvența cosinus devine [1, -0, -1, 0, 1, 0, -1, 0, …]. Acesta este un tip de filtru eficient trece-bandă centrat în frecvența f_s / 4 $. -$ f_B = \frac{f_s}{2} $, în acest caz secvența cosinus devine [1, -1, 1, -1, ...] și obținem un filtru trece-sus. Acum că știți toate acestea (sperăm că ați reținut și de la curs), aveți de făcut următoaree: - Generați filtrul low-pass hk(n) cu L = 65 eșantioane precum în laboratorul precedent. - Din această secvență obțineți a secvență corespunzătoare unui filtru trece-banda cu$ f_B = \frac{f_s}{4} $. - Obțineți o secvență corespunzătoare unui filtru trece-sus cu$ f_B = \frac{f_s}{2} $. - Generați trei sinusoide cu frecvențe diferite (ex:$ 0.1 f_s$, $0.25 f_s$, $0.4 f_s $și filtrați-le (folosind funcția //conv//) cu filtrele trece-sus și trece-bandă obținute mai sus. Plotați atât inputul și outputul in acelasi plot folosind stem pentru a observa efectele filtrelor. /*<hidden>*/ <note> A possible solution is this: <code matlab bandpass_highpass.m> clear; close all; %% Generate an ideal-filter sequence (in freq domain) N = 256; kc = N/16; HK = zeros(1, N); HK(1:kc) = 1; HK(end-kc+2:end) = 1; h = figure; fx = linspace(0,1,N); plot(fx, HK); xlabel('Frequency (x Fs)'); title('Ideal low-pass filter'); print(h, '-dpng', 'ideal_filter.png'); %% Obtain time-domain sequence hk = ifftshift(ifft(HK)); figure; plot(1:N, hk); xlabel('Index'); title('Time-domain sequence of ideal low-pass filter'); %% Generate windowed time-domain sequence % window = rect L = 65; w = ones(1, L); hw = hk(N/2+1-floor(L/2):N/2+1+floor(L/2)) .* w; h = figure; plot(1:L, hw); xlabel('Sample index'); title('Windowed filter sequence h(k) using rect(n)'); %% Obtain frequency-domain of windowed sequence HW = fft(hw); figure; fx = linspace(0,1,L); plot(fx, abs(HW)); xlabel('Frequency (x Fs)'); title('FFT of the time-domain sequence with w(n)=rect(n)'); %% Generate windowed time-domain sequence % window = blackman L = 65; w = blackman(L)'; hw = hk(N/2+1-floor(L/2):N/2+1+floor(L/2)) .* w; h = figure; plot(1:L, hw); xlabel('Sample index'); title('Windowed filter sequence h(k) using blackman(n)'); %% Obtain frequency-domain of windowed sequence HW = fft(hw); figure; fx = linspace(0,1,L); plot(fx, abs(HW)); xlabel('Frequency (x Fs)'); title('FFT of the time-domain sequence with w(n)=blackman(n)'); %% Filter a sinewave using the low-pass filter fs = 64000; N = 64; x = 0:(N-1); f = 3000; A = 1; xn = A*sin(2*pi*f/fs*x); yn = conv(xn, hw); h = figure; stem(x, xn, 'b'); hold on; stem(x, yn(1:N), 'r'); xlabel('Sample index'); legend('xn', 'yn'); title('low-pass filtering using hw'); %% Obtain band-pass from low-pas filter % Multiply with sequence cos(pi/2 * n) = % [0, -1, 0, 1, 0, -1, ...] hwb = hw; hwb(1:2:end) = 0; hwb(2:4:end) = -hwb(2:4:end); %% Obtain frequency-domain of band-pass sequence HWB = fft(hwb); figure; fx = linspace(0,1,L); plot(fx, abs(HWB)); xlabel('Frequency (x Fs)'); title('FFT of the time-domain band-pass sequence'); %% Filter sinewaves using the band-pass filter fs = 64000; N = 64; x = 0:(N-1); f1 = 3000; f2 = 15000; f3 = 30000; A = 1; xn1 = A*sin(2*pi*f1/fs*x); xn2 = A*sin(2*pi*f2/fs*x); xn3 = A*sin(2*pi*f3/fs*x); yn1 = conv(xn1, hwb); yn2 = conv(xn2, hwb); yn3 = conv(xn3, hwb); h = figure; subplot(2,1,1); stem(x, xn1, 'r'); hold on; stem(x, xn2, 'g'); stem(x, xn3, 'b'); xlabel('Sample index'); legend('xn 3kHz', 'xn 15kHz', 'xn 30kHz'); title('Input sinusoids of different frequencies'); subplot(2,1,2); stem(x, yn1(1:N), 'r'); hold on; stem(x, yn2(1:N), 'g'); stem(x, yn3(1:N), 'b'); xlabel('Sample index'); legend('yn 3kHz', 'yn 15kHz', 'yn 30kHz'); title('Applying band-pass filtering'); %% Obtain high-pass from low-pas filter % Multiply with sequence cos(pi * n) = % [1, -1, 1, -1, ...] hwh = hw; hwh(2:2:end) = -hwh(2:2:end); %% Obtain frequency-domain of high-pass sequence HWH = fft(hwh); figure; fx = linspace(0,1,L); plot(fx, abs(HWH)); xlabel('Frequency (x Fs)'); title('FFT of the time-domain high-pass sequence'); %% Filter sinewaves using the high-pass filter fs = 64000; N = 64; x = 0:(N-1); f1 = 3000; f2 = 15000; f3 = 30000; A = 1; xn1 = A*sin(2*pi*f1/fs*x); xn2 = A*sin(2*pi*f2/fs*x); xn3 = A*sin(2*pi*f3/fs*x); yn1 = conv(xn1, hwh); yn2 = conv(xn2, hwh); yn3 = conv(xn3, hwh); h = figure; subplot(2,1,1); stem(x, xn1, 'r'); hold on; stem(x, xn2, 'g'); stem(x, xn3, 'b'); xlabel('Sample index'); legend('xn 3kHz', 'xn 15kHz', 'xn 30kHz'); title('Input sinusoids of different frequencies'); subplot(2,1,2); stem(x, yn1(1:N), 'r'); hold on; stem(x, yn2(1:N), 'g'); stem(x, yn3(1:N), 'b'); xlabel('Sample index'); legend('yn 3kHz', 'yn 15kHz', 'yn 30kHz'); title('Applying high-pass filtering'); </code> </note> /*</hidden>*/ === Exercițiul 2 -- Proiectarea de FIR rapidă folosind MATLAB[2p] === Precum am văzut la curs, putem descrie un filtru (sistem liniar) cu feedback(IIR având termenii$ a_i $mai jos ) sau fără feedback (FIR) folosind o ecuație cu diferențe precum:$ y(n) = b_0 x(n) + b_1 x(n-1) + \ldots + b_q x(n-q) + a_1 y(n-1) + \ldots + a_p y(n-p) $Putem reprezenta întârzierile x(n-1) ca$ z^{-1} x(n)$, unde $z=e^{j2\pi} $. Apoi, obținem o ecuație care depinde doar de x(n) și y(n) și obținem funcția de transfer a filtrului$ H(z) = \frac{y(n)}{x(n)} $precum:$ H(z) = \frac{\sum_{k=0}^q b_q z^{-q}}{1 - \sum_{k=0}^p a_p z^{-p}} $În MATLAB, puteți folosi funcția //fir1// pentru a obține rapid elementele$ b_i$ ale unui filtru FIR trece-jos, trece-bandă sau trece-sus(ignorați ceoficienții $a_i$ deocamdată). Apoi puteți folosi funcția //filter// pentru a filtra orice secvență folosind coeficienții $a_i $dați de //fir1//. Pentru acest exercițiu, folosiți funcția //fir1// pentru a proiecta filtre FIR trece-jos, trece-bandă și trece-sus. Apoi folosiți funcția //filter// pentru a testa filtrele cu aceleași secvențe ca în exercițiul precedent. Puteți verifica designul filtrelor folosind tool-ul //fvtool//. /*<hidden>*/ <note> Possible solution: <code matlab fir_filter.m> clear; close all; %% Generate sinewaves fs = 64000; N = 64; x = 0:(N-1); f1 = 3000; f2 = 15000; f3 = 30000; A = 1; xn1 = A*sin(2*pi*f1/fs*x); xn2 = A*sin(2*pi*f2/fs*x); xn3 = A*sin(2*pi*f3/fs*x); %% Design band-pass filter with fir L = 64; wn = [0.2, 0.4]; b = fir1(L, wn); figure; fvtool(b,1); %% Filter sequence with band-pass yn1 = filter(b, 1, xn1); yn2 = filter(b, 1, xn2); yn3 = filter(b, 1, xn3); h = figure; subplot(2,1,1); stem(x, xn1, 'r'); hold on; stem(x, xn2, 'g'); stem(x, xn3, 'b'); xlabel('Sample index'); legend('xn 3kHz', 'xn 15kHz', 'xn 30kHz'); title('Input sinusoids of different frequencies'); subplot(2,1,2); stem(x, yn1(1:N), 'r'); hold on; stem(x, yn2(1:N), 'g'); stem(x, yn3(1:N), 'b'); xlabel('Sample index'); legend('yn 3kHz', 'yn 15kHz', 'yn 30kHz'); title('Applying band-pass filtering'); %% Design high-pass filter with fir L = 64; wn = 0.75; b = fir1(L, wn, 'high'); figure; fvtool(b,1); %% Filter sequence with band-pass yn1 = filter(b, 1, xn1); yn2 = filter(b, 1, xn2); yn3 = filter(b, 1, xn3); h = figure; subplot(2,1,1); stem(x, xn1, 'r'); hold on; stem(x, xn2, 'g'); stem(x, xn3, 'b'); xlabel('Sample index'); legend('xn 3kHz', 'xn 15kHz', 'xn 30kHz'); title('Input sinusoids of different frequencies'); subplot(2,1,2); stem(x, yn1(1:N), 'r'); hold on; stem(x, yn2(1:N), 'g'); stem(x, yn3(1:N), 'b'); xlabel('Sample index'); legend('yn 3kHz', 'yn 15kHz', 'yn 30kHz'); title('Applying high-pass filtering'); </code> </note> /*</hidden>*/ === Exercițiul 3 -- Proiectarea de IIR rapidă folosind MATLAB [2p] === Acum să proiectăm un filtru IIR folosind functia MATLAB //butter//. Aceasta funcție ne va returna atât coeficienții$ a_i$ cât și $b_i$.

Precum în exerecițiul precedent, procedați în felul următor:

Proiectați un filtru trece-jos folosind butter
Verificați filtrul folosind fvtool
Folosiți filtrul pentru a filtra secvențele de la exercițiul anterior cu funcția filter
Comparați performanța filtrelor IIR și FIR utilizând:
1. număr similar de coeficienți
2. performanță similară (scăderea în dB a frecvențelor în afara benzii de trecere)

Exercițiul 4 --- Proiectarea de filtre FIR si IIR folosind //fdtool// din MATLAB [2p]

Folosiți tool-ul MATLAB fdatool pentru a proiecta și analiza rapid performanțele filtrelor FIR și IIR:

Generați un filtru FIR trece-jos folosind o fereastră Kaiser, având Fs = 48000 Hz, Fc = 10000 Hz, cu N=10 coeficienți(ordin)
Generați un filtru IIR trece-jos (Butterworth sau Chebyshev type I) cu aceeași parametri.
Comparați amplitudinea răspunsului
Încercați sa creați un filtru FIR cu un răspuns similar cu al unui filtru IIR, prin cresterea numarului de coeficienti.

După proiectarea filtrelor precum ați dorit, le puteți salva (File→Generate MATLAB Code) și să le folosiți direct în alte scripturi MATLAB pentru a filtra diverse semnale.