0( ) cos(2 )X t A f t

Esempio: Oscillazione con ampiezza e fasi aleatorie

Tipico segnale presente in sistemi radiomobile, radar, sonar, etc.:

( )A R (0,2 )U indipendentiA e

0 0 0( ) cos(2 ) cos cos(2 ) sin sin(2 )X t A f t A f t A f t

2

2 2 222

( ) ( ), con { } 2a

A

af a e u a E A

1

2

cos

sin

X A

X A

2 21 2

2 1atan2( , )

A X X

X X

21

22

(0, )

(0, )

X

X

N

N 1 2, IIDX X

scomposizione del fasore con ampiezza e fase aleatorie in somma di due fasori in quadratura con ampiezze aleatorie:

1 0 2 0( ) cos(2 ) sin(2 )X t X f t X f t con e

È la trasformazione di coordinate polari-cartesiane, per cui si è visto che:

2 2 2{ } 2 2 { }iE A E X N.B.:

Tema 6: Analisi in Potenza di Processi Parametrici

- 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 s

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]

Oscillazione con ampiezza e fasi aleatorie

Analisi in Potenza di Processi Parametrici

1 0 2 0[ ] ( ) cos(2 ) sin(2 ), 0,1,2, , 1cX n X nT X F n X F n n N

0 0dove: cF f T

2 2{ } 2 1E A

0 50f KHz

Esempi di realizzazioni

M=3 realizzazioni

M=50 realizzazioni

0 10 20 30 40 50 60 70 80 90 100-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

t, sec

x(t)0 10 20 30 40 50 60 70 80 90 100

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

t, sec

x(t)

0( ) cos(2 )X t A f t

File.m: genoscampfasealea.mgenoscampfasealea.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)')

- Generare una realizzazione del processo tempo-discreto “oscillazione con

ampiezza e fase aleatorie” corrotto da rumore AWGN (tempo-discreto)

Parametri:

Segnale: come prima!

Rumore:

Oscillazione con ampiezza e fasi aleatorie + rumore AWGN

Generazione processi parametrici

1 0 2 0[ ] [ ] [ ] cos(2 ) sin(2 ) [ ]Y n X n W n X F n X F n W n

2 2 2

2 2 2 2

{ [ ]} { } 2 1 2

{ [ ]} { [ ]} W W

E X n E ASNR

E W n E W n

2[ ] 0, WW n N

210 1010 2

110Log ( ) 10 10

dB dBSNR SNR

dB WSNR SNR SNR

10, 5,0,5,10,50dBSNR dB

Analisi in Potenza di un Processo SSL

2{| ( ; ) | }( ) { ( ; )} lim T

X X T

E X fS f E S f

T

( ) ( )X XS f R FTTeorema di

Einstein-Wiener-Khintchine:

Definizione di densità spettrale di potenza (PSD) di un processo aleatorio tempo-continuo stazionario almeno in senso lato (SSL):

Per un processo aleatorio stazionario in senso lato si ha inoltre:

PSD dei segnali deterministici “realizzazioni”

( ) ( ) ( )XR E X t X t dove è la funzione di autocorrelazione (ACF)

2 ( ) 2

1

{| ( ; ) | } 1 | ( ) |( ) lim

iMT T

X iT

E X f X fS f

T M T

spettri di ampiezza delle realizzazioni

Misura sperimentale della PSD

interpretazione in termini di frequenza relativa

Date M realizzazioni del processo osservate su un intervallo di osservazione di durata finita T, misurare la PSD:

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

21( )

0( ) ( )

c

k njNi N

kT c i cnfN T

X f T x nT e

Si utilizza la FFT:

[ realizzazioni in tempo-discreto, Tc=tempo di campionamento ]

- 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]

Valutazione sperimentale della PSD

Esempio di risultati

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.50

5

10

15

20

25

f , Hz

dsp

2 2

20 0 0( ) cos(2 ) ( ) ( )

2 2XS f f f f f f FT

Effetto della finestratura:le divengono “discrete”

( ) 2 ( )Sinc

M=2000 realizzazioni

0 0.05 secf

0( ) cos(2 )X t A f t

100 secT

2 2{ } 2 1E A

20( ) cos(2 )XR f

Esempio di risultati

PSD misurata da tre realizzazioni

(la PSD é “phase blind”)

fluttuazioni di ampiezza da realizzazione a realizzazione

PSD del processo aleatorio

-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.20

5

10

15

20

25

30

f , Hz

dsp

-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.20

5

10

15

20

25

30

f , Hz

dsp

File.m: calcdsppa.mcalcdsppa.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');

File.m: calcdspcalcdspdetdet.m.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');

1 1

1 1ˆ( ) { ( ; )} lim ( ) ( ) ( )

M M

X i X ii iMt E X t x t t x t

M M

Valutazione della funzione valor medio ed ACF

Interpretazione in termini di frequenza relativa

Sono date M realizzazioni del processo, osservate su un intervallo di durata finita T:

[ con t nell’intervallo temporale osservato ]

1

1

1( , ) { ( ; ) ( ; )} lim ( ) ( )

1ˆ ( , ) ( ) ( )

M

X i iiM

M

X i ii

R t t E X t X t x t x tM

R t t x t x tM

Funzione valor medio valutata con M=2000 realizzazioni

Risultati per il p.a. “oscillazione

con ampiezza di Rayleigh e fase uniforme”

( ) 0X t

0 10 20 30 40 50 60 70 80 90 100-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

t

Funzione di autocorrelazione (ACF)

valutata con M=2000 realizzazioni

ACF calcolata per

t

0 10 20 30 40 50 60 70 80 90 100-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

1( , ) cos(2 0.05 )

2XR t t

Verifica di stazionarietà in autocorrelazione:

0 10 20 30 40 50 60 70 80 90 100-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

0 10 20 30 40 50 60 70 80 90 100-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

L’ACF, nei limiti dell’errore di valutazione sperimentale, non dipende da t

ACF calcolata per t

ACF calcolata per t

1( , ) cos(2 0.05 )

2XR t t

Interpretazione dei risultati di autocorrelazione:

Scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo

si nota che le variabili sono correlate positivamente

(t

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

2

2

( ) ( )( ) cos(2 0.05 19) 0.95

(0) (0)X X X

XX X X

c r

c r

Coefficiente di correlazione tra

X(t) e X(t+)=

Covarianza tra le due v.a. / prodotto delle deviazioni

standard

=

cX()/cX(0)

le variabili sono correlate negativamente

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

le variabili sono quasi incorrelate

Scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo

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!

scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

(t

Interpretazione dei risultati di autocorrelazione:

% 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)])

File.m: calcmediaeautocorrpa.mcalcmediaeautocorrpa.m

Filtraggio di un treno di oscillazioni in rumore AWGN

- 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]

Esempio: segnale in una comunicazione digitale, o in sistema radar/sonar/ecografico

- Filtrare una realizzazione del processo aleatorio prima definito mediante un filtro passa basso e poi mediante un filtro passa-banda

Progetto 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]

Filtraggio di un treno di oscillazioni in rumore AWGN

0 50 100 150 200 250 300 350 400 450 500-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Il treno di oscillazioni immerso in una realizzazione di rumore

Esempio di risultati

0 50 100 150 200 250 300 350 400 450 500-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Risultato del filtraggio passa-basso

Esempio di risultati

0 50 100 150 200 250 300 350 400 450 500-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Risultato del filtraggio passa-banda

Esempio di risultati

Esempio di risultati

0 50 100 150 200 250 300 350 400 450 500-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

0 50 100 150 200 250 300 350 400 450 500-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

0 50 100 150 200 250 300 350 400 450 500-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

% 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)])

File.m: genoscrum.mgenoscrum.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 passabanda di 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)

File.m: filtrpbasso.mfiltrpbasso.m