Laboratorul 07.

SNR, decibeli și DFT

Î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 utile:

Î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: \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\} \end{equation}

Exercițiul 1 ‐‐ unde cu zgomot [6p]

În acest exercițiu aveți dat un semnal cu zgomot, (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$.

Mai jos puteți vedea semnalul cu zgomot precum și spectrul pozitiv( folosind doar $k = \{0, 1, ..., N/2 - 1\}$).

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:

1. Încărcați semnalul zgomotos și plotați-l. Folosiți acest cod ca să încărcați semnalul:

load('noisy_signal.mat')

Acest semnal are $N = 128$ eșantioane și puteți considera că frecvența de eșantionare este $f_s = N = 128$ Hz.

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

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.

4. Care este echivalentul acestui SNR în decibeli(dB)?

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.

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 originalul zgomotos. Astfel, practic ați implementat un tip foarte simplu de filtru trece-jos.

Aici aveți o imagine cu semnalul fără zgomot.

Arată precum semnalul vostru reconstruit?

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

lab07_dtf.m
function fss = dft(s,K)
 
N = length(s);
fss = zeros(1, K);
for k=1:K
    for n=1:N
        fss(k) = fss(k) + s(n)*exp(-1j*2*pi*(k-1)*(n-1)/K);
    end
end

lab07_idtf.m
function s = idft(fss,N)
 
K = length(fss);
s = zeros(1, N);
for k=1:K
    for n=1:N
        s(n) = s(n) + fss(k)*exp(1j*2*pi*(k-1)*(n-1)/K);
    end
end
s = s / K;

O posibilă soluție:

noisy_waves.m
close all;
clear;
 
%% Load noisy signal
load('noisy_signal');
N = length(noisy_signal);
t = linspace(0, 1, N);
h = figure;
plot(t, noisy_signal);
xlabel('Time (s)');
title('Signal with noise');
print(h, '-dpng', 'noisy_signal.png');
 
%% Check fft of noisy signal
fss = fft(noisy_signal);
fs = 128;
fidx = 0:127;
h = figure;
stem(fidx, abs(fss));
xlabel('DFT frequency index');
title('Frequency spectrum of noisy signal');
print(h, '-dpng', 'noisy_signal_fft.png');
 
%% Plot only the positive spectrum
h = figure;
stem(fidx(1:N/2-1), abs(fss(1:N/2-1)));
xlabel('DFT frequency index');
title('Frequency spectrum of noisy signal');
print(h, '-dpng', 'noisy_signal_fftpos.png');
 
%% Plot spectrum centred at zero
h = figure;
fidx = linspace(-N/2, N/2-1, N);
stem(fidx, abs(fftshift(fss)));
xlabel('DFT frequency index');
title('zero-centred Frequency spectrum of noisy signal');
print(h, '-dpng', 'noisy_signal_fftcentered.png');
 
%% Compute snr
NS = 10;
psignal = (fss(1:NS)*fss(1:NS)' + fss(end-NS+2:end)*fss(end-NS+2:end)') / (2*NS);
pnoise = fss*fss' / N;
snr = psignal / pnoise
snr_db = 10*log10(snr)
 
%% Filter out frequencies above k=10
ffsf = zeros(1,N);
ffsf(1:11) = fss(1:11);
ffsf(end-9:end) = fss(end-9:end);
figure;
stem(0:127, abs(ffsf));
xlabel('DFT frequency index');
title('Frequency spectrum after removing high freqs');
 
%% Reconstruct signal
sr = ifft(ffsf);
figure;
plot(t, sr);
xlabel('Time (s)');
title('Reconstructed signal after removing high freqs');

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 (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:

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.

Pentru a încărca sunetul și frecvența de eșantionare folosiți următorul cod:

S = load('noisy_sound.mat');
s = S.noisy_sound;
fs = S.fs;

Pentru a asculta sunetul stocat în vectorul s, folosiți codul următor:

sound(s)

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ă.

O posibilă soluție:

noisy_sound.m
close all;
clear;
 
%% Load noisy signal
S = load('noisy_sound.mat');
s = S.noisy_sound;
fs = S.fs;
N = length(s);
t = linspace(0, N/fs, N);
h = figure;
plot(t, s);
xlabel('Time (s)');
title('Signal with noise');
print(h, '-dpng', 'noisy_sound.png');
 
%% Play noisy sound
sound(s);
pause(3);
 
%% Plot fft of noisy sound
ys = fft(s);
fidx = (fs/N)*linspace(0,N-1,N);
h = figure;
stem(fidx, abs(ys));
xlabel('Frequency (Hz)');
title('Frequency spectrum of noisy sound');
print(h, '-dpng', 'noisy_sound_fft.png');
 
%% Plot only the positive spectrum
fidx = (fs/N)*linspace(0,N/2-1,N/2);
h = figure;
stem(fidx, abs(ys(1:N/2)));
xlabel('Frequency (Hz)');
title('Positive frequency spectrum of noisy sound');
print(h, '-dpng', 'noisy_sound_fftpos.png');
 
%% Plot spectrum centred at zero
h = figure;
fidx = (fs/N)*linspace(-N/2, N/2-1, N);
stem(fidx, abs(fftshift(ys)));
xlabel('Frequency (Hz)');
title('zero-centred Frequency spectrum of noisy sound');
print(h, '-dpng', 'noisy_sound_fftcentered.png');
 
%% Compute snr
kratio = 500 / 8000;
nk = floor(N * kratio);
psignal = (ys(1:nk)*ys(1:nk)' + ys(end-nk+2:end)*ys(end-nk+2:end)') / (2*nk);
pnoise = ys(nk+1:end-nk)*ys(nk+1:end-nk)' / (N-2*nk);
snr = psignal / pnoise
snr_db = 10*log10(snr)
 
%% Filter out frequencies above f=500 Hz
sf = zeros(1,N);
sf(1:nk) = ys(1:nk);
sf(end-nk+2:end) = ys(end-nk+2:end);
 
%% Plot filtered spectrum
fidx = (fs/N)*linspace(0,N-1,N);
h = figure;
stem(fidx, abs(sf));
xlabel('Frequency (Hz)');
title('Frequency spectrum after removing high freqs');
 
%% Reconstruct signal
sr = ifft(sf);
h = figure;
plot(t, sr);
xlabel('Time (s)');
title('Reconstructed signal after removing high freqs');
print(h, '-dpng', 'reconstructed_sound.png');
 
%% Play reconstructed sound
sound(sr);

Exercițiul 3 ‐‐ transformata fourier discretă pe imagini [Bonus 2p]

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:

  1. încărcați o imagine folosind funcția imread din MATLAB: ex: img = imread('peppers.png')
  2. convertiți imaginea RGB în imagine gri și normalizați: img = double(rgb2gray(img))/255
  3. aplicați fft2 pe imagine și obțineți un spectru 2D (o matrice) S(k)
  4. 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.
  5. 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.
  6. obțineți din fiecare spectru o imagine: reconstructed_img = ifft2(S1)
  7. afișați cele 2 imagini cu funcția imshow()

Rezultatele ar trebui să arate în felul următor:

lab07_ex3.m
close all
clear all
 
% citim o imagine, o convertim in gri si normalizam
img = imread('peppers.png');
img = double(rgb2gray(img)) / 255;
 
h = figure, imshow(img);impixelinfo
title('imagine originala')
print(h, '-dpng', 'lab07_imagine_originala.png');
 
% aflam spectrul S(k) al imaginii
S = fft2(img);
 
% taiam din spectru frecventele joase 
% frecventele joase sunt doar in colturile spectrului(fara fftshift)
 
S1  = S;
val = 0;
d = 30;
S1(1:d,1:d)                 = val;
S1(1:d,end-d+2:end)         = val;
S1(end-d+2:end,1:d)         = val;
S1(end-d+2:end,end-d+2:end) = val;
 
% taiam frecventele inalte (pastram doar colturile, restul facem 0)
S2 = S;
S2(S1 ~= val) = 0;
 
% afisam cu spectrul trecut prin log ca sa atenuam (doar a vedea pe grafic)
% valorile mari
 
figure, imagesc(log(abs(S1)));
title('spectru filtru high-pass')
print(h, '-dpng', 'lab07_spectru_filtru_high-pass.png');
 
figure, imagesc(log(abs(S2)));
title('spectru filtru low-pass')
print(h, '-dpng', 'lab07_spectru_filtru_low-pass.png');
 
img2 = ifft2(S1);
h = figure, imshow(img2);
title('imagine filtrata cu high-pass')
print(h, '-dpng', 'lab07_imagine_high_pass.png');
 
img3 = ifft2(S2);
h = figure, imshow(img3);
title('imagine filtrata cu low-pass')
print(h, '-dpng', 'lab07_imagine_low_pass.png');

  • Data publicării: 12.11.2018
ps/laboratoare/07.txt · Last modified: 2018/11/17 13:51 by ionut.gorgos
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