现代数字信号处理仿真作业1.仿真题3.17仿真结果及图形：图 1 基于FFT的自相关函数计算图 3 周期图法和BT 法估计信号的功率谱图 2 基于式3.1.2的自相关函数的计算图 4 利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为：a1=-0.402637623107952-0.919787323662670i;a2=-0.013530139693503+0.024214641171318i;a3=-0.074241889634714-0.088834852915013i;a4=0.027881022353997-0.040734794506749i;a5=0.042128517350786+0.068932699075038i;a6=-0.0042799971761507 + 0.028686095385146i;a7=-0.048427890183189 - 0.019713457742372i;a8=0.0028768633718672 - 0.047990801912420ia9=0.023971346213842+ 0.046436389191530i;a10=0.026025963987732 + 0.046882756497113i;a11= -0.033929397784767 - 0.0053437929619510i;a12=0.0082735406293574 - 0.016133618316269i;a13=0.031893903622978 - 0.013709547028453i ;a14=0.0099274520678052 + 0.022233240051564i;a15=-0.0064643069578642 + 0.014130696335881i;a16=-0.061704614407581- 0.077423818476583i.仿真程序（3_17）:clear allclc%% 产生噪声序列N=32; %基于FFT的样本长度%N=256; %周期图法，BT法，AR模型功率谱估计的长度vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);%%产生复正弦信号f=[0.15 0.17 0.26]; %归一化频率SNR=[30 30 27]; %信噪比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)];%% 利用3.1.2估计Rr2=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算法中最接近单位圆的两个根为：-0.001156047541561 + 1.001503153449793i0.587376604261220 - 0.810845628739986i对应的归一化频率为：0.250183714447964-0.150223406926494相同信号的MUSIC谱估计结果如下图 5 对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')仿真结果及图形：单次ESPRIT算法中最接近单位元的两个特征值为：0.001826505974929 + 1.000690248438859i0.586994191014025 - 0.809491260856630i对应的归一化频率为：0.249709503383161-0.150146235268272仿真程序（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);仿真结果及图形：步长为0.05时失调参数为m1=0.0493；步长为0.005时失调参数为m2=0.0047。

图 6 步长为0.05时权向量的收敛曲线图 7 步长为0.005时权向量的收敛曲线图 8 步长分别为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=[0 0];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=e1_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 方程211(0)(1)(1)(0)0r r a r r σ⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可得到自相关矩阵247.048746.578346.578347.0487R ⎡⎤=⎢⎥⎣⎦相应的计算程序为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 方程212(0)(1)(2)1(1)(0)(1)0(2)(1)(0)0r r r r r r a r r r a σ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可得到自相关矩阵为347.048746.578346.112546.578347.048746.578346.112546.578347.0487R ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦相应的计算程序为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。