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
Top Related