典型时间序列模型分析.(9页).doc

《典型时间序列模型分析.(9页).doc》由会员分享，可在线阅读，更多相关《典型时间序列模型分析.(9页).doc（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、-第 1 页典型时间序列模型分析.-第 2 页实验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)的均值和方

2、差（3）画出理论的功率谱（4）估计X(n)的相关函数和功率谱【分析】给定二阶的AR 过程，可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：这是一个全极点的滤波器，具有无限长的冲激响应。对于功率谱，可以这样得到，可以看出，xP w完全由两个极点位置决定。对于 AR 模型的自相关函数，有下面的公式：这称为 Yule-Walker 方程，当相关长度大于p时，由递推式求出：这样，就可以求出理论的 AR 模型的自相关序列。1.1.产生样本函数，并画出波形产生样本函数，并画出波形2.题目中的 AR 过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按

3、照上面的方法进行描述。clear all;b=1;a=1 0.3 0.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 on;得到的输出序列波形为：2.2.估计均值和方差估计均值和方差-第 3 页可以首先计算出理论输出的均值和方

4、差，得到0 xm，对于方差可以先求出理论自相关输出，然后取零点的值。并且，带入有在最大值处输出的功率，也就是方差，为对实际数据进行估计，均值为 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,piGx=4*(abs(1./(

5、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.4.估计自相关函数和功率谱密度估计自相关函数和功率谱密度用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出最后的仿真图形。Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,coeff);m=-Mlag:Mlag;stem(m,Rx,r.);title(邹先雄自相关函

6、数);最终的值为可以看出，它和上面的理论输出值吻合程度很好。实际的功率谱密度可以用类似于上面的方法进行估计，window=hamming(20);%采用 hanmming 窗，长度为 20noverlap=10;%重叠的点数Nfft=512;%做 FFT 的点数Fs=1000;%采样频率，为 1000Hzb=1;a=1 0.3 0.5;%由描述的差分方程，得到系统传递函数h=impz(b,a,20);%得到系统的单位冲激函数，在20 点处已经可以认为值是0-第 4 页randn(state,0);w=normrnd(0,2,1,500);%产生题设的白噪声随机序列，标准差为2x=filter(

7、b,a,w);%通过线形系统，得到输出就是题目中要求的 2 阶 AR 过程Px,f=pwelch(x,window,noverlap,Nfft,Fs,onesided);%估计功率谱密度f=-fliplr(f)f(1:end);%构造一个对称的频率，范围是-Fs/2,Fs/2Py=-fliplr(Px)Px(1:end);%对称的功率谱plot(f,10*log10(Py),b);title(邹先雄实际的功率谱密度曲线);估计出来的功率谱密度为，将两幅图画在一起，可以看到拟合的情况比较好（两者相位刚好相反，但是基本波形相似）：代码如下：clear all;delta=2*pi/1000;w_m

8、in=-pi;w_max=pi;Fs=1000;w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是-pi,piGx=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);hold on;title(邹先雄理论和实际的功率谱密度曲线拟合);window=hamming(20);%采用 hanmming 窗，长度为 20noverlap=10;%重叠的点数Nfft=512;%做 FFT 的点数Fs=100

9、0;%采样频率，为 1000Hzb=1;a=1 0.3 0.5;%由描述的差分方程，得到系统传递函数h=impz(b,a,20);%得到系统的单位冲激函数，在 20 点处已经可以认为值是 0randn(state,0);w=normrnd(0,2,1,500);%产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w);%通过线形系统，得到输出就是题目中要求的 2 阶 AR 过程Px,f=pwelch(x,window,noverlap,Nfft,Fs,onesided);%估计功率谱密度f=-fliplr(f)f(1:end);%构造一个对称的频率，范围是-Fs/2,Fs/2Py

10、=-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(实际值,理论值);grid on;ARMA 模型分析模型分析设有 ARMA(2,2)模型，-第 5 页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)的均值

11、和方差（3）画出理论的功率谱（4）估计 X(n)的相关函数和功率谱【分析】给定(2,2)的ARMA 过程，也可以用递推公式得出最终的输出序列。或者按照一个白噪声通过线性系统的方式得到，这个系统的传递函数为：对于功率谱，可以这样得到，对于 ARMA 过程，当模型的所有极点均落在单位圆内时，才是一个渐进平稳的随机过程。这个过程的自相关函数不能简单地写成Yule-Walker 方程形式，它于模型的参数具有高度的非线性关系。1.产生样本函数，并画出波形产生样本函数，并画出波形题目中的 ARMA 过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。clear all;b=1

12、0.5-0.2;a=1 0.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);title(邹先雄输出的 AR 随机序列);得到的输出序列波形为：2.估计均值和方差估计均值和方差可以首先计算出理论输出的均值和方差，得到0 xm，对于方差可以先求出理论自相关输出，然后取零点的值。并且

13、，带入有在最大值处就是输出的功率，也就是方差，为对实际数据进行估计，均值为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,piNS=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);%分母-

14、第 6 页Gx=4*(abs(NS./DS).2);%计算出理论值Gx=Gx/max(Gx);f=w*Fs/(2*pi);%转化到模拟域上的频率plot(f,Gx,b);title(邹先雄理论的功率谱密度曲线);grid on;4.估计相关函数和功率谱密度曲线估计相关函数和功率谱密度曲线用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出仿真图形。%计算理论和实际的自相关函数序列Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,coeff);m=-Mlag:Mlag;stem(m,Rx,r.);title(邹先雄估计自相关函数);最终的值为实际的功率谱密度可以用

15、类似于上面的方法进行估计，window=hamming(20);%采用 hanmming 窗，长度为 20noverlap=10;%重叠的点数Nfft=512;%做 FFT 的点数Fs=1000;%采样频率，为 1000Hzb=1 0.5-0.2;a=1 0.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)阶 AR

16、MA 过程Px,f=pwelch(x,window,noverlap,Nfft,Fs,onesided);%估计功率谱密度f=-fliplr(f)f(1:end);%构造一个对称的频率，范围是-Fs/2,Fs/2Py=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;%得到数字域上的频率取样点

17、，范围是-pi,piNS=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(邹先雄理论和实际的功率谱密度曲线的拟合);hold on;window=hamming(20);%采用 hanmming 窗，长度为 20noverlap=10;%重叠的点数Nfft=512;%做 FFT 的点数-第 7 页Fs=1000;%采样频率，为

18、1000Hzb=1 0.5-0.2;a=1 0.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 过程Px,f=pwelch(x,window,noverlap,Nfft,Fs,onesided);%估计功率谱密度f=-fliplr(f)f(1:end);%构造一个对称的频率，范围是-Fs/2,Fs/2

19、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(实际值,理论值);grid on;3、实验内容1、熟悉实验原理，将实验原理上的程序应用 matlab 工具实现；2、设有 MA(2)模型，W(n)是零均值正态白噪声，方差为4。（1）用MATLAB 模拟产生X(n)的500 观测点的样本函数，并绘出波形（2）用产生的500 个观测点估计X(n)的均值和方差（3）画出理论的功率谱（4）估计 X(n)的相关函数和功率谱完

20、成 4 个问题的源代码如下clear all;%产生样本函数，并画出波形b=1-0.3 0.2;a=1;%由描述的差分方程，得到系统传递函数h=impz(b,a,10);%得到系统的单位冲激函数，在 10 点处已经可以认为值是 0randn(state,0);w=normrnd(0,2,1,500);%产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w);%通过线形系统，得到输出就是题目中要求的(2,2)阶 ARMA 过程figure(1);plot(x,r);title(邹先雄样本函数);Py_mean=mean(x)Py_var=var(x)%画出理论的功率谱密度曲线del

21、ta=2*pi/1000;w_min=-pi;w_max=pi;Fs=1000;w=w_min:delta:w_max;%得到数字域上的频率取样点，范围是-pi,piNS=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);%转化到模拟域上的频率-第 8 页figure(2);plot(f,Gx,b);title(邹先雄理论的功率谱密度曲线);%估计相关函数Mlag=20;%定义最大自相关长度Rx=xcorr(x,Mlag,coeff);m=-Mla

22、g:Mlag;figure(3);stem(m,Rx,r.);title(邹先雄估计相关函数);%画出估计的功率谱密度曲线window=hamming(20);%采用 hanmming 窗，长度为 20noverlap=10;%重叠的点数Nfft=512;%做 FFT 的点数Fs=1000;%采样频率，为 1000HzPx,f=pwelch(x,window,noverlap,Nfft,Fs,onesided);%估计功率谱密度f=-fliplr(f)f(1:end);%构造一个对称的频率，范围是-Fs/2,Fs/2Py=fliplr(Px)Px(1:end);%对称的功率谱figure(4)

23、;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,piNS=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(邹先雄实际和估计

24、两功率谱密度曲线的拟合);hold on;window=hamming(20);%采用 hanmming 窗，长度为 20noverlap=10;%重叠的点数Nfft=512;%做 FFT 的点数Fs=1000;%采样频率，为 1000Hzb=1-0.3 0.2;a=1;%由描述的差分方程，得到系统传递函数h=impz(b,a,10);%得到系统的单位冲激函数，在 10 点处已经可以认为值是 0randn(state,0);w=normrnd(0,2,1,500);%产生题设的白噪声随机序列，标准差为 2x=filter(b,a,w);%通过线形系统，得到输出就是题目中要求的(2,2)阶 AR

25、MA 过程Px,f=pwelch(x,window,noverlap,Nfft,Fs,onesided);%估计功率谱密度f=-fliplr(f)f(1:end);%构造一个对称的频率，范围是-Fs/2,Fs/2Py=fliplr(Px)Px(1:end);%对称的功率谱-第 9 页Py=10*log10(Py);Py=Py/max(Py);Py=-Py;Py=3*Py;Py=Py+4;%用来归一处理，使两者吻合plot(f,Py,b);legend(实际值,理论值);grid on;样本函数波形为：理论功率谱密度曲线为：估计相关函数波形为：估计功率谱密度曲线为：实际和估计两功率谱密度曲线的拟

26、合截图如下：附程序运行后得到的均值与方差的截图，其中 y_mean 为均值，大小为-0.1127；y_var 为方差，大小为 3.9324：4.实验总结：实验总结：通过实验，让我更加的了解了随机序列的均值、方差、功率谱密度以及自相关函数。通过软件的编程运行结果，加深了对书上理论知识的理解与掌握。首先第一个实验中两个随机序列的练习让我更容易着手于本次的实验。通过实验，我对matlab 有了更深的认识，在使用 matlab 的过程中，经常产生问题、发现问题并解决问题，这让我对 matlab 使用的更加熟练。这次统计信号实验并不是我第一次接触 matlab，以前对matlab 的应用让我有了一些基础，这次的实验更是让我学到了不少东西，从选题到做题，让我学到了以前没有接触过的 matlab 知识。比如要在一幅图中显示多条曲线时，可以使用 holdon 语句来实现，或者在画图函数中的参数进行设定（如：plot(f,(Py),r,f,Gy,b);），这在以前没有学过。