时间序列ARMA模型(shiyan2c.m)

3.设有ARMA（2，2）模型：

W(n)是零均值正态白噪声，方差为4。

（1）用MATLAB模拟产生X(n)的500观测点的样本函数，并绘出图形；

（2）用产生的500个观测点估计X(n)的均值和方差；

（3）画出X(n)的理论的功率谱；

（4）估计X(n)的自相关函数和功率谱。

%%AR(2.2)

clear all;

b=[1 0.5 -0.2];a=[1 0.3 -0.2];%ÓÉÃèÊöµÄ²î·Ö·½³Ì£¬µÃµ½ÏµÍ³´«µÝº¯Êý

h=impz(b,a,20);%µÃµ½ÏµÍ³µÄµ¥Î»³å¼¤º¯Êý£¬ÔÚ20µã´¦ÒѾ­¿ÉÒÔÈÏΪֵÊÇ0

randn('state',0);

w=normrnd(0,2,1,500);%²úÉúÌâÉèµÄ°×ÔëÉùËæ»úÐòÁУ¬±ê×¼²îΪ2

x=filter(b,a,w);%ͨ¹ýÏßÐÎϵͳ£¬µÃµ½Êä³ö¾ÍÊÇÌâÄ¿ÖÐÒªÇóµÄ2½×AR¹ý³Ì

subplot(3,2,[1,2])

plot(x,'r');ylabel('x(n)');title('²úÉúµÄARËæ»úÐòÁÐ');grid

%%

mean=mean(x)

var=var(x)

%%

R=xcorr(x,'coeff');

subplot(3,2,3)

plot(R)

title('ÀíÂÛ×ÔÏà¹Øº¯Êý');

%%

delta=2*pi/1000;

w_min=-pi;

w_max=pi;

Fs=1000;

w=w_min:delta:w_max;

NS=1 0.5*exp(-i*w)-0.2*exp(-2*i*w);

ND=1 0.3*exp(-i*w)-0.2*exp(-2*i*w);

Gx=4*(abs(NS./ND).^2);

%Gx=4*(abs(1./(1 0.3*exp(-i*w) 0.5*exp(-2*i*w))).^2);

Gx=Gx/max(Gx);

f=w*Fs/(2*pi);

subplot(3,2,4)

plot(f,Gx);

title('ÀíÂÛ¹¦ÂÊÆ×');

%%

Mlag=20;

Rx=xcorr(x,Mlag,'coeff');

m=-Mlag:Mlag;

subplot(3,2,5)

stem(m,Rx,'r');

title('¹À¼ÆX(n)µÄ×ÔÏà¹Øº¯Êý');

%%

window=hamming(20);

   noverlap=10;

   Nfft=512;

   Fs=1000;

   b=[1 0.5 -0.2];a=[1 0.3 -0.2];

   h=impz(b,a,10);

   randn('state',0);

   w=normrnd(0,2,1,500);

x=filter(b,a,w);

[Px,f]=pwelch(x,window,noverlap,Nfft,Fs,'onesided');

f=[-fliplr(f) f(1:end)];

Py=[-fliplr(Px) Px(1:end)];

subplot(3,2,6)

plot(f,10*log10(Py),'b');

title('¹À¼ÆX(n)µÄ¹¦ÂÊÆ×');