现代数字信号处理_仿真作业_万群.doc

作业（1） 已知二阶AR过程：，模型参数如下（分两组试验）：取三组参数： a1= -0.195 a2=0.95 a1= -1.5955 a2=0.95 a1= -1.9114 a2=0.95 （1）若要求，求=？ （2）由上述模型参数计算并画出这个二阶AR过程的真实自相关函数（k=0,1,2…255）和功率谱； （3）产生观测数据，，； （4）由观测数据估计并画出样本自相关函数，k=0,1,2…255； （5）由观测数据估计模型参数； （6）由估计的模型参数估计并画出这个二阶AR过程的功率谱； （7）由观测数据估计并画出周期图、相关图和BT法得到的功率谱估计 （8）对上述结果进行比较和分析，给出比较分析的结论。 程序： close all; clear; clc; a1 = -0.195; a2 = 0.95; r0 = 1; %a1= -1.5955; a2=0.95; r0 = 1; %a1= -1.9114; a2=0.95; r0 = 1; A = [1 -a1 -a2; 0 -(1+a2) 0; 0 -a1 -1]; b = [r0 a1*r0 a2*r0]; x = inv(A)*b; %x = [sigma_v r1 r2]; sigma_v=x(1); r1=x(2); r2=x(3); %求解Yule-Walker方程 p=2; r=zeros(1,256); r(1)=r0; r(2)=r1; r(3)= r2; for m=p+1:256-1 r(m+1)=-(a1*r(m)+a2*r(m-1)); %外推r(m) end plot(0:255,r,*-); axis( [0 255 -1.2 1.2] ); dw=2*pi/256; w=-pi:dw:pi-dw; Sw=zeros(1,length(w)); for i=1:length(Sw); Sw(i)=sigma_v/( abs(1+a1*exp(-j*w(i))+a2*exp(-j*w(i)*2)).^2 ); end Sw=10*log10(Sw); Sw=Sw-max(Sw); %Sar(w) figure(); subplot(2,1,1); plot(w,Sw); N=1000; xtemp=zeros(N+p,1); vn=wgn(N+p,1,10*log10(sigma_v)); %方差为sigma_v的高斯白噪声 for i=1:N %产生2阶AR信号 xtemp(i+2)=vn(i)-a1*xtemp(i+1)-a2*xtemp(i); end x=xtemp(p+1:length(xtemp)); Ln=length(x); x2=[0 0 x]; r_e=zeros(1,p+1); for m=1:p+1 index=p+1-(m-1); r_e(m)=x2(index:Ln+index-1)*x/Ln; %由观测数据估计r(m) end R_e=[r_e(1) r_e(2); r_e(2) r_e(1)]; b=[-r_e(2) -r_e(3)]; a_e=inv(R_e)*b; %a_e=[a1_e a2_e]; %由观测数据估计模型参数 sigma_e=r_e(1)+[r_e(2) r_e(3)]*a_e Sw_e=zeros(1,length(w)); for i=1:length(Sw_e); Sw_e(i)=sigma_v/( abs(1+[exp(-j*w(i)) exp(-j*w(i)*2)]*a_e).^2 ); end Sw_e=10*log10(Sw_e); Sw_e=Sw_e-max(Sw_e); subplot(2,1,2); plot(w,Sw_e); %Sar(w)估计 Sw_p=(abs(fft(x)).^2)/length(x); Sw_p=10*log10(Sw_p); Sw_p=Sw_p-max(Sw_p); figure(); dw=2*pi/length(x); w=-pi:dw:pi-dw; plot(w,Sw_p); %Sar(w)周期图谱估计 x_inv=x(length(x):-1:1); rm_e=conv(x,x_inv)./length(x); Sw_c=abs(fft(rm_e)); Sw_c=10*log10(Sw_c); Sw_c=Sw_c-max(Sw_c); figure(); dw=2*pi/length(rm_e); plot(-pi:dw: