现代数字信号处理仿真作业1.仿真题3.17仿真结果及图形：图基于FFT的自相关函数计算图图周期图法和BT法估计信号的功率谱图利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为：a1=--;a2=;a3=3i;a4=7;a5=68i;a6=7+6i;a7=9-2i;a8=2-0ia9=2+0i;a10=2+3i;a11=7-10i;a12=4-9i;a13=8-3i ;a14=2+4i;a15=2+1i;a16=3i.仿真程序（3_17）:clear allclc%%产生噪声序列N=32;%基于FFT的样本长度%N=256;%周期图法，BT法，AR模型功率谱估计的长度vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);%%产生复正弦信号f=[0.150.170.26];%归一化频率SNR=[303027];%信噪比A=10.^(SNR./20);%幅度signal=[A(1)*exp(1i*2*pi*f(1)*(0:N-1));%复正弦信号A(2)*exp(1i*2*pi*f(2)*(0:N-1));A(3)*exp(1i*2*pi*f(3)*(0:N-1))];%%产生观察样本un=sum(signal)+vn;%%利用3.1.1的FFT估计Uk=fft(un,2*N);Sk=(1/N)*abs(Uk).^2;r0=ifft(Sk);r1=[r0(N+2:2*N),r0(1:N)];%%r2=xcorr(un,N-1,'biased');%画图k=-N+1:N-1;figure(1)subplot(1,2,1)stem(k,real(r1))xlabel('m');ylabel('实部');subplot(1,2,2)stem(k,imag(r1))xlabel('m');ylabel('虚部');figure(2)subplot(1,2,1)stem(k,real(r2))xlabel('m');ylabel('实部');subplot(1,2,2)stem(k,imag(r2))xlabel('m');ylabel('虚部');%%周期图法NF=1024;Spr=fftshift((1/NF)*abs(fft(un,NF)).^2);kk=-0.5+(0:NF-1)*(1/(NF-1));Spr_norm=10*log10(abs(Spr)/max(abs(Spr)));%%BT法M=64;r3=xcorr(un,M,'biased');BT=fftshift(fft(r3,NF));BT_norm=10*log10(abs(BT)/max(abs(BT)));figure(3)subplot(1,2,1)plot(kk,Spr_norm)xlabel('w/2pi');ylabel('归一化功率谱/DB');title('周期图法')subplot(1,2,2)plot(kk,BT_norm)xlabel('w/2pi');ylabel('归一化功率谱/DB');title('BT法')%%LD迭代算法p=16;r0=xcorr(un,p,'biased');r4=r0(p+1:2*p+1);%计算自相关函数a(1,1)=-r4(2)/r4(1);sigma(1)=r4(1)-(abs(r4(2))^2)/r4(1);for m=2:p%LD迭代算法k(m)=-(r4(m+1)+sum(a(m-1,1:m-1).*r4(m:-1:2)))/sigma(m-1);a(m,m)=k(m);for i=1:m-1a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));endsigma(m)=sigma(m-1)*(1-abs(k(m))^2);endPar=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);%p阶AR模型的功率谱Par_norm=10*log10(abs(Par)/max(abs(Par)));figure(4)plot(kk,Par_norm)xlabel('w/2pi');ylabel('归一化功率谱/DB');title('16阶AR模型')2.仿真题3.20仿真结果及图形：单次Root-MUSIC算法中最接近单位圆的两个根为：+-对应的归一化频率为：相同信号的MUSIC谱估计结果如下图对3.20信号进行MUSIC谱估计的结果仿真程序（3_20）:clear allclc%%信号样本和高斯白噪声的产生N=1000;vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal=[exp(1i*0.5*pi*(0:N-1)+1i*2*pi*rand);%复正弦信号exp(-1i*0.3*pi*(0:N-1)+1i*2*pi*rand)];un=sum(signal)+vn;%%计算自相关矩阵M=8;for k=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=xs*xs'/(N-M);%%自相关矩阵的特征值分解[U,E]=svd(R);%%Root-MUSIC算法的实现G=U(:,3:M);Gr=G*G';co=zeros(2*M-1,1);for m=1:Mco(m:m+M-1)=co(m:m+M-1)+Gr(M:-1:1,m);endz=roots(co);ph=angle(z)/(2*pi);err=abs(abs(z)-1);%%计算MUSIC谱En=U(:,2+1:M);NF=2048;for n=1:NFAq=exp(-1i*2*pi*(-0.5+(n-1)/(NF-1))*(0:M-1)');Pmusic(n)=1/(Aq'*En*En'*Aq);endkk=-0.5+(0:NF-1)*(1/(NF-1));Pmusic_norm=10*log10(abs(Pmusic)/max(abs(Pmusic))); plot(kk,Pmusic_norm)xlabel('w/2*pi');ylabel('归一化功率谱/dB')3.仿真题3.21仿真结果及图形：单次ESPRIT算法中最接近单位元的两个特征值为：+-对应的归一化频率为：仿真程序（3_21）:clear allclc%%信号样本和高斯白噪声的产生N=1000;vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal=[exp(1i*0.5*pi*(0:N-1)+1i*2*pi*rand);%复正弦信号exp(-1i*0.3*pi*(0:N-1)+1i*2*pi*rand)];un=sum(signal)+vn;%%自相关矩阵的计算M=8;for k=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endRxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-M-1);Rxy=xs(:,1:end-1)*xs(:,2:end)'/(N-M-1);%%特征值分解[U,E]=svd(Rxx);ev=diag(E);emin=ev(end);Z=[zeros(M-1,1),eye(M-1);0,zeros(1,M-1)];Cxx=Rxx-emin*eye(M);Cxy=Rxy-emin*Z;%%广义特征值分解[U,E]=eig(Cxx,Cxy);z=diag(E);ph=angle(z)/(2*pi);err=abs(abs(z)-1);4.仿真题4.18仿真结果及图形：步长为0.05时失调参数为m1=0.0493；步长为0.005时失调参数为m2=0.0047。

图步长为0.05时权向量的收敛曲线图步长为0.005时权向量的收敛曲线图步长分别为0.05和0.005时100次独立实验的学习曲线仿真程序（4_18）:clear allclc%%产生100组独立样本序列data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;v=sqrt(sigma_v_2)*randn(data_len,1,trials);u0=[00];num=1;den=[1,a1,a2];Zi=filtic(num,den,u0);u=filter(num,den,v,Zi);%产生100组独立信号%%LMS迭代mu1=0.05;mu2=0.005;w1=zeros(2,data_len,trials);w2=w1;for m=1:100;temp=zeros(data_len,1);e1(:,:,m)=temp;e2(:,:,m)=temp;d1(:,:,m)=temp;d2(:,:,m)=temp;for n=3:data_len-1w1(:,n+1,m)=w1(:,n,m)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n,1,m));w2(:,n+1,m)=w2(:,n,m)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n,1,m));d1(n+1,1,m)=w1(:,n+1,m)'*u(n:-1:n-1,:,m);d2(n+1,1,m)=w2(:,n+1,m)'*u(n:-1:n-1,:,m);e1(n+1,1,m)=u(n+1,:,m)-d1(n+1,1,m);e2(n+1,1,m)=u(n+1,:,m)-d2(n+1,1,m);endendt=1:data_len;w1_mean=zeros(2,data_len);w2_mean=w1_mean;e1_mean=zeros(data_len,1);e2_mean=e 1_mean;for m=1:100w1_mean=w1_mean+w1(:,:,m);w2_mean=w2_mean+w2(:,:,m);e1_mean=e1_mean+e1(:,:,m).^2;e2_mean=e2_mean+e2(:,:,m).^2;endw1_mean=w1_mean/100;%100次独立实验权向量的均值w2_mean=w2_mean/100;e1_100trials_ave=e1_mean/100;%100次独立实验的学习曲线均值e2_100trials_ave=e2_mean/100;figure(1)plot(t,w1(1,:,90),t,w1(2,:,90),t,w1_mean(1,:),t,w1_mean(2,:))xlabel('迭代次数');ylabel('权向量')title('步长=0.05')figure(2)plot(t,w2(1,:,90),t,w2(2,:,90),t,w2_mean(1,:),t,w2_mean(2,:))xlabel('迭代次数');ylabel('权向量')title('步长=0.005')%%计算剩余误差和失调参数wopt=zeros(2,trials);Jmin=zeros(1,trials);sum_eig=zeros(trials,1);for m=1:trialsrm=xcorr(u(:,:,m),'biased');R=[rm(512),rm(513);rm(511),rm(512)];p=[rm(511);rm(510)];wopt(:,m)=R\p;[v,d]=eig(R);Jmin(m)=rm(512)-p'*wopt(:,m);sum_eig(m)=d(1,1)+d(2,2);endsJmin=sum(Jmin)/trials;Jex1=e1_100trials_ave-sJmin;%剩余均方误差mu1Jex2=e2_100trials_ave-sJmin;%剩余均方误差mu2sum_eig_100trials=sum(sum_eig)/100;Jexfin1=mu1*sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials));Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials));M1=Jexfin1/sJmin;%失调参数m1M2=Jexfin2/sJmin;%失调参数m2figure(3)plot(t,e1_100trials_ave,'*',t,e2_100trials_ave)xlabel('迭代次数');ylabel('均方误差')legend('u1=0.05','u2=0.005')axis([0,600,0,1])5.仿真题5.10仿真结果及图形：(1)M=2时，210.99,0.93627a σ=-=，求解Yule-Walker 方程可得到自相关矩阵相应的计算程序为r2=inv([1,-0.99;-0.99,1])*[0.93627;0];R2=[r2(1),r2(2);r2(2),r2(1)];%M=2(2)M=3时，2120.99,0,0.93627a a σ=-==，求解Yule-Walker 方程可得到自相关矩阵为相应的计算程序为r3=inv([1,-0.99,0;-0.99,1,0;0,-0.99,1])*[0.93627;0;0];R3=[r3(1),r3(2),r3(3);r3(2),r3(1),r3(2);r3(3),r3(2),r3(1)];%M=3(3)计算特征值扩展%%特征值分解eig_value_1=eig(R2);eig_value_2=eig(R3);%%特征值扩展eig_spread_1=max(eig_value_1)/min(eig_value_1);eig_spread_2=max(eig_value_2)/min(eig_value_2);M=2时特征值扩展是199.0000；M=3时特征值扩展是444.2790。