Show page

Differences

This shows you the differences between two versions of the page.

--- ps:labs:07 [2020/11/25 09:20]
ionut.gorgos
+++ ps:labs:07 [2022/11/15 11:18] (current)
ionut.gorgos
@@ Line 1: / Line 1: @@
 ===== Laboratorul 07. =====
 /*<hidden>*/
-==== SNR, decibeli și DFT ====
+==== DFT în detaliu: DFT leakage, zero-padding ====
+Prezentarea PowerPoint pentru acest laborator poate fi găsită aici: [[https://docs.google.com/presentation/d/1I0jGMdwlAporLfDDVjQGpfOmVjAoRhKS/edit?usp=sharing&ouid=110538702824281541719&rtpof=true&sd=true|aici]]
-Î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.
+În acest laborator vom continua să explorăm Transformata Fourier Discretă (DFT), urmărind efectul eșantionării în domeniul frecvență (apariția sinc-ului din cauza fenomenului de leakage) și metode de rezolvare a acestuia (zero-padding, ferestre, creșterea numărului de eșantioane).
-Materiale utile:
+=== Exercițiul 1 -- DFT leakage și zero-padding ===
-   * [[http://www.ece.rice.edu/~dhj/courses/elec241/col10040.pdf|D. Johnson's book]], secțiunile 5.4‐5.7
+[<color red>4p</color>]
-În exercițiile din acest laborator vom considera că DFT are același număr de componente ca și semnalul de intrare. Din curs știm formulele pentru DFT:
+În acest exercițiu veți face un exemplu cu 2 sinusoide pentru a vedea efectul fenomenului de leakage și a experimenta zero-padding. Pentru asta vom folosi următorul semnal:
-\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\} \\
-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\}
+$s(n) = A_1 \sin(2\pi f_1 n t_s) + A_2 \sin(2\pi f_2 n t_s)$,
-\end{equation}
-=== Exercițiul 1 ‐‐ unde cu zgomot [6p] ===
-Î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$.
+unde $f_1$ si $f_2$ sunt frecvențele celor două sinusoide care compun semnalul, $t_s = 1/f_s$ e perioada de eșantionare ($f_s = 1/t_s$ ) și $n=0:N-1$, unde $N$ reprezintă numărul de eșantioane.
-Mai jos puteți vedea semnalul cu zgomot precum și spectrul pozitiv( folosind doar $k = \{0, 1, ..., N/2 - 1\}$).
+Pentru a face asta urmăriți următorii pași:
+  - Creați si plotați semnalul, folosind $A_1=1$, $A_2=0.5$, $f_s = 8000$ Hz, $f_1 = 1000$ Hz, $f_2 = 2000$ Hz pentru $N=8$ eșantioane.
+  - Calculați DFT pentru acest semnal și plotați (cu //stem//) magnitudinea acesteia, ca în laboratoarele anterioare. Ar trebui să obțineti ceva de genul următor:
+{{:ps:labs:lab7_sinewaves_a.png?300|}} {{:ps:labs:lab7_sinewaves_a_fft.png?300|}}
+  - Apoi eliminați prima sinusoidă (ex.: făcând $f_1=0$) și verificați dacă aveți semnal doar la 2kHz. [<color red>1p</color>]
+  - Schimbați $f_2$ de la 2kHz la $f_2=2500$ Hz. Plotați spectrul (cu //stem//). Ce putem observa? Ar trebui să vedeți că toată energia de la frecvența de 2.5kHz a fost distribuită pe frecvențele adiacente. Acesta este fenomenul cunoscut ca "DFT leakage", și apare din cauza faptului că folosirea unui număr finit de eșantioane poate fi modelată ca înmulțirea unui semnal infinit eșantionat cu o funcție rectangulară (al cărui spectru este un sinc). Înainte, nu vedeam acest efect pentru că eșantioanele sinc-ului erau exact în punctele unde sinc-ul era 0. [<color red>1p</color>]
+  - Pentru a vedea mai bine efectul de leakage trebuie să creștem numărul de eșantioane folosite pentru DFT. Pentru asta adăugați zerouri semnalului vostru. De exemplu, adăugați la finalul semnalului 56 de zerouri ca să obțineți un total de $K=64$ eșantioane (din care doar $N=8$ sunt diferite de 0). Apoi calculați DFT pentru acest semnal. Ar trebui să vedeți forma sinc-ului mult mai clară și de asemenea că e centrată în jurul frecvenței semnalului (2.5kHZ). [<color red>1p</color>]
+  - Acum schimbați din nou frecvența la $f2=2000$ Hz, dar folosind în continuare zero-padding și plotați DFT (cu //stem//). Ar trebui să vedeți că într-adevăr sinc-ul era acolo, dar eșantioanele de la $1000, 3000, 4000 \ldots$ erau 0. [<color red>1p</color>]
-{{:ps:labs:noisy_signal.png?300|}}
+=== Exercițiul 2 -- DFT leakage și ferestre ===
-{{:ps:labs:noisy_signal_fftpos.png?300|}}
+[<color red>5p</color>]
-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:
+În acest exercițiu aveți dat un semnal({{:ps:labs:notes_signal.mat|click aici}}) care conține două note (două unde sinusoidale). Însă, una dintre ele este mult mai puternică decât cealaltă, așa că, a doua, cea mai slabă, nu e ușor de detectat din spectrul semnalului. În acest exercițiu vom încerca să folosim o funcție fereastră pentru a determina cele două note.
-. Încărcați semnalul zgomotos și plotați-l. Folosiți acest cod ca să încărcați semnalul:
+Să facem următoarele:
+  - Incărcați și plotați semnalul dat. Ar trebui să observați că se vor încărca variabilele //notes_signal// și //fs//, unde //fs// este frecvența de eșantionare (amintiți-vă ca aveți nevoie de ea pentru a înțelege rezultatul dat de DFT).
+  - Calculați DFT pentru semnal și plotați spectrul pozitiv al semnalului (cu //stem//), ca în laboratoarele precedente. Ar trebui să obțineți ceva precum aceasta:
-<code>
+{{:ps:labs:notes_signal_fftpos.png?300|}}
-load('noisy_signal.mat')
-</code>
-Acest semnal are $N = 128$ eșantioane și puteți considera că frecvența de eșantionare este $f_s = N = 128$ Hz.
+În acest moment probabil nu puteți spune care sunt cele două frecvențe ale semnalului, din cauza faptului că funcția sinc a primei sinusoide acoperă componenta celei de-a doua sinusoide.
+  - Folosind zero-padding în acest caz nu va ajuta prea mult (încercați). Așa că vom aplica semnalului o funcție fereastră (ex. //Hanning// sau //Hamming//; căutați aceste funcții în MATLAB folosind help). Ideea este să generăm o funcție fereastră pe care o vom înmulți cu semnalul original. Plotați semnalul după aplicarea funcției fereastră.
+  - Calculați DFT pentru semnalul rezultat după aplicarea funcției fereastră. Puteți spune, cel puțin aproximativ, care sunt cele două frecvente conținute de semnal?
-. 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$.
+Semnalul, fereastra //Hanning// și semnalul atenuat ar trebui să arate așa:
-. 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.
+{{:ps:labs:lab08_notes_signal.png?300|}} {{:ps:labs:lab08_hanning_window.png?300|}}
+{{:ps:labs:lab08_notes_signal_window.png?300|}}
-. Care este echivalentul acestui SNR în decibeli(dB)?
-. 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.
+=== Exercițiul 3 -- Tratarea DFT leakage prin creșterea numărului de eșantioane ===
+[<color red>1p</color>]
-. 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 originalul zgomotos. Astfel, practic ați implementat un tip foarte simplu de filtru trece-jos.
+Precum am discutat la curs, frecvențele date de sinc pot fi reduse prin creșterea numărului de eșantioane diferite de zero ale semnalului nostru. Așa că, atunci când este posibil, aceasta ne va ajuta să vizualizăm semnale foarte aproape în frecvență.
-Aici aveți o imagine cu semnalul fără zgomot.
+Pentru a vedea acest efect, să utilizăm aceleași note, dar cu un semnal mult mai lung {{:ps:labs:notes_signal_long.mat|click aici}}.
-{{:ps:labs:waves_signal.png?300|}}
-Arată precum semnalul vostru reconstruit?
+Procedați la fel ca înainte:
+  - Plotați semnalul și spectrul său pozitiv. Verificați dacă puteți distinge cele două frecvențe(ar trebui).
-<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>
+  - Aplicați funcția fereastră și verificați spectrul. Ar trebui sa fie mult mai clar.
+  - Ce note muzicale reprezintă aceste frecvențe? Puteți să redați acest sunet folosind funcția MATLAB //sound//.
-=== Exercițiul 2 ‐‐ sunet cu zgomot [4p] ===
-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:
-<code>
-S = load('noisy_sound.mat');
-s = S.noisy_sound;
-fs = S.fs;
-</code>
-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ă.
-=== Exercițiul 3 ‐‐ transformata fourier discretă pe imagini [Bonus 1p] ===
-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()
-<note>
-Pentru a utiliza funcția **rgb2gray** în Octave trebuie să instalați pachetul Octave **image** din Octave command prompt.
-<code>
-pkg install -forge image
-pkg load image
-</code>
-</note>
-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>*/

Laboratoare

Resurse

Tema

ps/labs/07.1606288822.txt.gz · Last modified: 2020/11/25 09:20 by ionut.gorgos

Show page Old revisions

Media Manager Back to top