Differences

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

Link to this comparison view

ps:laboratoare:08 [2018/11/18 14:36]
darius.necula
ps:laboratoare:08 [2020/10/07 18:31] (current)
ionut.gorgos
Line 1: Line 1:
 ===== Laboratorul 08. ===== ===== Laboratorul 08. =====
-<​hidden>​ +/*<​hidden>​*/ 
-==== DFT in detaliu: DFT leakage, zero-padding ====+==== DFT în detaliu: DFT leakage, zero-padding ====
  
  
-In acest laborator vom continua ​sa exploram transformata ​Fourier ​Discreta ​(DFT), ​urmarind ​efectul ​esantionarii in domeniul ​frecventa ​(aparitia ​sinc-ului din cauza fenomenului de leakage) ​si metode de rezolvare a acestuia (zero-padding,​ ferestre, ​cresterea numarului ​de esantioane).+Î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).
  
-=== Exercitiul ​1 -- DFT leakage ​si zero-padding [4p] ===+=== Exercițiul ​1 -- DFT leakage ​și zero-padding [4p] ===
  
-In acest exercitiu veti reface exemplul pe care l-am facut la curs cu 2 sinusoide pentru a vedea efectul fenomenului de leakage ​si a experimenta zero-padding. Pentru asta vom folosi ​urmatorul ​semnal:+În acest exercițiu veți reface exemplul pe care l-am făcut ​la curs cu 2 sinusoide pentru a vedea efectul fenomenului de leakage ​și a experimenta zero-padding. Pentru asta vom folosi ​următorul ​semnal:
  
 $s(n) = A_1 \sin(2\pi f_1 n t_s) + A_2 \sin(2\pi f_2 n t_s)$, $s(n) = A_1 \sin(2\pi f_1 n t_s) + A_2 \sin(2\pi f_2 n t_s)$,
  
  
-unde $f_1$ si $f_2$ sunt frecventele ​celor doua sinusoide care compun semnalul, ​si $t_s = 1/f_s$ e perioada de sampling ​($f_s = 1/t_s$).+unde $f_1$ si $f_2$ sunt frecvențele ​celor două sinusoide care compun semnalul, ​și $t_s = 1/f_s$ e perioada de eșantionare ​($f_s = 1/t_s$).
  
-Pentru a face asta urmariti urmatorii pasi+Pentru a face asta urmăriți următorii pași
-  - Creati ​si plotati ​semnalul, folosind $A_1=1$, $A_2=0.5$, $f_s = 8000$ Hz, $f_1 = 1000$ Hz, $f_2 = 2000$ Hz pentru $N=8$ esantioane+  - 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
-  - Calculati ​DFT pentru acest semnal ​si plotati ​magnitudinea acesteia, ca in laboratoarele anterioare. Ar trebui ​sa obtineti ​ceva de genul urmator<​code></​code>​{{:​ps:​laboratoare:​lab7_sinewaves_a.png?​300|}} {{:​ps:​laboratoare:​lab7_sinewaves_a_fft.png?​300|}}  +  - Calculați ​DFT pentru acest semnal ​și plotați ​magnitudinea acesteia, ca în laboratoarele anterioare. Ar trebui ​să obțineti ​ceva de genul următor 
-  - Apoi eliminati ​prima sinusoida ​(ex.: facand ​$f_1=0$) ​si verificati daca aveti semnal doar la 2kHz. +{{:​ps:​laboratoare:​lab7_sinewaves_a.png?​300|}} {{:​ps:​laboratoare:​lab7_sinewaves_a_fft.png?​300|}}  
-  - Schimbati ​$f_2$ de la 2kHz la $f_2=2500$ Hz. Plotati ​spectrul. Ce putem observa? Ar trebui ​sa vedeti ​ca toata energia de la frecventa ​de 2.5kHz a fost distribuita ​pe frecventele ​adiacente. ​Dupa cum am invatat ​la curs acesta este fenomenul cunoscut ca "DFT leakage", ​si apare din cauza faptului ​ca folosirea unui numar finit de esantioane ​poate fi modelata ​ca inmultirea ​unui semnal infinit ​esantionat ​cu o functie rectangulara ​(al carui spectru este un sinc). ​Inainte, nu vedeam acest efect pentru ​ca esantioanele ​sinc-ului erau exact in punctele unde sinc-ul era 0. +  - Apoi eliminați ​prima sinusoidă ​(ex.: făcând ​$f_1=0$) ​și verificați dacă aveți ​semnal doar la 2kHz. 
-  - Pentru a vedea mai bine efectul de leakage trebuie ​sa crestem numarul ​de sample-uri ​folosite pentru DFT. Pentru asta adaugati 0-uri semnalului vostru. De exemplu ​adaugati ​56 de zerouri ca sa obtineti ​un total de $K=64$ ​esantioane ​(din care doar $N=8$ sunt diferite de 0). Apoi calculati ​DFT pentru acest semnal. Ar trebui ​sa vedeti forma sinc-ului mult mai clar si de asemenea ​ca centrata in jurul frecventei ​semnalului (2.5kHZ). +  - Schimbați ​$f_2$ de la 2kHz la $f_2=2500$ Hz. Plotați ​spectrul. Ce putem observa? Ar trebui ​să vedeti ​că toată ​energia de la frecvența ​de 2.5kHz a fost distribuită ​pe frecvențele ​adiacente. ​După cum am învățat ​la cursacesta 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. 
-  - Acum schimbati din nou frecventa ​la $f2=2000$ Hz dar folosind ​in continuare zero-padding ​si plot-ati ​DFT. Ar trebui ​sa vedeti ca intradevar ​sinc-ul era acolo, dar esantioanele ​de la $1000, 3000, 4000 \ldots$ erau 0.+  - 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ăugati ​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ă vedeti forma sinc-ului mult mai clară și de asemenea ​că centrată în jurul frecvenței ​semnalului (2.5kHZ). 
 +  - Acum schimbati din nou frecvența ​la $f2=2000$ Hzdar folosind ​în continuare zero-padding ​și plotați ​DFT. Ar trebui ​să vedeți că într-adevăr ​sinc-ul era acolo, dar eșantioanele ​de la $1000, 3000, 4000 \ldots$ erau 0.
  
-/​*<​hidden>​*/​ +=== Exercițiul 2 -- DFT leakage și ferestre [5p] ===
-<​note>​ +
-Solution without zero-padding:​ +
-<code matlab lab7_sinewaves.m>​ +
-close all; +
-clear; +
-fs 8000; +
-f1 0000; +
-f2 2500; +
-A1 1; +
-A2 0.5; +
-8;+
  
-t = 0:(N-1)+În acest exercițiu aveți dat un semnal({{:ps:​laboratoare:​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.
-s1 = A1*sin(2*pi*(f1/​fs)*t);​ +
-s2 = A2*sin(2*pi*(f2/​fs)*t)+
-= s1 + s2;+
  
-%% Plot signals +Să facem următoarele
-h = figure; +  - 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). 
-plot(t, s1, '​r--'​);​ +  - Calculați ​DFT pentru semnal și plotați magnitudinea,​ ca în laboratoarele precedente. Ar trebui să obțineți ceva precum aceasta: {{:​ps:​laboratoare:​notes_signal_fftpos.png?​300|}}. Î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
-hold on; +  - 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), precum am discutat la curs. Ideea este să generăm o funcție fereastră pe care o vom înmulții cu semnalul original. Plotați semnalul după aplicarea funcției fereastră.
-plot(t, s2, '​b-.'​);​ +
-plot(t, s, '​k-'​);​ +
-xlabel('​Sample index'​);​ +
-ylabel('​Amplitude'​);​ +
-title('​Two sinewaves'​);​ +
-legend('​s1',​ '​s2',​ 's1 + s2'​);​ +
-print(h, '​-dpng',​ '​lab7_sinewaves_c.png'​);​ +
- +
- +
-%% Compute/​plot fft +
-ffs = fft(s); +
-fidx = (fs/​N)*linspace(0,​ N-1, N); +
-h = figure; +
-stem(fidx, abs(ffs));​ +
-title('​FFT of signal'​);​ +
-xlabel('​DFT frequency index'​);​ +
-print(h, '​-dpng',​ '​lab7_sinewaves_c_fft.png'​);​ +
-</​code>​ +
- +
-Code with zero-padding:​ +
-<code matlab sinewaves_zeropad.m>​ +
-close all; +
-clear; +
-fs = 8000; +
-f1 = 0; +
-f2 = 2000; +
-A1 = 1; +
-A2 = 0.5; +
-N = 8; +
-Z = 64; +
- +
-t = 0:(N-1); +
-tt = 0:(Z-1); +
-s1 = A1*sin(2*pi*(f1/​fs)*t);​ +
-s2 = A2*sin(2*pi*(f2/​fs)*t);​ +
- +
-%% Zero-pad the signals +
-s1 = [s1, zeros(1, Z-N)]; +
-s2 = [s2, zeros(1, Z-N)]; +
-s = s1 + s2; +
- +
- +
-%% Plot signals +
-h = figure; +
-plot(tt, s1, '​r--'​);​ +
-hold on; +
-plot(tt, s2, '​b-.'​);​ +
-plot(tt, s, '​k-'​);​ +
-xlabel('​Sample index'​);​ +
-ylabel('​Amplitude'​);​ +
-title('​Two sinewaves'​);​ +
-legend('​s1',​ '​s2',​ 's1 + s2'​);​ +
-print(h, '​-dpng',​ '​lab7_sinewaves_zeropad_b.png'​);​ +
- +
-%% Compute/​plot fft +
-ffs = fft(s); +
-fidx = (fs/​Z)*linspace(0,​ Z-1, Z); +
-h = figure; +
-stem(fidx, abs(ffs));​ +
-title('​FFT of signal'​);​ +
-xlabel('​DFT frequency index'​);​ +
-print(h, '​-dpng',​ '​lab7_sinewaves_zeropad_fft_b.png'​);​ +
-</​code>​ +
-</​note>​ +
-/​*</​hidden>​*/​ +
- +
- +
- +
- +
- +
- +
- +
- +
- +
- +
-=== Exercițiul 2 -- DFT leakage și ferestre [4p] === +
- +
-În acest exercițiu aveți dat un semnal({{:​ps:​laboratoare:​notes_signal.mat|click aici}}) care conține două note(două unde sinusoidale). Însă, una dintre ele este mult mai puternică decât cealaltă, asa că a doua, cea mai slabă nu e ușor de detectat din spectrul semnalului. În acest exercițiu vom încerca sa folosim o funcție fereastră pentru a determina cele doua note. +
- +
-Să facem urmatoarele+
-  - Incărcați și plotați semnalul dat. Ar trebui ​sa 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 DTF). +
-  - Calculați ​DTF pentru semnal și plotați magnitudinea,​ ca în laboratoarele precedente. Ar trebui să obțineți ceva precum aceasta: {{:​ps:​laboratoare:​notes_signal_fftpos.png?​300|}}. În acest moment probabil nu puteți spune care sunt cele doua frecvențe ale semnalului, din cauza faptului că funcția sinc a primei ​sinusoidale ​acoperă componenta celei de-a doua sinusoidale+
-  - Folosind zero-padding în acest caz nu va ajuta prea mult (încercați). ​Asa că vom aplica semnalului o funcție ​fereastra(ex. '​Hanning'​ sau '​Hamming'; ​cautați aceste funcții în Matlab folosind help ), precum am discutat la curs. Ideea este sa generăm o funcție fereastră pe care o vom înmulții 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 rezultat după aplicarea funcției fereastră. Puteți spune, cel puțin aproximativ,​ care sunt cele două frecvente conținute de semnal?
  
-Semnalul, fereastra ​henning si semnalul atenuat ar trebui ​sa arate asa:+Semnalul, fereastra ​hanning și semnalul atenuat ar trebui ​să arate așa:
  
 {{:​ps:​laboratoare:​lab08_notes_signal.png?​300|}} {{:​ps:​laboratoare:​lab08_hanning_window.png?​300|}} ​ {{:​ps:​laboratoare:​lab08_notes_signal.png?​300|}} {{:​ps:​laboratoare:​lab08_hanning_window.png?​300|}} ​
Line 136: Line 40:
  
  
-/​*<​hidden>​*/​ +=== Exercițiul 3 -- Tratarea DFT leakage prin creșterea numărului de eșantione [1p] ===
-<​note>​ +
-O posibilă soluție: +
-<code matlab notes_signal_short.m>​ +
-clear; +
-close all; +
-gen_signal ​0;+
  
-%% Generate or load the 2-note signal +Precum am discutat la cursfrecvențele date de sinc pot fi reduse prin creșterea numărului de eșantioane diferite de zero ale semnalului nostruAșa căatunci când este posibilaceasta ne va ajuta să vizualizăm semnale foarte aproape în frecvență.
-if gen_signal +
-    %% Generate the 2-note signal ​    +
-    fs = 1000; +
-    N = 64; +
-    x = 0:(N-1); +
-    f1 = 392; % (G4 - Sol major) +
-    A1 = 10; +
-    f2 = 440; % (A4 - La major) +
-    A2 = 1; +
-    s = zeros(1N); +
-    s1 = A1*sin(2*pi*f1/​fs*x);​ +
-    s2 = A2*sin(2*pi*f2/​fs*x);​ +
-    notes_signal = s1+s2; +
-    save('​notes_signal.mat''​notes_signal'​'​fs'​);​ +
-else +
-    %% Load signal +
-    fname = '​notes_signal.mat';​ +
-    load(fname);​ +
-end+
  
-%% Plot signal +Pentru a vedea acest efect, ​să utilizăm aceleași notedar cu un semnal mult mai lung {{:​ps:​laboratoare:​notes_signal_long.mat|click aici}}.
-N = length(notes_signal);​ +
-t = (N/​fs)*linspace(0,​ 1, N); +
-h = figure; +
-plot(t, notes_signal);​ +
-title('​Original signal'​);​ +
-xlabel('​Time (s)'​);​ +
-ylabel('​Amplitude'​);​ +
-print(h, '​-dpng',​ '​notes_signal.png'​);​ +
- +
-%% Plot fft of the signal +
-ys = fft(notes_signal);​ +
-fidx = linspace(0,​fs/​2-1,​N/​2);​ +
-h = figure; +
-stem(fidx, abs(ys(1:​N/​2)));​ +
-xlabel('​Frequency (Hz)'​);​ +
-title('​Frequency spectrum of notes signal'​);​ +
-print(h, '​-dpng',​ '​notes_signal_fftpos.png'​);​ +
- +
-% %% Try zero-padding +
-% Z = 128; +
-% s_zeropadding = zeros(1, Z); +
-% s_zeropadding(1:​N) = notes_signal;​ +
-%  +
-% %% Plot fft after zero-padding +
-% ys = fft(s_zeropadding);​ +
-% fidx = linspace(0,​fs/​2-1,​N/​2);​ +
-% h = figure; +
-% stem(fidx, abs(ys(1:​N/​2)));​ +
-% xlabel('​Frequency (Hz)'​);​ +
-% title('​Frequency spectrum of notes signal with zero-padding'​);​ +
-% print(h, '​-dpng',​ '​notes_signal_zp_fftpos.png'​);​ +
- +
-%% Try to use a window to limit the sinc's effect +
-w1 = hanning(N);​ +
-w1 = w1(:​)';​ +
-sw = notes_signal .* w1; +
- +
-%% Plot windowed signal +
-h = figure; +
-plot(t, sw); +
-title('​Windowed signal'​);​ +
-xlabel('​Time (s)'​);​ +
-ylabel('​Amplitude'​);​ +
-print(h, '​-dpng',​ '​notes_signal_window.png'​);​ +
- +
-%% Plot fft after window +
-ys = fft(sw); +
-fidx = linspace(0,​fs/​2-1,​N/​2);​ +
-h = figure; +
-stem(fidx, abs(ys(1:​N/​2)));​ +
-xlabel('​Frequency (Hz)'​);​ +
-title('​Frequency spectrum of notes signal with window'​);​ +
-print(h, '​-dpng',​ '​notes_signal_window_fftpos.png'​);​ +
-</​code>​ +
-</​note>​ +
-/​*</​hidden>​*/​ +
- +
-=== Exercise 3 -- Tratarea DFT leakage prin creșterea numărului de eșantione [2p] === +
- +
-Precum am discutat la curs, frecventele date de sinc pot fi reduse prin creșterea numărului de eșantioane diferite de zero ale semnalului nostru. Asa că, atunci când este posibil, aceasta ne va ajuta sa vizualizăm semnale foarte aproape în frecvență. +
- +
-Pentru a vedea acest efect, ​sa utilizăm aceleași note dar cu un semnal mult mai lung {{:​ps:​laboratoare:​notes_signal_long.mat|click aici}}.+
  
 Procedați la fel ca înainte: Procedați la fel ca înainte:
-  - Plotați semnalul și spectrul ​sau. Verificați dacă puteți distinge cele doua frecvențe(ar trebui).+  - Plotați semnalul și spectrul ​său. Verificați dacă puteți distinge cele două frecvențe(ar trebui).
   - Aplicați funcția fereastră și verificați spectrul. Ar trebui sa fie mult mai clar.   - 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'​.   - Ce note muzicale reprezintă aceste frecvențe? Puteți să redați acest sunet folosind funcția Matlab '​sound'​.
- 
-/​*<​hidden>​*/ ​ 
-<​note>​ 
-O posibilă soluție: 
-<code matlab notes_signal_long.m>​ 
-clear; 
-close all; 
-gen_signal = 0; 
- 
-%% Generate or load the 2-note signal 
-if gen_signal 
-    %% Generate the 2-note signal ​   ​ 
-    fs = 1000; 
-    N = 1024; 
-    x = 0:(N-1); 
-    f1 = 392; % (G4 - Sol major) 
-    A1 = 10; 
-    f2 = 440; % (A4 - La major) 
-    A2 = 1; 
-    s = zeros(1, N); 
-    s1 = A1*sin(2*pi*f1/​fs*x);​ 
-    s2 = A2*sin(2*pi*f2/​fs*x);​ 
-    notes_signal = s1+s2; 
-    save('​notes_signal_long.mat',​ '​notes_signal',​ '​fs'​);​ 
-else 
-    %% Load signal 
-    fname = '​notes_signal_long.mat';​ 
-    load(fname);​ 
-end 
- 
-%% Plot signal 
-N = length(notes_signal);​ 
-t = (N/​fs)*linspace(0,​ 1, N); 
-h = figure; 
-plot(t, notes_signal);​ 
-title('​Original signal'​);​ 
-xlabel('​Time (s)'); 
-ylabel('​Amplitude'​);​ 
-print(h, '​-dpng',​ '​notes_signal.png'​);​ 
- 
-%% Plot fft of the signal 
-ys = fft(notes_signal);​ 
-fidx = linspace(0,​fs/​2-1,​N/​2);​ 
-h = figure; 
-stem(fidx, abs(ys(1:​N/​2)));​ 
-xlabel('​Frequency (Hz)'​);​ 
-title('​Frequency spectrum of notes signal'​);​ 
-print(h, '​-dpng',​ '​notes_signal_fftpos.png'​);​ 
- 
-%% Try to use a window to limit the sinc's effect 
-w1 = hanning(N); 
-w1 = w1(:​)';​ 
-sw = notes_signal .* w1; 
- 
-%% Plot windowed signal 
-h = figure; 
-plot(t, sw); 
-title('​Windowed signal'​);​ 
-xlabel('​Time (s)'); 
-ylabel('​Amplitude'​);​ 
-print(h, '​-dpng',​ '​notes_signal_window.png'​);​ 
- 
-%% Plot fft after window 
-ys = fft(sw); 
-fidx = linspace(0,​fs/​2-1,​N/​2);​ 
-h = figure; 
-stem(fidx, abs(ys(1:​N/​2)));​ 
-xlabel('​Frequency(Hz)'​);​ 
-title('​Frequency spectrum of notes signal with window'​);​ 
-print(h, '​-dpng',​ '​notes_signal_window_fftpos.png'​);​ 
-</​code>​ 
-</​note>​ 
-</​hidden>​ 
  
ps/laboratoare/08.1542544570.txt.gz · Last modified: 2018/11/18 14:36 by darius.necula
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