現(xiàn)代數(shù)字信號處理仿真作業(yè)[共50頁]

1、實用標準文檔現(xiàn)代數(shù)字信號處理仿真作業(yè)1.仿真題3.17仿真結(jié)果及圖形：圖 1 基于FFT的自相關(guān)函數(shù)計算圖 Error! Main Document Only. 基于式3.1.2的自相關(guān)函數(shù)的計算 圖 3 周期圖法和BT法估計信號的功率譜圖 4 利用LD迭代對16階AR模型的功率譜估計16階AR模型的系數(shù)為：a1=-0.402637623107952-0.919787323662670i; a2=-0.013530139693503+0.024214641171318i;a3=-0.074241889634714-0.088834852915013i;a4=0.027881022353997-

2、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.0

3、053437929619510i;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% 產(chǎn)生噪聲序列N=32; %基于FFT的樣本長度%N=256; %周期圖法，BT

4、法，AR模型功率譜估計的長度vn=(randn(1,N)+1i*randn(1,N)/sqrt(2);%產(chǎn)生復正弦信號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);% 產(chǎn)生觀察樣本 un=sum(signal)+vn;% 利用3.1.1的FFT估計Uk=fft(un,2*N);Sk=(1/N)*abs(Uk).

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

6、l(虛部); % 周期圖法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(周期圖法)

7、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); %計算自相關(guān)函數(shù)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,

8、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仿真結(jié)果及圖形：單次Root-MUSIC算法中最接近單位圓的兩個根為：-0.001156047541561

9、+ 1.001503153449793i0.587376604261220 - 0.810845628739986i對應的歸一化頻率為：0.250183714447964-0.150223406926494相同信號的MUSIC譜估計結(jié)果如下圖 5 對3.20信號進行MUSIC譜估計的結(jié)果仿真程序（3_20）:clear allclc% 信號樣本和高斯白噪聲的產(chǎn)生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)+

10、1i*2*pi*rand);un=sum(signal)+vn;% 計算自相關(guān)矩陣M=8;for k=1:N-M xs(:,k)=un(k+M-1:-1:k).;endR=xs*xs/(N-M);% 自相關(guān)矩陣的特征值分解U,E=svd(R);% Root-MUSIC算法的實現(xiàn)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

11、=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仿真結(jié)果及圖形：單次ESPRIT算法中最接近單位元的兩個特征值為：0.001826505974929 + 1.000690248438859i0.

12、586994191014025 - 0.809491260856630i對應的歸一化頻率為：0.249709503383161-0.150146235268272仿真程序（3_21）:clear allclc% 信號樣本和高斯白噪聲的產(chǎn)生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;% 自相關(guān)矩陣的計算M=8;for k=1:N-M xs(:,

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

14、8仿真結(jié)果及圖形：步長為0.05時失調(diào)參數(shù)為m1=0.0493；步長為0.005時失調(diào)參數(shù)為m2=0.0047。圖 6 步長為0.05時權(quán)向量的收斂曲線圖 7 步長為0.005時權(quán)向量的收斂曲線圖 8 步長分別為0.05和0.005時100次獨立實驗的學習曲線仿真程序（4_18）:clear allclc% 產(chǎn)生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

15、,a1,a2;Zi=filtic(num,den,u0);u=filter(num,den,v,Zi); %產(chǎn)生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)

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

17、,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次獨立實驗權(quán)向量的均值w2_mean=w2_mean/100;e1_100trials_ave=e1_mean/100; %100次獨立實驗的學習曲線均值e2_100trials_ave=e2_mean/100;figure(1)plot(t,w1(1,:,

18、90),t,w1(2,:,90),t,w1_mean(1,:),t,w1_mean(2,:)xlabel(迭代次數(shù));ylabel(權(quán)向量)title(步長=0.05)figure(2)plot(t,w2(1,:,90),t,w2(2,:,90),t,w2_mean(1,:),t,w2_mean(2,:) xlabel(迭代次數(shù));ylabel(權(quán)向量)title(步長=0.005) % 計算剩余誤差和失調(diào)參數(shù)wopt=zeros(2,trials);Jmin=zeros(1,trials);sum_eig=zeros(trials,1);for m=1:trials rm=xcorr(u(:

19、,:,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*

20、sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials);Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials);M1=Jexfin1/sJmin; %失調(diào)參數(shù)m1M2=Jexfin2/sJmin; %失調(diào)參數(shù)m2figure(3)plot(t,e1_100trials_ave,*,t,e2_100trials_ave) xlabel(迭代次數(shù));ylabel(均方誤差)legend(u1=0.05,u2=0.005)axis(0,600,0,1)5.仿真題5.10仿真結(jié)果及圖形：

21、(1) M=2時， ，求解Yule-Walker方程可得到自相關(guān)矩陣相應的計算程序為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方程可得到自相關(guān)矩陣為相應的計算程序為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(R

22、2);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) 根據(jù)LMS算法均方誤差收斂特性，M=2時步長因子應在區(qū)間（0,0.0213）中，M=3時，步長因子應在區(qū)間（0,0.0142）之間，因此題中的步長因子不合理。故在仿真中，M=2時采用步長因子0.001，M=3時采用步長因子0.0006.圖 9 500次獨立實驗M=2步

23、長為0.001時權(quán)向量收斂曲線圖 10 500次獨立實驗M=3步長為0.0006時權(quán)向量收斂曲線圖 11 500次獨立實驗M=2步長為0.001時的學習曲線圖 12 500次獨立實驗M=3步長為0.0006時的學習曲線仿真程序（5_10_4）:clear allclc% 產(chǎn)生系統(tǒng)輸入白噪聲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

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

25、(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(迭代次數(shù)n); ylabel(抽頭權(quán)值);title(M=2,步長0.001的權(quán)向量收斂曲線)fig

26、ure(2)plot(t,e_mean)xlabel(迭代次數(shù)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(

27、:,2),t,w_mean(:,3)xlabel(迭代次數(shù)n); ylabel(抽頭權(quán)值);title(M=3,步長0.0006的權(quán)向量收斂曲線)figure(2)plot(t,e_mean)xlabel(迭代次數(shù)n); ylabel(MSE);title(M=2,步長0.0006的學習曲線)6.仿真題6.13仿真結(jié)果及圖形：濾波器抽頭個數(shù)為4時圖 13 M=4時的MVDR譜圖 14 M=4時基于奇異值分解的MVDR譜從上面兩圖可以看出，M=4時并沒有將3個頻點分辨出來，增加濾波器階數(shù)可以解決此問題，因此當M=20時仿真結(jié)果如下兩圖所示：圖 15 M=20時的MVDR譜圖 16 M=20時基于

28、奇異值分解的MVDR譜仿真程序（6_13）:clear allclc% 產(chǎn)生觀測信號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;end

29、Un=Un.;% 構(gòu)建矩陣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;% 構(gòu)建向量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);end

30、Pmvdr=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(歸一化頻

31、譜/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仿真結(jié)果及圖形：圖 17 單次實驗估計權(quán)值以及500次獨立實驗的估計權(quán)值圖 18 500次獨立實驗的學習曲線仿真程序（6_15）:clear allclc% 產(chǎn)生AR模型的輸入信號a1=0.99;sigma=0.995;N=1000;trials=500;vn=sqrt(sigma)*randn(N

32、,1,trials);nume=1;deno=1 a1;u0=zeros(length(deno)-1,1);xic=filtic(nume,deno,u0);un=filter(nume,deno,vn,xic); %產(chǎn)生500組獨立的信號% 產(chǎn)生期望信號和觀測數(shù)據(jù)矩陣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算法求最優(yōu)權(quán)向量delta=

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

34、.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次獨立實驗權(quán)向量的均值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(迭代次數(shù)n);ylabel(權(quán)值)figure(2)pl

35、ot(1:L,MSE_mean)xlabel(迭代次數(shù)n);ylabel(MSE)8.仿真題6.15仿真結(jié)果及圖形：圖 19 100次獨立實驗權(quán)值變化曲線圖 20 單次實驗權(quán)值變化曲線圖 21 100次獨立實驗MSE取對數(shù)的變化曲線仿真程序（7_14）:clear allclc% 產(chǎn)生白噪聲N=2000;gv=0.0332;for m=1:100 v(m,:)=randn(1,N)*sqrt(gv);end% 產(chǎn)生AR模型序列a=1.6 -1.46 0.616 -0.1525;u1=zeros(1,N,100);for m=1:100 % 產(chǎn)生100組獨立實驗信號 for i=1:(N-4)

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

37、=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,

38、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次獨立實驗的權(quán)向量均值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(迭代次數(shù));ylabel(權(quán)值)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(迭代次數(shù));ylabel(權(quán)值)figure(3)semilogy(t,MSE)xlabel(迭代次數(shù));ylabel(對數(shù)MSE)9.仿真題8.16仿真結(jié)果及圖形：單次RootMUSIC算法得到的DOA估計為：-9.9938 39.9904單次ESPRIT算法得到的DOA估計為：-9.9716 39.7