典型时间序列模型分析
- 格式:doc
- 大小:254.37 KB
- 文档页数:18
实验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)用MA TLAB 模拟产生X(n)的500 观测点的样本函数,并绘出波形 (2)用产生的500 个观测点估计X(n)的均值和方差 (3)画出理论的功率谱
(4)估计X(n)的相关函数和功率谱
【分析】给定二阶的AR 过程,可以用递推公式得出最终的输出序列。或者按照一个白噪声 通过线性系统的方式得到,这个系统的传递函数为:
1
2
1
()10.30.5H z z z --=
++
这是一个全极点的滤波器,具有无限长的冲激响应。 对于功率谱,可以这样得到,
()()
2
2
12
12exp 11x w
z jw P w a z a z σ--==++
可以看出,
()
x P w 完全由两个极点位置决定。
对于 AR 模型的自相关函数,有下面的公式:
这称为 Yule-Walker 方程,当相关长度大于p 时,由递推式求出:
这样,就可以求出理论的 AR 模型的自相关序列。
1.产生样本函数,并画出波形
2.题目中的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 on;
得到的输出序列波形为:
2.估计均值和方差
可以首先计算出理论输出的均值和方差,得到
x
m
,对于方差可以先求出理论自相
关输出,然后取零点的值。
并且,,带入有
在最大值处输出的功率,也就是方差,为
对实际数据进行估计,均值为mean(x)=-0.0703,而方差为var(x)=5.2795,两者合理论值吻合得比较好。
程序及运行结果图如下,其中y_mean表示均值,y_var表示方差。
3.画出理论的功率谱密度曲线
理论的功率谱为,
用下面的语句产生:
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);
title('邹先雄——理论功率谱密度曲线');
grid on;
得到的图形为:
可以看出,这个系统是带通系统。
4.估计自相关函数和功率谱密度
用实际数据估计自相关函数和功率谱的方法前面已经讨论过,在这里仅给出最后的仿真图形。
Mlag=20; % 定义最大自相关长度
Rx=xcorr(x,Mlag,'coeff');
m=-Mlag:Mlag;
stem(m,Rx,'r.');
title('邹先雄——自相关函数');
最终的值为
可以看出,它和上面的理论输出值吻合程度很好。实际的功率谱密度可以用类似于上面的方法进行估计,
window=hamming(20); % 采用hanmming 窗,长度为20
noverlap=10; % 重叠的点数
Nfft=512; % 做FFT 的点数
Fs=1000; % 采样频率,为1000Hz
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 过程
[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('邹先雄——实际的功率谱密度曲线');
估计出来的功率谱密度为,