This shows you the differences between two versions of the page.
ps:laboratoare:08 [2019/11/16 18:39] ionut.gorgos |
ps:laboratoare:08 [2020/10/07 18:31] (current) ionut.gorgos |
||
---|---|---|---|
Line 6: | Line 6: | ||
Î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). | Î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 și zero-padding [4p] === | + | === Exercițiul 1 -- DFT leakage și zero-padding [4p] === |
Î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: | Î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: | ||
Line 23: | Line 23: | ||
- 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ă e centrată în jurul frecvenței semnalului (2.5kHZ). | - 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ă e centrată în jurul frecvenței semnalului (2.5kHZ). | ||
- Acum schimbati din nou frecvența la $f2=2000$ Hz, dar 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. | - Acum schimbati din nou frecvența la $f2=2000$ Hz, dar 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> | ||
- | <note> | ||
- | Solution without zero-padding: | ||
- | <code matlab lab7_sinewaves.m> | ||
- | close all; | ||
- | clear; | ||
- | fs = 8000; | ||
- | f1 = 0000; | ||
- | f2 = 2500; | ||
- | A1 = 1; | ||
- | A2 = 0.5; | ||
- | N = 8; | ||
- | |||
- | t = 0:(N-1); | ||
- | s1 = A1*sin(2*pi*(f1/fs)*t); | ||
- | s2 = A2*sin(2*pi*(f2/fs)*t); | ||
- | s = s1 + s2; | ||
- | |||
- | %% Plot signals | ||
- | h = figure; | ||
- | plot(t, s1, 'r--'); | ||
- | hold on; | ||
- | 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 [5p] === | === Exercițiul 2 -- DFT leakage și ferestre [5p] === | ||
- | Î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. | + | Î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. |
- | Să facem urmatoarele: | + | Să facem următoarele: |
- | - 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). | + | - 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 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. | + | - 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. |
- | - 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ă. | + | - 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ă. |
- 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 129: | 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 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ță. |
- | 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(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.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 note, dar 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 [1p] === | + | |
- | + | ||
- | 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> | ||