[ ]h n[ ]W n [ ] [ 1] [ ]X n X n W n

Processi Autoregressivi AR(1)

AR(1): eq. alle differenze ricorsiva di

ordine 1

MATLAB:

w=randn(N,1)

x=filter(b,a,w) 22

1

1j f

j fH e

e

White

Gaussian Noise

(WGN)

Filtro IIR di ordine 1

[ ] [ ]nh n u n

ACF e PSD di un processo AR(1)

2 22 2

2 2

1( )

1 1 2 cos 2j f W

X W j fS e

e f

22

21W

X

2| |

2[ ]

1mW

XR m

f

passa-basso2( )j fXS e 0 :

0.5

0.1

0.2

N=64 N=1024

Periodogramma del processo di ingresso WGN

22 1ˆ ( ) [ ] ,02

[ ] [ ]

[ ] [ ] [ ], 02

j fkW

fN

NS e W k k

N

W k FFT w n

NW k W N k W k k

2 2( )j fW WS e

- Generare e graficare tre realizzazioni di un p.a. AR(1) Gaussiano.

Parametri:

Intervallo temporale discreto di osservazione: N=128

Potenza del processo:

Parametro AR(1):

[istruzioni utili: randn, for, filter, subplot, stem]

- Rappresentare i grafici della funzione di autocorrelazione e della densità spettrale di potenza.

Generazione processo AR(1)

0.1, 0.5, 0.9 e 0.9

2 1X

Esempio di risultati

=0.1

Esempio di risultati

=0.5

Esempio di risultati

=0.9

Esempio di risultati

=-0.9

Confronto …..

=0.1 =0.9

Funzione di autocorrelazione (1/2)

= 0.5

0 5 10 15 20 25 30 35 40 45 50

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

m

Rx(

m)

Funzione di autocorrelazione, a=0.5

Funzione di autocorrelazione (2/2)

= - 0.9= 0.9

0 20 40 60 80 100 120

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

m

Rx(

m)

Funzione di autocorrelazione, a=0.9

0 20 40 60 80 100 120

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

m

Rx(

m)

Funzione di autocorrelazione, a=-0.9

Densità spettrale di potenza (1/2)

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

0.5

1

1.5

2

2.5

frequenza

Sx(

f)Densità Spettrale di Potenza

a=0.1

a=0.2a=0.5

Densità spettrale di potenza (2/2)

= - 0.9= 0.9

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

2

4

6

8

10

12

14

16

18

frequenza

Sx(

f)

Densità Spettrale di Potenza

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

2

4

6

8

10

12

14

16

18

frequenza

Sx(

f)

Densità Spettrale di Potenza

Esempio di risultati (1/2)

% generazione realizzazioni p.a. AR(1) Gaussiano di potenza unitaria% function [x] = genproar1(Nc,nr,a)

% in: numero di campioni, Nc;% numero di realizzazioni, nr;% coefficiente di reazione, a;

% out: matrice di realizzazioni, una per riga;% uscita su video di grafico delle nr realizzazioni.

sigmaw=sqrt(1-a^2); % calcola la deviazione standard per il processo di forzamento;w=sigmaw*randn(nr,Nc); % genera nr realiz. di durata Nc del processo di forzamento;

for i=1:nr;x(i,:)=filter(1,[1 -a],w(i,:),0); % filtraggio IIR del primo ordine di ogni realiz.;end

figure; % grafica le nr realizzazioni del p.a. AR(1);for i=1:nr; subplot(nr,1,i) stem(x(i,:)) xlabel(‘n’) ylabel(‘x(n)’) axis([0 Nc -2.5 2.5])end

Esempio di risultati (2/2)

% Funzione di autocorrelazionem=1:Nc;sigmax=sqrt(sigmaw^2/(1-a^2));Rx=sigmax^2*a.^m;figure(2);stem(Rx);axis([0 Nc -sigmax^2 sigmax^2])grid onxlabel('m')ylabel ('Rx(m)')title('Funzione di autocorrelazione')

% Densità spettrale di potenzaaa=[0.1 0.2 0.5]for n=1:3 Nf=1024; Df=1/Nf; f=-0.5:Df:0.5-Df; Sx=sigmaw^2./(1+aa(n)^2-2*aa(n)*cos(2*pi*f)); plot(f,Sx); hold on axis([-0.5 0.5 0 max(Sx)]) grid on xlabel('frequenza') ylabel ('Sx(f)') title('Densità Spettrale di Potenza')end;

- Generare N campioni di una realizzazione di un p.a. AR(1) Gaussiano,

calcolare e graficare il periodogramma “mediato” su M segmenti dati,

e confrontarlo con la PSD teorica

- Parametri:

Intervallo temporale discreto di osservazione: N=2048

Potenza del processo:

Parametro AR(1):

Numero di segmenti: M=16 (ciascuno di L=N/M=128 campioni)

[istruzioni utili: randn, filter, fft, fftshift, abs, mean]

Periodogramma mediato di processo AR(1)

0.5

2 1X

Periodogramma: esempio di risultati (1/3)

= 0.5

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-40

-35

-30

-25

-20

-15

-10

-5

0

5

10

frequenza

DS

P s

timat

a

periodogramma

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-22

-20

-18

-16

-14

-12

-10

-8

-6

frequenza

DS

P s

timat

a

periodogramma mediato su L blocchi

Periodogramma: esempio di risultati (2/3)

= 0.9

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-50

-40

-30

-20

-10

0

10

20

frequenza

DS

P s

timat

a

periodogramma

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-30

-25

-20

-15

-10

-5

0

5

frequenza

DS

P s

timat

a

periodogramma mediato su L blocchi

Periodogramma: esempio di risultati (3/3)

% Stima spettrale AR(1)clear allclose allclc N=2048;M=16;L=N/M;sigmax=1;a=0.5;sigmaw=sigmax*sqrt(1-a^2);w=sigmaw*randn(1,N);x=filter(1,[1 -a],w,0); Nfft=1024;Df=1/Nfft;f=[-0.5:Df:0.5-Df]; % periodogrammaTDFx=fftshift(fft(x,Nfft));modx=abs(TDFx);modx2=modx.^2;period=modx2./N;figure(1);plot(f,10*Log10(period));xlabel('frequenza');ylabel('DSP stimata');grid on;title('periodogramma');

% periodogramma mediatofor i=1:M TDFxmed(i,:)=fftshift(fft(x(((i-1)*L+1):i*L),Nfft)); modxmed=abs(TDFxmed); modx2med=modxmed.^2; period_med=modx2med./N;end;period_fin=mean(period_med,1);figure(2);plot(f,10*Log10(period_fin));xlabel('frequenza');ylabel('DSP stimata');grid on;title('periodogramma mediato su L blocchi');