典型时间序列模型分析doc.docx
《典型时间序列模型分析doc.docx》由会员分享，可在线阅读，更多相关《典型时间序列模型分析doc.docx（25页珍藏版）》请在冰豆网上搜索。
典型时间序列模型分析doc
实验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过程，可以用递推公式得出最终的输出序列。
或者按照一个白噪声
通过线性系统的方式得到，这个系统的传递函数为：
这是一个全极点的滤波器，具有无限长的冲激响应。
对于功率谱，可以这样得到，
可以看出，
完全由两个极点位置决定。
对于AR模型的自相关函数，有下面的公式：
这称为Yule-Walker方程，当相关长度大于p时，由递推式求出：
这样，就可以求出理论的AR模型的自相关序列。
1.产生样本函数，并画出波形
2.题目中的AR过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。
clearall;
b=[1];a=[10.30.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随机序列'）;
gridon;
得到的输出序列波形为：
2.估计均值和方差
可以首先计算出理论输出的均值和方差，得到
，对于方差可以先求出理论自相
关输出，然后取零点的值。
并且，
，带入有
在最大值处输出的功率，也就是方差，为
对实际数据进行估计，均值为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（'邹先雄——理论功率谱密度曲线'）;
gridon;
得到的图形为：
可以看出，这个系统是带通系统。
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=[10.30.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（'邹先雄——实际的功率谱密度曲线'）;
估计出来的功率谱密度为，
将两幅图画在一起，可以看到拟合的情况比较好（两者相位刚好相反，但是基本波形相似）：
代码如下：
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）））.^2）;%计算出理论值
Gx=Gx/max（Gx）;%归一化处理
f=w*Fs/（2*pi）;%转化到模拟域上的频率结束
plot（f,Gx,'r'）;
holdon;
title（'邹先雄——理论和实际的功率谱密度曲线拟合'）;
window=hamming（20）;%采用hanmming窗，长度为20
noverlap=10;%重叠的点数
Nfft=512;%做FFT的点数
Fs=1000;%采样频率，为1000Hz
b=[1];a=[10.30.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）];%对称的功率谱
Py=-10*log10（Py）;
Py=Py/max（Py）;
Py=-Py;Py=3*Py;Py=Py+2.6;%用来归一处理，使两者吻合
plot（f,Py,'b'）;
legend（'实际值','理论值'）;
gridon;
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过程，也可以用递推公式得出最终的输出序列。
或者按照一
个白噪声通过线性系统的方式得到，这个系统的传递函数为：
对于功率谱，可以这样得到，
对于ARMA过程，当模型的所有极点均落在单位圆内时，才是一个渐进平稳的随机过程。
这个过程的自相关函数不能简单地写成Yule-Walker方程形式，它于模型的参数具有高度的
非线性关系。
1.产生样本函数，并画出波形
题目中的ARMA过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照
上面的方法进行描述。
clearall;
b=[10.5-0.2];a=[10.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'）;
title（'邹先雄——输出的AR随机序列'）;
得到的输出序列波形为：
2.估计均值和方差
可以首先计算出理论输出的均值和方差，得到
，对于方差可以先求出理论自相
关输出，然后取零点的值。
并且，
，带入有
在最大值处就是输出的功率，也就是方差，为
对实际数据进行估计，均值为mean（x）=-0.0547，而方差为var（x）=3.8，两者和理论值吻
合的比较好。
附代码及运行结果截图如下：
3.画出理论的功率谱密度曲线
理论的功率谱为，
用下面的语句产生：
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'）;
title（'邹先雄——理论的功率谱密度曲线'）;
gridon;
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=[10.5-0.2];a=[10.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过程
[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（'邹先雄——实际的功率谱密度曲线'）;
估计出来的功率谱密度为
把两幅图画在一起，可以得到下面的图形，可以看出两者的吻合度比较高。
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,'r'）;
title（'邹先雄——理论和实际的功率谱密度曲线的拟合'）;
holdon;
window=hamming（20）;%采用hanmming窗，长度为20
noverlap=10;%重叠的点数
Nfft=512;%做FFT的点数
Fs=1000;%采样频率，为1000Hz
b=[10.5-0.2];a=[10.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过程
[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;
3、实验内容
1、熟悉实验原理，将实验原理上的程序应用matlab工具实现；
2、设有MA
（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点处已经可以认为值是0
randn（'state',0）;
w=normrnd（0,2,1,500）;%产生题设的白噪声随机序列，标准差为2
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）.^2）;%计算出理论值
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（'邹先雄——估计相关函数'）;
%画出估计的功率谱密度曲线
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（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）.^2）;%计算出理论值
Gx=Gx/max（Gx）;f=w*Fs/（2*pi）;%转化到模拟域上的频率
figure（5）;
plot（f,Gx,'r'）;
title（'邹先雄——实际和估计两功率谱密度曲线的拟合'）;
holdon;
window=hamming（20）;%采用hanmming窗，长度为20
noverlap=10;%重叠的点数
Nfft=512;%做FFT的点数
Fs=1000;%采样频率，为1000Hz
b=[1-0.30.2];a=[1];%由描述的差分方程，得到系统传递函数
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过程
[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;
样本函数波形为：
理论功率谱密度曲线为：
估计相关函数波形为：
估计功率谱密度曲线为：
实际和估计两功率谱密度曲线的拟合截图如下：
附程序运行后得到的均值与方差的截图，其中y_mean为均值，大小为-0.1127；y_var为方差，大小为3.9324：
4.实验总结：
通过实验，让我更加的了解了随机序列的均值、方差、功率谱密度以及自相关函数。
通过软件的编程运行结果，加深了对书上理论知识的理解与掌握。
首先第一个实验中两个随机序列的练习让我更容易着手于本次的实验。
通过实验，我对matlab有了更深的认识，在使用matlab的过程中，经常产生问题、发现问题并解决问题，这让我对matlab使用的更加熟练。
这次统计信号实验并不是我第一次接触matlab，以前对matlab的应用让我有了一些基础，这次的实验更是让我学到了不少东西，从选题到做题，让我学到了以前没有接触过的matlab知识。
比如要在一幅图中显示多条曲线时，可以使用holdon语句来实现，或者在画图函数中的参数进行设定（如：
plot（f,（Py）,'r',f,Gy,'b'）;），这在以前没有学过。
升级会员 