Esempio : Oscillazione con ampiezza e fasi aleatorie

1. Tema 6:Analisi in Potenza di Processi Parametrici Esempio: Oscillazione con ampiezza e fasi aleatorie Tipico segnale presente in sistemi radiomobile, radar, sonar, etc.: scomposizione del fasore con ampiezza e fase aleatorie in somma di due fasoriin quadratura con ampiezze aleatorie: con e È la trasformazione di coordinate polari-cartesiane, per cui si è visto che: N.B.:

2. Analisi in Potenza di Processi Parametrici Oscillazione con ampiezza e fasi aleatorie • Generare M realizzazioni di un p.a. parametrico costituito da una oscillazione con ampiezza distribuita secondo Rayleigh e fase uniforme, utilizzando il metodo della scomposizione: • Parametri:Intervallo di campionamento: Tc=1 ms Intervallo temporale di osservazione: T=NTc=0.1 ms Valore quadratico medio dell’ampiezza dell’oscillazione: Frequenza dell’oscillazione: • [istruzioni utili:randn, repmat, for, plot, hold on]

3. Esempi di realizzazioni M=50 realizzazioni M=3 realizzazioni

4. File.m: genoscampfasealea.m % generazione realizzazioni oscillazione con ampiezza di Rayleigh e fase % uniforme (metodo della scomposizione in due fasori Gaussiani in quadratura) function [x] = genoscfasealea(n,aqm,f0,tfin) % IN: numero di realizzazioni, n;% valor quadratico medio della ampiezza della oscillazione, aqm;% frequenza della oscillazione, f0 (il tempo di campionamento e' unitario);% ampiezza intervallo temporale, tfin (a partire da zero); % OUT: matrice di realizzazioni, una per riga;% uscita su video di grafico delle n realizzazioni. sigma=sqrt(aqm/2); % calcola la dev. standard per le ampiezze dei due fasori x1=sigma*randn(1,n); % genera n realizz. della v.a. Gaussiana ampiezza fasore 1;x2=sigma*randn(1,n); % genera n realizz. della v.a. Gaussiana ampiezza fasore 2;t=[1:tfin]; % calcola i tempi su cui campionare; x1=repmat(x1.',1,tfin); % organizza tempi e ampiezze fasorix2=repmat(x2.',1,tfin); % per calcolo efficientet=repmat(t,n,1); % delle realizzazioni; x=x1.*cos(2*pi*f0*t)-x2.*sin(2*pi*f0*t); % calcola tutte le realizzazioni; figure; % grafica le realizzazioni;for i=1:n;plot(x(i,:))hold onendxlabel('t, sec')ylabel('x(t)')

5. Generazione processi parametrici Oscillazione con ampiezza e fasi aleatorie + rumore AWGN • Generare una realizzazione del processo tempo-discreto “oscillazione con ampiezza e fase aleatorie” corrotto da rumore AWGN (tempo-discreto) • Parametri:Segnale: come prima! Rumore:

6. Analisi in Potenza di un Processo SSL Definizione di densità spettrale di potenza (PSD) di un processo aleatorio tempo-continuo stazionario almeno in senso lato (SSL): PSD dei segnali deterministici “realizzazioni” Per un processo aleatorio stazionario in senso lato si ha inoltre: Teorema di Einstein-Wiener-Khintchine: dove è la funzione di autocorrelazione (ACF)

7. Misura sperimentale della PSD Date M realizzazioni del processo osservate su un intervallo di osservazione di durata finita T, misurare la PSD: interpretazione in termini di frequenza relativa spettri di ampiezza delle realizzazioni Si utilizza la FFT: [ realizzazioni in tempo-discreto, Tc=tempo di campionamento ] Nota bene: M finito causa errori aleatori sulla valutazione della PSD T finito causa errori deterministici sulla valutazione della PSD originati dalla “finestratura” temporale (leakage): T=NTc Il campionamento può causare ovviamente errori di aliasing

8. Valutazione sperimentale della PSD • Disponendo di M=2000 realizzazioni del p.a. “oscillazione con ampiezza di Rayleigh e fase uniforme§” valutare sperimentalmente la PSD del processo (effettuare una FFT con zero-padding a 1024 campioni) • [istruzioni utili:fft, abs, mean, fftshift, plot] • Confrontare il risultato con PSD “deterministica” di tre realizzazioni[istruzioni utili:fft, abs, plot, fftshift, hold on]

9. M=2000 realizzazioni Effetto della finestratura:le divengono “discrete” Esempio di risultati

10. Esempio di risultati PSD del processo aleatorio PSD misurata da tre realizzazioni (la PSD é “phase blind”) fluttuazioni di ampiezza da realizzazione a realizzazione

11. File.m: calcdsppa.m % calcolo funzione densita' spettrale di potenza di un p.a.function [dsp,f] = calcdsppa(x) % in: matrice di realizzazioni del p.a. (una per riga, % con tempo di campionamento assunto unitario), x; % out: funzione densita' spettrale di potenza, dsp;% ascisse frequenziali, f;% uscita su video di grafico di densita' spettrale di potenza. tfin=size(x,2); % finestra temporale disponibile n=size(x,1); % numero di realizzazioni disponibilitc=1; % tempo di campionamento unitariozeropad=1024; % zero-padding a 1024 campioni if tfin>zeropad disp('aumentare lo zero-padding!') % (per tc=1)end spettri=tc*fft(x,zeropad,2); % calcola gli spettri complessi delle realizzazionidsprealizzazioni=abs(spettri).^2/(tfin-1); % calcola le dsp deterministiche delle realizzazionidsp=mean(dsprealizzazioni); % calcola la dsp del processo dsp=fftshift(dsp); % centra la componente a frequenza zeroncamp=zeropad; % numero di campioni dopo lo zero-paddingdeltaf=1/tc/ncamp; % passo di campionamento in frequenzaf=[-ncamp/2:ncamp/2-1]*deltaf; % scala delle frequenze figure; % grafico dsp del processoplot(f,dsp);xlabel('f , Hz');ylabel('dsp');

12. File.m: calcdspdet.m % calcolo densita' spettrali di potenza deterministiche di 3 realizzazioni di un p.a.function calcdspdet(x) % in: matrice di realizzazioni del p.a. (una per riga, % con tempo di campionamento assunto unitario), x; % out: uscita su video di grafico di densita' spettrali di potenza deterministiche. tfin=size(x,2); % finestra temporale disponibiletc=1; % tempo di campionamento unitariozeropad=1024; % zero-padding a 1024 campioni if 3>size(x,1) disp('non sono disponibili le realizzazioni richieste!') endif tfin>zeropad disp('aumentare lo zero-padding!') % (per tc=1)end spettri=tc*fft(x(2:4,:),zeropad,2); % calcola gli spettri complessi delle 3 realizzazionidsprealizzazioni=abs(spettri).^2/(tfin-1); % calcola le dsp deterministiche delle 3 realizzazioni ncamp=zeropad; % numero di campioni dopo lo zero-paddingdeltaf=1/tc/ncamp; % passo di campionamento in frequenzaf=[-ncamp/2:ncamp/2-1]*deltaf; % scala delle frequenze figure; % grafico dsp delle 3 realizzazioniplot(f,fftshift(dsprealizzazioni(1,:)),'g');hold onplot(f,fftshift(dsprealizzazioni(2,:)),'k');hold onplot(f,fftshift(dsprealizzazioni(3,:)),'r');xlabel('f , Hz');ylabel('dsp');

13. Valutazione della funzione valor medio ed ACF Sono date M realizzazioni del processo, osservate su un intervallo di durata finita T: [ con t nell’intervallo temporale osservato ] Interpretazione in termini di frequenza relativa

14. Risultati per il p.a. “oscillazione con ampiezza di Rayleigh e fase uniforme” Funzione valor medio valutata con M=2000 realizzazioni t

15. Funzione di autocorrelazione (ACF) valutata con M=2000 realizzazioni ACFcalcolata per t=1 t

16. Verifica di stazionarietà in autocorrelazione: ACFcalcolata per t=5 t ACFcalcolata per t=1 t L’ACF, nei limiti dell’errore di valutazione sperimentale, non dipende da t

17. Interpretazione dei risultati di autocorrelazione: Scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo t Coefficiente di correlazione tra X(t) e X(t+t) = Covarianza tra le due v.a. / prodotto delle deviazioni standard = cX(t)/cX(0) t=19 (t=1) t=19 si nota che le variabili sono correlate positivamente

18. Interpretazione dei risultati di autocorrelazione: Scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo t t=29 t=24 le variabili sono correlate negativamente le variabili sono quasi incorrelate

19. Interpretazione dei risultati di autocorrelazione: Se l’ampiezza non é più aleatoria ma costante (ad es. di valore 1): il processo di segnale non è più Gaussiano, si vede anche dallo scatterplot! (t=1) t=19 scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo t

20. File.m: calcmediaeautocorrpa.m % calcolo funzione valor medio e funzione di autocorrelazione di un p.a.function [eta,r] = calcmediaeautocorrpa(x,t1) % in: matrice di realizzazioni del p.a. (una per riga, % con tempo di campionamento assunto unitario), x;% istante di riferimento per la funzione di autocorr. (intero non nullo), t1; % out: funzione valor medio, eta;% funzione di autocorr. (in funzione del tempo di ritardo, positivo), r;% uscita su video di grafico di valor medio e autocorrelazione. tfin=size(x,2); eta=mean(x); % calcola la funzione valor medio; xt1=repmat(x(:,t1),1,tfin); % organizza realizz. di x(t1) per calcolo efficiente% della funzione di autocorrelazione; r=mean(xt1.*x); % calcola la funzione di autocorrelazione; r=r(t1+1:tfin); % seleziona i valori con tempo di ritardo positivo; figure;plot(eta)axis([0,tfin,min(min(x)),max(max(x))])figure;plot(r)v=axis;axis([0,tfin,v(3),v(4)])

21. Filtraggio di un treno di oscillazioni in rumore AWGN Esempio: segnale in una comunicazione digitale, o in sistema radar/sonar/ecografico • Parametri del processo aleatorio costituito da un teno di oscillazioni immerso in rumore addtivo Gaussiano bianco in banda:Intervallo di campionamento: 1 sDurata intervallo di osservazione: 500 sTempo di inizio oscillazione: 100 sDurata temporale oscillazione: 100 sAmpiezza oscillazione: 1Frequenza oscillazione: 0.05 HzPotenza rumore: 0.1 [il rumore é bianco in banda di Nyquist] • [istruzioni utili: randn]

22. Filtraggio di un treno di oscillazioni in rumore AWGN • Filtrare una realizzazione del processo aleatorio prima definito mediante un filtro passa basso e poi mediante un filtro passa-bandaProgetto dei filtri: FIR, 200 prese (finestra di Hamming) tempo di campionamento: unitario freq. taglio passabasso: 0.06 Hz ... freq. taglio inf. e sup. passabanda: 0.04 Hz, 0.06 Hz ... [istruzioni utili:FDATOOL, save, load, filter]

23. Esempio di risultati Il treno di oscillazioni immerso in una realizzazione di rumore

24. Esempio di risultati Risultato del filtraggio passa-basso

25. Esempio di risultati Risultato del filtraggio passa-banda

26. Esempio di risultati

27. File.m: genoscrum.m % generazione realizzazione spezzone di oscillazione immersa in rumore Gaussianofunction [x] = genoscrum() % out: realizzazione;% uscita su video di grafico della realizzazione. tfin=500; % ampiezza intervallo temporale (tempo di campionamento unitario);tin=100; % istante di inizio oscillazione;tdur=100; % durata temporale oscillazione;a=1; % ampiezza della oscillazione;f0=0.05; % frequenza della oscillazione;pot=0.1; % potenza del rumore; t=[tin:tin+tdur-1]; % calcola tempi su cui campionare lo spezzone;x=sqrt(pot)*randn(1,tfin); % genera una realizzazione di rumore Gaussiano bianco% (in banda di Nyquist); x(tin:tin+tdur-1)=x(tin:tin+tdur-1)+a*sin(2*pi*f0.*t); % somma lo spezzone; plot(x)axis([1 tfin -a-3*sqrt(pot) a+3*sqrt(pot)])

28. File.m: filtrpbasso.m % filtraggio passabasso di spezzone di oscillazione immersa in rumorefunction [y] = filtrpbasso(x) % in: realizzazione di spezzone di oscillazione immersa in rumore;% out: realizzazione filtrata;% uscita su video di grafico della realizzazione filtrata. load cpbasso cpbasso;% carica coefficienti filtro FIR a 200 prese con finestra di Hamming, % freq. taglio 0.06 Hz, tempo di campionamento unitario (da FDATOOL); y=filter(cpbasso,1,x); % esegue il filtraggio plot(y) % filtraggio passabandadi spezzone di oscillazione immersa in rumorefunction [y] = filtrpbanda(x) % in: realizzazione di spezzone di oscillazione immersa in rumore;% out: realizzazione filtrata;% uscita su video di grafico della realizzazione filtrata. load cpbanda cpbanda;% carica coefficienti filtro FIR a 200 prese con finestra di Hamming, % freq. taglio 0.04 - 0.06 Hz, tempo di campionamento unitario (da FDATOOL); y=filter(cpbanda,1,x); % esegue il filtraggio plot(y)