Î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

[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 să 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 să 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ță 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 Python ș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 np.convolve din NumPy sau signal.convolve din scipy.signal. Obțineți aceleași rezultate? Care sunt diferențele?
6. [BONUS] Teorema Convoluției. Convoluția în timp [ $h(k) \ast x(n)$ ] este echivalentă cu înmulțirea element cu element a spectrelor în frecvență [ $H(m) \cdot X(m)$ ] și invers. Pentru a testa asta faceți următorii pași:

• Folosiți secvența h(k) = [0.1, 0.2, 0.2, 0.2, 0.1], dar adăugați încă 59 de zerouri la finalul ei, precum am făcut în laboratoarele precedente.
• Aplicați FFT peste această secvență și afișați spectrul ei (cu stem).
• Aplicați FFT peste sinusoida de la subpunctul 3 și afișați spectrul ei (cu stem).
• Înmulțiți element cu element cele două spectre și afișați rezultatul cu stem. S-a modificat ceva?
• Reveniți în domeniul timp folosind ifft și np.real.
• Comparați rezultatul de mai sus cu cel al funcției np.convolve din NumPy sau signal.convolve, setând parametrul mode='full' și păstrând doar primele N elemente. Obțineți aceleași rezultate? Care sunt diferențele?

[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 Python) 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ă dreptunghiulară 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 dreptunghiulară (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ă (hk(n)), dar înmulțiți-o cu o fereastră precum Blackman (np.blackman în Python). 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 np.convolve din NumPy sau signal.convolve din scipy.signal). Plotați intrarea și ieșirea în aceeași figură folosind stem pentru a observa efectele filtrului.

[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 Python, puteți folosi funcția signal.firwin din SciPy, pentru a obține coeficienții ($b_i$) ai unui filtru FIR trece-jos, trece-bandă sau trece-sus(ignorați coeficienții $a_i$ deocamdată). Apoi puteți folosi funcția signal.lfilter din SciPy, pentru a filtra orice secvență folosind coeficienții $b_i$ dați de signal.firwin, iar $a_0 = 1$.

Pentru acest exercițiu, trebuie să proiectați filtre FIR trece-jos (folosiți cutoff = 0.1), trece-bandă (folosiți cutoff = [0.2, 0.5]) și trece-sus (folosiți cutoff = 0.75). Puteți alege numtaps = 65. Apoi folosiți funcția signal.lfilter pentru a testa filtrele cu niște sinusoide ca în Exercițiul 1, cu $f=3$ kHz, $15$ kHz, $30$ kHz. Afișați cu stem, în subplot-uri sinusoidele inițiale și pe cele filtrate. Pentru răspunsul în frecvență, vă puteți folosi de următorul cod:

freq, H = signal.freqz(b_low, 1, fs=fs)
plt.figure()
plt.plot(freq, 20 * np.log10(abs(H)))
plt.xlabel('Frequency [Hz], from 0 to fs/2')
plt.ylabel('Amplitude [dB]')
plt.title('Digital Filter Frequency Response')
plt.grid()
plt.show()

, unde b_low sunt coeficienții $b_i$ ai unui filtru FIR trece-jos.

[BONUS]

În acest exercițiu, încercăm să urmăm pașii de la exercițiul 3, dar în MATLAB. Pentru a accesa MATLAB, puteți folosi MATLAB instalat pe PC-urile din laborator sau să folosiți MATLAB Online, din browser. Pentru MATLAB Online, trebuie să vă creați un cont cu adresa de e-mail de la facultate pe site-ul acesta.

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

De asemenea puteți folosi tool-ul MATLAB fdatool pentru a proiecta și analiza rapid performanțele filtrelor FIR și IIR:

1. 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)
2. Generați un filtru IIR trece-jos (Butterworth sau Chebyshev type I) cu aceiași parametri.
3. Comparați amplitudinea răspunsului.
4. Încercați să 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.