This shows you the differences between two versions of the page.
ps:laboratoare:11 [2017/12/10 21:29] ionut.gorgos |
ps:laboratoare:11 [2020/10/07 18:34] (current) ionut.gorgos |
||
---|---|---|---|
Line 1: | Line 1: | ||
===== Laboratorul 11. ===== | ===== Laboratorul 11. ===== | ||
+ | /*<hidden>*/ | ||
==== Filtre FIR trece bandă și trece-sus, filtre IIR ==== | ==== 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ă). | + | Î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: | + | Materiale utile: |
- | * Vedeți slide-urile din R. Lyons [[http://cs.curs.pub.ro/2015/pluginfile.php/27567/mod_resource/content/0/iir_lyons_part1.pdf|aici]] | + | * Vedeți slide-urile din R. Lyons {{:ps:laboratoare:iir_lyons_part1.pdf|aici}} |
=== Exercițiul 1 -- Filtre FIR trece bandă și trece-sus [4p] === | === 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 valorile unei sinusoide de o anumită frecvență (centrul frecvențelor pentru filtrul trece-bandă). | + | Similar cu al doilea exercițiu din laboratorul 9, 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 să 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: | 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) | + | - Generați filtrul în timp h(k) precum în laboratorul 9 (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. | - Î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: | 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}{4}$, în acest caz secvența cosinus devine [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. | - $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: | 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 low-pass hk(n) cu L = 65 eșantioane precum în laboratorul precedent. | + | - Generați filtrul trece-jos hk(n) cu L = 65 eșantioane precum în laboratorul 9. |
- | - Din această secvență obțineți a secvență corespunzătoare unui filtru trece-banda cu $f_B = \frac{f_s}{4}$. | + | - 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}$. | - 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. | + | - Generați trei sinusoide cu frecvențe diferite (ex: $0.1 f_s$, $0.25 f_s$, $0.4 f_s$, unde $f_s = 64000$ și $N = 64$) ș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 cât și output-ul în același plot folosind //stem// pentru a observa efectele filtrelor. |
- | <hidden> | ||
- | <note> | ||
- | O soluție posibilă: | ||
- | <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] === | === 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: | + | 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)$ | $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)$ | ||
Line 212: | Line 39: | ||
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//. | 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> | ||
- | O soluție posibilă:: | ||
- | <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] === | === 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$. | + | Acum să proiectăm un filtru IIR folosind funcția din 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: | Precum în exerecițiul precedent, procedați în felul următor: | ||
- Proiectați un filtru trece-jos folosind //butter// | - Proiectați un filtru trece-jos folosind //butter// | ||
Line 304: | Line 52: | ||
- | === Exercițiul 4 --- Proiectarea de filtre FIR si IIR folosind //fdtool// din MATLAB [2p] === | + | === 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: | Folosiți tool-ul MATLAB //fdatool// pentru a proiecta și analiza rapid performanțele filtrelor FIR și IIR: | ||
Line 310: | Line 58: | ||
- Generați un filtru IIR trece-jos (Butterworth sau Chebyshev type I) cu aceeași parametri. | - Generați un filtru IIR trece-jos (Butterworth sau Chebyshev type I) cu aceeași parametri. | ||
- Comparați amplitudinea răspunsului | - 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. | + | - Încercați sa creați un filtru FIR cu un răspuns similar cu al unui filtru IIR, prin creșterea numărului de coeficienți. |
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. | 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. | ||
+ | /*</hidden>*/ | ||
+ | |||
+ | * Responsabil: [[ionutgorgos@gmail.com|Ionuț Gorgos]] | ||
+ | |||
+ | * Data publicării: 9.12.2019 | ||