Tema 6: Analisi in Potenza di Processi Parametrici
Esempio: Oscillazione con ampiezza e fasi aleatorie
X (t )  A cos(2 f 0t  )
Tipico segnale presente in sistemi
radiomobile, radar, sonar, etc.:
A R ( )  U (0,2 ) A e  indipendenti
scomposizione del fasore con ampiezza
e fase aleatorie in somma di due fasori
in quadratura con ampiezze aleatorie:
f A (a) 
a

e
2

a2
2 2
u (a ), con E{ A2 }  2 2
X (t )  A cos(2 f0t  )  A cos  cos(2 f 0t )  Asin  sin(2 f 0t )
X (t )  X1 cos(2 f0t )  X 2 sin(2 f0t )
con
 A  X 2  X 2
 X 1  A cos 
1
2
e 

 X 2  A sin 
  atan2( X 2 , X 1 )
È la trasformazione di coordinate polari-cartesiane, per cui si è visto che:
X 1  N (0, 2 )
X 2  N (0, 2 )
X1, X 2 IID
N.B.: E{ A2 }  2 2  2 E{ X i2 }
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:
X [n]  X (nTc )  X 1 cos(2 F0n)  X 2 sin(2 F0n), n  0,1, 2,
, N 1
dove: F0  f 0Tc
Parametri:
Intervallo di campionamento: Tc=1 ms
Intervallo temporale di osservazione: T=NTc=0.1 ms
Valore quadratico medio dell’ampiezza dell’oscillazione: E{ A2 }  2 2  1
Frequenza dell’oscillazione: f0  50 KHz
[istruzioni utili: randn, repmat, for, plot, hold on]
Esempi di realizzazioni
X (t )  A cos(2 f 0t  )
2.5
2
1.5
1
M=50 realizzazioni
0
2.5
-0.5
2
-1
1.5
-1.5
1
-2
-2.5
0.5
0
10
20
30
40
50
t, sec
60
70
80
90
100
x(t)
x(t)
0.5
0
-0.5
M=3 realizzazioni
-1
-1.5
-2
-2.5
0
10
20
30
40
50
t, sec
60
70
80
90
100
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);
x2=sigma*randn(1,n);
t=[1:tfin];
% genera n realizz. della v.a. Gaussiana ampiezza fasore 1;
% genera n realizz. della v.a. Gaussiana ampiezza fasore 2;
% calcola i tempi su cui campionare;
x1=repmat(x1.',1,tfin); % organizza tempi e ampiezze fasori
x2=repmat(x2.',1,tfin); % per calcolo efficiente
t=repmat(t,n,1);
% delle realizzazioni;
x=x1.*cos(2*pi*f0*t)-x2.*sin(2*pi*f0*t);
figure;
for i=1:n;
plot(x(i,:))
hold on
end
xlabel('t, sec')
ylabel('x(t)')
% calcola tutte le realizzazioni;
% grafica le realizzazioni;
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)
Y [n]  X [n]  W [n]  X 1 cos(2 F0n)  X 2 sin(2 F0n)  W [n]
Parametri:
Segnale: come prima!
Rumore: W [n]  N  0, W2 
SNRdB  10, 5,0,5,10,50 dB
E{ X 2 [n]} E{ A2 } 2  2 1 2
SNR 

 2  2
2
2
E{W [n]} E{W [n]}  W  W
SNRdB  10Log10 ( SNR)  SNR  10
SNRdB
10
  
2
W
1

2
10

SNRdB
10
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):
E{| X T (; f ) |2 }
S X ( f )  E{S X (; f )}  lim
T 
T
PSD dei segnali deterministici “realizzazioni”
Per un processo aleatorio stazionario in senso lato si ha inoltre:
Teorema di
Einstein-Wiener-Khintchine:
S X ( f )  FT RX ( )
dove RX ( )  E  X (t ) X (t   ) è la funzione di autocorrelazione (ACF)
Misura sperimentale della PSD
Date M realizzazioni del processo osservate su un intervallo di
osservazione di durata finita T, misurare la PSD:
E{| X T (; f ) |2 } 1
S X ( f )  lim

T 
T
M
| X T(i ) ( f ) |2
i1 T
M
interpretazione in termini
di frequenza relativa
Si utilizza la FFT:
[ realizzazioni in tempo-discreto,
Tc=tempo di campionamento ]
(i )
T
X (f)
spettri di ampiezza
delle realizzazioni
 Tc  n0 xi (nTc ) e
N 1
f
k
N Tc
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
 j 2
kn
N
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]
Esempio di risultati
X (t )  A cos(2 f 0t  )
E{ A2 }  2 2  1
25
f 0  0.05 sec
20
T  100 sec
15
dsp
M=2000 realizzazioni
RX ( )   2 cos(2 f 0 )
10
5
0
-0.5
-0.4
-0.3
-0.2
-0.1
0
f , Hz
0.1
0.2
S X ( f )  FT  2 cos(2 f 0 ) 
0.3
2
2
0.4
Effetto della finestratura:
le  () divengono Sinc2 ()
“discrete”
0.5
 ( f  f0 ) 
2
2
 ( f  f0 )
Esempio di risultati
PSD del processo aleatorio
PSD misurata
da tre realizzazioni
30
25
(la PSD é “phase blind”)
15
30
10
25
5
20
0
-0.2
-0.15
-0.1
-0.05
0
f , Hz
0.05
0.1
0.15
0.2
dsp
dsp
20
15
10
5
fluttuazioni
di ampiezza
da realizzazione a realizzazione
0
-0.2
-0.15
-0.1
-0.05
0
f , Hz
0.05
0.1
0.15
0.2
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);
n=size(x,1);
tc=1;
zeropad=1024;
%
%
%
%
finestra temporale disponibile
numero di realizzazioni disponibili
tempo di campionamento unitario
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 realizzazioni
dsprealizzazioni=abs(spettri).^2/(tfin-1); % calcola le dsp deterministiche delle realizzazioni
dsp=mean(dsprealizzazioni);
% calcola la dsp del processo
dsp=fftshift(dsp);
ncamp=zeropad;
deltaf=1/tc/ncamp;
f=[-ncamp/2:ncamp/2-1]*deltaf;
%
%
%
%
centra la componente a frequenza zero
numero di campioni dopo lo zero-padding
passo di campionamento in frequenza
scala delle frequenze
figure;
plot(f,dsp);
xlabel('f , Hz');
ylabel('dsp');
% grafico dsp del processo
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);
tc=1;
zeropad=1024;
% finestra temporale disponibile
% tempo di campionamento unitario
% zero-padding a 1024 campioni
if 3>size(x,1) disp('non sono disponibili le realizzazioni richieste!')
end
if 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 realizzazioni
dsprealizzazioni=abs(spettri).^2/(tfin-1); % calcola le dsp deterministiche delle 3 realizzazioni
ncamp=zeropad;
deltaf=1/tc/ncamp;
f=[-ncamp/2:ncamp/2-1]*deltaf;
% numero di campioni dopo lo zero-padding
% passo di campionamento in frequenza
% scala delle frequenze
figure;
plot(f,fftshift(dsprealizzazioni(1,:)),'g');
hold on
plot(f,fftshift(dsprealizzazioni(2,:)),'k');
hold on
plot(f,fftshift(dsprealizzazioni(3,:)),'r');
xlabel('f , Hz');
ylabel('dsp');
% grafico dsp delle 3 realizzazioni
Valutazione della funzione valor medio ed ACF
Sono date M realizzazioni del processo, osservate su un
intervallo di durata finita T:
1
 X (t )  E{ X (; t )}  lim
M  M
1
 i1 xi (t )  ˆX (t )  M
M
Interpretazione in termini
di frequenza relativa
1
RX (t , t   )  E{ X (; t ) X (; t   )}  lim
M  M
1
M
ˆ
 RX (t , t   ) 
x (t ) xi (t   )

i 1 i
M

M
x (t )
i 1 i
[ con t nell’intervallo
temporale osservato ]

M
x (t ) xi (t   )
i 1 i
Risultati per il p.a. “oscillazione
con ampiezza di Rayleigh e fase uniforme”
Funzione valor medio valutata con M=2000 realizzazioni
2.5
2
1.5
1
 X (t )  0
0.5
0
-0.5
-1
-1.5
-2
-2.5
0
10
20
30
40
50
t
60
70
80
90
100
Funzione di autocorrelazione (ACF)
valutata con M=2000 realizzazioni
0.5
0.4
ACF calcolata per
0.3
0.2
t1
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5
0
10
20
30
40
50

60
70
80
90
1
RX (t , t   )  cos(2  0.05  )
2
100
Verifica di stazionarietà in autocorrelazione:
1
RX (t , t   )  cos(2  0.05  )
2
0.5
ACF calcolata per t5
0.4
0.3
0.2
0.5
0.1
0.4
0
0.3
0.2
-0.1
0.1
-0.2
0
-0.3
-0.1
-0.4
-0.5
-0.2
0
10
20
30
40
50

60
70
80
ACF calcolata per t1
90
100
-0.3
-0.4
-0.5
0
10
20
30
40
50

60
70
80
90
L’ACF, nei limiti dell’errore di valutazione sperimentale, non dipende da t
100
Interpretazione dei risultati di autocorrelazione:
Scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo 
2.5
2
1.5
Coefficiente di
correlazione tra
X(t) e X(t+)
=
19 (t1
19
1
0.5
0
Covarianza tra le
due v.a. / prodotto
delle deviazioni
standard
-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
=
si nota che le variabili sono correlate positivamente
cX ( ) rX ( )   X 2
 X ( ) 

 cos(2  0.05 19)  0.95
cX (0) rX (0)   X 2
cX()/cX(0)
Interpretazione dei risultati di autocorrelazione:
Scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo 
2.5
2
29
1.5
1
2.5
0.5
1.5
-0.5
1
-1
0.5
-1.5
0
-2
-2.5
-2.5
24
2
0
-0.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
-1
-1.5
le variabili sono correlate negativamente
-2
-2.5
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
le variabili sono quasi incorrelate
2.5
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!
1
0.8
0.6
19 (t1
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
scatterplot (M=2000 realizz.) di coppia di v.a. estratte con ritardo 
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)])
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 s
Durata intervallo di osservazione: 500 s
Tempo di inizio oscillazione: 100 s
Durata temporale oscillazione: 100 s
Ampiezza oscillazione: 1
Frequenza oscillazione: 0.05 Hz
Potenza rumore: 0.1 [il rumore é bianco in banda di Nyquist]
[istruzioni utili: randn]
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-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]
Esempio di risultati
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0
50
100
150
200
250
300
350
400
450
500
Il treno di oscillazioni immerso in una realizzazione di rumore
Esempio di risultati
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0
50
100
150
200
250
300
350
400
450
500
Risultato del filtraggio passa-basso
Esempio di risultati
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0
50
100
150
200
250
300
350
400
450
500
Risultato del filtraggio passa-banda
Esempio di risultati
2
1.5
1
0.5
0
-0.5
-1
2
-1.5
1.5
-2
0
50
100
150
200
250
300
350
400
450
1
500
0.5
0
-0.5
2
-1
1.5
-1.5
1
-2
0
50
100
150
200
250
300
350
400
450
0.5
500
0
-0.5
-1
-1.5
-2
0
50
100
150
200
250
300
350
400
450
500
File.m: genoscrum.m
% generazione realizzazione spezzone di oscillazione immersa in rumore Gaussiano
function [x] = genoscrum()
% out: realizzazione;
%
uscita su video di grafico della realizzazione.
tfin=500;
tin=100;
tdur=100;
a=1;
f0=0.05;
pot=0.1;
%
%
%
%
%
%
ampiezza intervallo temporale (tempo di campionamento unitario);
istante di inizio oscillazione;
durata temporale oscillazione;
ampiezza della oscillazione;
frequenza della oscillazione;
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: filtrpbasso.m
% filtraggio passabasso di spezzone di oscillazione immersa in rumore
function [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 rumore
function [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)