实验5-典型时间序列模型分析

２４实验5典型时间序列模型分析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过程，可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：121()10.30.5Hzzz−−=++这是一个全极点的滤波器，具有无限长的冲激响应。对于功率谱，可以这样得到，221212exp()1()1xwzjPazazωωσ−−==++可以看出，()xPω完全由两个极点位置决定。对于AR模型的自相关函数，有下面的公式：２５211(0)(1)()(1)(0)(1)0()(1)(0)0xxxwxxxpxxxrrrparrrparprprσ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥−⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−⎢⎥⎣⎦⎣⎦⎣⎦#####这称为Yule-Walker方程，当相关长度大于p时，由递推式求出：1()(1)()0xprrararp+−++−=这样，就可以求出理论的AR模型的自相关序列。1.产生样本函数，并画出波形题目中的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随机序列');grid得到的输出序列波形为：050100150200250300350400450500-8-6-4-202468x(n)产生的AR随机序列2.估计均值和方差可以首先计算出理论输出的均值和方差，得到0xm=，对于方差可以先求出理论自相关输出，然后取零点的值。()()*()*()xwrmhmhmrm=−并且，()4()wrmmδ=，带入有()4[()*()]xrmhmhm=−可以采用上面介绍的方法，对式中的卷积进行计算。计算出的卷积输出图形为：２６0510152025303540-1-0.500.511.5nh(n)*h(-n)在最大值处就是输出的功率，也就是方差，为2(0)5.6xxrσ==对实际数据进行估计，均值为mean(x)=-0.0703，而方差为var(x)=5.2795，两者和理论值吻合的比较好。3.画出理论的功率谱密度曲线理论的功率谱为，22()()()4()jjxwPPHeHeωωωω==用下面的语句产生：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))).^2);%计算出理论值Gx=Gx/max(Gx);%归一化处理f=w*Fs/(2*pi);%转化到模拟域上的频率plot(f,Gx,’b’),gridon;得到的图形为：-500-400-300-200-1000100200300400500-12-10-8-6-4-20f(Hz)Gx(w)dB可以看出，这个系统是带通系统。4.估计自相关函数和功率谱密度用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出最后的仿真２７图形。%计算理论和实际的自相关函数序列Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,'coeff');m=-Mlag:Mlag;stem(m,Rx,'r.');最终的值为-20-15-10-505101520-0.500.51mRx(m)可以看出，它和上面的理论输出值吻合程度很好。实际的功率谱密度可以用类似于上面的方法进行估计，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(2:end)];%构造一个对称的频率，范围是[-Fs/2,Fs/2]Py=[-fliplr(Py)Py(2:end)];%对称的功率谱plot(f,10*log10(Py),’b’);估计出来的功率谱密度为，-500-400-300-200-1000100200300400500-14-12-10-8-6-4-20f(Hz)Gx(w)dB将两幅图画在一起，可以看到拟合的情况比较好：２８-500-400-300-200-1000100200300400500-12-10-8-6-4-20f(Hz)Gx(w)dB理论值实际值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过程，也可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：121210.50.2()10.30.2zzHzzz−−−−+−=+−对于功率谱，可以这样得到，22exp()()()xwzjPHzωωσ==对于ARMA过程，当模型的所有极点均落在单位圆内时，才是一个渐进平稳的随机过程。这个过程的自相关函数不能简单地写成Yule-Walker方程形式，它于模型的参数具有高度的非线性关系。1.产生样本函数，并画出波形题目中的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’);得到的输出序列波形为：050100150200250300350400450500-6-4-20246nX(n)2.估计均值和方差可以首先计算出理论输出的均值和方差，得到0xm=，对于方差可以先求出理论自相关输出，然后取零点的值。()()*()*()xwrmhmhmrm=−并且，()4()wrmmδ=，带入有()4[()*()]xrmhmhm=−可以采用上面介绍的方法，对式中的卷积进行计算。计算出的卷积输出图形为：0510152025303540-0.200.20.40.60.811.2nh(n)*h(-n)３０-500-400-300-200-10001002003004005000.20.30.40.50.60.70.80.91f(Hz)Gx(w)dB在最大值处就是输出的功率，也就是方差，为2(0)4.2xxrσ==对实际数据进行估计，均值为mean(x)=-0.0547，而方差为var(x)=3.8，两者和理论值吻合的比较好。3.画出理论的功率谱密度曲线理论的功率谱为，22()()()4()jjxwPPHeHeωωωω==用下面的语句产生：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).^2);%计算出理论值Gx=Gx/max(Gx);f=w*Fs/(2*pi);%转化到模拟域上的频率plot(f,Gx,’b’),gridon;4.估计相关函数和功率谱密度曲线用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出仿真图形。%计算理论和实际的自相关函数序列Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,'coeff');m=-Mlag:Mlag;stem(m,Rx,'r.');最终的值为３１-20-15-10-505101520-0.200.20.40.60.811.2nRx(n)实际的功率谱密度可以用类似于上面的方法进行估计，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(2:end)];%构造一个对称的频率，范围是[-Fs/2,Fs/2]Py=[fliplr(Py)Py(2:end)];%对称的功率谱plot(f,10*log10(Py),’b’);估计出来的功率谱密度为，-500-400-300-200-1000100200300400500-6-5-4-3-2-10f(Hz)Px(w)dB把两幅图画在一起，可以得到下面的图形，可以看出两者的吻合度比较高。-500-400-300-200-1000100200300400500-6-5-4-3-2-10f(Hz)Gx(w)dB理论值真实值３２3、实验内容1、熟悉实验原理，将实验原理上的程序应用matlab工具实现；2、设有MA(2)模型，x(n)=W(n)-0.3W(n-1)+0.2W(n-2)W(n)是零均值正态白噪声，方差为4。（1）用MATLAB模拟产生X(n)的500观测点的样本函数，并绘出波形（2）用产生的500个观测点估计X(n)的均值和方差（3）画出理论的功率谱（4）估计X(n)的相关函数和功率谱4、实验要求（1）用MATLAB编写程序。（2）写出详细试验报告(要有自己对实验结果的结论)。