Studiul transformatelor aplicate secventelor numerice - lucrarea de laborator

studiul transformatelor aplicate secventelor numerice - lucrarea de laborator

Obiectivele lucrarii

1) Asimilarea functiilor MATLAB pentru realizarea transformatelor aplicate secventelor 1D si 2D;

2 Studiul algoritmilor pentru calculul transformatelor

3) Studiul proprietatilor transformatei Fourier discrete si a catorva aplicatii ale transformatelor;

4) Studiul interactiv al transformatelor aplicate secventelor numerice utilizand mediul DIDACTICIEL.

Desfasurarea lucrarii

Transformata Fourier discreta a unei secvente numerice 1D

Calculati si vizualizati spectrul secventei 1D

cu .

N=8; x=ones(1, N);

X=fft(x);

subplot(311);

stem(x);

title('Reprezentarea secventei discrete x[n]');

subplot(312);

stem(abs(X));

title('Reprezentarea modulului transformatei Fourier');

subplot(313);

stem(angle(X));

title('Reprezentarea fazei transformatei Fourier');

Densitatea spectrala de putere si transformata Fourier

Reprezentati densitatea spectrala de putere a unui semnal format din doua sinusoide de 50 Hz si 120 Hz si un zgomot aditiv, alb, gaussian, de medie nula si dispersie 4.

t=0:0.001:0.8;

x=sin(2*pi*50*t)+sin(2*pi*120*t)+2*randn(1, length(t));

subplot(211); plot(x);

title('Semnalul x(t)');

grid;

X=fft(x, 512);

Px=X.*conj(X)/512;

f=1000*(0:255)/512;

subplot(212); plot(f, Px(1:256));

title('Densitatea spectrala de putere');

grid

Proprietati ale transformatei Fourier

Sa ve verifice proprietatile de translatie (deplasare, decalaj) in timp si de modulatie ale transformatei Fourier. Se va considera semnalul:

figure

w = -pi:2*pi/255:pi;

wo = 0.4*pi; D = 10;

num = [1 2 3 4 5 6 7 8 9];

h1 = freqz(num, 1, w);

h2 = freqz([zeros(1, D) num], 1, w);

subplot(2, 2, 1)

plot(w/(2*pi), abs(h1)); grid

title('Spectrul de amplitudine al secventei initiale')

subplot(2, 2, 2)

plot(w/(2*pi), abs(h2)); grid

title('Spectrul de amplitudine al secventei decalate')

subplot(2, 2, 3)

plot(w/(2*pi), angle(h1)); grid

title('Spectrul de faza al secventei initiale')

subplot(2, 2, 4)

plot(w/(2*pi), angle(h2)); grid

title('Spectrul de faza al secventei decalate')

figure

w = -pi:2*pi/255:pi;

wo = 0.4*pi; num1 = [1 3 5 7 9 11 13 15 17];

L = length(num1);

h1 = freqz(num1, 1, w);

n = 0:L-1;

num2 = exp(wo*i*n).*num1;

h2 = freqz(num2, 1, w);

subplot(2, 2, 1)

plot(w/(2*pi), abs(h1)); grid

title('Spectrul de amplitudine al secventei initiale')

subplot(2, 2, 2)

plot(w/(2*pi), abs(h2)); grid

title('Spectrul de amplitudine al secventei modulate')

subplot(2, 2, 3)

plot(w/(2*pi), angle(h1)); grid

title('Spectrul de faza al secventei initiale')

subplot(2, 2, 4)

plot(w/(2*pi), angle(h2)); grid

title('Spectrul de faza al secventei modulate')

Influenta lungimii semnalului asupra spectrului estimat  cu TFD

Sa se studieze influenta lungimii semnalului asupra spectrului sau, estimat prin intermediul TFD. Se va considera un semnal sinusoidal de frecventa 10 Hz, esantionat la 100 Hz. Sa se studieze apoi efectul procedeului de zero paddding si cazurile in care numarul de perioade este intreg sau nu.

t = 0:0.01:0.5-0.01;

x = cos(20*pi*t);

N = length(x) ; X = fft(x, N);

fp = (0:N-1)/N/0.01; fp = fp - 1/2/0.01;

Referat: Contabilitatea tertilor Referat: Structurarea elementelor de activ si pasiv dupa criterii financiare

stem(fp, fftshift(abs(X)));

axis([ -1/(2*0.01) 1/(2*0.01) 0 25]);

xlabel('frecventa (Hz)')

ylabel('spectru de amplitudine')

X = fft(x, 20*N); N = length(X);

fp = (0:N-1)/N/0.01;

fp = fp - 1/2/0.01;

plot(fp, abs(fftshift(X)));

axis([ -1/(2*0.01) 1/(2*0.01) 0 25]);

xlabel('frecventa (Hz)')

ylabel('spectru de amplitudine')

t = 0:0.01:0.5-0.01-1/20;

x = cos(20*pi*t); N = length(x);

X = fft(x, N);

fp = (0:N-1)/N/0.01;

fp = fp - 1/2/0.01;

plot(fp, abs(fftshift(X)));

axis([ -1/(2*0.01) 1/(2*0.01) 0 25]);

xlabel('frecventa (Hz)'); ylabel('spectru de amplitudine')

Influenta tipului ferestrei asupra spectrului estimat cu TFD

a Sa se genereze o sinusoida de amplitudine 1V, de frecventa 100 Hz, esantionata la 256 Hz, in 32 de puncte. Sa se calculeze TFD in 1024 de puncte considerand succesiv o fereastra dreptunghiulara, apoi triunghiulara, Hamming, Hanning si Blackman. Sa se concluzioneze asupra rezolutiilor spectrale si dinamice.

t=(1:32); f1=50; Fe=256; Nfft=1024 ;

y1=sin(2*pi*f1/Fe*t); sig=y1.*boxcar(32)' ;

y_rect=abs(fftshift((fft(sig, Nfft))));

sig=y1.*triang(32)';

y_tria=abs(fftshift((fft(sig, Nfft))));

sig=y1.*hamming(32)' ;

y_hamm=abs(fftshift((fft(sig, Nfft))));

sig=y1.*hanning(32)' ;

y_hann=abs(fftshift((fft(sig, Nfft))));

sig=y1.*blackman(32)';

y_blac=abs(fftshift((fft(sig, Nfft))));

f=[-Fe/2:Fe/Nfft:(Fe/2-Fe/Nfft)];

subplot(511)

semilogy(f(513:1024), y_rect(513:1024)); axis([0 128 1e-3 100]); grid

legend('fereastra rectangulara', -1)

subplot(512); semilogy(f(513:1024), y_tria(513:1024));

axis([0 128 1e-3 100]); grid

legend('fereastra triangulaire', -1)

subplot(513); semilogy(f(513:1024), y_hamm(513:1024));

grid; ylabel('amplitudine spectrala')

axis([0 128 1e-4 100]); legend('fereastra Hamming', -1)

subplot(514); semilogy(f(513:1024), y_hann(513:1024));

axis([0 128 1e-3 100]); grid ; legend('fereastra Hanning', -1)

subplot(515); semilogy(f(513:1024), y_blac(513:1024));

axis([0 128 1e-3 100]); grid

legend('fereastra Blackman', -1)

b Sa se genereze un semnal compus din suma a doua sinusoide de frecvente 100 si 94 Hz si de amplitudine 1V. Se va considera aceeasi lungime a semnalului de 32 puncte, iar TFD va fi calculata tot in 1024 de puncte.

t=(1:32);

f1=94; f2=100 ;

Fe=256; Nfft=1024 ;

y1= sin(2*pi*f1/Fe*t)+ sin(2*pi*f2/Fe*t);

sig=y1.*boxcar(32)' ;

y_rect=abs(fftshift((fft(sig, Nfft))));

sig=y1.*triang(32)';

y_tria=abs(fftshift((fft(sig, Nfft))));

sig=y1.*hanning(32)' ;

y_hann=abs(fftshift((fft(sig, Nfft))));

sig=y1.*hamming(32)' ;

y_hamm=abs(fftshift((fft(sig, Nfft))));

sig=y1.*blackman(32)';

y_blac=abs(fftshift((fft(sig, Nfft))));

f=[-Fe/2:Fe/Nfft:(Fe/2-Fe/Nfft)];

subplot(511)

plot(f(513:1024), y_rect(513:1024));

grid; axis([0 128 0 20])

legend('fereastra rectangulara', -1)

subplot(512);

plot(f(513:1024), y_tria(513:1024)); grid;

axis([0 128 0 7])

legend('fereastra triangulaire', -1)

subplot(513);

plot(f(513:1024), y_hann(513:1024));

grid;

axis([0 128 0 7])

legend('fereastra Hanning', -1)

ylabel('amplitudine spectrala')

subplot(514);

plot(f(513:1024), y_hamm(513:1024));

grid;

axis([0 128 0 7])

legend('fereastra Hamming ', -1)

subplot(515);

plot(f(513:1024), y_blac(513:1024));

grid;

axis([0 128 0 6])

legend('fereastra Blackman', -1)

xlabel('frecventa (Hz)')

c Sa se genereze semnalul urmator in 32 de puncte:

Reprezentati spectrul sau utilizand diverse ponderari si interpretati rezultatele.

t=(1:32); f1=74; f2=100 ; Fe=256; Nfft=1024 ;

y1= 0.1*sin(2*pi*f1/Fe*t)+ sin(2*pi*f2/Fe*t);

sig=y1.*boxcar(32)' ; y_rect=abs(fftshift(fft(sig, Nfft)));

sig=y1.*triang(32)'; y_tria=abs(fftshift(fft(sig, Nfft)));

sig=y1.*hanning(32)' ; y_hann=abs(fftshift(fft(sig, Nfft)));

sig=y1.*hamming(32)' ; y_hamm=abs(fftshift(fft(sig, Nfft)));

sig=y1.*blackman(32)' ; y_blac=abs(fftshift(fft(sig, Nfft)));

f=[-Fe/2:Fe/Nfft:(Fe/2-Fe/Nfft)];

subplot(511)

plot(f(513:1024), y_rect(513:1024)); grid; axis([0 128 0 20])

legend('fereastra rectangulara', -1)

subplot(512); plot(f(513:1024), y_tria(513:1024)); grid;

axis([0 128 0 10])

legend('fereastra triangulaire', -1)

subplot(513); plot(f(513:1024), y_hann(513:1024)); grid;

axis([0 128 0 10])

legend('fereastra Hanning', -1)

ylabel('amplitudine spectrala')

subplot(514); plot(f(513:1024), y_hamm(513:1024)); grid;

axis([0 128 0 10])

legend('fereastra Hamming ', -1)

subplot(515); plot(f(513:1024), y_blac(513:1024)); grid;

axis([0 128 0 8])

legend('fereastra Blackman', -1) ; xlabel('frecventa (Hz)')

Transformata cosinus discreta 1D

Generati o secventa sinusoidala de 10 Hz, esantionata cu 1000 Hz. Calculati si vizuzalizati coeficientii descompunerii in DCT.

t=0:0.001:1;

x=sin(2*pi*10*t);

y=dct(x);

stem(y(1:40));

Compresia de imagini cu transformata cosinus discreta 2D

Sa se studieze efectul compresiei utilizand TCD asupra unei imagini simple. Se va utiliza pentru aceasta imaginea dots", care se gaseste in fisierul "imdemos" al toolbox-ului "images" si functia dct2 din MATLAB.

a) Sa se reprezinte variatia erorii medii patratice in functie de numarul de coeficienti retinuti.

b) Sa se realizeze compresia conservand numai cei mai importanti 1024 coeficienti ai transformatei.

a)

load imdemos dots; [N1, N2]=size(dots); sl=dots; s=zeros(N1, N2);

s(sl)=ones(1, length(find(sl)));

y=dct2(s); [ly, cy]=size(y); yv=y(:);

[yvs, idxs]=sort(abs(yv)); yvs=flipud(yvs);

idxs=flipud(idxs); cont=0;

for r=2.^[0:14]

yvr=yv; yr=y;

yvr(idxs(r+1:N1*N2))=zeros(N1*N2-r, 1);

for k=1:N2

yr(:, k)=yvr((k-1)*N1+1:k*N1, 1);

end

sr=idct2(yr); cont=cont+1;

eqm(cont)=(sum(sum(s-sr).^2))/(N1*N2);

end

plot(eqm);

set(gca, 'XTickLabels', 2.^(str2num(get(gca, 'XTickLabels'))-1))

xlabel('Numar de coeficienti retinuti'); ylabel('Eroare')

title('Variatia erorii medii patratice')

b)

load imdemos dots

r=1024; [N1, N2]=size(dots); sl=dots; s=zeros(N1, N2);

s(sl)=ones(1, length(find(sl)));

y=dct2(s); [ly, cy]=size(y); yv=y(:);

[yvs, idxs]=sort(abs(yv)); yvs=flipud(yvs);

idxs=flipud(idxs); cont=0;

yv(idxs(r+1:N1*N2))=zeros(N1*N2-r, 1);

for k=1:N2

y(:, k)=yv((k-1)*N1+1:k*N1, 1);

end

sr=idct2(y);

subplot(121); imagesc(s)

colormap(gray); title('Imagine initiala')

subplot(122); imagesc(sr)

colormap(gray); title('Imagine reconstruita')

Compresia de imagini utilizand TKL

Sa se realizeze compresia aceleeasi imagini utilizand transformarea Karhunen-Loève, aplicata vectorilor obtinuti prin impartirea imaginii in blocuri de 4 4 pixeli.

a) Sa se reprezinte variatia erorii medii patratice in functie de numarul de coeficienti retinuti.

b) Sa se realizeze compresia imaginii conservand numai componenta cea mai importanta a fiecarui vector din spatiul transformat.

a)

load imdemos dots; [N1, N2]=size(dots); sl=dots; s=zeros(N1, N2);

s(sl)=ones(1, length(find(sl)));

x=im2col(s, [4 4], 'distinct'); xr=x'; sigx=cov(xr);

[vectp, valp]=eig(sigx);

vlpr=diag(valp); [ld, idx]=sort(vlpr);

idx=flipud(idx); ld=flipud(ld);

for r=1:16

vectps='=[';

for k=1:r

no=num2str(k);

eval(['vectps=[vectps ''vectp(:, idx(', no, ')) '']; '])

end

vectps=[vectps ']''; '];

eval(['kmat', vectps]);

kmatc=zeros(16, 16); kmatc(1:r, :)=kmat; sc=kmat*x;

srcol=kmat'*sc;

sr=col2im(srcol, [4 4], [N1 N2], 'distinct');

eqm(r)=(sum(sum(s-sr).^2))/(N1*N2);

end

plot(eqm); xlabel('Numar de coeficienti retinuti');

ylabel('Eroare')

title('Variatia erorii medii patratice')

b)

load imdemos dots

r=1; [N1, N2]=size(dots); sl=dots; s=zeros(N1, N2);

s(sl)=ones(1, length(find(sl)));

x=im2col(s, [4 4], 'distinct'); xr=x'; sigx=cov(xr);

[vectp, valp]=eig(sigx);

vlpr=diag(valp); [ld, idx]=sort(vlpr);

idx=flipud(idx); ld=flipud(ld);

vectps='=[';

for k=1:r

no=num2str(k);

eval(['vectps=[vectps ''vectp(:, idx(', no, ')) '']; '])

end

vectps=[vectps ']''; '];

eval(['kmat', vectps]); kmatc=zeros(16, 16); kmatc(1:r, :)=kmat;

sc=kmat*x; srcol=kmat'*sc;

sr=col2im(srcol, [4 4], [N1 N2], 'distinct');

subplot(121); imagesc(s);

colormap(gray); title('Imagine initiala')

subplot(122); imagesc(sr);

colormap(gray); title('Imagine reconstruita')

Studiul transformatei Fourier  discrete utilizand mediul DIDACTICIEL

1) Se lanseaza DIDACTICIEL-ul prin introducerea comenzii:

didact

2) Se studiaza interactiv transformata Fourier discreta cu ajutorul meniurilor definite in:

Fourier analysis

Tema

Sa se verifice relatia lui Parceval:

in cazul secventei numerice: