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

1、作业（1）已知二阶AR过程：，模型参数如下（分两组试验）：取三组参数： a1= -0.195a2=0.95 a1= -1.5955a2=0.95 a1= -1.9114a2=0.95（1）若要求，求=？（2）由上述模型参数计算并画出这个二阶AR过程的真实自相关函数（k=0,1,2255）和功率谱；（3）产生观测数据，；（4）由观测数据估计并画出样本自相关函数，k=0,1,2255；（5）由观测数据估计模型参数；（6）由估计的模型参数估计并画出这个二阶AR过程的功率谱； （7）由观测数据估计并画出周期图、相关图和BT法得到的功率谱估计（8）对上述结果进行比较和分析，给出比较分析的结论。程序：cl

2、ose 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)=

3、r2;for m=p+1:256-1 r(m+1)=-(a1*r(m)+a2*r(m-1); %外推r(m)endplot(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 );endSw=10*log10(Sw); Sw=Sw-max(Sw); %Sar(w)figure(); subplot(2,1

4、,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);endx=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(

5、m)endR_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_eSw_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 );endSw_e=10*log10(Sw_e); Sw_e=Sw_e-max(Sw_e);subplot(2,1,2);

6、 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/le

7、ngth(rm_e);plot(-pi:dw:pi-dw,fftshift(Sw_c); %Sar(w)相关图谱估计M=64;index=ceil(length(rm_e)/2);rM_e=rm_e(index-M:index+M);Sw_bt=abs(fft(rM_e);Sw_bt=10*log10(Sw_bt); Sw_bt=Sw_bt-max(Sw_bt);dw=2*pi/length(rM_e);figure();plot(-pi:dw:pi-dw,fftshift(Sw_bt); %Sar(w)BT法谱估计作业（2）已知：，其中为零均值、方差等于的白高斯噪声，信噪比，即，要求：1)

8、产生并画出序列，；2) 估计并画出，；3) 估计并画出周期图；4) 构造4阶自相关矩阵，并计算奇异值SVD分解：，是特征向量，5) 令，估计并画出MUSIC伪谱；6) 估计并画出MVDR功率谱；7) 用ESPRIT方法估计中的两个角频率；8) 对上述功率谱估计和频率估计结果进行比较分析，给出比较分析的结论。程序：c close all; clear; clc;N=1000; Pv=-10;%-10dBn=1:N;vn=wgn(1,N,Pv,'complex');xn=exp(j*pi*(0.12.*n-0.61)+exp(j*pi*(0.31.*n-0.79)+vn;for m

9、=0:500 rm_e(m+1)=xn(m+1:N)*xn(1:N-m)'/(N-m); %r(m)估计end plot(1:500,rm_e(1:500); figure();rm_temp=conj(rm_e(length(rm_e):-1:2);R4_e=rm_e(1:4); rm_e(2)' rm_e(1:3); conj(rm_e(3:-1:2),rm_e(1:2); conj(rm_e(4:-1:1);Rxx_e=R4_e;Rxy_e=R4_e(2:4,:);conj(rm_e(5:-1:2);rm_e=rm_temp rm_e;Pw_bt=abs(fft(rm_e

10、); %周期图谱估计Pw_bt=10*log10(Pw_bt); Pw_bt=Pw_bt-max(Pw_bt);dw=2*pi/length(Pw_bt);plot(-pi:dw:pi-dw,fftshift(Pw_bt);V D=eig(R4_e); %MUSIC伪谱估计G=V(:,1) V(:,2);dw=2*pi/2048; w=-pi:dw:pi-dw;for i=1:length(w) aw=1;exp(-j*w(i);exp(-j*2*w(i);exp(-j*3*w(i); Pw_music(i)=1/abs(aw'*G*G'*aw);endPw_music=10*

11、log10(Pw_music); Pw_music=Pw_music-max(Pw_music);figure(); plot(w,Pw_music);R4_einv=inv(R4_e); %MVDR谱估计for i=1:length(w) aw=1;exp(-j*w(i);exp(-j*2*w(i);exp(-j*3*w(i); Pw_mvdr(i)=1/abs(aw'*R4_einv*aw);endPw_mvdr=10*log10(Pw_mvdr);Pw_mvdr=Pw_mvdr-max(Pw_mvdr);figure(); plot(w,Pw_mvdr);d=eig(Rxx_e)

12、; %用ESPRIT方法估计两个角频率Cxx=Rxx_e-min(d)*eye(4,4);Z=zeros(3,1) eye(3,3);zeros(1,4);Cxy=Rxy_e-min(d)*Z;d=eig(Cxx,Cxy);d2=abs(abs(d)-1);w1,i1=min(d2);d2(i1)=100;w2,i2=min(d2);w1=angle(d(i1)/pi %估计角频率1w2=angle(d(i2)/pi %估计角频率2作业（3）用LMS算法设计一个二阶的自适应一步预测横向滤波器，滤波器的输入为二阶AR过程：，模型参数如下（分两组试验）：分别在步长和，的情况下，分析设计这个二阶的自

13、适应一步预测横向滤波器的LMS 算法：1) 画出由200次独立实验的平均得到的平均学习曲线：；2) 给出学习结束时平均的滤波器系数和，并与和比较；3) 比较平均的和；4) 比较平均的失调量和；5) 分析LMS算法的计算量；6) 思考：有收敛得更快、更准确的方法吗？程序：close all; clear; clc;N=512; p=2; a1=-0.975; a2=0.95; sigma_v=0.0731;%a1=-1.5955; a2=0.95; sigma_v=0.0322;u=0.05; %u=0.005;J_aver=zeros(1,N); w_aver=zeros(2,1); for

14、k=1:200 xtemp=zeros(N+p,1); vn=wgn(N+p,1,10*log10(sigma_v); %方差为sigma_v的高斯白噪声 for i=1:N xtemp(i+2)=vn(i)-a1*xtemp(i+1)-a2*xtemp(i); end x=xtemp(p+1:length(xtemp); %产生2阶AR信号 x=zeros(p,1);x; w_e=zeros(p,1); for i=1:N %LMS迭代 dn=x(i+p); un=x(i+p-1:-1:i); en(i)=dn-w_e'*un; w_e=w_e+u*un*en(i); end J_a

15、ver=J_aver+en.2; w_aver=w_aver+w_e;endJ_aver=J_aver./200;plot(J_aver); axis(0 500 0 0.6); %200次独立实验的平均学习曲线w_aver=w_aver./200 Jmin=sigma_vJex=J_aver(N)-JminA=1 a1 a2; a1 (1+a2) 0; a2 a1 1;b=sigma_v 0 0'r=inv(A)*b;R=r(1) r(2);r(2) r(1);d=eig(R);EJex=Jmin*(1+u*sum(d)/(2-u*sum(d)miu=Jex/JminEmiu=EJe

16、x/Jmin作业4：¡ 二阶AR模型: l u(n)+a1*u(n-1)+a2*u(n-2)=v(n)，n=1,2,3,120;¡ 针对如下两组参数l a1=-0.975; a2=0.95; v2=0.0731l a1=-1.5955; a2=0.95; v2=0.0322¡ 分别利用LMS算法 (step=0.05) 和Kalman算法 (Jmin0.005), 由u(n-1)和u(n-2)的线性组合, 即w(1, n)*u(n-1)+w(2,n)*u(n-2), 预测u(n), n=1,2,3,N2, N2=100;¡ 分别画出两种算法的J(n)曲线

17、(300次平均)¡ 分别给出w(1, n=N2)和w(2, n=N2)程序：close all; clear; clc;N=500; p=2; a1=-0.975; a2=0.95; sigma_v=0.0731;u=0.05; Jmin=0.005;Jlms=zeros(1,N); Jklm=zeros(1,N);for k=1:300 xtemp=zeros(N+p,1); vn=wgn(N+p,1,10*log10(sigma_v); %方差为sigma_v的高斯白噪声 for i=1:N xtemp(i+2)=vn(i)-a1*xtemp(i+1)-a2*xtemp(i);

18、end x=xtemp(p+1:length(xtemp); %产生2阶AR信号 x=zeros(2,1);x; wlms=zeros(2,1); %LMS 迭代 for i=1:N dn=x(i+p); un=x(i:i+p-1); en(i)=dn-wlms'*un; wlms=wlms+2*u*un*en(i); end Jlms=Jlms+en.2; wklm=zeros(2,1); %Kalman 迭代 Pn=eye(2,2); for i=1:N zn=x(i+p); un=x(i:i+p-1); Kn=Pn*un/(un'*Pn*un+Jmin); an=zn-u

19、n'*wklm; wklm=wklm+Kn*an; Pn=Pn-Kn*un'*Pn; en(i)=zn-wklm'*un; end Jklm=Jklm+en.2; endJlms=Jlms./300; Jklm=Jklm./300;plot(Jlms); hold on; plot(Jklm,'r');作业5：观测序列为二谐波与噪声叠加，为零均值高斯白噪声，令，由估计功率谱（用最大的两个奇异值和奇异向量计算）做三种情况令，利用FLP估计功率谱Case1:用RLS算法，求，Case2:用LMS算法，求，之后用求功率谱，要求画出100次平均的结果，计算即第2

20、00次的。并由估计功率谱。程序：close all; clear; clc;N=25; Pv=-10;%-10dB M=8;n=1:N; vn=wgn(1,N,Pv,'complex');un=exp(j*pi*(n-1)+exp(j*pi*(0.79.*n+1.2)+vn;AHf=zeros(M,N-M);bfH=un(M+1:N);AHb=zeros(M,N-M);bbH=un(1:N-M);for i=0:M-1 AHf(i+1,:)=un(M-i:N-1-i); AHb(i+1,:)=un(2+i:N-M+1+i);endAH=AHf conj(AHb); %FBLP输

21、入数据矩阵bH=bfH conj(bbH); Y,S,X = svd(AH'); %奇异值分解s=svd(AH');wls=X(:,1)*X(:,1)'*AH*bH'/s(1)2+X(:,2)*X(:,2)'*AH*bH'/s(2)2;ar=-conj(wls); %模型参数估计dw=2*pi/512;w=-pi:dw:pi-dw;Sw_e=zeros(1,length(w);n=1:8;for i=1:length(Sw_e); ewi=exp(-j*w(i)*n); Sw_e(i)=0.1/(abs(1+ewi*ar).2);endSw_e=

22、10*log10(Sw_e); Sw_e=Sw_e-max(Sw_e);plot(-pi:dw:pi-dw,Sw_e); title('基于FBLP的最小二乘奇异值分解谱估计');%RLSclear;N=200; Pv=-10;%-10dBM=8; deta=0.004; landa=1;vn=wgn(N,1,Pv,'complex'); n=1:N; n=n'xn=exp(j*pi*(n-1)+exp(j*pi*(0.79.*n+1.2)+vn;xn=zeros(M,1);xn;Pn=eye(M,M)/deta; wrls=zeros(M,1);for

23、 k=1:N %RLS迭代 un=xn(M+k-1:-1:k); dn=xn(M+k); kn=(1/landa)*Pn*un/(1+(1/landa)*un'*Pn*un); en=dn-wrls'*un; wrls=wrls+kn*conj(en); Pn=(1/landa)*Pn-(1/landa)*kn*un'*Pn; endar=-conj(wrls);dw=2*pi/512; w=-pi:dw:pi-dw;Sw_e=zeros(1,length(w); n=1:8;for i=1:length(Sw_e); ewi=exp(-j*w(i)*n); Sw_e(i)=0.1/(abs(1+ewi*ar).2);endSw_e=10*log10(Sw_e); Sw_e=Sw_e-max(Sw_e);fig