实验1典型时间序列模型分析1、实验目的熟悉三种典型的时间序列模型：AR 模型，MA 模型与ARMA 模型，学会运用Matlab 工具对对上述三种模型进行统计特性分析，通过对2阶模型的仿真分析，探讨几种模型的适用范围， 并且通过实验分析理论分析与实验结果之间的差异。

2、实验原理AR 模型分析：设有AR(2)模型,X( n)=-0.3X( n-1)-0.5X( n-2)+W( n)其中：W(n)是零均值正态白噪声，方差为4。

(1 )用MATLAB 模拟产生X(n)的500观测点的样本函数，并绘出波形 (2) 用产生的500个观测点估计X(n)的均值和方差 (3) 画出理论的功率谱(4) 估计X(n)的相关函数和功率谱【分析】给定二阶的 AR 过程，可以用递推公式得出最终的输出序列。

或者按照一个白噪声 通过线性系统的方式得到，这个系统的传递函数为：这是一个全极点的滤波器，具有无限长的冲激响应。

对于功率谱，可以这样得到，可以看出，FX w完全由两个极点位置决定。

对于AR 模型的自相关函数，有下面的公式:\(0)打⑴匚⑴…^(0)■ 1'G 2W 0JAP) 人9-1)…凉0) _这称为Yule-Walker 方程，当相关长度大于p 时，由递推式求出:r (r) +-1) + -■ + (7r - JJ )= 0这样，就可以求出理论的AR 模型的自相关序列。

H(z)二1 1 0.3z ，P x w +W1 1 a 才 a 2z^1. 产生样本函数，并画出波形2. 题目中的AR过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。

clear all;b=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数h=impz(b,a,20); % 得到系统的单位冲激函数，在20点处已经可以认为值是0randn('state',0);w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2阶AR过程plot(x,'r');ylabel('x(n)');title(' 邹先雄——产生的AR随机序列');grid on;得到的输出序列波形为：邹先雄——产生的AR随机序列2. 估计均值和方差可以首先计算出理论输出的均值和方差，得到m x =0，对于方差可以先求出理论自相关输出，然后取零点的值。

并且， 」，带入有在最大值处输出的功率，也就是方差，为a; =r (0) = 56对实际数据进行估计，均值为 mean(x)=-0.0703，而方差为var(x)=5.2795， 合得比较好。

程序及运行结果图如下，其中y_mean 表示均值，y_var 表示方差。

>> clear all;b=Lll; a=[l 0.3 0.B]; %由摒迷的差分方程，僚到索猊倍谨硒数 l^i^z(b jaj 20); %得到系统的单f 立冲數函数，在20点雉已经可以认为值是0 r andnC st ate-J , 0).2, 1,500); %产生题设的白囁声随机序別，标准差为2 x=filt e r(b jaj w); %通过线形紊统，课到输出就是题目中要卡的2盼AR 过程 plot 甌，r ，); yl abel C K Cn));生的AR 龍机序貝N ;grid an;yl _jn&an=ine3n (x) y2_var=var (x) yl_jnean =-0.0703 y2_var = 5.27953.画出理论的功率谱密度曲线理论的功率谱为，£92恥训丹冲)f "|H (严)「用下面的语句产生：delta=2*pi/1000; w_min=-pi; w_max=pi; Fs=1000;w=w_min:delta:w_max; %得到数字域上的频率取样点，范围是[-pi,pi]Gx=4*(abs(1丿(1+0.3*exp(-i*w)+0.5*exp(-2*i*w)))A2); % 计算出理论值Gx=Gx/max(Gx); % 归一化处理f=w*Fs/(2*pi); % 转化到模拟域上的频率plot(f,Gx); title(' 邹先雄一一理论功率谱密度曲线');grid on;得到的图形为：两者合理论值吻邹先雄一理论功率谱密度曲线可以看出，这个系统是带通系统。

4•估计自相关函数和功率谱密度用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出最后的仿真图形。

Mlag=20; % 定义最大自相关长度Rx=xcorr(x,Mlag,'coeff');m=-Mlag:Mlag;stem(m,Rx,'r.');title(' 邹先雄-------- 自相关函数’)；最终的值为邹先雄—自相关函数可以看出，它和上面的理论输出值吻合程度很好。

实际的功率谱密度可以用类似于上面的方法进行估计，window=hamming(20); % 采用hanmming 窗，长度为20noverlap=10; % 重叠的点数Nfft=512; % 做FFT 的点数Fs=1000; % 采样频率，为1000Hzb=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数h=impz(b,a,20); % 得到系统的单位冲激函数，在20点处已经可以认为值是0randn('state',0);w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2阶AR过程[Px,f]=pwelch(x,window,noverlap,Nfft,Fs, 'onesided'); % 估计功率谱密度f=[-fliplr(f) f(1:end)]; % 构造一个对称的频率，范围是[-Fs/2, Fs/2]Py=[-fliplr(Px) Px(1:end)]; % 对称的功率谱plot(f,10*log10(Py),'b');title(' 邹先雄一一实际的功率谱密度曲线');估计出来的功率谱密度为，邹先雄一实际的功率谱密度曲线将两幅图画在一起，可以看到拟合的情况比较好(两者相位刚好相反，但是基本波形相似)代码如下：clear all;delta=2*pi/1000;w_mi n=-pi;w_max=pi;Fs=1000;w=w_mi n:delta:w_max; % 得到数字域上的频率取样点，范围是[-pi,pi]Gx=4*(abs(1丿(1+0.3*exp(-i*w)+0.5*exp(-2*i*w)))A2); % 计算出理论值Gx=Gx/max(Gx); % 归一化处理f=w*Fs/(2*pi); % 转化到模拟域上的频率结束plot(f,Gx,'r');hold on;title(' 邹先雄一一理论和实际的功率谱密度曲线拟合');window=hamming(20); % 采用hanmming 窗，长度为20noverlap=10; %重叠的点数Nfft=512; % 做FFT 的点数Fs=1000; % 采样频率，为1000Hzb=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数h=impz(b,a,20); % 得到系统的单位冲激函数，在20点处已经可以认为值是0randn ('state',0);w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2阶AR过程[Px,f]=pwelch(x,wi ndow, no verlap,Nfft,Fs, 'on esided'); % 估计功率谱密度f=[-fliplr(f) f(1:e nd)]; % 构造一个对称的频率，范围是[-Fs/2, Fs/2] Py=[-fliplr(Px) Px(1:e nd)]; % 对称的功率谱Py=-10*log10(Py);Py=Py/max(Py);Py=-Py;Py=3*Py;Py=Py+2.6；% 用来归一处理，使两者吻合plot(f,Py,'b');legend('实际值','理论值');grid on;ARMA模型分析设有ARMA(2,2)模型，X( n)+0.3X( n-1)-0.2X( n-2)=W( n)+0.5W( n-1)-0.2W( n-2)W(n)是零均值正态白噪声，方差为4。

(1 )用MATLAB模拟产生X(n)的500观测点的样本函数，并绘出波形(2)用产生的500个观测点估计X(n)的均值和方差(3)画出理论的功率谱(4)估计X(n)的相关函数和功率谱【分析】给定(2,2)的ARMA过程，也可以用递推公式得出最终的输出序列。

或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：I 1+0円厂1+03^-02尹对于功率谱，可以这样得到,对于ARMA过程，当模型的所有极点均落在单位圆内时，才是一个渐进平稳的随机过程。

这个过程的自相关函数不能简单地写成Yule-Walker方程形式，它于模型的参数具有高度的非线性关系。

1. 产生样本函数，并画出波形题目中的ARMA过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。

clear all;b=[1 0.5 -0.2]; a=[1 0.3 -0.2]; % 由描述的差分方程，得到系统传递函数h=impz(b,a,10); % 得到系统的单位冲激函数，在10点处已经可以认为值是0randn ('state',0);w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的(2,2)阶ARMA过程plot(x,'r');title(' 邹先雄一一输出的AR随机序列');得到的输出序列波形为：邹先雄一输出的AR随机序列2. 估计均值和方差可以首先计算出理论输出的均值和方差，得到m x =° ，对于方差可以先求出理论自相关输出，然后取零点的值。