现代数字信号处理仿真作业Word格式文档下载.docx

1、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（mylabel（实部subplot（1,2,2）stem（k,imag（r1）虚部figure（2）stem（k,real（r2）stem（k,imag（r2）% 周期图法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（ab

2、s（Spr）;% BT法M=64;r3=xcorr（un,M,BT=fftshift（fft（r3,NF）;BT_norm=10*log10（abs（BT）/max（abs（BT）;figure（3）plot（kk,Spr_norm）w/2pi归一化功率谱/DBtitle（周期图法）plot（kk,BT_norm）BT法% LD迭代算法p=16;r0=xcorr（un,p,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）

3、+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）16阶AR模型2.仿真题3.20单次Root-MUSIC算

4、法中最接近单位圆的两个根为：对应的归一化频率为：相同信号的MUSIC谱估计结果如下图 5 对3.20信号进行MUSIC谱估计的结果仿真程序（3_20）:% 信号样本和高斯白噪声的产生N=1000;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:k）.;R=xs*xs/（N-M）;% 自相关矩阵的特征值分解U,E=svd（R）;% Root-MUSIC算法的实现G=

5、U（:,3:M）;Gr=G*Gco=zeros（2*M-1,1）;for m=1:M co（m:m+M-1）=co（m:m+M-1）+Gr（M:1,m）;z=roots（co）;ph=angle（z）/（2*pi）;err=abs（abs（z）-1）;% 计算MUSIC谱En=U（:,2+1: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）;Pmusic_norm=10*log10（abs（Pmusic）/max（abs（Pmusic）;plot（kk,Pmusic_n

6、orm）w/2*pi归一化功率谱/dB3.仿真题3.21单次ESPRIT算法中最接近单位元的两个特征值为：仿真程序（3_21）:% 自相关矩阵的计算Rxx=xs（:,1:end-1）*xs（:end-1）/（N-M-1）;Rxy=xs（:,2:end）% 特征值分解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）;4.仿真题4.18步长为0.05时失调参

7、数为m1=0.0493；步长为0.005时失调参数为m2=0.0047。图 6 步长为0.05时权向量的收敛曲线图 7 步长为0.005时权向量的收敛曲线图 8 步长分别为0.05和0.005时100次独立实验的学习曲线仿真程序（4_18）:% 产生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=f

8、ilter（num,den,v,Zi）; %产生100组独立信号% LMS迭代mu1=0.05;mu2=0.005;w1=zeros（2,data_len,trials）; w2=w1;100; temp=zeros（data_len,1）; e1（:,:,m）=temp;e2（:d1（:d2（: for n=3:data_len-1 w1（:,n+1,m）=w1（:,n,m）+mu1*u（n-1:n-2,:,m）*conj（e1（n,1,m）; w2（:,n+1,m）=w2（:,n,m）+mu2*u（n-1:,m）*conj（e2（n,1,m）; d1（n+1,1,m）=w1（:,n+1,m

9、）*u（n:n-1,:,m）; d2（n+1,1,m）=w2（: 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）;t=1:w1_mean=zeros（2,data_len）; w2_mean=w1_mean;e1_mean=zeros（data_len,1）;e2_mean=e1_mean;100 w1_mean=w1_mean+w1（: w2_mean=w2_mean+w2（: e1_mean=e1_mean+e1（:,m）.2; e2_mean=e2_mean+e2（:w1_mean=w1_mea

10、n/100; %100次独立实验权向量的均值w2_mean=w2_mean/100;e1_100trials_ave=e1_mean/100; %100次独立实验的学习曲线均值e2_100trials_ave=e2_mean/100;plot（t,w1（1,:,90）,t,w1（2,:,90）,t,w1_mean（1,:）,t,w1_mean（2,:）迭代次数权向量步长=0.05plot（t,w2（1,:,90）,t,w2（2,:,90）,t,w2_mean（1,:）,t,w2_mean（2,:） 步长=0.005% 计算剩余误差和失调参数wopt=zeros（2,trials）;Jmin=z

11、eros（1,trials）;sum_eig=zeros（trials,1）;trials rm=xcorr（u（:,m）, 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（: sum_eig（m）=d（1,1）+d（2,2）;sJmin=sum（Jmin）/trials;Jex1=e1_100trials_ave-sJmin; %剩余均方误差mu1Jex2=e2_100trials_ave-sJmin; %剩余均方误差mu2sum_ei

12、g_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; %失调参数m2plot（t,e1_100trials_ave,*,t,e2_100trials_ave） 均方误差legend（u1=0.05,u2=0.005axis（0,600,0,1）5.仿真题5.

13、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时，可得到自相关矩阵为r3=inv（1,-0.99,0;-0.99,1,0;0,-0.99,1）*0.93627;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）;% 特征值扩展ei

14、g_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步

15、长为0.0006时权向量收敛曲线图 11 500次独立实验M=2步长为0.001时的学习曲线图 12 500次独立实验M=3步长为0.0006时的学习曲线仿真程序（5_10_4）:% 产生系统输入白噪声L=10000;sigma_v1_2=0.93627;500v（:,m）=sqrt（sigma_v1_2）*randn（L,1）;% 生成500组独立的AR模型信号a1=-0.99; u（1,1,m）=v（1,m）; for k=2:L u（k,1,m）=-a1*u（k-1,1,m）+v（k,m）;% LMS迭代算法M=2;%M=3;mu=0.001;%mu=0.0006;w=zeros（L,M

16、,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: w（k,:,m）=w（k-1,:,m）+mu*e（k-1,m）*uu; uu=u（k-1,1,m） uu（1: dd=（w（k,: e（k,m）=u（k,m）-dd;% M=2e_mean=zeros（10000,1）;w_mean=zeros（10000,2）; w_mean=w_mean+w（: e_mean=e_mean+e

17、（:w_mean=w_mean/500;e_mean=e_mean/500;L;plot（t,w（:,1,100）,t,w（:,2,100）,t,w_mean（:,1）,t,w_mean（:,2）迭代次数n ylabel（抽头权值M=2,步长0.001的权向量收敛曲线plot（t,e_mean）MSEM=2,步长0.001的学习曲线% M=3w_mean=zeros（10000,3）;,2,100）,t,w（:,3,100）,t,w_mean（:,2）,t,w_mean（:,3）M=3,步长0.0006的权向量收敛曲线M=2,步长0.0006的学习曲线6.仿真题6.13滤波器抽头个数为4时图

18、13 M=4时的MVDR谱图 14 M=4时基于奇异值分解的MVDR谱从上面两图可以看出，M=4时并没有将3个频点分辨出来，增加滤波器阶数可以解决此问题，因此当M=20时仿真结果如下两图所示：图 15 M=20时的MVDR谱图 16 M=20时基于奇异值分解的MVDR谱仿真程序（6_13）:% 产生观测信号M=4;%M=20;f=0.1 0.25 0.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

19、）*randn（size（t）;Un=vn;length（f） s=Am（k）*exp（1i*2*pi*N*f（k）.*t+1i*phi（k）; Un=Un+s;Un=Un.% 构建矩阵A=zeros（M,N-M+1）;N-M+1 A（:,n）=Un（M+n-1:n）;U,S,V=svd（Ainvphi=V*inv（S*S）*V% 构建向量P=1024;f=linspace（-0.5,0.5,P）;omega=2*pi*f;a=zeros（M,P）;P for m=1: a（m,k）=exp（-1i*omega（k）*（m-1）;% 计算MVDR谱Pmvdr=zeros（size（omega）

20、; Pmvdr（k）=1/（a（:,k）*invphi*a（:,k）;Pmvdr=abs（Pmvdr/max（abs（Pmvdr）;Pmvdr=10*log10（Pmvdr）;P-1）*（1/（P-1）;plot（kk,Pmvdr）% 基于习题6.11的奇异值分解的MVDR方法 Sw=zeros（1,M）; Sw（i）=（a（:）*V（:,i）/S（i,i）; P_svd（k）=1/sum（abs（Sw）.2）;P_svd=abs（P_svd/max（abs（P_svd）;P_svd=10*log10（P_svd）;归一化频谱/dBM=4的MVDR谱%title（M=20的MVDR谱plot（

21、kk,P_svd）M=4的基于SVD的MVDR频谱M=20的基于SVD的MVDR频谱7.仿真题6.15图 17 单次实验估计权值以及500次独立实验的估计权值图 18 500次独立实验的学习曲线仿真程序（6_15）:% 产生AR模型的输入信号a1=0.99;sigma=0.995;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;b=un（n0+1:N,:L=length（b）;A=zeros（M,L,trials）; un1=zeros（M-1,1）.,un（:,m）. for k=1:,k,m）=un1（M-1+k:k）;% RLS算法求最优权向量delta=0.004;lam