This is an old revision of the document!

Laboratorul 08.

DFT in 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).

Exercitiul 1 -- DFT leakage si 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:

$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$).

Pentru a face asta urmariti urmatorii pasi:

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.
Calculati DFT pentru acest semnal si plotati magnitudinea acesteia, ca in laboratoarele anterioare. Ar trebui sa obtineti ceva de genul urmator:
Apoi eliminati prima sinusoida (ex.: facand $f_1=0$) si verificati daca aveti semnal doar la 2kHz.
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.
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 e centrata in jurul frecventei semnalului (2.5kHZ).
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.

Solution without zero-padding:

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 with zero-padding:

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');

Exercițiul 2 -- DFT leakage și ferestre [5p]

În acest exercițiu aveți dat un semnal(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: . Î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?

Semnalul, fereastra henning si semnalul atenuat ar trebui sa arate asa:

O posibilă soluție:

notes_signal_short.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 = 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
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');

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 click aici.

Procedați la fel ca înainte:

Plotați semnalul și spectrul sau. Verificați dacă puteți distinge cele doua frecvențe(ar trebui).
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'.

O posibilă soluție:

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');