Estimare spectrala si analiza timp-frecventa - lucrarea de laborator

estimare spectrala si analiza timp-frecventa - lucrarea de laborator

Obiectivele lucrarii

1 Studiul estimatorilor spectrali clasici, parametrici si neparametrici

2 Studiul metodelor de analiza timp-frecventa a proceselor nestationare

3) Utilizarea functiilor MATLAB pentru implementarea tehnicilor de estimare spectrala si de analiza timp-frecventa;

4) Studiul interactiv al metodelor de estimare spectrala si de analizatimp-frecventa utilizand mediul DIDACTICIEL.

Desfasurarea lucrarii

Estimatorul spectral simplu

Fie un proces alb, gaussian, centrat () si de dispersie , definit pe 2048 esantioane. Sa se realizeze analiza spectrala a acestui proces utilizand estimatorul spectral simplu.

a) Sa se determine deplasarea si dispersia estimatorului.

b) Sa se refaca aceeasi analiza pentru  512, 1024 si 4096 esantioane.

c) Sa se verifice ca acest estimator este deplasat si ca dispersia sa nu depinde de durata observatiei (estimator inconsistent).

longueur_P=2048;

P=randn(1, longueur_P)*sqrt(0.25);

[spectre_simple, Frecventa]=periodo_simple(P, longueur_P);

semilogy(Frecventa, spectre_simple, 'b');      

xlabel('Frecventa normalizata');

ylabel('Amplitudine (dB)');

title('Periodograma simpla');

biais=mean(spectre_simple)

variance=std(spectre_simple)^2                            

Estimatorul spectral mediat

In cazul aceluiasi proces aleator sa se realizeze analiza spectrala folosind estimatorul spectral mediat. Sa se studieze influenta numarului de secvente folosite pentru mediere. Sa se determine deplasarea si dispersia estimatorului.

longueur_sequence=256;   

[spectre_moyenne, Frecventa]=periodo_moyenne(P, longueur_sequence);

semilogy(Frecventa, spectre_moyenne, 'b');

xlabel('frecventa normalizata');

ylabel('Amplitudine (dB)');

hold on;

longueur_sequence=128;   

[spectre_moyenne, Frecventa]=periodo_moyenne(P, longueur_sequence);

semilogy(Frecventa, spectre_moyenne, '--r');

title('Periodograma mediata');

legend('mediere pe 8 ferestre', 'mediere pe 16 ferestre');

Sa se verifice ca acest estimator este tot deplasat, dar consistent. Ce se poate spune despre variatia rezolutiei frecventiale in raport cu numarul de sectiuni mediate

Rezultatele se vor trece intr-un tabel de tipul:

Estimatorul spectral modificat

Sa se realizeze analiza spectrala a aceluiasi proces aleator folosind estimatorul spectral modificat. Sa se studieze influenta numarului de secvente folosite pentru mediere si a diferitelor ferestre utilizate pentru netezire. Sa se determine deplasarea si dispersia estimatorului.

fenetre=[boxcar(longueur_sequence)];

fenetre=fenetre*sqrt(longueur_sequence/ sum(fenetre.*fenetre));

fenetre=fenetre.';

[spectre_modifie, Frecventa]=periodo_modifie(P, longueur_sequence, fenetre);

semilogy(Frecventa, spectre_modifie, 'b');

xlabel('frecventa normalizata');

ylabel('Amplitudine (dB)');

hold on;

title('Periodograma modificata');

Rezolutia dinamica

Sa se realizeze analiza spectrala a unui semnal compus din doua sinusoide afectate de zgomot. Parametrii componentelor semnalului sunt urmatorii:

prima sinusoida:                a doua sinusoida: zgomot alb gaussian:

frecventa = 25 Hz frecventa = 50 Hz medie = 0

amplitudine = 1                 amplitudine = 0.01 deviatie standard = 0.031

Frecventa de esantionare se considera 200 Hz. Lungimea a semnalului va lua valorile urmatoare : =.

Sa se testeze cei trei estimatori spectrali studiati anterior pentru: , , , si ferestrele rectangulara, Hamming si Blackman, fiind numarul de subsecvente in care este impartita secventa initiala.

longueur_signal=512;

t=0:1/200:(longueur_signal/200)-1/200;

y1=sin(2*pi*25*t); y2=0.01 * sin(2*pi*50*t);

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

bruit=randn(1, longueur_signal); signal=y1+y2+bruit*0.031;

[spectre_simple, Frecventa]=periodo_simple(signal, longueur_signal);

semilogy(Frecventa, spectre_simple, 'b'); title('512 puncte');

longueur_sequence=256;

[spectre_moyenne, Frecventa]=periodo_moyenne(signal, longueur_sequence);

semilogy(Frecventa, spectre_moyenne, 'b');

title('K=512, L=2');

longueur_sequence=256;

fenetre=[hamming(longueur_sequence)];

fenetre=fenetre*sqrt(longueur_sequence/ sum(fenetre.*fenetre));

fenetre=fenetre.';      

[spectre_modifie, Frecventa]=periodo_modifie(signal, longueur_sequence, fenetre);

semilogy(Frecventa, spectre_modifie, 'b'); title('Hamming');

xlabel('frecventa normalizata'); ylabel('Amplitudine')

Proces aleator AR

Sa se genereze un zgomot alb, gaussian, centrat, de dispersie 0.25.

Sa se filtreze apoi acest semnal cu ajutorul unui filtru, avand caracteristica de transfer .

Sa se identifice procesul AR rezultat si sa se regaseasca valorile parametrilor (, , ) corespunzator datelor generate. Se va presupune cunoscut ordinul procesului ().

bruit=randn(1, 2048)*0.5;

b=1; a=[1 0.5 0.75];

y=filter(b, a, bruit); y=y.';

[coef, sigma2]=parametre_AR(y, 2);

title('Polii si zerourile unui proces A.R. de ordinul 2');

ordre=2; dsp=0;

for nu=0:200

denominateur=1;

for ind=1:ordre

coef(ind+1, 1);

coefm=coef(ind+1, 1)*exp(-2*pi*j*ind*nu/200);

denominateur=denominateur+coefm ;

end;

dsp=[dsp sigma2/(abs(denominateur).^2)];

end;

nu=0:1/200:0.5;

plot(nu, dsp(1, 1:length(nu)));

title('Raspunsul in frecventa al filtrului');

xlabel('frecventa normalizata');

ylabel('Amplitudine');

Analiza spectrala de inalta rezolutie parametrica

Sa se genereze o sinusoida pura de frecventa 50 Hz esantionata la 200 Hz (definita pe 4096 esantioane). Reprezentati densitatea sa spectrala de putere considerand acest semnal ca un proces AR, al carui ordin se va determina folosind criteriul lui Akaike.

longueur=4096; t=0:1/200:(longueur/200)-1/200; y=sin(2*pi*50*t); y=y.';

dsp=[]; FPE=zeros(1, 10);

for ordre=1:10

[coef, sigma2]=parametre_AR(y, ordre);

FPE(1, ordre)=((longueur+ordre)/(longueur-ordre))*sigma2;

[mini, position]=min(FPE);

ordre_selectionne=position;        

end

[coef, sigma2]=parametre_AR(y, ordre_selectionne);

for nu=0:200

denominateur=1;

for ind=1:ordre_selectionne

coef(ind+1, 1);

coefm= coef(ind+1, 1)*exp(-2*pi*j*ind*nu/200);

denominateur=denominateur+coefm ;

end;

dsp=[dsp sigma2/(abs(denominateur).^2)];

end;

nu=0:1/200:0.5; dsp=dsp/max(dsp);

plot(nu, dsp(1, 1:length(nu))); title('DSP a semnalului');

xlabel('frecventa normalizata'); ylabel('Amplitudine');

Spectrograma unui semnal chirp

Un semnal MLF are urmatoarea expresie analitica:

unde , fiind durata semnalului.

a) Sa se genereze un astfel de semnal pentru .

b) Sa se reprezinte semnalul in domeniile temporal si frecvential.

c) Sa se reprezinte acelasi semnal in planul timp frecventa cu ajutorul spectrogramei.

a)

f1=2000;

f2=8000;

pulselength=0.025;

Fe=20000;   

t=(0:1/Fe:pulselength);

beta=(f2-f1)/(2*pulselength);

chirp=sin(2*pi*(f1+beta*t).*t);

chirp2=vco(sawtooth((2*pi/pulselength)*t, 1), [f1/Fe, f2/Fe]*Fe, Fe);

figure(1); clf;

subplot(211);

plot(t, chirp);

xlabel('Timp');

ylabel('Amplitudine');

title('Analiza unei modulatii liniare de frecventa')

b)

C=fftshift(abs(fft(chirp)).^2);

lc=length(chirp);

mc=lc/2;

freq=(-mc:1:mc-1)*Fe/lc;

subplot(212);

plot(freq, C);

xlabel('Frecventa [Hz]');

ylabel('Densitatea spectrala de putere')

c)

Wsize=32;

[Cspec, F, T] = specgram(chirp, 128, Fe, Wsize);

figure(2); clf;

imagesc(1000*T, F/1000, abs(Cspec));

xlabel('Timp [ms]');

ylabel('Frecventa [kHz]');

title('Analiza unei modulatii liniare de frecventa')

Ce concluzie se poate trage in privinta informatiei aduse de cele 3 reprezentari

Compromisul rezolutie spectrala rezolutie temorala

Sa se genereze un semnal compus dintr o sinusoida si un Dirac.

Sa se analizeze acest semnal cu ajutorul spectrogramei pentru diverse dimensiuni ale ferestrei de analiza.

Sa se interpreteze rezultatele obtinute.

Fe=10000; f1=1000; f2=4000;

T=0.015; t=[0.005:1/Fe:T];

delta=0.005*Fe;

sig=[zeros(1, delta), sin(2*pi*f1*t), zeros(1, delta),

5*ones(1, 1), zeros(1, 2*delta)];

subplot(311);

plot(0:1000*(1/Fe):1000*(length(sig)/Fe-1/Fe), sig);

title('Semnal temporal');

xlabel('Timp [ms]'); ylabel('Amplitudine');

[S, F, T] = specgram(sig, 128, Fe, 64);

subplot(312); imagesc(T*1000, F/1000, abs(S));

xlabel('Timp [ms]'); ylabel('Frecventa [kHz]');

Title('Spectrograma cu o fereastra de 64 puncte'),

[S, F, T] = specgram(sig, 128, Fe, 16);

subplot(313); imagesc(T*1000, F/1000, abs(S));

xlabel('Timp [ms]'); ylabel('Frecventa [kHz]');

title('Spectrograma cu o fereastra de 16 puncte')

Scalograma

Scalograma este o distributie de energie care se obtine ca modulul la patrat al transformatei wavelet continue.

a Undisoara Morlet este definita prin relatia:

Aceasta nu verifica proprietatea de medie nula. Pentru a aproxima aceasta conditie, conservand un numar mic de oscilatii Morlet a propus compromisul:

Sa se vizualizeze alura acestei undisoare pentru diversi parametri.

b) Puneti in evidenta variatiile rezolutiilor temporala si frecventiala in functie de pozitia in planul timp frecventa. Se va utiliza un semnal compus dintr un Dirac si doua sinusoide trunchiate, de frecvente diferite.

a

[onde, ts]= ond_mor(129, 10, 100);

figure(1); clf;

subplot(211); plot(ts, real(onde));

hold on; plot(ts, abs(onde), 'r');

grid; title('Undisoara Morlet')

b)

Fe=64000;

f1=500; f2=8000; T=0.01;

t=[0:1/Fe:T];

delta=0.005*Fe;

sig1=[zeros(1, delta), sin(2*pi*f1*t), zeros(1, 2*delta)];

sig2=[zeros(1, delta), sin(2*pi*f2*t), zeros(1, 2*delta)];

sig=sig1+sig2; sig(4*delta+1:4*delta+5)=10*ones(1, 5);

figure(1); subplot(212); plot(sig);

title('Semnal studiat: suma a doua sinusoide si un Dirac')

NB_OCT=8; NB_VOIES=1;

[TOnd_SIG,   freqs] = morlet( sig, Fe, length(sig), Fe/2,

NB_OCT, NB_VOIES);

vec_Timp=[0:1/Fe:length(sig)/Fe]*1000; vec_freqs=[0:1:7];

figure(2); clf; colormap(gray);

imagesc(vec_Timp, vec_freqs, 2*log10(abs(TOnd_SIG)));

set(gca, 'YTickLabels', 250*2.^([7:-1:0]));

title('Scalograma unui Dirac si doua sinusoide');

xlabel('Timp [ms]'); ylabel('Frecventa [Hz]')

Distributia Wigner-Ville

a Sa se calculeze distributia Wigner-Ville a unui semnal MLF.

b) Sa se genereze un semnal compus din 3 sinusoide trunchiate: doua debutand la momentul , de frecvente si diferite, a treia debutand la momentul , de frecventa .

Sa se reprezinte distributia Wigner-Ville a acestui semnal si sa se concluzioneze asupra interferentelor temporale si frecventiale.

a)

f1=2000; f2=8000; T=0.008; Fe=20000; t=(0:1/Fe:T);

beta=(f2-f1)/(2*T); delta=0.001*Fe;

x=[zeros(1, delta), sin(2*pi*(f1+beta*t).*t),

zeros(1, delta)]; [W, Wr, abscisse]=wigner(x, Fe);

max(max(real(W)))

max(max(imag(W)))

sum(sum(W))

sum(x.^2)

b)

f1=1000; deltaf=3000; f2=f1+deltaf;

Fe=round(3*f2); T=0.01; t=[0:1/Fe:T];

deltat=0.01*Fe; marge=20;

sig1=[zeros(1, marge), sin(2*pi*f1*t),

zeros(1, T*Fe+deltat+marge)];

sig2=[zeros(1, marge), sin(2*pi*f2*t), zeros(1, T*Fe+deltat+marge)];

sig3=[zeros(1, marge+T*Fe+deltat), sin(2*pi*f1*t), zeros(1, marge)]; sig=sig1+sig2+sig3;

[W, Wr, fwv, abscisse]=wigner(sig, Fe); [N1, N2]=size(W);

figure(2); clf; zoom on ;

subplot(311); milieuf=round(((f1+f2)/2)/(Fe/2)*N1);

plot(abscisse, Wr(milieuf, :)); xlabel('Timp (ms)');

title('Studiul interferentelor')

subplot(312); milieut=round(marge+T*Fe+deltat);

plot(fwv, Wr(:, milieut)); xlabel('Frecventa [kHz]')

subplot(313); repere=round(marge+T/2*Fe);

plot(fwv, Wr(:, repere)); xlabel('Frecventa [kHz]')

Analiza multirezolutie

Sa se studieze cu ajutorul analizei multirezolutie un semnal format dintr o sinusoida cu zgomot. Pentru analiza se va folosi undisoara lui Haar.

Sa se realizeze sinteza aceluiasi semnal plecand de la coeficientii descompunerii si sa se arate ca este posibila filtrarea semnalului anuland o parte dintre acestia.

clear; Fe=5000; taille=1024;

Timp=[1:1:taille]/Fe; numpts=length(Timp);

data=sin(2*pi*25*Timp)+rand(1, taille);

[approx, detail]=analysehaar(data);

newsig=synthesehaar(approx, detail);

selection=detail; selection(8:10, :)=zeros(3, numpts);

sig_debruite=synthesehaar(approx, selection);

figure(1); clf; zoom on

for plt = 1:size(detail, 1),

subplot(10, 2, 2*(plt-1)+1), plot(Timp, detail(plt, 1:numpts));

grid;

if (plt == 1),

title('Detaliu');

end;

if (plt == 10);

xlabel('Timp [s]');

end

end

for plt = 1:size(approx, 1),

subplot(10, 2, 2*plt), plot(Timp, approx(plt, 1:numpts)); grid

if (plt == 1),

title('Aproximare');

end;

if (plt == 10);

xlabel('Timp [s]');

end

end

figure(2); clf; zoom on;

subplot(311); plot(Timp, data); title('Semnal original');

subplot(312); plot(Timp, newsig); title('Semnal reconstruit');

subplot(313); plot(Timp, sig_debruite);

title('Semnal filtrat'); xlabel('Timp [s]')

Studiul operatorilor de estimare spectrala si analiza timp frecventa utilizand mediul DIDACTICIEL

1) Se lanseaza DIDACTICIEL-ul prin introducerea comenzii:

didact

2) Se studiaza interactiv operatorii de estimare spectrala si analiza timp frecventa cu ajutorul meniurilor definite in:

Fourier analysisParametric Models

Fourier analysisEstimators

Deep StudyTime-Frequency Analysis

Tema

Sa se realizeze analiza spectrala folosind algoritmul MUSIC a unui semnal compus din doua sinusoide afectate de zgomot. Parametrii componentelor semnalului sunt urmatorii:

prima sinusoida:                     a doua sinusoida: zgomot alb gaussian:

frecventa = 25 Hz frecventa = 30 Hz medie = 0

amplitudine = 1                      amplitudine = 1 deviatie standard=0.25

Frecventa de esantionare se considera 100 Hz. Se considera ca se dispune de 32 esantioane ale semnalului util.

(Indicatie: matricea de autocorelatie se estimeaza folosind subsecvente ale semnalului, obtinute prin parcurgerea acestuia cu o fereastra glisanta de lungime 16, iar dimensiunea subspatiului semnal se considera egala cu 4).