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

实验报告

姓名： 学号： 班级： 组号：

实验名称： 课程名称： 实验室名称： 实验日期：

离散时间随机过程建模 统计信号处理基础

2012.10.10

一、实验目的、要求

本实验的目的是在了解了Matlab编程语言的编程和调试的基础上，利用Matlab本身自带的函数来验证随机信号建模，并掌握子函数的编写方法。计算机根据理论模型生成随机数，学生需要根据观测的数据编程来计算随机过程的参数。本实验主要是为了让学生在充分理解不同的随机过程建模的理论方法的基础上，用计算机来认识理论和仿真模型之间的差异。

要求包括以下几个部分：

1．要求完成实验的内容所要求的各项功能，编制完整的Matlab程序，并在程序中注释说明各段程序的功能。

2．要填写完整的实验报告，报告应包含程序、图形和结论。要求记录在实验过程中碰到的问题，以及解决的方法和途径。

3．实验报告是现场用Word填写并打印完成。个人或组必须在报告上署名。

二、实验环境

验所要求的设备： 每组包含完整的计算机 1 台；

可共用的打印机1台，A4纸张若干；

计算机上安装的软件包括： Matlab 6.5以上（应包含Signal Processing Toolbox, Filter Design Toolbox）； Word 2000以上；

三、实验原理

实验内容包括2个，

实验1．本实验主要是采用FIR最小二乘逆滤波器来实现反卷积。假定观测的数据y(n)是由信号x(n)通过脉冲响应为

专业.整理

. WORD格式.资料 .

cos(0.2[n25])exp{0.01[n25]2};0n50 g(n);其它0的滤波器而生成的。如果从y(n)中恢复的信号x(n)是一组脉冲序列，

x(n)x(k)(nnk)

k110其中x(k)和nk的取值为

nk 25 40 0.8 55 0.7 65 0.5 85 0.7 95 0.2 110 0.9 130 0.5 140 0.6 155 0.3 x(k) 1 a. 根据上面的关系，画出观测数据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)= 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('幅值'); 专业.整理

. WORD格式.资料 .

图形

b. 用教材中给出的spike.m函数来设计长度N50的最小二乘逆滤波器hN(n)，并确定最佳的延迟。 程err1=10;N=50;n1=0; 序 for n0=0:200; [h,err]=spike(g,n0,N); if err. WORD格式.资料 .

图形 最佳延迟为37

ˆ(n)hN(n)y(n)，图中的峰c. 用估计的hN(n)来滤波y(n)，并画出滤波器的输出x值的位置和幅度是否与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('幅值'); subplot(2,1,2),plot(n2,x1);axis([0,200,0,1]);title('逆滤波器输出y');xlabel('n');ylabel('幅值'); 图形 专业.整理

. WORD格式.资料 .

d. 如果观测数据中还包含噪声，即观测数据为y(n)x(n)g(n)v(n)，其中v(n)是

222方差为v的高斯白噪声，分别取v0.0001,v0.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]; 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); 专业.整理

. WORD格式.资料 .

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. 如果g(n)的测量也包含噪声，即g(n)g(n)w(n)，而w(n)是在[0.005,0.005]间均匀分布的白噪声，重复b和c中的计算分析。评论这时获得的结果。 程r=0.001/12.*rand(1,length(g)); 序 g1=g+r; y2=conv(x,g1); H1=[];err0=1; for n0=0:200; [h,err]=spike(g1,n0,N); if err. WORD格式.资料 .

err0=err;H1=h;N0=n0;end end figure(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. 根据教材上给出的方法，编写一个给定自相关序列rx(k)，采用修改的Yule-Walker

专业.整理

. WORD格式.资料 .

方程方法来求解ARMA(p,q)模型参数的程序

程序 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]; if q == 0 rt = []; be = zp2tf(rt,[],1); else rt = roots(pd); be = zp2tf(rt(q+1:end,1),[],1); end 图形

b. 让单位方差的高斯白噪声通过下列滤波器

10.9z10.18z2H(z)

11.978z12.853z21.877z30.904z4得到观测数据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) 专业.整理

. WORD格式.资料 .

plot(pw);title('理率谱');xlabel('n');ylabel('幅值'); 图形

c. 用a中编制的程序根据观测数据来求解ARMA(4,2)模型的参数，把计算结果与理论模型的系数相比，有什么结论。重复10次不同的样本实现，并计算10次的模型参数再取平均，与理想的系数相比，平均是否有效果。

程序 sa=zeros(5,1);sb=zeros(1,3); for i=1:10 v=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; end b1=sb./10 a1=sa./10 err1=sum((b-B).*(b-B))+sum((a-A').*(a-A')) err10=sum((b1-B).*(b1-B))+sum((a1-A').*(a1-A')) 图形 A= 1 -1.978 2.853 -1.877 0.904 B= 1 -0.9 0.18 b = 1.0000 -0.5106 0.3082 a = 1.0000 -1.61 2.7162 -1.7549 0.8683 b1 = 1.0000 -0.98 0.4565 a1 = 1.0000 -2.2047 3.2952 -2.3071 1.1074 专业.整理

. WORD格式.资料 .

平均效果不佳

指导教师评语:

成绩： 指导教师签名：

专业.整理 批阅日期：