典型时间序列模型分析.docx
- 文档编号:7783697
- 上传时间:2023-01-26
- 格式:DOCX
- 页数:26
- 大小:188.04KB
典型时间序列模型分析.docx
《典型时间序列模型分析.docx》由会员分享,可在线阅读,更多相关《典型时间序列模型分析.docx(26页珍藏版)》请在冰豆网上搜索。
典型时间序列模型分析
典型时间序列模型分析实验1
、实验目的1工具对Matlab模型,学会运用模型,MA模型与ARMA熟悉三种典型的时间序列模型:
AR
阶模型的仿真分析,探讨几种模型的适用范围,2,通过对对上述三种模型进行统计特性分析并且通过实验分析理论分析与实验结果之间的差异。
2、实验原理
模型分析:
AR
AR
(2)模型,设有
X(n)=-0.3X(n-1)-0.5X(n-2)+W(n)
4。
其中:
W(n)是零均值正态白噪声,方差为
500观测点的样本函数,并绘出波形模拟产生X(n)的
(1)用MATLAB
X(n)的均值和方差)用产生的500个观测点估计(23)画出理论的功率谱(X(n)的相关函数和功率谱(4)估计过程,可以用递推公式得出最终的输出序列。
或者按照一个白噪声AR【分析】给定二阶的通过线性系统的方式得到,这个系统的传递函数为:
1?
z)H(
2?
?
1z0.5z?
0.3?
1
这是一个全极点的滤波器,具有无限长的冲激响应。
对于功率谱,可以这样得到,21?
?
2?
?
wPwx2?
?
1za1?
z?
a21?
?
jwz?
exp
?
?
wPx完全由两个极点位置决定。
可以看出,对于AR模型的自相关函数,有下面的公式:
p时,由递推式求出:
Yule-Walker方程,当相关长度大于这称为
AR模型的自相关序列。
这样,就可以求出理论的
1
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.估计均值和方差
m?
0,对于方差可以先求出理论自相可以首先计算出理论输出的均值和方差,得到x关输出,然后取零点的值。
2
,带入有并且,
在最大值处输出的功率,也就是方差,为
,两者合理论值吻,而方差为var(x)=5.2795mean(x)=-0.0703对实际数据进行估计,均值为合得比较好。
y_var表示方差。
程序及运行结果图如下,其中y_mean表示均值,
画出理论的功率谱密度曲线3.理论的功率谱为,
用下面的语句产生:
delta=2*pi/1000;
w_min=-pi;
w_max=pi;
Fs=1000;
3
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('邹先雄——自相关函数');
最终的值为
4
实际的功率谱密度可以用类似于上面的方可以看出,它和上面的理论输出值吻合程度很好。
法进行估计,20hanmming窗,长度为window=hamming(20);%采用重叠的点数noverlap=10;%
的点数Nfft=512;%做FFT
1000Hz
采样频率,为Fs=1000;%
由描述的差分方程,得到系统传递函数b=[1];a=[10.30.5];%
020得到系统的单位冲激函数,在点处已经可以认为值是h=impz(b,a,20);%
randn('state',0);
2
w=normrnd(0,2,1,500);%产生题设的白噪声随机序列,标准差为过程2阶AR通过线形系统,得到输出就是题目中要求的x=filter(b,a,w);%
估计功率谱密度[Px,f]=pwelch(x,window,noverlap,Nfft,Fs,'onesided');%
[-Fs/2,Fs/2]f=[-fliplr(f)f(1:
end)];%构造一个对称的频率,范围是对称的功率谱Py=[-fliplr(Px)Px(1:
end)];%
plot(f,10*log10(Py),'b');
');
title('邹先雄——实际的功率谱密度曲线估计出来的功率谱密度为,5
:
(两者相位刚好相反,但是基本波形相似)将两幅图画在一起,可以看到拟合的情况比较好代码如下:
clearall;
delta=2*pi/1000;
w_min=-pi;
w_max=pi;
Fs=1000;
[-pi,pi]得到数字域上的频率取样点,范围是w=w_min:
delta:
w_max;%
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('邹先雄——理论和实际的功率谱密度曲线拟合20窗,长度为window=hamming(20);%采用hanmming
重叠的点数noverlap=10;%
FFT的点数做Nfft=512;%
1000Hz
Fs=1000;%采样频率,为b=[1];a=[10.30.5];%由描述的差分方程,得到系统传递函数0点处已经可以认为值是h=impz(b,a,20);%得到系统的单位冲激函数,在20
randn('state',0);
2
w=normrnd(0,2,1,500);%产生题设的白噪声随机序列,标准差为过程2阶AR通过线形系统,得到输出就是题目中要求的x=filter(b,a,w);%
估计功率谱密度[Px,f]=pwelch(x,window,noverlap,Nfft,Fs,'onesided');%
6
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过程,也可以用递推公式得出最终的输出序列。
或者按照一【分析】给定
个白噪声通过线性系统的方式得到,这个系统的传递函数为:
7
对于功率谱,可以这样得到,
对于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');
');AR随机序列邹先雄——输出的title('得到的输出序列波形为:
2.估计均值和方差
8
m?
0,对于方差可以先求出理论自相可以首先计算出理论输出的均值和方差,得到x关输出,然后取零点的值。
并且,,带入有
在最大值处就是输出的功率,也就是方差,为
对实际数据进行估计,均值为mean(x)=-0.0547,而方差为var(x)=3.8,两者和理论值吻
合的比较好。
附代码及运行结果截图如下:
画出理论的功率谱密度曲线3.
理论的功率谱为,
用下面的语句产生:
delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;
[-pi,pi]w=w_min:
delta:
w_max;%得到数字域上的频率取样点,范围是分子NS=1+0.5*exp(-i*w)-0.2*exp(-2*i*w);%
9
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('邹先雄——估计自相关函数最终的值为
10
实际的功率谱密度可以用类似于上面的方法进行估计,窗,长度为20window=hamming(20);%采用hanmmingnoverlap=10;%重叠的点数FFT的点数Nfft=512;%做1000HzFs=1000;%采样频率,为由描述的差分方程,得到系统传递函数b=[10.5-0.2];a=[10.3-0.2];%
010点处已经可以认为值是h=impz(b,a,10);%得到系统的单位冲激函数,在randn('state',0);
产生题设的白噪声随机序列,标准差为2w=normrnd(0,2,1,500);%
过程阶通过线形系统,得到输出就是题目中要求的(2,2)ARMAx=filter(b,a,w);%
估计功率谱密度[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('邹先雄——实际的功率谱密度曲线估计出来的功率谱密度为11
把两幅图画在一起,可以得到下面的图形,可以看出两者的吻合度比较高。
delta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;
[-pi,pi]w=w_min:
delta:
w_max;%得到数字域上的频率取样点,范围是分子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;
窗,长度为20window=hamming(20);%采用hanmmingnoverlap=10;%重叠的点数的点数Nfft=512;%做FFT1000HzFs=1000;%采样频率,为由描述的差分方程,得到系统传递函数b=[10.5-0.2];a=[10.3-0.2];%
点处已经可以认为值是010h=impz(b,a,10);%得到系统的单位冲激函数,在randn('state',0);
产生题设的白噪声随机序列,标准差为2w=normrnd(0,2,1,500);%
过程阶ARMA(2,2)x=filter(b,a,w);%通过线形系统,得到输出就是题目中要求的[Px,f]=pwelch(x,window,noverlap,Nfft,Fs,'onesided');%估计功率谱密度[-Fs/2,Fs/2]f=[-fliplr(f)f(1:
end)];%构造一个对称的频率,范围是Py=[fliplr(Px)Px(1:
end)];%对称的功率谱Py=10*log10(Py);
Py=Py/max(Py);
12
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
13
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('邹先雄——实际和估计两功率谱密度曲线的拟合');
14
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;
样本函数波形为:
理论功率谱密度曲线为:
15
估计相关函数波形为:
16
估计功率谱密度曲线为:
实际和估计两功率谱密度曲线的拟合截图如下:
17
附程序运行后得到的均值与方差的截图,其中y_mean为均值,大小为-0.1127;y_var为方差,大小为3.9324:
4.实验总结:
通过实验,让我更加的了解了随机序列的均值、方差、功率谱密度以及自相关函数。
通过软件的编程运行结果,加深了对书上理论知识的理解与掌握。
我对通过实验,首先第一个实验中两个随机序列的练习让我更容易着手于本次的实验。
的过程中,经常产生问题、发现问题并解决问题,有了更深的认识,在使用matlabmatlab,以前对matlab使用的更加熟练。
这次统计信号实验并不是我第一次接触这让我对matlab让从选题到做题,matlab的应用让我有了一些基础,这次的实验更是让我学到了不少东西,hold比如要在一幅图中显示多条曲线时,可以使用我学到了以前没有接触过的matlab知识。
,这在以前plot(f,(Py),'r',f,Gy,'b');语句来实现,或者在画图函数中的参数进行设定(如:
)on没有学过。
18
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 典型 时间 序列 模型 分析