离散时间随机过程建模

-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函数）b.用教材中给出的spike.m函数来设计长度50N=的最小二乘逆滤波器()Nhn，并确定最佳的延迟。c.用估计的()Nhn来滤波()yn，并画出滤波器的输出ˆ()()()Nxnhnyn=∗，图中的峰值的位置和幅度是否与()xn中的结果一致。d.如果观测数据中还包含噪声，即观测数据为()()()()ynxngnvn=∗+，其中()vn是方差为2vσ的高斯白噪声，分别取220.0001,0.001vvσσ==，重复b和c中的计算分析。评论这时获得的结果。e.如果()gn的测量也包含噪声，即()()()gngnwn=+，而()wn是在[0.005,0.005]−间均匀分布的白噪声，重复b和c中的计算分析。评论这时获得的结果。实验二、本实验是用计算机编程来求解ARMA过程的模型参数。-2-a.根据教材上给出的方法，编写一个给定自相关序列()xrk，采用修改的Yule-Walker方程方法来求解ARMA(,)pq模型参数的程序b.让单位方差的高斯白噪声通过下列滤波器12123410.90.18()11.9782.8531.8770.904zzHzzzzz−−−−−−−+=−+−+得到观测数据()xn的100个样本，画出()xn的理论功率谱。c.用a中编制的程序根据观测数据来求解ARMA(4,2)模型的参数，把计算结果与理论模型的系数相比，有什么结论。重复10次不同的样本实现，并计算10次的模型参数再取平均，与理想的系数相比，平均是否有效果。closeall;clear;clc;N=50;n=0:N;m=0:N+155-1;g=cos(0.2*(n-0.25)).*exp(-0.01*(n-25).^2);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.3;%x=[zeros(1,24)1zeros(1,14)0.8zeros(1,14)0.7zeros(1,9)0.5...%zeros(1,19)0.7zeros(1,9)0.2zeros(1,14)0.9zeros(1,19)0.5zeros(1,9)0.6...%zeros(1,14)0.3];%%%%%%%%%%%%%a%%%%%%%%%%%%%y=conv(x,g);plot(m,y,'b-');holdon;plot(n,g,'r--');gridon;legend('y(n)','h(n)');xlabel('自变量：n');ylabel('函数值：y(n)');title('h(n)和y(n)信号');%%%%%%%%%%%%%b%%%%%%%%%%%%%Err=[];min=1;forn0=0:N[herr]=spike(g,n0,N);Err=[Errerr];ifminerrmin=err;index=n0;-3-hn=h;endendfigure;plot(n,Err);set(gca,'XTickMode','manual','XTick',[index]);set(gca,'YTickMode','manual','YTick',[min0.40.60.8]);gridon;xlabel('延迟时间：n0');ylabel('误差：\xi(N)');title('延迟时间和对应误差分布图');%%%%%%%%%%%%%c%%%%%%%%%%%%%xn=conv(hn,y);%nn=0:length(xn)-1;nn=0:154+index;mm=[0:154]+index;figure;%plot(nn,xn,'b',mm,x,'ro');plot(nn,xn(nn+1),'b',mm,x,'r.');legend('估计信号','真实信号');xlabel('自变量：n');ylabel('信号：x(n)');title('估计信号与真实信号分布图');%%%%%%%%%%%%%d%%%%%%%%%%%%%fori=-4:-3yd=conv(x,g)+sqrt(10^i)*randn(size(m));xn=conv(hn,yd);%nn=0:length(xn)-1;nn=0:154+index;mm=[0:154]+index;figure;%plot(nn,xn,'b',mm,x,'ro');plot(nn,xn(nn+1),'b',mm,x,'r.');legend('估计信号','真实信号');xlabel('自变量：n');ylabel('信号：x(n)');title('加噪声后估计信号与真实信号分布图');end%%%%%%%%%%%%%e%%%%%%%%%%%%%gn=g+0.01*(rand(size(n))-0.5);Err=[];min=1;forn0=0:N-4-[herr]=spike(gn,n0,N);Err=[Errerr];ifminerrmin=err;index=n0;hn=h;endendfori=-4:-3ye=conv(x,g)+sqrt(10^i)*randn(size(m));xn=conv(hn,ye);%nn=0:length(xn)-1;nn=0:154+index;mm=[0:154]+index;figure;%plot(nn,xn,'b',mm,x,'ro');plot(nn,xn(nn+1),'b',mm,x,'r.');legend('估计信号','真实信号');xlabel('自变量：n');ylabel('信号：x(n)');title('加噪声后估计信号与真实信号分布图');endcloseall;clear;clc;N=1000;cnt=20;n=0:N-1;b=[1-0.90.18];a=[1-1.9782.853-1.8770.904];%%%%%%%%%%%%%b%%%%%%%%%%%%%noise=randn(size(n));x=filter(b,a,noise);[Hw]=freqz(b,a);plot(w/pi,abs(H).^2);gridon;xlabel('频率:\pi');ylabel('X(e^{j\omega})');title('X(n)的理论功率谱');%%%%%%%%%%%%%c%%%%%%%%%%%%%M=armax(x',[42]);aa=M.a;bb=M.c;[HHw]=freqz(bb,aa);figure;-5-plot(w/pi,abs(HH).^2,'b-');holdon;plot(w/pi,abs(H).^2,'r--');gridon;legend('估计值','理论值');xlabel('频率:\pi');ylabel('X(e^{j\omega})');title('X(n)的估计功率谱');fori=0:cnt-2noise=randn(size(n));x=filter(b,a,noise);MM=armax(x',[42]);ac=MM.a;bc=MM.c;bb=bb+bc;aa=aa+ac;endbb=bb/cnt;aa=aa/cnt;[HHHw]=freqz(bb,aa);figure;plot(w/pi,abs(HHH).^2,'b-');holdon;plot(w/pi,abs(H).^2,'r--');gridon;legend('估计平均值','理论值');xlabel('频率:\pi');ylabel('X(e^{j\omega})');title('X(n)的估计功率谱');function[bq,ap]=armahat(rx,p,q)N=length(rx);Rx=convm(rx((N+1)/2:N),p+1);Rq=Rx(q+p-1:(N+1)/2,2:p+1);ap=[1;-Rq\Rx(p+q-1:(N+1)/2,1)];Rxa=xcorr(ap);Rxx=convm(rx,2*p+1);Rxb=Rxx*Rxa;Rxb=Rxb((length(Rxb)-1)/2-q+1:(length(Rxb)-1)/2+q+1);root=roots(Rxb);M=length(root);t=1;fork=1:Mif(abs(root(k))1)xzeros(t)=root(k);t=t+1;endend-6-[bq,aa]=zp2tf(xzeros',0,1);functionX=convm(x,p)N=length(x)+2*p-2;x=x(:);xpad=[zeros(p-1,1);x;zeros(p-1,1)];fori=1:pX(:,i)=xpad(p-i+1:N-i+1);end;function[h,err]=spike(g,n0,n)g=g(:);m=length(g);ifm+n-1=n0error('Delaytoolarge!')endG=convm(g,n);d=zeros(m+n-1,1);d(n0+1)=1;h=G\d;err=1-G(n0+1,:)*h;