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

1、现代数字信号处理仿真作业1.仿真题 3.17仿真结果及图形：图 1 基于 FFT 的自相关函数计算图 3 周期图法和 BT 法估计信号的功率谱图 Error! Main Document Only. 基于式 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.0278810

2、22353997-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.033929397784

3、767 - 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=25

4、6; %周期图法， 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/

5、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)

6、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(归一化功率谱/

7、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

8、,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.0011

9、56047541561 + 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(-1

10、i*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 谱

11、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 算法中最接近单位元的两个特征值为：0.001826505974929 + 1.

12、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;

13、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

14、(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

15、_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

16、-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

17、;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;figu

18、re(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

19、:trialsrm=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

20、;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.仿

21、真题 5.10仿真结果及图形：(1) M=2 时， ，求解 Yule-Walker 方程210.9,.367a21(0)0ra可得到自相关矩阵 247.86.5734R相应的计算程序为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 方程210.9,0.9367a212()()102rra可得到自相关矩阵为 347.086.5734.58.120R相应的计算程序为r3=inv(1,-0.99,0;-0.99,1,0;0,-0.99,1)*0.93627;0;

22、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 时步长因子应在区

23、间（ 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.9

24、3627;for m=1:500v(:,m)=sqrt(sigma_v1_2)*randn(L,1);end% 生成500 组独立的 AR模型信号a1=-0.99;for m=1:500u(1,1,m)=v(1,m);for k=2:Lu(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:500e(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);

25、dd=(w(2,:,m)*uu);e(2,m)=u(3,m)-dd;for k=3:Lw(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:500w_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(

26、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:500w_mean=w_mean+w(:,:,m);e_mean=e_mean+e(:,m).2;endw

27、_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

28、 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)

29、;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+1A(:,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);

30、for k=1:Pfor m=1:Ma(m,k)=exp(-1i*omega(k)*(m-1);endend% 计算MVDR 谱Pmvdr=zeros(size(omega);for k=1:PPmvdr(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:PSw=zeros(1,M);for i=1:MSw(i)=(a(:,

31、k)*V(:,i)/S(i,i);endP_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 单次实验

32、估计权值以及 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

33、=length(b);A=zeros(M,L,trials);for m=1:trialsun1=zeros(M-1,1).,un(:,:,m).;for k=1:LA(:,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:trialsepsilon=zeros(L,1);P1=eye(M)/delta;for k=1:LPIn=P1*A(:,k,m); denok=lambda+A(:,k,m)*PIn; kn=PIn/denok;epsilon(k)=b

34、(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;endMSE(m,:)=abs(epsilon).2;endMSE_mean=zeros(1,L);w_mean=zeros(2,N);for m=1:trialsw_mean=w_mean+w(:,:,m);MSE_mean=MSE_mean+MSE(m,:);endw_mean=w_mean/trials; %500次独立实验权向量的均值MSE_mean=MSE_mean/trials; %

35、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=

36、1:100v(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:100for i=5:N2U(:,i,m)=u1(1,i-1,m);u1(1,i-2,m

37、);u1(1,i-3,m);u1(1,i-4,m);endendW_esti=zeros(4,N2,100);for m=1:100P_esti=1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;for l=1:N2-1P_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;ende

38、ndw=zeros(4,N2);e=zeros(1,N2,100);d=zeros(1,N2,100);MSE=zeros(1,N2);for m=1:100w=w+W_esti(:,:,m);for n=5:N2d(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);endMSE=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,:

39、),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 算

40、法得到的 DOA 估计为：-9.9716 39.7934图 22 MUSIC 算法实现 DOA 估计图 23 MVDR 算法实现 DOA 估计仿真程序（8_16）:clear allclc% 产生阵列接收信号N=100;M=10;K=2;theta=-10;40*pi/180;SNR=10;20;Am=sqrt(10.(SNR/10);S=Am*ones(1,N);S(2,:)=S(2,:).*exp(1i*2*pi*rand(1,N);for a=1:Mfor b=1:KA(a,b)=exp(-1i*(a-1)*pi*sin(theta(b);endendV=zeros(M,N);for m

41、=1:Mv=wgn(1,N,0,complex);v=v-mean(v);v=v/std(v);V(m,:)=v;endX=A*S+V;% 估计空间自相关矩阵R=zeros(M,M);for i=1:NR=R+X(:,i)*X(:,i);endR=R/N;% MUSIC算法 VR,D=eig(R);B,IX=sort(diag(D);G=VR(:,IX(M-K:-1:1);Pmusic=;for n=-pi/2:pi/180:pi/2a=exp(-1i*0:M-1*pi*sin(n);Pmusic=Pmusic,1/(a*G*G*a);end% MVDR算法Pmvdr=;for n=-pi/2

42、:pi/180:pi/2a=exp(-1i*0:M-1*pi*sin(n);Pmvdr=Pmvdr,1/(a*inv(R)*a);end% RootMUSIC算法syms zpz=z.(0:M-1);pz1=(z(-1).(0:M-1);fz=z(M-1)*pz1*G*G*pz;a=sym2poly(fz);r=roots(a);r1=abs(r);for i=1:2*KY,I(i)=min(abs(r1-1);r1(I(i)=inf;endfor i=1:2*Ktheta_esti_music(i)=asin(-angle(r(I(i)/pi)*180/pi;end% ESPRIT算法S=VR(:,IX(M:-1:M-K+1);S1=S(1:M-1,:);S2=S(2:M,:);fai=S1/S2;