Laboratorul 09.

Convoluția, filtre FIR și metoda de proiectare 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.

Materiale utile:

  • Vedeți slide-urile de la R. Lyons aici

Exercițiul 1 -- Convoluția [4p]

Precum am spus la curs, putem defini operația de convoluție dintre două secvențe h(k) și x(n) după cum urmează:

$y(n) = h(k) \ast x(n) = \sum_{k=0}^{M-1} h(k) \cdot x(n-k)$,

unde $M$ este lungimea secvenței h(k). De reținut că această operație definește un singur element de ieșire y(n). Pentru următorul element de ieșire, y(n+1), trebuie sa shiftăm secvența h(k) astfel încât să se potrivească cu elementele [x(n+1), x(n), …, x(n-M+2)]. Vedeți slide-urile 5 și 6 aici.

În general presupunem că secvența h(k) are un număr $M$ finit și, în general, mic de elemente.

Rezolvați următoarele exerciții:

  1. 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)
  2. 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).
  3. 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.
  4. Î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?
  5. Încercați operațiile de mai sus folosind funcția Matlab/Octave conv. Obțineți aceleași rezultate? Care sunt diferențele?

Pentru a plota secvențe de lungimi diferite în același plot, doar afișați primele N elemente ale secvenței, renunțând la restul.

Exercițiul 2 -- Filtre FIR [4p]

Î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.

Încercați să urmați acești pași:

  1. 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.
  1. 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 să aplicați funcția ifftshift pe rezultat pentru a obține o funcție sinc simetrică, adică folosiți ceva precum:
    hk = ifftshift(ifft(HK));
  2. 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.
  3. 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.
  4. 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?.
  5. Î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.

Exercițiul 3 -- Proiectarea rapidă de filtre FIR folosind MATLAB/Octave [2p]

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 \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} \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(z)}{X(z)}$ precum: $H(z) = \frac{\sum_{k=0}^q b_k \cdot z^{-k}}{1 - \sum_{k=1}^p a_k \cdot z^{-k}}$

În MATLAB/Octave, 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 niște sinusoide ca în Exercițiul 1, cu $f=3$ kHz, $15$ kHz, $30$ kHz . Puteți verifica designul filtrelor folosind freqz.

freqz(b,1);

Pentru a putea folosi funcțiile fir1 și freqz trebuie să instalați pachetul signal din Octave.

pkg install -forge -nodeps signal

, iar apoi să-l încărcați:

pkg load signal

Pentru răspunsul la impuls puteți folosi:

stem(b);

Ar trebui să vedeți sinc-ul digital.

Pentru răspunsul în frecvență, vă puteți folosi de următorul cod:

[H, freq] = freqz(b,1,N, fs);
figure;
plot(freq, abs(H));
xlabel('Frequency (Hz)');
ylabel('Frequency response');
ps/labs/09.txt · Last modified: 2021/12/12 21:50 by ionut.gorgos
CC Attribution-Share Alike 3.0 Unported
www.chimeric.de Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0