This is an old revision of the document!
Laboratorul 11.
Filtre FIR trece bandă și trece-sus, filtre IIR
În acest laborator ne vom îmbunătății cunoștințele despre filtrele FIR, experimentând cu filtre trece bandă și trece-sus și vom învăța despre implementarea filtrelor IIR( filtre cu răspuns infinit la impuls), care pot fi mult mai eficiente decât filtrele FIR (ca și număr de operații de înmulțire pentru performanță 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-bandă și trece-sus. Precum am văzut la curs, putem folosi același principiu pentru a crea filtre trece-jos, trece-bandă sau trece-sus. Tot ce trebuie sa facem este să înmulțim coeficienții filtrelor (adică secvența h(k)) cu valorile unei sinusoide de o anumită frecvență (centrul frecvențelor pentru filtrul trece-bandă).
Pentru a crea un filtru trece-bandă cu frecvența centrală $f_B$ ar trebui să 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țiți 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ătoarele:
- Generați filtrul trece-jos hk(n) cu L = 65 eșantioane precum în laboratorul precedent.
- Din această secvență obțineți o secvență corespunzătoare unui filtru trece-bandă 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 input-ul și output-ul în acelasi plot folosind stem pentru a observa efectele filtrelor.
- 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');
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.
- 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');
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:
- număr similar de coeficienți
- performanță similară (scăderea în dB a frecvențelor în afara benzii de trecere)
Exercițiul 4 --- Proiectarea de filtre FIR si IIR folosind fdatool 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.
