现代信号处理大作业-1217.docx
现代信号处理大作业姚苏洋11303491711.1 .使用matlab实现LD迭代算法1.2 1.evinson-Durbin算法功率谱估计大致可以分为经典谱估计和现代功率谱估计，经典谱估计方法存在着以下几点缺陷：A) .数据加窗或自相关加窗，都隐含着假定在窗外未观测到的数据或自相关系数为零，该假设不切实际。B) .要性能好往往需要较长的数据，但实际数据长度有限.窗函数容易造成谱的模糊。采用AR模型的现代谱估计方法可以克服这些不足。其中LD递推算法可以在计算机上方便实现。1.D递推算法具体计算步骤如下：(1)Yule-Walker方程的矩阵形式所示：Ho)L(-1)心(-2)a(-)12噎Mo)噎)rt(-+l)ak.0噎(Z)MkT)H"2)噎(0)_au0系数矩阵尺I=EV,为Hermitian矩阵，对角线上元素相同，即为ToPIieZ矩阵。(2)P-I阶Yule-Walker方程为:K(O)%K(P-I)4(-1)Rx(-p)K(O)&(2-P)RX(P-2)RA(O)1Vp1=-2-叫OO其中，(3=EejW)为误差功率。写成联立方程：P-IX*R(m-&)=<A:=0Y-0,m二pw=01,-,/?-1取共枕得：TAR：(*幻=<Jt=O0,n=,m=O1,p-l(4)变量替换,并利用R<)=R：(。得:PTE*iR(m-k)=<k=02pm=p-0,/?=0,p-2(5)表示成矩阵：-&(。)&Q)K(T)4RX(I-P)-RXQ-P)-pP-1VpP2-00(6)求解得：R(P-I)RX(P-2)6(0)_I.1.ap.k=aP-U+KPdiP_k，k=O,、pp=M+KQ；)(8)O=G+"3KP(10)(U)*=p-+K/bK；%二%I1-区1(3) 当k=l时，即一阶递推为:RX(O)J(T)T1&RX(O)N求解可得:o = LRA(12)(13)(14)(15)(4) 对于P2时，递推为:册,。三1，%=Qp"+Kp0g,=1l-Cp2r*p=ap,p=-一'"=R(p)+EaPMRX(P-Qp-A=I矩阵RX已知，可得到各阶AR模型系数为：4二一焉，P1=(l-kl)2)rtx(O)i(Q+£-(/)G伏-。ak(k)=J=(16)PkTPk-ak(Z)=ak(Z)÷ak(k)ak_x(k-i)i=1,2,A-1(17)Pk=(IT4(幻以T(18)1.2实验结果假设P=5,使用matlab求得递推结果矩阵A如图1所示。画出每次迭代矩阵A的对应值如图2所示。A=图2多次迭代输出1.000000001.00000-2.0000001.0000-0.6667-2.00000.333301.0000-0.7917-1.25000.5833-0.3750图1LD算法得到递推结果A2.令信号(f) = ,ej2flt ej2tej2afy0<<T4T4<t<T2由三个不同频率的复正弦信号首尾相连而成。T2<t<T其中，工=4月),人=2工),力=工)。(1)试求(t)的WV分布，并画出三维WV分布图(2)指出并分析其WV分布的信号项和交叉项。2.1 原理分析已知WV分布公式如(19)所示：+OOW(z,/)=z(+*(r-2Jr(19)-OO(20)Mf)=e"""(/)+eyw()÷ej2h,u3(t)求得X的WV分布为:÷xW(Zj)=J=2而(如用)*(/1(。+2雨(妙。2)*%(/)+2昉(。-3)*，3")+4点8$(姐-。2)6(-色记)*(712")+ 4%cos(6¾-(j)(<-码)*U式r)+4/rcos(牡外)6("；')*U2i(Z)=火力*卬。+对九)*%")+的力)*(/向+2侬2乃(/；九)须，-写&他2(。+ 2cos(iiW-A)*U)+2cos优工)收/-A*%*)其中：U(r)=%(f)*%(f)CZ2(Z)=M2(Z)*M2(r)U3«)="3«)*%。)C12(r)=M1(0*w2(0U3(f)=%")*%3U23(t)=u2(t)*uy(t)2.2 实验结果使用Matlab时频分析工具箱，可以很方便地帮助我们画出WV分布以及信号项，交叉项。将信号看作三个分段信号的叠加，分别用信号1,2,3表示，其结果分别如图3至图9所示。图3信号1信号项图4信号2信号项3.一非平稳信号由两个高斯信号叠加而成：/、“4z.(a.a,.2.2(0=IIexp(-3(/TO)+/%)分别求出z(t)的WV分布及模糊函数，画出二者的波形图，指出并分析其信号项和交叉项。3.1 原理分析根据WV分布的定义式，可以得到z(r)的WV分布为;(26)(27)(28)所以z(力的WV分布的信号项为:II7Zrr一("%/一左)2()=2e°(24)交叉项为：WzJtJ)=2/"告一",(25)叫w*J)=2Re叫&J)-(t(-lm)2ffm)2r.f.r->c/r=2eacos4(f-fmytd+2fdt其中乙=空，=咛万a=。，力=()所以z(t)的WV分布为：叱。J)=%*j)+叱B(FJ)-(r-)2-T-<-o)2-(m)-(/-A,)22_x>,、=2e。+2eos4-(7-fm)td+2ftlt(1)由模糊函数的定义：4A(r,v)=+：”*("；)"2如/<22可以计算得出z】(t)的模糊函数为:+OOAz(r,v)= z(r + -)z,t(r-)ey2xnWr¢0f÷*( a -(z÷)2÷<¾(f÷)=L1/2Ct y+<o -(-fo)*r + j例 r+2rw e 4dt-<o(29)1/2.7TV Ct 2 b* 2f ÷oo flt(l-F-y-)* L>* ÷ j(-lt?)J e a 4 0rdt(Qf Y 7 r2 u2+(ftr-2zrrwu)二6七©2Er 2 一二-4产 +y(<r-2,0t>)=e 4 c模糊函数信号项为:2-T22+j(r-2f0)4(工，射)=e4 '(30)交叉项为:AcmSSG，)= AI Zl (工,V)-(2>+<w=e a.eJ<m+2tm+dt,ny(31)其中，q=?，/=制为，C=0,0=（）所以可得z(t)的模糊函数为：A(Z, u) = Aaiito(, V)+ Aenfsx(, v)_ pK一 曝八 W2 飞-i(2w)2-(D2-(Wmr+2y)- CVZI VZC(32)3.2实验结果4003200m O图10信号IWV分布信号项图11信号2WV分布交叉项如图12至图17所示，分别为WV分布信号项，交叉项，模糊函数，信号项，模糊函数交叉项。图16信号1信号项图17信号2信号项附：实验代码1LD算法%PNsequencegenerationfbconnection=1000010011;n=length(fbconnection);N=2n-1;register=zeros(1,n-l)1;mseq(1)=register(n);fori=2:Nnewregister(1)=mod(sum(fbconnection.*register),2);forj=2:nnewregister(j)=register(j-l);end;register=newregister;mseq(i)=register(n);end%initializationrxx=abs(xcorr(mseq(1:6);P=5;Rxx_0=rxx(6);Rxx_ii=rxx(1:5);Rxx_jj=rxx(7:11);Rxx=zeros(p,p);fori=1:pforj=1:pifi=jRxx(i,j)=Rxx_0;elseifi<jRxx(i,j)=Rxx_ii(l+j-i);elseifi>jRxx(i,j)=Rxx_jj(l+i-j);endendend%LDalgorithmP=p-i；A=zeros(p,p+1);forloop=1:pifloop=1P=1；a_pO=1;a_pl=-Rxx(2,1)/Rxx(1,1);sigma=Rxx(1,1)+a_pl*Rxx(2,1);A(loop,loop)=a_pO;A(loop,loop+l)=a_pl;elsepre_sum=0;fork=1:loop-1pre_sum=pre_sum+A(loop-1,k)*Rxx(1,loop-k);enddelta_p=Rxx(loop,loop)+pre_sum;K=-delta_p/sigma;sigma=sigma*(l-abs(K).2);fork=1:loop+1a_pk=A(loop-1,k)+K*conj(A(loop-1,loop-k+2);A(loop,k)=a_pk;endendfigure;stem(A(loop,:);end2WV分布wl=400;w2=200;w3=100;t=1:1024;xl=exp(1j*2*pi*wl*t(1:256)/1024);x2=exp(Ij*2*pi*w2*t(257:512)/1024);x3=exp(Ij*2*pi*w3*t(513:1024)/1024);sig_l=zeros(1,1024);sig_2=sig_l;sig_3=sig_l;sig_l(1:256)=xl;%sig1autosig_temp=sig_l.,;tfrwv(sig_temp);sig_2(257:512)=x2;%sig2auto%sig_temp=sig_2.,;%tfrwv(sig_temp);sig_3(513:end)=x3;%sig3auto%sig_temp=sig_3.,;%tfrwv(sig_temp);%totalsigautosig=sig_l+sig_2+sig_3;sig=sig.,;%tfrwv(sig);%sig1sig2cross%sig_temp=sig_l,sig_2.,;%tfrwv(sig_temp);%sig1sig3cross%sigtemp=sig1,sig3.'%tfrwv(sig_temp);%sig2sig3crosssig_temp=sig_2,sig_3).,;tfrwv(sig_temp);3模糊函数t=1:1024;a=0.5;t=64;wn=100;sig_l=(api).0.25*exp(-0.25*(t-t).2);sig_2=(api).0.25*exp(+1j*wn*t);%sig=sig.,;sig_l=sig_l.,;sig_2=sig_2.'%tfrwv(sig_l);%tfrwv(sig_2);sig=sig_l.,sig_2tfrwv(sig);%ambifuncsigcross%ans=ambifunb(sig);%T,F=meshgrid(l:2047,1:2048);%mesh(TzF,abs(ans);%ambifuncsig_lauto%ans=ambifunb(sig_l);%T,F=meshgrid(l:1023,1:1024);%mesh(T,F,abs(ans);%ambifuncsig2autoans=ambifunb(sig_2);T,F=meshgrid(1:1023)z1:1024);mesh(T,F,abs(ans);
