Differences
This shows you the differences between two versions of the page.
|
ps:labs_python:08-10 [2026/09/25 21:35] marian_gabriel.dinu created |
ps:labs_python:08-10 [2026/09/25 21:46] (current) marian_gabriel.dinu |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| - | ===== Laboratorul 08. ===== | + | ===== Secțiunea 3 - RTL-SDR, Convoluție și Tipuri de filtre ===== |
| + | ==== Partea I - RTL-SDR ==== | ||
| /*<hidden>*/ | /*<hidden>*/ | ||
| Line 450: | Line 451: | ||
| </file> | </file> | ||
| + | |||
| + | /*</hidden>*/ | ||
| + | |||
| + | ===== 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 {{:ps:labs:fir_lyons.pdf|aici}} | ||
| + | |||
| + | |||
| + | === 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ă: | ||
| + | |||
| + | $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 {{:ps:labs:fir_lyons.pdf|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: | ||
| + | - 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 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). | ||
| + | - 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. | ||
| + | - Î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 //np.convolve// din NumPy sau //signal.convolve// din //scipy.signal//. Obțineți aceleași rezultate? Care sunt diferențele? | ||
| + | - [<color red>BONUS</color>] 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? | ||
| + | |||
| + | <note tip> | ||
| + | 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. | ||
| + | </note> | ||
| + | |||
| + | === Exercițiul 2 -- Filtre FIR === | ||
| + | [<color red>4p</color>] | ||
| + | |||
| + | Î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: | ||
| + | - 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: | ||
| + | {{: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 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: <code python>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ă dreptunghiulară centrată în punctul maxim al funcției sinc. Plotați secvența. | ||
| + | - 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. | ||
| + | - 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?. | ||
| + | - Î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. | ||
| + | |||
| + | === Exercițiul 3 -- Proiectarea rapidă de filtre FIR folosind Python === | ||
| + | [<color red>2p</color>] | ||
| + | |||
| + | 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: | ||
| + | |||
| + | <code python> | ||
| + | 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() | ||
| + | </code> | ||
| + | , unde b_low sunt coeficienții $b_i$ ai unui filtru FIR trece-jos. | ||
| + | |||
| + | === Exercițiul 4 -- Proiectarea rapidă de filtre FIR/IIR folosind MATLAB === | ||
| + | [<color red>BONUS</color>] | ||
| + | |||
| + | Î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 [[https://www.mathworks.com/mwaccount/register|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: | ||
| + | - 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 aceiași parametri. | ||
| + | - Comparați amplitudinea răspunsului. | ||
| + | - Î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. | ||
| + | |||
| + | ===== Laboratorul 10. ===== | ||
| + | |||
| + | ==== Filtre FIR trece bandă și trece-sus, filtre IIR ==== | ||
| + | /*<hidden>*/ | ||
| + | === Exercițiul 1 -- Filtre FIR trece bandă și trece-sus === | ||
| + | [<color red>8p</color>] | ||
| + | |||
| + | Pentru acest exercițiu 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_{ideal}$) 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_{ideal}$ precum în laboratorul 9, ex.2 (creați filtrul în frecvență, treceți în domeniul timp, înmulțiți cu o fereastră precum Blackman sau alta) | ||
| + | * Înmulțiți $h_{ideal}$ element cu element cu secvența $\cos(\frac{2\pi f_B n}{f_s})$, unde $f_B$ este frecvența din centrul benzii dorite, iar $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. | ||
| + | |||
| + | <note tip>Spectrul 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.</note> | ||
| + | |||
| + | 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 o secvență de filtru ideal trece-jos $H_{ideal}$ 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 să treacă, pe când totul mai sus trebuie sa fie oprit (folosiți un dreptunghi care se oprește 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//). 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: {{:ps:labs:9.png?300|}} | ||
| + | - Acum aplicați inversa DFT (în practică inversa FFT) pentru a obține secvența corespunzătoare în domeniul timp $h_{ideal}$. | ||
| + | - Trunchiați secvența $h_{ideal}$ 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 $h_{ideal}$ cu o fereastră dreptunghiulară centrată în punctul maxim al funcției sinc. Plotați secvența. | ||
| + | - Aplicați DFT($fft$) pe secvența trunchiată înmulțită cu fereastra dreptunghiulară (care conține doar 1) și plotați spectrul. 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. | ||
| + | - Folosiți aceeași secvență trunchiată ca mai sus, dar înmulțiți-o cu o fereastră precum $Blackman$. Aplicați din nou $fft$ și plotați spectrul. | ||
| + | - Din această secvență puteți obține o secvență corespunzătoare unui filtru trece-bandă cu $f_B = \frac{f_s}{4}$. Afișați spectrul. | ||
| + | - Obțineți o secvență corespunzătoare unui filtru trece-sus cu $f_B = \frac{f_s}{2}$. Afișați spectrul. | ||
| + | - Generați trei sinusoide cu frecvențe diferite (ex: $ f = 3 kHz$, $15 kHz$, $30 kHz$, cu $f_s = 64000$ și $N = 64$) și filtrați-le (folosind funcția //np.convolve// din NumPy sau //signal.convolve// din //scipy.signal//) 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. | ||
| + | |||
| + | === Exercițiul 2 -- Proiectarea rapidă de filtre IIR folosind Python === | ||
| + | [<color red>2p</color>] | ||
| + | |||
| + | 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.butter// din SciPy, pentru a obține coeficienții ($b_i$ și $a_i$) ai unui filtru IIR trece-jos, trece-bandă sau trece-sus. Apoi puteți folosi funcția //signal.lfilter// din SciPy, pentru a filtra orice secvență folosind coeficienții $b_i$ și $a_i$ dați de //signal.butter//. | ||
| + | |||
| + | Pentru acest exercițiu, trebuie să proiectați filtre IIR 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 = 5//. 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: | ||
| + | |||
| + | <code python> | ||
| + | freq, H = signal.freqz(b_low, a_low, 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() | ||
| + | </code> | ||
| + | , unde b_low sunt coeficienții $b_i$ și a_low coeficienții $a_i$ ai unui filtru IIR trece-jos. Observați diferențe față de cele FIR? | ||
| /*</hidden>*/ | /*</hidden>*/ | ||
