实验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过程，可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：
H(z)=
1
1+0.3z−1+0.5z−2
这是一个全极点的滤波器，具有无限长的冲激响应。 对于功率谱，可以这样得到，
2
Px(ω)=σw
11+a1z−1+a2z−2
2
z=exp(jω)
可以看出，Px(ω)完全由两个极点位置决定。 对于AR模型的自相关函数，有下面的公式：
２４
2
rx(1)\"rx(p)⎤⎡1⎤⎡σw⎤⎡rx(0)
⎢r(1)⎥⎢a⎥⎢⎥−(0)(1)rrpxxx⎢⎥⎢1⎥=⎢0⎥
⎢#⎥⎢#⎥⎢#⎥##⎢⎥⎢⎥⎢⎥()(1)(0)rprpr\"−⎢ap⎦⎥⎣0⎦xx⎣x⎦⎣
这称为Yule-Walker方程，当相关长度大于p时，由递推式求出：
rx(r)+a1(r−1)+\"+ap(r−p)=0
这样，就可以求出理论的AR模型的自相关序列。 1. 产生样本函数，并画出波形
题目中的AR过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。
clear all;
b=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数
h=impz(b,a,20); % 得到系统的单位冲激函数，在20点处已经可以认为值是0 randn('state',0);
w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为2 x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2阶AR过程 plot(x,'r');ylabel('x(n)');title('产生的AR随机序列');grid
得到的输出序列波形为：
产生的AR随机序列82)n(x0-2-4-6-8050100150200250300350400450500
2.
估计均值和方差
可以首先计算出理论输出的均值和方差，得到mx=0，对于方差可以先求出理论自相关输出，然后取零点的值。
rx(m)=h(m)*h(−m)*rw(m)
并且，rw(m)=4δ(m)，带入有
rx(m)=4[h(m)*h(−m)]
可以采用上面介绍的方法，对式中的卷积进行计算。计算出的卷积输出图形为：
２５
1.51)n-(h*)n(h0.50-0.5-105101520n25303540
在最大值处就是输出的功率，也就是方差，为
σx2=rx(0)=5.6
对实际数据进行估计，均值为mean(x)= -0.0703，而方差为var(x)= 5.2795，两者和理论值吻合的比较好。
画出理论的功率谱密度曲线 理论的功率谱为，
3.
Px(ω)=Pw(ω)H(ejω)=4H(ejω)
用下面的语句产生：
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’),grid on;
22
得到的图形为：
0-2-4Bd )w(xG-6-8-10-12-500-400-300-200-1000f (Hz)100200300400500
4.
可以看出，这个系统是带通系统。 估计自相关函数和功率谱密度
用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出最后的仿真
２６
图形。
% 计算理论和实际的自相关函数序列 Mlag=20; % 定义最大自相关长度 Rx=xcorr(x,Mlag,'coeff'); m=-Mlag:Mlag; stem(m,Rx,'r.');
最终的值为
10.5)m(xR0-0.5-20-15-10-50m5101520
可以看出，它和上面的理论输出值吻合程度很好。 实际的功率谱密度可以用类似于上面的方法进行估计，
window=hamming(20); % 采用hanmming窗，长度为20 noverlap=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’);
估计出来的功率谱密度为，
0-2-4Bd )w(xG-6-8-10-12-14-500-400-300-200-1000f (Hz)100200300400500
将两幅图画在一起，可以看到拟合的情况比较好：
２７
0理论值实际值-2-4Bd )w(xG-6-8-10-12-500-400-300-200-1000f (Hz)100200300400500
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过程，也可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：
1+0.5z−1−0.2z−2
H(z)=−1−2
1+0.3z−0.2z
对于功率谱，可以这样得到，
2
Px(ω)=σwH(z)z=exp(jω)
2
对于ARMA过程，当模型的所有极点均落在单位圆内时，才是一个渐**稳的随机过程。这个过程的自相关函数不能简单地写成Yule-Walker方程形式，它于模型的参数具有高度的非线性关系。
产生样本函数，并画出波形
题目中的ARMA过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。 1.
２８
clear all;
b=[1 0.5 -0.2]; a=[1 0.3 -0.2]; % 由描述的差分方程，得到系统传递函数 h=impz(b,a,10); % 得到系统的单位冲激函数，在10点处已经可以认为值是0 randn(‘state’,0);
w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为2
x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的(2,2)阶ARMA过程 plot(x,’r’);
得到的输出序列波形为：
2)n(X0-2-4-6050100150200250n300350400450500
2. 估计均值和方差
可以首先计算出理论输出的均值和方差，得到mx=0，对于方差可以先求出理论自相
关输出，然后取零点的值。
rx(m)=h(m)*h(−m)*rw(m)
并且，rw(m)=4δ(m)，带入有
rx(m)=4[h(m)*h(−m)]
可以采用上面介绍的方法，对式中的卷积进行计算。计算出的卷积输出图形为：
1.210.8)n-(h*)n(h0.60.40.20-0.205101520n25303540
２９
在最大值处就是输出的功率，也就是方差，为
σx2=rx(0)=4.2
对实际数据进行估计，均值为mean(x)= -0.0547，而方差为var(x)=3.8，两者和理论值吻合的比较好。
画出理论的功率谱密度曲线 理论的功率谱为，
3.
Px(ω)=Pw(ω)H(ejω)=4H(ejω)
用下面的语句产生：
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’),grid on;
22
10.90.80.7Bd )w(xG0.60.54.
0.4
0.3
0.2-500-400-300-200-1000100200300400500f (Hz)
估计相关函数和功率谱密度曲线
用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出仿真图形。
% 计算理论和实际的自相关函数序列 Mlag=20; % 定义最大自相关长度 Rx=xcorr(x,Mlag,'coeff'); m=-Mlag:Mlag; stem(m,Rx,'r.');
最终的值为
３０
1.210.80.6)n(xR0.40.20-0.2-20-15-10-50n5101520
实际的功率谱密度可以用类似于上面的方法进行估计，
window=hamming(20); % 采用hanmming窗，长度为20 noverlap=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’);
估计出来的功率谱密度为，
0-1-2Bd )w(xP-3-4-5-6-500-400-300-200-1000f (Hz)100200300400500
把两幅图画在一起，可以得到下面的图形，可以看出两者的吻合度比较高。
0理论值真实值-1-2Bd )w(xG-3-4-5-6-500-400-300-200-1000100f (Hz)200300400500
３１
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）写出详细试验报告(要有自己对实验结果的结论)。
３２
