搜索
实验报告姓名：实验名称：离散时间随机过程建模学号：课程名称：统计信号处理基础班级：实验室名称：组号：实验日期：2012.10.10一、实验目的、要求本实验的目的是在了解了Matlab编程语言的编程和调试的基础上，利用Matlab本身自带的函数来验证随机信号建模，并掌握子函数的编写方法。计算机根据理论模型生成随机数，学生需要根据观测的数据编程来计算随机过程的参数。本实验主要是为了让学生在充分理解不同的随机过程建模的理论方法的基础上，用计算机来认识理论和仿真模型之间的差异。要求包括以下几个部分：1．要求独立完成实验的内容所要求的各项功能，编制完整的Matlab程序，并在程序中注释说明各段程序的功能。2．要填写完整的实验报告，报告应包含程序、图形和结论。要求记录在实验过程中碰到的问题，以及解决的方法和途径。3．实验报告是现场用Word填写并打印完成。个人或组必须在报告上署名。二、实验环境验所要求的设备：每组包含完整的计算机1台；可共用的打印机1台，A4纸张若干；计算机上安装的软件包括：Matlab6.5以上（应包含SignalProcessingToolbox,FilterDesignToolbox）；Word2000以上；三、实验原理实验内容包括2个，实验1．本实验主要是采用FIR最小二乘逆滤波器来实现反卷积。假定观测的数据()yn是由信号()xn通过脉冲响应为2cos(0.2[25])exp{0.01[25]};050()0;nnngn其它的滤波器而生成的。如果从()yn中恢复的信号()xn是一组脉冲序列，101()()()kkxnxknn其中()kxkn和的取值为kn254055658595110130140155()xk10.80.70.50.70.20.90.50.60.3a.根据上面的关系，画出观测数据()()()ynxngn，并看看是否能通过()yn的峰值来确定()xn的幅度和位置。（需要调用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)=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,1,3),plot(y);title('滤波器输出');xlabel('n');ylabel('幅值');图形b.用教材中给出的spike.m函数来设计长度50N的最小二乘逆滤波器()Nhn，并确定最佳的延迟。程序err1=10;N=50;n1=0;forn0=0:200;[h,err]=spike(g,n0,N);iferrerr1err1=err;H=h;n1=n0;endendfigure(2)plot(H);title('逆滤波器冲击响应');xlabel('n');ylabel('幅值');图形最佳延迟为37c.用估计的()Nhn来滤波()yn，并画出滤波器的输出ˆ()()()Nxnhnyn，图中的峰值的位置和幅度是否与()xn中的结果一致。程序x1=conv(y,H);n2=[-37:length(x1)-38];figure(3)subplot(2,1,1),plot(x);title('输入序列x(n)');xlabel('n');ylabel('幅值');subplot(2,1,2),plot(n2,x1);axis([0,200,0,1]);title('逆滤波器输出y');xlabel('n');ylabel('幅值');图形d.如果观测数据中还包含噪声，即观测数据为()()()()ynxngnvn，其中()vn是方差为2v的高斯白噪声，分别取220.0001,0.001vv，重复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];figure(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('幅值');图形0.001时程序v=0.001;y1=y(1:205)+v.*randn(1,205);x2=conv(y1,H);n3=[-n1:length(x2)-n1-1];figure(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('幅值');图形e.如果()gn的测量也包含噪声，即()()()gngnwn，而()wn是在[0.005,0.005]间均匀分布的白噪声，重复b和c中的计算分析。评论这时获得的结果。程序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);iferrerr0err0=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(n)');xlabel('n');ylabel('幅值');subplot(2,1,2),plot(n2,x3);axis([0,200,0,1]);title('逆滤波器输出y2');xlabel('n');ylabel('幅值');图形实验2.本实验是用计算机编程来求解ARMA过程的模型参数。a.根据教材上给出的方法，编写一个给定自相关序列()xrk，采用修改的Yule-Walker方程方法来求解ARMA(,)pq模型参数的程序程序function[be,ae]=arma1(r,p,q)r1=r(1,p+2*q+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.让单位方差的高斯白噪声通过下列滤波器12123410.90.18()11.9782.8531.8770.904zzHzzzzz得到观测数据()xn的100个样本，画出()xn的理论功率谱。程序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)plot(pw);title('理论功率谱');xlabel('n');ylabel('幅值');图形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((b1-B).*(b1-B))+sum((a1-A').*(a1-A'))图形A=1-1.9782.853-1.8770.904B=1-0.90.18b=1.0000-0.51060.3082a=1.0000-1.89612.7162-1.75490.8683b1=1.0000-0.98980.4565a1=1.0000-2.20473.2952-2.30711.1074平均效果不佳指导教师评语:成绩：指导教师签名：批阅日期：