离散时间随机过程建模实验报告

1、实验报告姓名:学号:班级:组号:实验名称：离散时间随机过程建模课程名称：统计信号处理基础实验室名称：实验日期：2012.10.10一、实验目的、要求本实验的目的是在了解了Matlab编程语言的编程和调试的基础上，利用Matlab本身自带的函数来验证随机信号建模，并掌握子函数的编写方法。计算机根据理论模型生成随机数，学生需要根据观测的数据编程来计算随机过程的参数。本实验主要是为了让学生在充分理解不同的随机过程建模的理论方法的基础上，用计算机来认识理论和仿真模型之间的差异。要求包括以下几个部分：1 .要求独立完成实验的内容所要求的各项功能，编制完整的Matlab程序，并在程序中注释说明各段程序的功

2、能。2 .要填写完整的实验报告，报告应包含程序、图形和结论。要求记录在实验过程中碰到的问题，以及解决的方法和途径。3 .实验报告是现场用Word填写并打印完成。个人或组必须在报告上署名。二、实验环境验所要求的设备：每组包含完整的计算机1台；可共用的打印机1台，A4纸张若干；计算机上安装的软件包括：Matlab6.5以上(应包含SignalProcessingToolbox,FilterDesignToolbox);Word2000以上；三、实验原理实验内容包括2个，实验1.本实验主要是采用FIR最小二乘逆滤波器来实现反卷积。假定观测的数据y(n)是由信号x(n)通过脉冲响应为g(n)cos(0

3、.2n25)exp0.01n252;0n500;其它的滤波器而生成的。如果从y(n)中恢复的信号x(n)是一组脉冲序列,10x(n)x(k)(nnJk1其中x(k)和nk的取值为nk254055658595110130140155x(k)10.80.70.50.70.20.90.50.60.3a.根据上面的关系，画出观测数据y(n)x(n)g(n),并看看是否能通过y(n)的峰值来确定x(n)的幅度和位置。(需要调用conv函数)程序n=0:50;g=cos(n-25)/5).*exp(-(n-25).*(n-25)/100);g(51:200)=0;x=zeros(200,1);x(25)=

4、1;x(40)=0.8;x(55)=0.7;x(65)=0.5;x(85)=0.7;x(95)=0.2;x(110)=0.9;x(130)=0.5;x(140)=0.6;x(155)=0.2;y=conv(x,g);figure(1)subplot(3,1,1),plot(g);title('滤波器冲击响应);xlabel('n');ylabel('响应幅值');subplot(3,1,2),plot(x);title('输入序列x');xlabel('n');ylabel('幅值')；subplot(3,

5、1,3),plot(y);title('滤波器输出);xlabel('n');ylabel('幅值');整理文本SS'滤浪盥冲击口向应2040 G0 001001201401C0130200n 输入序列支20406。301001201401601802008J曾n滤波器输出50100115020025030。350400b.用教材中给出的spike.m函数来设计长度N50的最小二乘逆滤波器hN(n),并确定最佳的延迟。程序err1=10;N=50;n1=0;forn0=0:200;h,err=spike(g,n0,N);iferr<err1

6、err1=err;H=h;n1=n0;endendfigure(2)plot(H);title('逆滤波器7响应);xlabel('n');ylabel('幅值')；最佳延迟为37c.用估计的hN(n)来滤波y(n),并画出滤波器的输出及(n)hN(n)y(n),图中的峰值的位置和幅度是否与x(n)中的结果一致。x1=conv(y,H);n2=-37:length(x1)-38;figure(3)subplot(2,1,1),plot(x);title('输入序列x(n)');xlabel('n');ylabel('

7、;幅值')；subplot(2,1,2),plot(n2,x1);axis(0,200,0,1);title('逆滤波器输出y');xlabel('n');ylabel('幅值');输入序列d.如果观测数据中还包含噪声，即观测数据为%n)x(n)g(n)v(n),其中v(n)是方差为2的高斯白噪声，分别取:0,0001,20.001,重复b和c中的计算分析。评论这时获得的结果。0.0001时v=0.0001;y1=y(1:205)+v.*randn(1,205);x2=conv(y1,H);n3=-n1:length(x2)-n1-1;f

8、igure(4)subplot(2,1,1),plot(x);title('输入序列 x(n)');xlabel('n');ylabel('幅值);subplot(2,1,2),plot(n3,x2);axis(0,200,0,5);title('逆滤波器输出y1');xlabel('n');ylabel('幅值);输入序列乂 (n)逆漉波器输出V10.001时v=0.001;y1=y(1:205)+v.*randn(1,205);x2=conv(y1,H);n3=-n1:length(x2)-n1-1;figur

9、e(4)subplot(2,1,1),plot(x);Mle('输入序列x(n)');xlabel('n');ylabel('幅值')；subplot(2,1,2),plot(n3,x2);axis(0,200,0,5);title('逆滤波器输出y1');xlabel('n');ylabel('幅值')；输入序列乂(n)逆漉浪器输出v1e.如果g(n)的测量也包含噪声，即g%n)g(n)w(n),而w(n)是在0.005,0.005间均匀分布的白噪声，重复b和c中的计算分析。评论这时获得的结果。程

10、序r=0.001/12.*rand(1,length(g);g1=g+r;y2=conv(x,g1);H1=;err0=1;forn0=0:200;h,err=spike(g1,n0,N);iferr<err0err0=err;H1=h;N0=n0;endendfigure(5)plot(H1);title('逆滤波器冲击响应);xlabel('n');ylabel('幅值')；x3=conv(y2,H1);n2=-37:length(x3)-n1-1;figure(6)subplot(2,1,1),plot(x);title('输入序列x

11、(n)');xlabel('n');ylabel('幅值')；subplot(2,1,2),plot(n2,x3);axis(0,200,0,1);title('逆滤波器输出y2');xlabel('n');ylabel('幅值')；实验2.本实验是用计算机编程来求解ARMAi程的模型参数。a.根据教材上给出的方法，编写一个给定自相关序列rx(k),采用修改的Yule-Walker方程方法来求解ARMA(p,q)模型参数的程序程序functionbe,ae=arma1(r,p,q)r1=r(1,p+2*q+

12、1:end-1);r2=r(1,p+2*q+1:-1:2*q+2);R1=toeplitz(r1,r2);r3=r(1,p+2*q+2:end)'ae=1;-inv(R1)*r3;R2=toeplitz(r(1,p+q+1:2*q+p+1),r(1,p+q+1:-1:q+1);c=R2*ae;d=conv(c,flipud(ae);dc=d(p+2:end,1);pd=flipud(dc);d(p+1,1);dc;ifq=0rt=;be=zp2tf(rt,口,1);elsert=roots(pd);be=zp2tf(rt(q+1:end,1),口,1);end图形b.让单位方差的高斯白

13、噪声通过下列滤波器12»、10.9z0.18zH(z)123411.978z2.853z1.877z0.904z得到观测数据x(n)的100个样本，画出x(n)的理论功率谱。程序A=1,-1.978,2.853,-1.877,0.904;B=1,-0.9,0.18;v=randn(1,100);x=filter(B,A,v);p=length(A)-1;q=length(B)-1;rx=xcorr(x,p+q,'biased');pw=abs(fft(rx)/100;figure(1)c.用a中编制的程序根据观测数据来求解ARMA(4,2)模型的参数，把计算结果与理论模型的系数相比，有什么结论。重复10次不同的样本实现，并计算10次的模型参数再取平均，与理想的系数相比，平均是否有效果。程序sa=zeros(5,1);sb=zeros(1,3);fori=1:10v=randn(1,100);x=filter(B,A,v);rx=xcorr(x,length(A)+length(B)-2,'biased');b,a=arma1(rx,p,q);sa=sa+a;sb=sb+b;endb1=sb./10a1=sa./10err1=sum(b-B).*(b-B)+sum(a-A).*(a-A')err10=sum