《现代信号处理大作业1.docx》由会员分享，可在线阅读，更多相关《现代信号处理大作业1.docx（27页珍藏版）》请在三一办公上搜索。
1、现代信号处理大作业1现代信号处理大作业姚苏洋 现代信号处理大作业 姚苏洋1130349171 1.使用matlab实现LD迭代算法 1.1 Levinson-Durbin算法 功率谱估计大致可以分为经典谱估计和现代功率谱估计，经典谱估计方法存在着以下几点缺陷： A). 数据加窗或自相关加窗，都隐含着假定在窗外未观测到的数据或自相关系数为零，该假设不切实际。 B). 要性能好往往需要较长的数据，但实际数据长度有限 C). 窗函数容易造成谱的模糊。采用AR模型的现代谱估计方法可以克服这些不足。其中LD递推算法可以在计算机上方便实现。 LD递推算法具体计算步骤如下： Yule-Walker方程的矩阵
2、形式(1)所示： rxx(-2)rxx(0)rxx(-1)r(1)r(0)rxx(-1)xxxxMMMrxx(k)rxx(k-1)rxx(k-2)rxx(-k)1s2aLrx(-k+1)k,1=0MM (1) OMLrxx(0)ak,k0LH系数矩阵Rxx=Rxx,为Hermitian矩阵，对角线上元素相同，即为Topliez矩阵。 P-1阶Yule-Walker方程为： 21Rx(-1)LRx(1-p)spRx(0)-1R(1)a,1R(0)R(2-p)0p-1xxx= (2) MMMLa,p-1R(p-1)R(p-2)LR(0)0p-1xxx22其中，sp-1=Eep-1(l)为误差功率。
3、 写成联立方程： 2sp-1,m=0 (3) aR(m-k)=p-1,kxk=00,m=1,L,p-1p-1 取共轭得： ak=0p-1*p-1,k2sp-1,m=0 (4) R(m-k)=0,m=1,L,p-1*x1 现代信号处理大作业姚苏洋 *变量替换，并利用Rx(l)=Rx(l)得： ak=0p-1*p-1,p-1-k2sp-1,m=p-1 (5) Rx(m-k)=0,m=0,L,p-2表示成矩阵： Rx(-1)LRx(1-p)a*0Rx(0)p-1,p-10R(1)*R(0)R(2-p)a,p-2xxxp-1 (6) =MMLM2R(p-1)R(p-2)LR(0)1xxxsp-1求解得
4、： ap.k=ap-1,k+Kpa*p-1,p-k,k=0,L,p (7) 22*sp=sp+KD-1pp (8) 20=D2p+Kpsp-1 (9) Kp=ap,p (10) 22*22sp=sp-1+Kp-Kpsp-1=sp-11-Kp (11) 2 当k=1时，即一阶递推为： Rx(0)Rx(-1)1s12R(1)R(0)a= xx1,10求解可得： a1,0=1, a1,1=-Rx(1)Rx(0)(12) s12=Rx(0)+a1,1Rx(1) 对于p2时，递推为: ap,01, ap,k=ap-1,k+KaKp=ap,p=-Dpp-1*pp-1,p-k,22sp=sp-11-Kp(1
5、3) 2s2p-1Dp=Rx(p)+ap-1,kRx(p-k) (14) k=1矩阵Rx已知，可得到各阶AR模型系数为： a1(1)=-rxx(1)2, r1=(1-a1(1)rxx(0) (15)rxx(0)2 现代信号处理大作业姚苏洋 ak(k)=-Dkrk-1=-rxx(k)+ak-1(l)rxx(k-l)l=1k-1rk-1 (16) *ak(i)=ak-1(i)+ak(k)ak-1(k-i)2i=1,2,L,k-1 (17) rk=(1-ak(k)rk-1 (18) 1.2实验结果 假设p=5,使用matlab求得递推结果矩阵A如图1所示。画出每次迭代矩阵A的对应值如图2所示。 图1
6、 LD算法得到递推结果A 10.90.80.70.60.50.40.30.20.1011.522.533.544.5510.50-0.5-1-1.5-211.522.533.544.55110.50.500-0.5-0.5-1-1-1.5-211.522.533.544.55-1.511.522.533.544.55图2 多次迭代输出 ej2pf1tj2pf2t2.令信号x(t)=eej2pf3t0tT/4T/4tT/2由三个不同频率的复正弦信号首尾相连而成。T/2tT其中，f1=4f0,f2=2f0,f3=f0。 3 现代信号处理大作业姚苏洋 试求x(t)的WV分布，并画出三维WV分布图 指
7、出并分析其WV分布的信号项和交叉项。 2.1 原理分析 已知WV分布公式如(19)所示： W(z,f)=t*t-j2pft(19) z(t+)z(t-)edt22-+x(t)=ej2pf1tu1(t)+ej2pf2tu2(t)+ej2pf3tu3(t) (20) 求得x(t)的wv分布为： +W(z,f)=-+jw1(t+)2x(t+)x*(t-)e-jwtdt22tjw2(t+)2ttt =-(eu1(t)+eu2(t)+ejw3(t+)2tu3(t)(-e-jw1(t-)2tu1(t)-e-jw2(t-)2tu2(t)-e-jw3(t-)2tu3(t)ejwtdt (21) =2pd(w-
8、w1)*U1(t)+2pd(w-w2)*U2(t)+2pd(w-w3)*U3(t)+4pcos(w1-w2)d(w- +4pcos(w1-w3)d(w-w1+w2)*U12(t)2w1+w3w+w3) *U13(t)+4pcos(w2-w3)d(w-2)*U23(t)22f+f =d(f-f1)*U1(t)+d(f-f2)*U2(t)+d(f-f3)*U3(t)+2cos2p(f1-f2)td(f-12)*U12(t)2f+ff+f +2cos(f1-f3)td(f-13) *U13(t)+2cos(f2-f3)td(f-23)*U23(t)22 其中： U1(t)=u1(t)*u1(t)U2
9、(t)=u2(t)*u2(t)U3(t)=u3(t)*u3(t)U12(t)=u1(t)*u2(t)U13(t)=u1(t)*u13(t)U23(t)=u2(t)*u3(t)2.2 实验结果 使用Matlab时频分析工具箱，可以很方便地帮助我们画出WV分布以及信号项，交叉项。将信号看作三个分段信号的叠加，分别用信号1,2,3表示，其结果分别如图3至图9所示。 (22) 4 现代信号处理大作业姚苏洋 图3 信号1信号项图4 信号2信号项 图5 信号3信号项图6 x(t) WV分布 图7 信号1 信号2 交叉项图8 信号1 信号3 交叉项 图9 信号2 信号3 交叉项 3.一非平稳信号由两个高斯信
10、号叠加而成： aaz(t)=exp(-(t-t0)2+jw0t) 2p分别求出z(t)的WV分布及模糊函数，画出二者的波形图，指出并分析其信号项和交叉项。 3.1原理分析 5 1/4现代信号处理大作业姚苏洋 根据WV分布的定义式，可以得到z(t)的WV分布为; Wz(t,f)=+-z(t+)z*(t-)e-j2pftdt221/4ttta = -p+e+-(t+-t0)2+jw0(t+)222atapdt1/4e-(t-t0)2-jw0(t-)222-jwtattedta = pa1/2-e-a(t-t0)2-t2-j(w-w0)t4a1/22+-t2-j(w-w)ta (23) = (t-t
11、0)pe-a4-e0dt1/2-1 = (w-wa-a(t-t0)20)2ape4pae = 2e-a(t-t10)2-a(w-w0)2 =2e-a(t-t0)2-4p2a(f-f0)2 所以z(t)的WV分布的信号项为: Wt,f)=2e-a(t-t0)2-4p2a(f-f0)2 auto( 交叉项为： 24p2Wz,z(t,f)=2e-a(t-t0)-a(f-f20)e-j4p2(f-f0)2td-j2pfdt Wcross(t,f)=2ReWz,z(t,f) = 2e-a(t-tm)22-4p(f-fm)2 acos4p2(f-fm)2td+2pfdt其中tt0+t0m=2,ff0+f0
12、m=2,td=0,fd=0 所以z(t)的WV分布为： Wz(t,f)=Wauto(t,f)+Wcross(t,f) =2e-a(t-t0)2-4p2(f-f-4p2a0)2 +2e-a(t-tm)a(f-fm)2cos4p2(f-f2m)td+2pfdt 由模糊函数的定义： +A(t,v)=z(t+t)z*(t-t)ej2ptvdt -22可以计算得出z1(t)的模糊函数为： 6 (24) (25) (26) (27) (28) 现代信号处理大作业姚苏洋 +Az(t,v)=*j2ptvz(t+)z(t-)edt22-+tt =-ap1/21/4e-a(t+-t0)2+jw0(t+)222tt
13、ap1/4e-a2(t-t0)2-jw0(t-)22ttej2ptvdta = pa = pa = pa+-+ee-a(t-t0)2-t2+jw0t+j2ptu4adtdt (29) 1/2-a(t-t0-jpu2a2p22)-t-u+j(w0t-t0u)a4a-1/2tp-ae4a2-p22u+j(w0t-2pt0u)a =e-t2-u2+j(w0t-2pt0u)4ap2模糊函数信号项为： Aauto(t,v)= e交叉项为： -t2-u2+j(w0t-2pt0u)4aap2 (30) Across(t,v)=Az1,z1(t,v) = e其中，tm=-1a(2pu+wd)2-a4(t-td
14、)2ej(wmt+2ptmu+wdtm) (31) t0+t0w0+w0w=td=0，wd=0 m，22所以可得z(t)的模糊函数为： Az(t,v)=Aauto(t,v)+Across(t,v) =e3.2 实验结果 -t2-u24aap2ej(w0t-2pt0u)+e1a-(2pu)2-(t)2a4 (32) ej(wmt+2ptmu)图10 信号1wv分布信号项图11 信号2 wv分布交叉项 如图12至图17所示，分别为wv分布信号项，交叉项，模糊函数，信号项，模糊函数交叉项。 7 现代信号处理大作业姚苏洋 图12 WV分布交叉项图14 模糊函数 图15 模糊函数交叉项 图16 信号1
15、信号项 图17 信号2 信号项 附：实验代码 1 LD算法 % PN sequence generation fbconnection = 1 0 0 0 0 1 0 0 1 1 ; n=length(fbconnection); N=2n-1; register=zeros(1,n-1) 1; mseq(1)=register(n); fori=2:N 8 现代信号处理大作业姚苏洋 newregister(1)=mod(sum(fbconnection.*register),2); for j=2:n newregister(j)=register(j-1); end; register=n
16、ewregister; mseq(i)=register(n); end % initialization rxx = 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:p for j =1:p ifi = j Rxx(i,j) = Rxx_0; elseifij Rxx(i,j) = Rxx_jj(1+i-j); end end end % LD algorithm p = p-1; A = zeros(p,p+1); for
17、 loop = 1:p if loop =1 p = 1; a_p0 = 1; a_p1 = -Rxx(2,1)/Rxx(1,1); sigma = Rxx(1,1)+a_p1*Rxx(2,1); A(loop,loop) = a_p0; A(loop,loop+1) = a_p1; else pre_sum = 0; for k = 1:loop-1 pre_sum = pre_sum + A(loop-1,k)*Rxx(1,loop-k); end delta_p = Rxx(loop,loop)+pre_sum; K = -delta_p/sigma; sigma = sigma*(1-
18、abs(K).2); 9 现代信号处理大作业姚苏洋 for k =1:loop+1 a_pk = A(loop-1,k)+K*conj(A(loop-1,loop-k+2); A(loop,k) = a_pk; end end figure; stem(A(loop,:); end 2 WV分布 w1=400; w2=200; w3=100; t = 1:1024; x1 = exp(1j*2*pi*w1*t(1:256)/1024); x2 = exp(1j*2*pi*w2*t(257:512)/1024); x3 = exp(1j*2*pi*w3*t(513:1024)/1024); si
19、g_1 = zeros(1,1024); sig_2 = sig_1; sig_3 = sig_1; sig_1(1:256) = x1; % sig 1 auto sig_temp = sig_1.; tfrwv(sig_temp); sig_2(257:512) = x2; % sig 2 auto % sig_temp = sig_2.; % tfrwv(sig_temp); sig_3(513:end) = x3; % sig 3 auto % sig_temp = sig_3.; % tfrwv(sig_temp); % total sig auto sig = sig_1+sig_
20、2+sig_3; sig = sig.; % tfrwv(sig); % sig 1 sig 2 cross % sig_temp = sig_1,sig_2.; % tfrwv(sig_temp); % sig 1 sig 3 cross % sig_temp = sig_1,sig_3.; 10 现代信号处理大作业姚苏洋 % tfrwv(sig_temp); % sig2 sig 3 cross sig_temp = sig_2,sig_3.; tfrwv(sig_temp); 3 模糊函数 t = 1:1024; a = 0.5; t0 = 64; wn = 100; sig_1 = (
21、a/pi).0.25*exp(-0.25*(t-t0).2); sig_2 = (a/pi).0.25*exp(+1j*wn*t); % sig = sig.; sig_1 = sig_1.; sig_2 = sig_2.; % tfrwv(sig_1); % tfrwv(sig_2); sig = sig_1.,sig_2.; tfrwv(sig); % ambifunc sig cross % ans = ambifunb(sig); % T,F = meshgrid(1:2047,1:2048); % mesh(T,F,abs(ans); % ambifunc sig_1 auto % ans = ambifunb(sig_1); % T,F = meshgrid(1:1023,1:1024); % mesh(T,F,abs(ans); % ambifunc sig 2 auto ans = ambifunb(sig_2); T,F = meshgrid(1:1023,1:1024); mesh(T,F,abs(ans); 11
