典型时间序列模型分析

实验1典型时间序列模型分析1、实验目的熟悉三种典型的时间序列模型：AR模型，MA模型与ARMA模型，学会运用Matlab工具对对上述三种模型进彳T统计特性分析，通过对2阶模型的仿真分析，探讨几种模型的适用范围，并且通过实验分析理论分析与实验结果之间的差异。2、实验原理AR模型分析：设有AR(2)模型,X(n)=-0.3X(n-1)-0.5X(n-2)+W(n)其中：W⑻是零均值正态白噪声，方差为4。(1)用MATLAB模拟产生X(n)的500观测点的样本函数，并绘出波形(2)用产生的500个观测点估计X(n)的均值和方差(3)画出理论的功率谱(4)估计X(n)的相关函数和功率谱【分析】给定二阶的AR过程，可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：H(z)=H(z)=110.3z,0.5z?七0)匚⑴O七0)匚⑴O73-1)这称为Yule-Walker方程，当相关长度大于这是一个全极点的滤波器，具有无限长的冲激响应。对于功率谱，可以这样得到，1121alza2z可以看出，巳(亚)完全由两个极点位置决定。对于AR模型的自相关函数，有下面的公式：r(r)+-1)+--+(7r-JJ)=0这样，就可以求出理论的AR模型的自相关序列。1,产生样本函数，弁画出波形2.题目中的AR过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。clearall;b=[1];a=[10.30.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随机序列’)；gridon;得到的输出序列波形为：邹先雄——产生的AR随机序列o1111n1111nI*0501001602002503003504004505002.估计均值和方差可以首先计算出理论输出的均值和方差，得到mx=0,对于方差可以先求出理论自相关输出，然后取零点的值。

并且，9如)=4双用)，带入有rx(m)=在最大值处输出的功率，也就是方差，为<7；=r(0>=$6mean(x)=-0.0703,而方差为mean(x)=-0.0703,而方差为var(x)=5.2795,y_mean表示均值，y_var表示方差。两者合理论值吻程序及运行结果图如下，其中>>clearall;b=LH;a=[l0.3O.B];%由描逑的差分方程，得到系缢传递函数反imp工通叫2。)；%得到系统的单位冲数函数，在前点处已经可以认为值是口randn('stat&!1,0),^nomrndCO,^1,500);%产生题设的白噪声随机序列，标;隹差为2x=filter(bjajw);%通过建形系缢，筹到输出就是题目中要求的2阶AR过程plot(xJ;ylabelktn)");十这工“'产生的AR随机序列」；gridan;yl-mean二mean(x)y2_var=var(k)yl_jnean=-0.0703y2_var=6.2795»l.画出理论的功率谱密度曲线理论的功率谱为，用下面的语句产生：delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;

邹先雄一理论功率谱密度曲线w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是邹先雄一理论功率谱密度曲线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('邹先雄一一理论功率谱密度曲线');gridon;得到的图形为10.90.80.70.60.50.40.30.20.1$0。400-300-200-1000100200300400500可以看出，这个系统是带通系统。.估计自相关函数和功率谱密度用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出最后的仿真图形。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=[10.30.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('邹先雄一一实际的功率谱密度曲线’);估计出来的功率谱密度为，将两幅图画在一起，可以看到拟合的情况比较好(两者相位刚好相反，但是基本波形相似)代码如下：clearall;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,'r');holdon;title('邹先雄一一理论和实际的功率谱密度曲线拟合,);window=hamming(20);%采用hanmming窗，长度为20noverlap=10;%重叠的点数Nfft=512;%做FFT的点数Fs=1000;%采样频率，为1000Hzb=[1];a=[10.30.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)];%对称的功率谱Py=-10*log10(Py);Py=Py/max(Py);Py=-Py;Py=3*Py;Py=Py+2.6;%用来归一处理，使两者吻合plot(f,Py,'b');legend('实际值’,‘理论值');gridon;邹先雄——理论和实际的功率谱密度曲线拟合-02■Q4IIL」JIIL」I--500400-300-200-1000100200300400500ARMA模型分析设有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过程，也可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：目l+o&T-oz'产一1+0一3二1-0后对于功率谱，可以这样得到,对于ARMA过程，当模型的所有极点均落在单位圆内时，才是一个渐进平稳的随机过程。这个过程的自相关函数不能简单地写成Yule-Walker方程形式，它于模型的参数具有高度的非线性关系。.产生样本函数，弁画出波形题目中的ARMA过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。clearall;b=[10.5-0.2];a=[10.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随机序列.估计均值和方差mx=0mx=0,对于方差可以先求出理论自相/(附)=并且，工⑺…⑺，带入有匕(削)=4{/r(ffl)*/f(-ffl)]在最大值处就是输出的功率，也就是方差，为大/所41对实际数据进行彳t计，均值为mean(x)=-0.0547,而方差为var(x)=3.8,两者和理论值吻合的比较好。附代码及运行结果截图如下：»clearall;b=[l0.5-0.2];a=[l0,3-0.2];%由描述的差分方程，得到系统传递函数h=iJnpz(bJaJ10);%得到系统的单位冲激函数，在1。点处已经可以认为值是0randn(statej0);^normrndCOj2jL500);%产生题设的白噪声随机序列，标准差为2x^ilterCb^.w);K通过线形系统，得到输出就是题目中要求的⑵幻阶皿过程plot⑸'/),Py_mean=mean(x)Py_var=var(x)Pygmean=-0,1488Py^var=3.7986»3.画出理论的功率谱密度曲线理论的功率谱为，以的=尸.3)凶(*)(二4|单/)『用下面的语句产生：delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是[-pi,pi]NS=1+0.5*exp(-i*w)-0.2*exp(-2*i*w);%分子

DS=1+0.3*exp(-i*w)-0.2*exp(-2*i*w);%分母Gx=4*(abs(NS./DS).A2);%计算出理论值Gx=Gx/max(Gx);f=w*Fs/(2*pi);%转化到模拟域上的频率plot(f,Gx,'b');title('邹先雄一一理论的功率谱密度曲线');gridon;OO2^a邹先雄一理论的功率谱密度曲线OO2^a邹先雄一理论的功率谱密度曲线10.90.80.70.60.50.40.3OO1三O504.估计相关函数和功率谱密度曲线4.估计相关函数和功率谱密度曲线用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出仿真图形。%计算理论和实际的自相关函数序列Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,'coeff);m=-Mlag:Mlag;stem(m,Rx,'r.');title('邹先雄一一估计自相关函数,);最终的值为030.60.40.2邹先雄一估计自相关函数030.60.40.2邹先雄一估计自相关函数JTill-02L1JL0510162005101620实际的功率谱密度可以用类似于上面的方法进行估计,window=hamming(20);%采用hanmming窗，长度为20noverlap=10;%重叠的点数Nfft=512;%做FFT的点数Fs=1000;%采样频率，为1000Hzb=[10.5-0.2];a=[10.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过程[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('邹先雄一一实际的功率谱密度曲线');估计出来的功率谱密度为邹先雄一实际的功率谱密度曲线*28--々9l।।j।।।।।-500400-300-200-1000100200300400500把两幅图画在一起，可以得到下面的图形，可以看出两者的吻合度比较高。delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是[-pi,pi]NS=1+0.5*exp(-i*w)-0.2*exp(-2*i*w);%分子DS=1+0.3*exp(-i*w)-0.2*exp(-2*i*w);%分母Gx=4*(abs(NS./DS).A2);%计算出理论值Gx=Gx/max(Gx);f=w*Fs/(2*pi);%转化到模拟域上的频率plot(f,Gx,'r');title('邹先雄一一理论和实际的功率谱密度曲线的拟合,);holdon;window=hamming(20);%采用hanmming窗，长度为20noverlap=10;%重叠的点数Nfft=512;%做FFT的点数Fs=1000;%采样频率，为1000Hzb=[10.5-0.2];a=[10.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过程[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)];%对称的功率谱Py=10*log10(Py);Py=Py/max(Py);

111111111111||111||产q■■■■\Aj——卖标值..理论值一——/\_LK-\।…一.DJ..JL」!\-一,!一…00]T＞=1==libih|iPF1hHhbibII卜==一一="==卜二一hkHHV1a-一hhhHgBHNHh-T1邹先雄——理论和实际的功率诣密度曲线的拟合legend('实际值’,'理论值');gridon;1.210.806040.20-0.2-0.4-600-400-300-200-10001002003004005003、实验内容1、熟悉实验原理，将实验原理上的程序应用matlab工具实现；2、设有MA(2)模型，x(n)=W(n)-03W(n-1)+0.2W(n-2)W(n)是零均值正态白噪声，方差为4。(1)用MATLAB模拟产生X(n)的500观测点的样本函数，并绘出波形(2)用产生的500个观测点估计X(n)的均值和方差(3)画出理论的功率谱(4)估计X(n)的相关函数和功率谱完成4个问题的源代码如下clearall;%产生样本函数，并画出波形b=[1-0.30.2];a=[1];%由描述的差分方程，得到系统传递函数h=impz(b,a,10);%得到系统的单位冲激函数，在10点处已经可以认为值是0randn('state',0);w=normrnd(0,2,1,500);%产生题设的白噪声随机序列，标准差为x=filter(b,a,w);%通过线形系统，得到输出就是题目中要求的(2,2)阶ARMA过程figure(1);plot(x,'r');title('-先雄样本函数');Py_mean=mean(x)Py_var=var(x)%画出理论的功率谱密度曲线delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是[-pi,pi]NS=1-0.3*exp(-i*w)+0.2*exp(-2*i*w);%分子DS=1;%分母Gx=4*(abs(NS./DS).A2);%计算出理论值Gx=Gx/max(Gx);f=w*Fs/(2*pi);%转化到模拟域上的频率figure(2);plot(f,Gx,'b');title('邹先雄一一理论的功率谱密度曲线,);%估计相关函数Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,'coeff);m=-Mlag:Mlag;figure(3);stem(m,Rx,'r.');title('邹先雄一一估at相关函数,);%画出估计的功率谱密度曲线window=hamming(20);%采用hanmming窗，长度为20noverlap=10;%重叠的点数Nfft=512;%做FFT的点数Fs=1000;%采样频率，为1000Hz[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)];%对称的功率谱figure(4);plot(f,10*log10(Py),'b');title('邹先雄一一估计的功率谱密度曲线');%对实际和估计两功率谱密度曲线进行拟合delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是[-pi,pi]NS=1-0.3*exp(-i*w)+0.2*exp(-2*i*w);%分子DS=1;%分母Gx=4*(abs(NS./DS).A2);%计算出理论值Gx=Gx/max(Gx);f=w*Fs/(2*pi);%转化到模拟域上的频率figure(5);plot(f,Gx,'r');title('邹先雄一一实际和估计两功率谱密度曲线的拟合');holdon;window=hamming(20);%采用hanmming窗，长度为20noverlap=10;%重叠的点数Nfft=512;%做FFT的点数Fs=1000;%采样频率，为1000Hzb=[1-0.30.2];a=[1];%由描述的差分方程，得到系统传递函数h=impz(b,a,10);%得到系统的单位冲激函数，在10点处已经可以认为值是0randn('state',0);w=normrnd(0,2,1,500);%产生题设的白噪声随机序列，标准差为2x=filter(b,a,w);%通过线形系统，得到输出就是题目中要求的(2,2)阶ARMA过程[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)];%对称的功率谱Py=10*log10(Py);Py=Py/max(Py);Py=-Py;Py=3*Py;Py=Py+4;%用来归一处理，使两者吻合plot(f,Py,'b');legend('实际值’,'理论值');gridon;样本函数波形为：邹先雄^一^样本函数理论功率谱密度曲线为:估计相关函数波形为:邹先雄一估计相关函数040.2-02O4111L-20-15dO-5101520估计功率谱密度曲线为:实际和估计两功率谱密度曲线的拟合截图如下:19邹先雄——实际和估计两功率谱密度曲线的拟合1080.60.4020r：1实C理j市值电值z—1VaIV—\/\/!iiII1i■u.z-5()04()0-300200-100()100200300400500

附