This shows you the differences between two versions of the page.
ps:labs:09 [2021/11/15 22:10] darius.necula |
ps:labs:09 [2022/12/06 15:19] (current) ionut.gorgos |
||
---|---|---|---|
Line 1: | Line 1: | ||
===== Laboratorul 09. ===== | ===== Laboratorul 09. ===== | ||
- | <hidden> | + | /*<hidden>*/ |
==== Convoluția, filtre FIR și metoda de proiectare folosind ferestre ==== | ==== Convoluția, filtre FIR și metoda de proiectare folosind ferestre ==== | ||
+ | Prezentarea PowerPoint pentru acest laborator poate fi găsită aici: [[https://docs.google.com/presentation/d/16NWgSnRxFygNBcrHc7IvEfMR_MzZxyu7/edit?usp=sharing&ouid=110538702824281541719&rtpof=true&sd=true|aici]] | ||
În acest laborator vom face câteva exerciții pentru a ne familiariza cu operația de convoluție precum și cu filtre cu răspunsul finit la impuls (finite impulse response - FIR) și cu metode de proiectare a acestora prin metode folosind ferestre. | În acest laborator vom face câteva exerciții pentru a ne familiariza cu operația de convoluție precum și cu filtre cu răspunsul finit la impuls (finite impulse response - FIR) și cu metode de proiectare a acestora prin metode folosind ferestre. | ||
Line 9: | Line 10: | ||
- | === Exercițiul 1 -- Convoluția [4p] === | + | === Exercițiul 1 -- Convoluția === |
+ | [<color red>4p</color>] | ||
Precum am spus la curs, putem defini operația de convoluție dintre două secvențe h(k) și x(n) după cum urmează: | Precum am spus la curs, putem defini operația de convoluție dintre două secvențe h(k) și x(n) după cum urmează: | ||
Line 18: | Line 20: | ||
În general presupunem că secvența h(k) are un număr $M$ finit și, în general, mic de elemente. | În general presupunem că secvența h(k) are un număr $M$ finit și, în general, mic de elemente. | ||
- | <hidden> In our context, these elements will usually be the //taps// of a filter, and the convolution operation will be used to filter an input sequence x(n), as we shall see in the following exercises. /*</hidden>*/ | + | <hidden> In our context, these elements will usually be the //taps// of a filter, and the convolution operation will be used to filter an input sequence x(n), as we shall see in the following exercises. </hidden> |
Rezolvați următoarele exerciții: | Rezolvați următoarele exerciții: | ||
- Dacă x(n) are N elemente și h(k) are M elemente, câte elemente are secvența obținută prin convoluție $h(k) \ast x(n)$? (presupunând că efectuăm convoluția doar pe elementele diferite de zero, adică atunci când elementele celei mai scurte secvențe se suprapun în totalitate cu elementele celeilalte) | - Dacă x(n) are N elemente și h(k) are M elemente, câte elemente are secvența obținută prin convoluție $h(k) \ast x(n)$? (presupunând că efectuăm convoluția doar pe elementele diferite de zero, adică atunci când elementele celei mai scurte secvențe se suprapun în totalitate cu elementele celeilalte) | ||
- Fie x(n) secvența [1, 3, 5, 7, 5, 4, 2] și h(k) secvența [0.1, 0.3, 0.1]. Secvența obținută prin convoluție $y(n) = h(k) \ast x(n), n \in \{1,2,3,4,5\}$ ar trebui sa aibă 5 elemente. Scrieți fiecare dintre aceste elemente ca un produs scalar. Observați că trebuie să inversați ordinea secvenței x(n) (doar partea care se înmulțeste cu filtrul) înainte să o înmulțiți cu h(k) (alternativ, puteți inversa o singură dată elementele filtrului). | - Fie x(n) secvența [1, 3, 5, 7, 5, 4, 2] și h(k) secvența [0.1, 0.3, 0.1]. Secvența obținută prin convoluție $y(n) = h(k) \ast x(n), n \in \{1,2,3,4,5\}$ ar trebui sa aibă 5 elemente. Scrieți fiecare dintre aceste elemente ca un produs scalar. Observați că trebuie să inversați ordinea secvenței x(n) (doar partea care se înmulțeste cu filtrul) înainte să o înmulțiți cu h(k) (alternativ, puteți inversa o singură dată elementele filtrului). | ||
- | - Fie x(n) o secvența de $N=64$ elemente corespunzătoare unei sinusoide de frecvență $f=3$ kHz, eșantionată cu $f_s=64$ kHz. Fie h(k) secvența [0.1, 0.2, 0.2, 0.2, 0.1]. Generați aceste secvențe în Matlab/Octave și implementați convoluția pentru a obține elementele $y(n) = h(k) \ast x(n)$. Plotați inputul x(n) și ieșirea y(n); folosiți funcția //stem// în locul funcției //plot// pentru acest exercițiu. | + | - Fie x(n) o secvență de $N=64$ elemente corespunzătoare unei sinusoide de frecvență $f=3$ kHz, eșantionată cu $f_s=64$ kHz. Fie h(k) secvența [0.1, 0.2, 0.2, 0.2, 0.1]. Generați aceste secvențe în MATLAB și implementați convoluția pentru a obține elementele $y(n) = h(k) \ast x(n)$. Plotați inputul x(n) și ieșirea y(n); folosiți funcția //stem// în locul funcției //plot// pentru acest exercițiu. |
- Înlocuiți x(n) cu secvența de impuls cu 9 elemente [0, 0, 0, 0, 1, 0, 0, 0, 0] și efectuați convoluția cu aceeași secvență h(k) ca mai sus. Ce obțineți ca y(n)? Cum se numește aceasta? | - Înlocuiți x(n) cu secvența de impuls cu 9 elemente [0, 0, 0, 0, 1, 0, 0, 0, 0] și efectuați convoluția cu aceeași secvență h(k) ca mai sus. Ce obțineți ca y(n)? Cum se numește aceasta? | ||
- | - Încercați operațiile de mai sus folosind funcția Matlab/Octave //conv//. Obțineți aceleași rezultate? Care sunt diferențele? | + | - Încercați operațiile de mai sus folosind funcția MATLAB //conv//. Obțineți aceleași rezultate? Care sunt diferențele? |
<note tip> | <note tip> | ||
Line 31: | Line 33: | ||
</note> | </note> | ||
- | === Exercițiul 2 -- Filtre FIR [4p] === | + | === Exercițiul 2 -- Filtre FIR === |
+ | [<color red>4p</color>] | ||
- | În acest exercițiu vom creea secvența x(n) pentru un filtru trece-jos plecând de la un filtru ideal în domeniul frecvență și, folosind IFFT, vom obține secvența x(n). Vom folosi de asemenea diferite ferestre pentru a le compara performanța în proiectarea de filtre trece-jos. | + | În acest exercițiu vom crea secvența x(n) pentru un filtru trece-jos plecând de la un filtru ideal în domeniul frecvență și, folosind IFFT, vom obține secvența x(n). Vom folosi de asemenea diferite ferestre pentru a le compara performanța în proiectarea de filtre trece-jos. |
- | Urmați următorii pași: | + | Încercați să urmați acești pași: |
- | - Generați o secvență de filtru ideal trece-jos HK având N = 256 elemente, reprezentând spectrul de frecvență al unui filtru trece-jos. Folosiți o frecvență de cut-off de fs/16. Adică totul înainte de fs/16 trebuie sa treacă, pe când totul mai sus trebuie sa fie oprit (folosiți un dreptunghi care se opreste la fs/16). Observați că trebuie să generați un spectru simetric pentru a obține o secvență reală la următorul pas. Plotați această secvență (folosind //plot//). Ar trebui să obțineți ceva precum: {{:ps:labs:9.png?300|}} | + | - Generați o secvență de filtru ideal trece-jos HK având //N = 256// elemente, reprezentând spectrul de frecvență al unui filtru trece-jos. Folosiți o frecvență de cut-off de $f_s/16$. Adică totul înainte de $f_s/16$ trebuie sa treacă, pe când totul după $f_s/16$ trebuie sa fie oprit (folosiți un dreptunghi care se oprește la $f_s/16$). Observați că trebuie să generați un spectru simetric pentru a obține o secvență reală la următorul pas. Plotați această secvență (folosind //plot//). Notați axa frecvențelor (axa x) ca o funcție de $F_s$, adică de la //0// la //1//. Ar trebui să obțineți ceva precum: |
- | * Țineți minte că acest spectru poate fi văzut ca ieșirea din DFT(FFT), adică primul element corespunde frecvenței 0, pe când următoarele N/2-1 corespund frecvențelor pozitive, iar ultimele N/2 componente reprezintă frecvențele negative. | + | {{:ps:labs:9.png?300|}} |
+ | * Țineți minte că acest spectru poate fi văzut ca ieșirea din DFT(FFT), adică primul element corespunde frecvenței //0//, pe când următoarele //N/2-1// corespund frecvențelor pozitive, iar ultimele //N/2// componente reprezintă frecvențele negative. | ||
- | - Acum aplicați inversa DFT(în practică inversa FFT, //ifft// în Matlab/Octave) pentru a obține secvența corespunzătoare în domeniul timp hk(n). Rețineți: trebuie sa aplicați funcția //ifftshift// pe rezultat pentru a obține o funcție sinc simetrică, adică folosiți ceva precum: <code matlab>hk = ifftshift(ifft(HK));</code> | + | - Acum aplicați inversa DFT(în practică inversa FFT, //ifft// în MATLAB) pentru a obține secvența corespunzătoare în domeniul timp hk(n). Rețineți: trebuie să aplicați funcția //ifftshift// pe rezultat pentru a obține o funcție sinc simetrică, adică folosiți ceva precum: <code matlab>hk = ifftshift(ifft(HK));</code> |
- | - Trunchiați secvența hk(n) prin selectarea a doar L=65 de eșantioane din centru(32 din stânga maximului funcției sinc, maximul funcției, și 32 de eșantioane din dreapta). Aceasta corespunde multiplicării secvenței hk(n) cu o fereastră rectangulară centrată în punctul maxim al funcției sinc. Plotați secvența. | + | - Trunchiați secvența hk(n) prin selectarea a doar //L=65// de eșantioane din centru(32 din stânga maximului funcției sinc, maximul funcției, și 32 de eșantioane din dreapta). Aceasta corespunde multiplicării secvenței hk(n) cu o fereastră rectangulară centrată în punctul maxim al funcției sinc. Plotați secvența. |
- | - Aplicați DFT(//fft//) pe secvența trunchiată înmulțită cu fereastra rectangulară (care conține doar 1) și plotați spectrul (cu //plot//). Rețineți: este important aici, precum și la primul plot pentru filtru trece-jos ideal, să notăm axa frecvențelor (axa x) ca o funcție de Fs, adică de la 0 la 1. Vedeți diferențe față de filtrul ideal trece-jos? Acestea sunt efectele ferestrei dreptunghiulare. | + | - Aplicați DFT(//fft//) pe secvența trunchiată, adică cea înmulțită cu fereastra rectangulară (care conține doar 1) și plotați spectrul (cu //plot//). Rețineți: este important aici, precum și la primul plot pentru filtru trece-jos ideal, să notăm axa frecvențelor (axa x) ca o funcție de $F_s$, adică de la //0// la //1//. Vedeți diferențe față de filtrul ideal trece-jos? Acestea sunt efectele ferestrei dreptunghiulare. |
- | - Folosiți aceeași secvență trunchiată ca mai sus, dar înmulțiți-o cu o fereastră precum //Blackman// (//blackman// în MATLAB/OCTAVE). Efectuați din nou DFT și plotați spectrul (cu //plot//). Arată mai bine?. | + | - Folosiți aceeași secvență trunchiată (hk(n)), dar înmulțiți-o cu o fereastră precum //Blackman// (//blackman// în MATLAB). Efectuați din nou DFT și plotați spectrul (cu //plot//). Arată mai bine?. |
- | - În final, folosiți ca intrare sinusoida din Exercițiul 1 ca x(n) și filtrați-o printr-o convoluție cu secvența obținută mai sus după folosirea ferestrei Blackman(folosiți funcția //conv// din Matlab/Octave). Plotați intrarea și ieșirea în aceeași figură folosind //stem// pentru a observa efectele filtrului. | + | - În final, folosiți ca intrare sinusoida din Exercițiul 1 ca x(n) și filtrați-o printr-o convoluție cu secvența obținută mai sus după folosirea ferestrei Blackman(folosiți funcția //conv// din MATLAB). Plotați intrarea și ieșirea în aceeași figură folosind //stem// pentru a observa efectele filtrului. |
- | === Exercițiul 3 -- Proiectarea de FIR rapidă folosind MATLAB [2p] === | + | === Exercițiul 3 -- Proiectarea rapidă de filtre FIR folosind MATLAB === |
+ | [<color red>2p</color>] | ||
- | 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: | + | 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 \cdot x(n) + b_1 \cdot x(n-1) + \ldots + b_q \cdot x(n-q) + a_1 \cdot y(n-1) + \ldots + a_p \cdot y(n-p)$ |
- | Putem reprezenta întârzierile x(n-1) ca $z^{-1} x(n)$, unde $z=e^{j2\pi}$. | + | Putem reprezenta întârzierile $x(n-1)$ ca $z^{-1} \cdot 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: | + | 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(z)}{X(z)}$ precum: |
- | $H(z) = \frac{\sum_{k=0}^q b_q z^{-q}}{1 - \sum_{k=0}^p a_p z^{-p}}$ | + | $H(z) = \frac{\sum_{k=0}^q b_k \cdot z^{-k}}{1 - \sum_{k=1}^p a_k \cdot z^{-k}}$ |
- | Î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//. | + | Î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 coeficienții $a_i$ deocamdată). Apoi puteți folosi funcția //filter// pentru a filtra orice secvență folosind coeficienții $b_i$ dați de //fir1//, iar $a_0=1$. |
- | + | ||
- | 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 niște sinusoide ca în Exercițiul 1, cu $f=3$ kHz, $15$ kHz, $30$ kHz . Puteți verifica designul filtrelor folosind //freqz//. | ||
+ | <code matlab>freqz(b,1);</code> | ||
+ | <hidden> | ||
+ | <note warning>Pentru a putea folosi funcțiile //fir1// și //freqz// trebuie să instalați pachetul signal din Octave. | ||
+ | <code matlab> | ||
+ | pkg install -forge -nodeps signal | ||
+ | </code> | ||
+ | , iar apoi să-l încărcați: | ||
+ | <code matlab> | ||
+ | pkg load signal | ||
+ | </code> | ||
+ | </note> | ||
</hidden> | </hidden> | ||
+ | |||
+ | Pentru răspunsul la impuls puteți folosi: | ||
+ | <code matlab>stem(b);</code> | ||
+ | Ar trebui să vedeți sinc-ul digital. | ||
+ | |||
+ | Pentru răspunsul în frecvență, vă puteți folosi de următorul cod: | ||
+ | |||
+ | <code matlab> | ||
+ | [H, freq] = freqz(b,1,N, fs); | ||
+ | figure; | ||
+ | plot(freq, abs(H)); | ||
+ | xlabel('Frequency (Hz)'); | ||
+ | ylabel('Frequency response'); | ||
+ | </code> | ||
+ | |||
+ | /*</hidden>*/ | ||