| 姓名： | 实验名称： | 离散时间随机过程建模 | |
| 学号： | 课程名称： | 统计信号处理基础 | |
| 班级： | 实验室名称： | ||
| 组号： | 实验日期： | 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最小二乘逆滤波器来实现反卷积。假定观测的数据是由信号通过脉冲响应为
的滤波器而生成的。如果从中恢复的信号是一组脉冲序列，
其中的取值为
| 25 | 40 | 55 | 65 | 85 | 95 | 110 | 130 | 140 | 155 | |
| 1 | 0.8 | 0.7 | 0.5 | 0.7 | 0.2 | 0.9 | 0.5 | 0.6 | 0.3 |
| 程序 | 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('幅值'); |
| 图形 |
| 程序 | err1=10;N=50;n1=0; for n0=0:200; [h,err]=spike(g,n0,N); if err end figure(2) plot(H);title('逆滤波器冲击响应');xlabel('n');ylabel('幅值'); |
| 图形 |
c. 用估计的来滤波，并画出滤波器的输出，图中的峰值的位置和幅度是否与中的结果一致。
| 程序 | 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('幅值'); |
| 图形 |
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('幅值'); |
| 图形 |
| 程序 | 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('幅值'); |
| 图形 |
| 程序 | 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 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('幅值'); |
| 图形 |
a. 根据教材上给出的方法，编写一个给定自相关序列，采用修改的Yule-Walker方程方法来求解模型参数的程序
| 程序 | 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 |
| 图形 |
得到观测数据的100个样本，画出的理率谱。
| 程序 | 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('幅值'); |
| 图形 |
| 程序 | 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 |
指导教师评语:
成绩： 指导教师签名：
批阅日期：下载本文
