现代数字信号处理仿真作业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）A2）;kk=-0.5+（0:NF-1）*（1/（NF-1）;Spr_norm=10*log10（abs（Spr）/ma

2、x（abs（Spr）;%B祛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/2piylabel（ title（ 周期图法 ）subplot（1,2,2） plot（kk,BT_norm）BT 法 ） %LDt代算法归一化功率谱 /DBp=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）F2）/r4（1

3、）;for m=2:p %LD1代算法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）F2）;Par=sigma（p）./fftshift（abs（fft（1,a（p,:）,NF）.A2）;Par_norm=10*log10（abs（Par）/max（abs（Par）;figure（4）plot（kk,Par_norm） 归一化功率谱 /DBti

4、tle（ 16 阶 AR模型）%p阶AR模型的功率谱2.仿真题 3.20单次Root-MUSIC算法中最接近单位圆的两个根为:+对应的归一化频率为：相同信号的MUSIC谱估计结果如下图对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）;%计算自相关矩阵M=8;for k=1:N-Mxs（:,k）=un（k+M-1:k）.;R=xs*xs/（N-M）;%自相关矩阵的特征

5、值分解U,E=svd（R）;%Root-MUSIC 算法的实现G=U（:,3:M）;Gr=G*Gco=zeros（2*M-1,1）;for m=1:Mco（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:NFAq=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）/ma

6、x（abs（Pmusic）;plot（kk,Pmusic_norm）w/2*pi 归一化功率谱 /dB3.仿真题 3.21单次ESPRITT法中最接近单位元的两个特征值为：-仿真程序（ 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,Cx

7、y）;z=diag（E）;4.仿真题 4.18步长为 0.05时失调参数为 m1=0.0493；步长为 0.005时失调参数为 m2=0.0047。图步长为 0.05时权向量的收敛曲线图步长为 0.005时权向量的收敛曲线图步长分别为 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=00;num=1;d

8、en=1,a1,a2;Zi=filtic（num,den,u0）;u=filter（num,den,v,Zi）; %产生 100 组独立信号%LM链代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-1w1（:,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（

9、n,1,m）;d1（n+1,1,m）=w1（:,n+1,m）*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=e 1_mean;100w1_mean=w1_mean+w1（:w2_mean=w2_mean+w2（:e1_mean=e1_mean+e1（:,m）.A2;e2_mean=e2_

10、mean+e2（:w1_mean=w1_mean/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%计算剩余误差和失调参数w

11、opt=zeros（2,trials）;Jmin=zeros（1,trials）;sum_eig=zeros（trials,1）;trialsrm=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）; endsJmin=sum（Jmin）/trials;Jex1=e1_100trials_ave-sJmin; %余均方误差 mu1Jex2=e2_100trials_ave-s

12、Jmin; %余均方误差 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; %调参数 mlM2=Jexfin2/sJmin; %调参数 m2plot（t,e1_100trials_ave, * ,t,e2_100trials_ave）xlabel（迭代次数均方误差） legend（ u

13、1=0.05 , u2=0.005axis（0,600,0,1）5.仿真题5.10（1）M=2 时，a =-0.99,a2 =0.93627 ,求解 Yule-Walker方程可得到自相关矩阵 相应的计算程序为 r2=inv（1,-0.99;-0.99,1）*0.93627;0;R2=r2（1）,r2;r2（2）,r2（1）; %M=2（2）M=3 时，a1 = -0.99,a2 =0户2 =0.93627 ,求解 Yule-Walker方程可得到自相关矩阵为 相应的计算程序为 r3=inv（1,-0.99,0;-0.99,1,0;0,-0.99,1）*0.93627;0;R3=r3（1）,r

14、3（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。根据LMS算法均方误差收敛特性，M=2时步长因子应在区间（0,0.0213）中，M=3时, 步长因子应在区间

15、（0,0.0142）之间，因此题中的步长因子不合理。故在仿真中， M=2时采用步长因子0.001, M=3时采用步长因子0.0006.图 500 次独立实验 M=2 步长为 0.001 时权向量收敛曲线图 500 次独立实验 M=3 步长为 0.0006时权向量收敛曲线图 500 次独立实验 M=2 步长为 0.001 时的学习曲线图 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组独立的A

16、R模型信号a1=-0.99;500 u（1,1,m）=v（1,m）;for k=2:Lu（k,1,m）=-a1*u（k-1,1,m）+v（k,m）;%LM链代算法M=2;%M=3;mu=0.001;%mu=0.0006;w=zeros（L,M,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:d

17、d=（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（:,m2;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 抽头权值 M=2, 步长 0.001 的权向量收敛曲线 figure（2） plot（t,e_mean）MSEM=2, 步长 0.001 的学习曲线 %M=3w_mean

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

19、.250.27;sigma=1;Am=sqrt（sigma*10.A（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;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+1A（:,n）=Un（M+n-1:n）;U,S,V=svd（Ainvphi=V*inv（S*S）*V%构建向量P=102

20、4;f=linspace（-0.5,0.5,P）;omega=2*pi*f;a=zeros（M,P）;Pa（m,k）=exp（-1i*omega（k）*（m-1）;%刻算 MVDRlfPmvdr=zeros（size（omega）;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的奇异值分解的 MVD昉法Sw=zeros（1,M）;Sw（i）=（a（:）*V（:,i）/S（i,i）;P_svd（k）

21、=1/sum（abs（Sw）-2）;P_svd=abs（P_svd/max（abs（P_svd）;P_svd=10*log10（P_svd）; 归一化频谱 /dBM=4 的 MVDRif）%title（M=20 的 MVDRfplot（kk,P_svd）M=4 的基于 SVD的 MVDR谱M=20 的基于 SVD 的 MVD砌谱7.仿真题 6.15图单次实验估计权值以及 500 次独立实验的估计权值图 500 次独立实验的学习曲线仿真程序（ 6_15） :%花生AR模型的输入信号a1=0.99;sigma=0.995;trials=500;vn=sqrt（sigma）*randn（N,1,trials）;nume=1;deno=1a1;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