Differences
This shows you the differences between two versions of the page.
|
ps:labs:06 [2020/11/11 13:21] ionut.gorgos |
ps:labs:06 [2022/11/08 14:58] (current) ionut.gorgos |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| ===== Laboratorul 06. ===== | ===== Laboratorul 06. ===== | ||
| /*<hidden>*/ | /*<hidden>*/ | ||
| - | ==== Semnale digitale - procesarea simplă şi eşantionarea ==== | + | ==== SNR, decibeli și DFT ==== |
| + | Prezentarea PowerPoint pentru acest laborator poate fi găsită aici: [[https://docs.google.com/presentation/d/1s85BCzqZ3ziRwyqL1ecqilqaMDBMRxq-/edit?usp=share_link&ouid=110538702824281541719&rtpof=true&sd=true|aici]] | ||
| - | În acest laborator vom începe să experimentăm efectele semnalelor eşantionate şi procesarea semnalelor digitale. | + | În acest laborator vom rezolva câteva exerciții legate de calcularea raportului semnal-zgomot(eng. signal-to-noise ratio - SNR), vom exprima acest raport în decibeli și vom căpăta experiență cu transformata Fourier discretă (DFT). Deși nu am discutat prea mult despre ea, în toate exercițiile care folosesc DFT vom utiliza funcția //fft// (Fast Fourier Transform) și inversa acesteia //ifft//. |
| - | Materiale ajutătoare: | + | Materiale utile: |
| - | - [[http://www.ece.rice.edu/~dhj/courses/elec241/col10040.pdf|D. Johnson's book]] | + | * {{:ps:labs:fundamentals-of-electrical-engineering-i-9.72.pdf|Cartea lui Don Johnson}}, secțiunile 5.4‐5.7 |
| - | - Secţiunile 5.1-5.3 | + | |
| + | În exercițiile din acest laborator vom considera că DFT are același număr de componente ca și semnalul de intrare. Formulele pentru DFT și IDFT: | ||
| + | \begin{equation} | ||
| + | DFT: S(k) = \sum^{N-1}_{n = 0}{s(n)e^{\frac{-j 2 \pi n k}{K}}}, k \in \{0, ..., K-1\} \\ | ||
| - | === Exerciţiul 1 -- Procesarea simplă [4p] === | + | IDFT\ (inversa\ DFT): s(n) = \frac{1}{K}\sum^{K-1}_{k = 0}{S(k)e^{\frac{j 2 \pi n k}{K}}}, n \in \{0, ..., N-1\} |
| - | Vom încerca să implementăm un sistem de procesare digitală simplă, care funcţionează după cum urmează: | + | \end{equation} |
| - | * luaţi 5 eşantioane (samples) din semnalul de intrare, $d_1, d_2, \ldots, d_5$. | + | === Exercițiul 1 ‐‐ unde cu zgomot === |
| - | * faceţi o medie a celor 5 eşantioane (samples), $a = 0.2 \cdot (d_1 + d_2 + d_3 + d_4 + d_5)$ | + | [<color red>6p</color>] |
| - | * puneţi la ieşire valoarea curentă $y(n) = a$ | + | |
| - | Urmăriţi aceşti paşi: | + | În acest exercițiu aveți dat un semnal cu zgomot, ({{:ps:labs:noisy_signal.mat|click aici}}), pentru care știți că informația de interes se află în primele 10 componente DFT, corespunzătoare, de fapt cu $k \in \{ -9, ..., 9\}$, adică componenta continuă / directă (DC: $k = 0$) și componentele frecvențelor pozitive $k = 1, ..., 9$ și negative $k = −1, ... , −9$. |
| - | - Creaţi câteva secvenţe digitale care reprezintă sinusoide de diferite frecvenţe (1, 2, 10, 20, 100 Hz), la aceeaşi frecvenţă de eşantionare (folosiţi acelaşi număr de eşantioane, e.g. $N=128$). | + | |
| - | - Implementaţi sistemul de procesare menţionat mai sus | + | |
| - | - Creaţi secvenţele corespunzătoare de ieşire. | + | |
| - | - Plotaţi toate secvenţele de intrare/ieşire | + | |
| - | - Ce fel de sistem de procesare este acesta ? | + | |
| - | - Cum ați putea implementa sistemul de procesare menționat mai sus, folosind funcția 'filter'? (opțional) | + | |
| - | <note tip> | + | |
| - | Pentru fiecare eşantion de intrare ($x(n)$) trebuie să produceţi o valoare de ieşire ($y(n)$). Pentru primele patru eşantioane nu aveţi suficiente valori de intrare diferite de zero, şi pentru simplitate ori setaţi ieşirea pe zero, ori faceţi media dintre 1, 2, 3, 4 eşantioane, în timp ce consideraţi ca celelalte să fie zero. Dar de la al 5-lea eşantion veţi avea câte 5 eşantioane la îndemână (actualul $x(n)$ şi anterioarele patru eşantioane). | + | |
| - | </note> | + | |
| - | <note> Pentru crearea timpilor de eșantionare puteți folosi 0:perioada_esantionare:timp_maxim sau funcția linspace(0,timp_maxim,numar_esantioane)</note> | + | Mai jos puteți vedea semnalul cu zgomot precum și spectrul pozitiv( folosind doar $k = \{0, 1, ..., N/2 - 1\}$). |
| - | === Exercițiul 2 -- Subeşantionare (Bandpass sampling) [3p] === | + | {{:ps:labs:noisy_signal.png?300|}} |
| + | {{:ps:labs:noisy_signal_fftpos.png?300|}} | ||
| - | În acest exercițiu vom experimenta subeșantionarea, metodă care reduce frecvența de eșantionare, necesară pentru a obține întreg spectru al semnalului. | + | Vrem să calculam SNR în domeniul frecvență și apoi să eliminăm zgomotul din semnal folosind DFT și păstrând doar frecvențele de interes. Pentru aceasta va trebui să urmați pașii următori: |
| - | Un semnal real cu o frecvență purtătoare (centrată) la 3kHz este eşantionată la 8 kHz. Lăţimea de bandă a semnalului este 1 kHz, cum puteţi vedea mai jos: | + | 1. Încărcați semnalul zgomotos și plotați-l. Folosiți acest cod ca să încărcați semnalul: |
| - | {{:ps:labs:bandpass.png?500|}} | + | <code> |
| - | + | load('noisy_signal.mat') | |
| - | Task-uri/întrebări: | + | |
| - | - De ce avem energie în ambele frecvenţe $f_c$ şi $-f_c$ ? | + | |
| - | - Folosiţi Octave pentru a calcula transformata Fourier Discretă a unei sinusoide de $3kHz$ eşantionată la $8kHz$. Folosiți $N=64$ pentru fft. Folosiți acest cod pentru DFT: <code matlab> | + | |
| - | figure; | + | |
| - | fx = zeros(1, N); | + | |
| - | fidx = (fs/N) * linspace(0,N-1,N); | + | |
| - | spectrum = fft(x1, N); | + | |
| - | stem(fidx, abs(spectrum)); | + | |
| - | xlabel('Frequency (Hz)'); | + | |
| - | ylabel('Amplitude'); | + | |
| - | title('Spectrum of x1'); | + | |
| </code> | </code> | ||
| - | - Plotaţi (cu 'stem') rezultatele şi observaţi frecvenţele. Acum folosiţi "fftshift" pe rezultatul obţinut de la "fft" şi plotaţi din nou rezultatele. Folosiți acest cod: <code matlab> | ||
| - | figure; | ||
| - | fidx = (fs/N)*linspace(-N/2, N/2-1, N); | ||
| - | stem(fidx, abs(fftshift(spectrum))); | ||
| - | xlabel('Frequency (Hz)'); | ||
| - | ylabel('Amplitude'); | ||
| - | title('Zero-centred frequency spectrum of x1'); | ||
| - | </code> | ||
| - | - De ce într-un caz avem semnal la $3 kHz$ şi $5 kHz$ în timp ce, în celălalt caz pare să avem la $-3 kHz$ şi $3 kHz$ ? | ||
| - | - Care este semnificaţia semnalului de $5 kHz$ ? Când plotaţi rezultatele de la "fft" luaţi în considerare că, în mod implicit output-ul transformatei Discrete Fourier, deci şi "fft", arată frecvenţele de la $0$ la $f_s$ (frecvenţa de eşantionare). | ||
| - | - Putem să reducem frecvenţa de eşantionare în timp ce încă detectăm spectrul dorit de frecvenţe ? (vedeţi bandpass sampling). Cât de mult putem reduce frecvenţa de eşantionare? Ce se întâmplă dacă reducem frecvenţa de eşantionare sub 2B? | ||
| - | === Exercițiul 3 -- Reconstruire de semnal [3p] === | + | Acest semnal are $N = 128$ eșantioane și puteți considera că frecvența de eșantionare este $f_s = N = 128$ Hz. |
| - | Pentru reconstruirea unui semnal ideea este: | + | 2. Calculați DFT pentru semnalul zgomotos folosind formula de mai sus, apoi comparați rezultatul cu cel al funcției fft din MATLAB. Plotați spectrul, atât în forma dată de //fft//, cât și în forma centrată în zero, formată aplicând funcția //fftshift// pe rezultatul de la //fft//. Ne amintim că funcția //fft// conține coeficienții pentru frecvențele pozitive între $0$ și $N/2 - 1$, iar următorii coeficienți aparțin de fapt spectrului negativ, de la $-N/2$ la $-1$. [<color red>1p</color>] |
| - | - dacă se face sampling pe un semnal cu $f_s > 2 \cdot W$ (unde $W$ este frecvenţa maximă a semnalului, i.e. respectiv Nyquist), atunci puteţi să reconstruiţi din semnalul digital semnalul analog (exceptând o eroare de cuantizare) pe baza unui low pass filter (cu un sinc). | + | |
| - | Vom simula un semnal sinus analog, pe care îl vom eșantiona, apoi îl vom reconstrui. Pentru asta ar trebui să faceţi următorii paşi: | + | 3. Calculați raportul semnal-zgomot în domeniul frecvență, ca raportul dintre puterea semnalului și puterea zgomotului. Pentru aceasta vom considera că semnalul util se află doar în componentele de frecvență $k \in \{-9, ..., 9\}$ pe când zgomotul se află în toate componentele de frecvență. Formula pentru calcularea puterii unui semnal este: $power = \frac{1}{N}\sum^{N-1}_{k = 0}{|S(k)|^2}$, unde N este numărul total de componente considerate în semnal. [<color red>1p</color>] |
| - | * simulați un semnal sinus 'analog'. Pentru aceasta generați un semnal sinus de frecventa 1Hz cu un număr mare de puncte (i.e. 5000 de puncte). Pentru timp folosiți: <code matlab> t = linspace(0, 5, N_analog); </code> (i.e. sinusoida face 5 perioade). | + | |
| - | * construiți semnalul digital, adică varianta eșantionată a semnalului analog, folosind un număr mai mic de puncte (i.e. de 100 de ore mai puține) din semnalul analog. | + | 4. Care este echivalentul acestui //SNR// în decibeli(//dB//)? [<color red>1p</color>] |
| - | * faceți din semnalul digital un semnal continuu (analog), prin crearea unui semnal de aceeași lungime ca cel analog inițial, folosind valorile din semnalul digital pentru valorile de la momentele eșantionate și zero în rest. | + | |
| - | * faceţi convoluţie între un sinc (de frecvenţă fs, dar centrat pe zero, adică un semnal simetric) şi semnalul continuu creat anterior. | + | 5. Din ieșirea funcției //fft// eliminați frecvențele din afara zonei de interes (faceți-le zero). Aveți grijă la indicii frecvențelor pozitive și negative și luați în considerare că primul element corespunde componentei directe, care nu are corespondent negativ. [<color red>1p</color>] |
| + | |||
| + | 6. Pentru a reconstrui semnalul corespunzător acestui spectru, care are componente diferite de zero doar pentru frecvențele de interes, calculați transformata Fourier discretă inversă (//IDFT//) folosind formula de mai sus, iar apoi comparați rezultatul cu cel al funcției //ifft// din MATLAB. Plotați acest semnal și comparați-l cu semnalul original zgomotos. Astfel, practic ați implementat un tip foarte simplu de filtru trece-jos. [<color red>2p</color>] | ||
| + | |||
| + | Aici aveți o imagine cu semnalul fără zgomot. | ||
| + | |||
| + | {{:ps:labs:reconstructed.png?300|}} | ||
| + | |||
| + | Arată precum semnalul vostru reconstruit? | ||
| + | |||
| + | <note> Pentru a putea reface perfect semnalul discret original avem nevoie să calculam un spectru de lungime egală cu semnalul, adică K trebuie să fie egal cu N. </note> | ||
| + | |||
| + | === Exercițiul 2 ‐‐ sunet cu zgomot === | ||
| + | [<color red>4p</color>] | ||
| + | |||
| + | Acest exercițiu este similar cu cel precedent. Acum aveți dat o înregistrare a unei secvențe vorbite ({{:ps:labs:noisy_sound.mat|click aici}}) care a fost amestecată cu zgomot de frecvență înaltă. Trebuie sa aplicați pașii precedenți pentru a elimina zgomotul și pentru a determina secvența vorbită. | ||
| + | |||
| + | |||
| + | Aici aveți ploturile sunetului cu zgomot precum și spectrul său pozitiv: | ||
| + | |||
| + | {{:ps:labs:noisy_sound.png?300|}} | ||
| + | {{:ps:labs:noisy_sound_fftpos.png?300|}} | ||
| + | |||
| + | Pentru acest sunet știm că semnalul de interes se află în intervalul 0-500 Hz și că sunetul a fost înregistrat la o frecvență de eșantionare de $f_s = 8000$ Hz. Observați că trebuie să convertiți indecșii DFT (adică $k = 0, 1, ..., N-1$) în frecvențe știind că fiecare element DFT corespunde unei frecvențe $\frac{k \cdot f_s}{N}$. Folosind această informație ar trebui să determinați ce frecvențe să păstrați, astfel încât să filtrați zgomotul, precum în exercițiul anterior. | ||
| + | |||
| + | Calculați SNR și filtrați zgomotul, asemănător ca la exercițiul anterior. Pentru calculul puterii zgomotului trebuie să luați în considerare că zgomotul nu se mai află în toate componentele de frecvență ca la exercițiul anterior. | ||
| + | |||
| + | Pentru a încărca sunetul și frecvența de eșantionare folosiți următorul cod: | ||
| - | Pentru a face convolutie cu un semnal sinc, folositi codul urmator: | ||
| - | <note> | ||
| <code> | <code> | ||
| - | N_sinc = N_analog; | + | S = load('noisy_sound.mat'); |
| - | t_sinc = linspace(-0.2, 0.2, N_sinc/10); | + | s = S.noisy_sound; |
| - | sincvec = sinc(samples_digital*t_sinc); | + | fs = S.fs; |
| - | s_cont_filtered = conv(s_cont, sincvec); | + | |
| </code> | </code> | ||
| - | </note> | ||
| - | , unde samples_digital, reprezintă numărul de eșantioane pentru semnalul digital. | ||
| - | <note warning> | ||
| - | Trebuie să folosiţi funcţia stem, nu plot pentru reprezentare ! | ||
| - | </note> | ||
| + | Pentru a asculta sunetul stocat în vectorul //s//, folosiți codul următor: | ||
| + | <code> | ||
| + | sound(s) | ||
| + | </code> | ||
| + | |||
| + | Ascultați sunetul înainte și după eliminarea zgomotului. Dacă ați eliminat zgomotul într-un mod corect, atunci veți putea determina secvența vorbită. | ||
| + | |||
| + | **Atenție** la redarea semnalului înaintea procesării pentru că zgomotul înalt poate fi foarte neplăcut. Așa că setați volumul foarte jos pentru această etapă. După filtrare puteți mări volumul pentru a auzi secvența vorbită. | ||
| + | <note important>În acest exercițiu folosiți pentru calculul spectrului funcția fft, respectiv ifft pentru reconstruirea semnalului.</note> | ||
| + | |||
| + | === Exercițiul 3 ‐‐ transformata fourier discretă pe imagini === | ||
| + | [Bonus <color red>1p</color>] | ||
| + | |||
| + | Până acum am descompus semnale 1D în serii de sinusoide cu transformata Fourier. Însă putem aplica același procedeu și pentru semnale 2D, descompunându-le în sinusoide 2D. Putem să ne gândim la o sinusoidă 2D ca la un val. Cel mai comun semnal 2D este o imagine cu nivele de gri. | ||
| + | |||
| + | Putem aplica transformata Fourier discretă pe imagini cu funcția fft2 din MATLAB. | ||
| + | în acest exercițiu veți face următoarele: | ||
| + | - încărcați o imagine folosind funcția imread din MATLAB: ex: img = imread('peppers.png') | ||
| + | - convertiți imaginea RGB în imagine gri și normalizați: img = double(rgb2gray(img))/255 | ||
| + | - aplicați fft2 pe imagine și obțineți un spectru 2D (o matrice) S(k) | ||
| + | - tăiați frecvențele joase din spectru obținând un nou spectru S1(k). Frecvențele joase 2D sunt cele de frecvență mică pe ambele axe, adică, pe rezultatul dat de fft sunt în colturile matricii. | ||
| + | - tăiați frecvențele înalte din spectru obținând un nou spectru S2(k). Frecvențele înalte sunt formate din tot spectrul (toată matricea) mai puțin colțurile. | ||
| + | - obțineți din fiecare spectru o imagine: reconstructed_img = ifft2(S1) | ||
| + | - afișați cele 2 imagini cu funcția imshow() | ||
| + | |||
| + | Folosiți această imagine: | ||
| + | {{:ps:labs:peppers.png?direct&200|peppers.png}} | ||
| + | |||
| + | Rezultatele ar trebui să arate în felul următor: | ||
| + | |||
| + | {{:ps:labs:lab07_imagine_originala.png?300|}} | ||
| + | {{:ps:labs:lab07_spectru_filtru_low-pass.png?300|}} | ||
| + | {{:ps:labs:lab07_spectru_filtru_high-pass.png?300|}} | ||
| + | {{:ps:labs:lab07_imagine_low_pass.png?300|}} | ||
| + | {{:ps:labs:lab07_imagine_high_pass.png?300|}} | ||
| + | |||
| + | /*</hidden>*/ | ||
