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

现代数字信号处理仿真作业1.仿真题3.17仿真结果及图形：图 1 基于FFT的自相关函数计算图 Error! Main Document Only. 基于式3.1.2的自相关函数的计算 图 3 周期图法和BT法估计信号的功率谱图 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-1 a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i); end sigma(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-M xs(:,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:M co(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:NF Aq=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算法中最接近单位元的两个特征值为：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-M xs(:,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。图 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-1 w1(:,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:100 w1_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:trials rm=xcorr(u(:,:,m),biased); R=rm(512),rm(513);rm(511),rm(512); p=rm(511);rm(510); wopt(:,m)=Rp; 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时， ，求解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时， ，求解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。(3) 根据LMS算法均方误差收敛特性，M=2时步长因子应在区间（0,0.0213）中，M=3时，步长因子应在区间（0,0.0142）之间，因此题中的步长因子不合理。故在仿真中，M=2时采用步长因子0.001，M=3时采用步长因子0.0006.图 9 500次独立实验M=2步长为0.001时权向量收敛曲线图 10 500次独立实验M=3步长为0.0006时权向量收敛曲线图 11 500次独立实验M=2步长为0.001时的学习曲线图 12 500次独立实验M=3步长为0.0006时的学习曲线仿真程序（5_10_4）:clear allclc% 产生系统输入白噪声L=10000;sigma_v1_2=0.93627;for m=1:500v(:,m)=sqrt(sigma_v1_2)*randn(L,1);end% 生成500组独立的AR模型信号a1=-0.99;for m=1:500 u(1,1,m)=v(1,m); for k=2:L u(k,1,m)=-a1*u(k-1,1,m)+v(k,m); endend% LMS迭代算法M=2;%M=3;mu=0.001;%mu=0.0006;w=zeros(L,M,500);for m=1:500 e(1,m)=u(1,m); uu=zeros(1,M); w(2,:,m)=w(1,:,m)+mu*e(1,m)*uu; uu=u(1,m) uu(1:M-1); dd=(w(2,:,m)*uu); e(2,m)=u(3,m)-dd; for k=3:L w(k,:,m)=w(k-1,:,m)+mu*e(k-1,m)*uu; uu=u(k-1,1,m) uu(1:M-1); dd=(w(k,:,m)*uu); e(k,m)=u(k,m)-dd; endend% M=2e_mean=zeros(10000,1);w_mean=zeros(10000,2);for m=1:500 w_mean=w_mean+w(:,:,m); e_mean=e_mean+e(:,m).2;endw_mean=w_mean/500;e_mean=e_mean/500;t=1:L;figure(1)plot(t,w(:,1,100),t,w(:,2,100),t,w_mean(:,1),t,w_mean(:,2)xlabel(迭代次数n); ylabel(抽头权值);title(M=2,步长0.001的权向量收敛曲线)figure(2)plot(t,e_mean)xlabel(迭代次数n); ylabel(MSE);title(M=2,步长0.001的学习曲线) % M=3e_mean=zeros(10000,1);w_mean=zeros(10000,3);for m=1:500 w_mean=w_mean+w(:,:,m); e_mean=e_mean+e(:,m).2;endw_mean=w_mean/500;e_mean=e_mean/500;t=1:L;figure(1)plot(t,w(:,1,100),t,w(:,2,100),t,w(:,3,100),t,w_mean(:,1),t,w_mean(:,2),t,w_mean(:,3)xlabel(迭代次数n); ylabel(抽头权值);title(M=3,步长0.0006的权向量收敛曲线)figure(2)plot(t,e_mean)xlabel(迭代次数n); ylabel(MSE);title(M=2,步长0.0006的学习曲线)6.仿真题6.13仿真结果及图形：滤波器抽头个数为4时图 13 M=4时的MVDR谱图 14 M=4时基于奇异值分解的MVDR谱从上面两图可以看出，M=4时并没有将3个频点分辨出来，增加滤波器阶数可以解决此问题，因此当M=20时仿真结果如下两图所示：图 15 M=20时的MVDR谱图 16 M=20时基于奇异值分解的MVDR谱仿真程序（6_13）:clear allclc% 产生观测信号M=4;%M=20;N=1000;f=0.1 0.25 0.27;SNR=30 30 27;sigma=1;Am=sqrt(sigma*10.(SNR/10);t=linspace(0,1,N);phi=2*pi*rand(size(f);vn=sqrt(sigma/2)*randn(size(t)+1i*sqrt(sigma/2)*randn(size(t);Un=vn;for k=1:length(f) s=Am(k)*exp(1i*2*pi*N*f(k).*t+1i*phi(k); Un=Un+s;endUn=Un.;% 构建矩阵A=zeros(M,N-M+1);for n=1:N-M+1 A(:,n)=Un(M+n-1:-1:n);endU,S,V=svd(A);invphi=V*inv(S*S)*V;% 构建向量P=1024;f=linspace(-0.5,0.5,P);omega=2*pi*f;a=zeros(M,P);for k=1:P for m=1:M a(m,k)=exp(-1i*omega(k)*(m-1); endend% 计算MVDR谱Pmvdr=zeros(size(omega);for k=1:P Pmvdr(k)=1/(a(:,k)*invphi*a(:,k);endPmvdr=abs(Pmvdr/max(abs(Pmvdr);Pmvdr=10*log10(Pmvdr);kk=-0.5+(0:P-1)*(1/(P-1);figure(1)plot(kk,Pmvdr)% 基于习题6.11的奇异值分解的MVDR方法for k=1:P Sw=zeros(1,M); for i=1:M Sw(i)=(a(:,k)*V(:,i)/S(i,i); end P_svd(k)=1/sum(abs(Sw).2);endP_svd=abs(P_svd/max(abs(P_svd);P_svd=10*log10(P_svd);xlabel(w/2*pi);ylabel(归一化频谱/dB)title(M=4的MVDR谱)%title(M=20的MVDR谱)figure(2)plot(kk,P_svd)xlabel(w/2*pi);ylabel(归一化频谱/dB)title(M=4的基于SVD的MVDR频谱)%title(M=20的基于SVD的MVDR频谱)7.仿真题6.15仿真结果及图形：图 17 单次实验估计权值以及500次独立实验的估计权值图 18 500次独立实验的学习曲线仿真程序（6_15）:clear allclc% 产生AR模型的输入信号a1=0.99;sigma=0.995;N=1000;trials=500;vn=sqrt(sigma)*randn(N,1,trials);nume=1;deno=1 a1;u0=zeros(length(deno)-1,1);xic=filtic(nume,deno,u0);un=filter(nume,deno,vn,xic); %产生500组独立的信号% 产生期望信号和观测数据矩阵n0=1;M=2;b=un(n0+1:N,:,:);L=length(b);A=zeros(M,L,trials);for m=1:trials un1=zeros(M-1,1).,un(:,:,m).; for k=1:L A(:,k,m)=un1(M-1+k:-1:k); endend% RLS算法求最优权向量delta=0.004;lambda=0.98;w=zeros(M,L+1,trials);for m=1:trials epsilon=zeros(L,1); P1=eye(M)/delta; for k=1:L PIn=P1*A(:,k,m); denok=lambda+A(:,k,m)*PIn; kn=PIn/denok; epsilon(k)=b(k,1,m)-w(:,k,m)*A(:,k,m); w(:,k+1,m)=w(:,k,m)+kn*conj(epsilon(k); P1=P1/lambda-kn*A(:,k,m)*P1/lambda; end MSE(m,:)=abs(epsilon).2;endMSE_mean=zeros(1,L);w_mean=zeros(2,N);for m=1:trials w_mean=w_mean+w(:,:,m); MSE_mean=MSE_mean+MSE(m,:);endw_mean=w_mean/trials; %500次独立实验权向量的均值MSE_mean=MSE_mean/trials; %500次独立实验的MSEt=1:1000;figure(1)plot(t,w(1,:,100),t,w(2,:,100),t,w_mean(1,:),t,w_mean(2,:)xlabel(迭代次数n);ylabel(权值)figure(2)plot(1:L,MSE_mean)xlabel(迭代次数n);ylabel(MSE)8.仿真题6.15仿真结果及图形：图 19 100次独立实验权值变化曲线图 20 单次实验权值变化曲线图 21 100次独立实验MSE取对数的变化曲线仿真程序（7_14）:clear allclc% 产生白噪声N=2000;gv=0.0332;for m=1:100 v(m,:)=randn(1,N)*sqrt(gv);end% 产生AR模型序列a=1.6 -1.46 0.616 -0.1525;u1=zeros(1,N,100);for m=1:100 % 产生100组独立实验信号 for i=1:(N-4) u1(1,i+4,m)=a(1)*u1(1,i+3,m)+a(2)*u1(1,i+2,m)+a(3)*u1(1,i+1,m)+a(4)*u1(1,i,m)+v(m,i+4); endend% 卡尔曼滤波N2=2000;Jmin=0.005;for m=1:100 for i=5:N2 U(:,i,m)=u1(1,i-1,m);u1(1,i-2,m);u1(1,i-3,m);u1(1,i-4,m); endendW_esti=zeros(4,N2,100);for m=1:100 P_esti=1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1; for l=1:N2-1 P_pre=P_esti; A=(U(:,l,m)*P_pre*U(:,l,m)+Jmin; K=P_pre*U(:,l,m)/A; alpha(l)=u1(1,l,m)-(U(:,l,m)*W_esti(:,l,m); W_esti(:,l+1,m)=W_esti(:,l,m)+K*alpha(l); P_esti=P_pre-K*(U(:,l,m)*P_pre; endendw=zeros(4,N2);e=zeros(1,N2,100);d=zeros(1,N2,100);MSE=zeros(1,N2);for m=1:100 w=w+W_esti(:,:,m); for n=5:N2 d(1,n,m)=W_esti(:,n,m)*u1(1,n-1:-1:n-4,m); e(1,n,m)=u1(1,n,m)-d(1,n,m); end MSE=MSE+e(:,:,m).2;endw=w/100; %100次独立实验的权向量均值MSE=MSE/100; % 100次独立实验的均方误差t=1:N2;figure(1)plot(t,w(1,:),t,w(2,:),t,w(3,:),t,w(4,:)legend(w1,w2,w3,w4)xlabel(迭代次数);ylabel(权值)figure(2)plot(t,W_esti(1,:,50),t,W_esti(2,:,50),t,W_esti(3,:,50),t,W_esti(4,:,50)legend(w1,w2,w3,w4);xlabel(迭代次数);ylabel(权值)figure(3)semilogy(t,MSE)xlabel(迭代次数);ylabel(对数MSE)9.仿真题8.16仿真结果及图形：单次RootMUSIC算法得到的DOA估计为：-9.9938 39.9904单次ESPRIT算法得到的DOA估计为：-9.9716 3