现代信号处理大作业王成志1

现代信号处理大作业姓名：王成志学号： 1140349078一. LD 迭代算法的 matlab 实现1.1 Levinson-Durbin 算法介绍功率谱估计大致可以分为经典谱估计和现代功率谱估计，经典谱估计方法存在着以下三点缺陷：（1 ）数据加窗或自相关加窗，都隐含着假定在窗外未观测到的数据或自相关系数为零，该假设不切实际。（2 ）要性能好往往需要较长的数据，但实际数据长度有限（3 ）窗函数容易造成谱的模糊。采用 AR 模型的现代谱估计方法可以克服这些不足。其中 LD 递推算法可以在计算机上方便实现。LD 递推算法具体计算步骤如下：（1 ） Yule-Walker 方程的矩阵形式 (1)所示：01)0()2()1()0)()()() 2, kkxxxx xxxx arkrrkrrr系数矩阵 ,为 Hermitian 矩阵，对角线上元素相同，即为 Topliez 矩阵。HR（2 ） P-1 阶 Yule-Walker 方程为：211(0)(1)(1),20,()()(0)xxx ppxxxRaRp 其中， 为误差功率。2211Eel写成联立方程：21 1,0 ,0()p pkxk maR取共轭得：21*1,0 ,()0p pkxk变量替换，并利用 得：*()xxRll21* 1,0 ,()0p pkxk maR 表示成矩阵：*1210()()(),12()()(0)xxxppxxxaRpR 求解得：*.1,1,pkkpkaKa2p210p,Ka22*211ppppK（3 ） 当 k=1 时，即一阶递推为：0)0(21,1aRxx求解可得：)1()0(,21,1, xxxRa（4 ） 对于 时，递推为:2p, , 0,pa*1,1, kpkpkK 122ppK21,pK1,)()(kxx kRaR矩阵 Rx 已知，可得到各阶 AR 模型系数为：)0( ,)0(1)( 211 xx rara11)()()( kklxxkk lr1,2)()()(*11 kikaiaikkk )(kkka1.2 实验结果(1) 输入 p=3，rr = 70,60,50,40 时，求得 AR 模型估计参数为 :a =1.0000 -0.8571 0 01.0000 -0.5275 -0.3846 01.0000 -0.7572 -0.6996 0.5972各阶求得的方差为：sigma = 18.5714 15.8242 10.18013阶时，a 3 (1)= -0.7572 a3 (2)= -0.6996 a3 (3)= -0.5972(2) 输入p=5，rr = 30,45,26,33,47,43时，AR模型估计参数为：a =1.0000 -1.5000 0 0 0 01.0000 0.2800 -1.1867 0 0 01.0000 0.8227 -1.3147 -0.4573 0 01.0000 1.9708 1.9858 -2.5226 -2.5105 01.0000 1.0869 1.0977 -1.8235 -1.8166 0.3521各阶求得的方差为：sigma =37.5000 15.3067 12.1054 64.1881 56.23165阶时， a5 (1)= 1.0869 a5(2)= 1.0977 a5(3)= -1.8235 a5(4)= -1.8166 a5(5)= 0.3521二 一维平稳信号由两个高斯信号叠加而成，其中12 2411()exp()exp()zttjttjt12,t，分别求出 的 WV 分布及其模糊函数，画出二者的波12zt形图，指出并分析其信号项和交叉项。(1) 信号 的 WV 分布求取公式为()zt2(,)()2jfzWtfztted在这里我们使用两种方法来对 WV 分布进行求取，分别为手工求法和 matlab 求法。A. 求解过程由二次叠加原理知，Wigner-Ville 分布的自项和交叉项分别为12(,)(,)(,)autozzWwtWtw1221, ,crszz式中 411exp()tjt1422()ztjwt由 Wigner-Ville 分布的定义知111(,)()()2jzWtwztted2 211xp()()4jwtd 在积分公式 中2eexpACBABCd 令 ， 和 ，则得4A1Bjw21t1 2211(,)2expzWttw类似地，我们有 2 22(,)expzttw故 2211(,)auto i iiwt又由交叉项定义知 12, 12()()()jwzWtwztted2221212112exp()4tjttjtd 于是，积分公式 中的参数为epexpACBAxBCd 4A1212Btjw2112Cttjt由此得 12 2,()expexpz mmmdWttwjwt故 Wigner-Ville 分布的交叉项 12,(,)R()crosztwtw2214expcosmmmdt wt式中 和12mtt12分别为两个谐波信号时延的平均值和频率的平均值，而 和12dtt，分别是两个谐波信号时延差与频率差。12dw将上述自项和交叉项相加就可以得到 WV 分布(,)(,)(,)zautocrosWtWtw22222111exp4expcosi i mmmdit twwtB. 程序实现由上面步骤可以发现，手工求取 WV 分布非常繁琐。在这里，我使用matlab 来求取 WV 分布并画图。为了进行比较，我还同时画出了 WV 分布信号项的图和交叉项的图。程序运行结果为当参数值 时，WV 分布表达式：12120,3,twWz=502186294043021/1407374883553280*exp(-10*t2)*pi(1/2)*5(1/2)*2(1/2)*(exp(50*t-505/8+5/2*i+i*t+1/2*w-i*w-1/10*w2)+exp(40*t-202/5-1/10*w2+2/5*w)+exp(-1/10*w2-909/10+3/5*w+60*t)+exp(50*t-505/8-5/2*i-i*t+1/2*w+i*w-1/10*w2)得到的结果图为：WV 分布图 WV 分布信号项图WV 分布交叉项图当设置如上参数时，WV 分布图中信号项和交叉项混叠在一起。我们只能把信号项图和交叉项图分别显示出来才能看清楚。当参数值 时，WV 分布表达式：121235,3,twWz=93950452832135/492581209243648*exp(-35*t2)*35(1/2)*pi(1/2)*(exp(210*t-11034/35-1/35*w2+6/35*w)+exp(140*t-4904/35-1/35*w2+4/35*w)+exp(175*t-3065/14+5/2*i+i*t+1/7*w-i*w-1/35*w2)+exp(175*t-3065/14-5/2*i-i*t+1/7*w+i*w-1/35*w2)程序中，Wz 中原来的分量 用 代替。f2wf得到的结果图为：WV 分布图 WV 分布信号项图WV 分布交叉项图当设置如上参数时,从图 4 可以看到, 信号项部分和交叉项部分被分了开来，结合上一组参数的经验，我们很容易对哪部分是信号项和哪部分是交叉项做出判断。而图 5 和图 6 的分布图正好证实了我们的判断。当我们对其他参数进行改变时，WV 分布图也会发生变化。(2) 信号 的模糊函数求取公式为()zt2,)()jtvzAvztted和求取 WV 分布类似，我同样使用两种方法来求取模糊函数，分别为手工求法和 matlab 求法。A 求解过程信号的 模糊函数的信号项和交叉项分别为()zt1221, ,()()()autozzcrsA由模糊函数公式知11122221 11(,)()()expexp4jtzAzttdjttjtdt又由积分公式 代入得：22 ACBABCd1 2211(,)expexp4zAjt类似地，我们有 2 222(,)4z jt故可以得到： 22121(,)expexpexp4autoAjtjt对于交叉项，我们有： 12, 12 21 22221211expexp4jtz jtzttdjtjtedtt t 12jt 于是，积分公式 中的参数为22expexpACBABxCd 12214ddmtjtj由此得 12 22,()expexp44z ddmdAtjt类似地有： 21 22, 1()eez ddmdmtjt所以模糊函数的交叉项1221, , 222(,)()()1expexp4444croszzddddmmAAt tjt 所以，该信号的模糊函数为：22122 2(,)(,)(,)1expexpexp4444expzautocrosddddmmAAjtjtt tjt 式中 和 分别为两个谐波信号时延的平均值和频率12mtt12的平均值，而 和 ，分别是两个谐波信号时延差与频率差。dtdB 程序实现在这一部分，我使用 matlab 来求取模糊函数并画图。为了进行比较，我还同时画出模糊函数信号项的图和交叉项的图。程序运行结果为当参数值 时，模糊函数表达式：1212,twAz=5081767996463981/18014398509481984*exp(-1/4*t2)/pi(1/2)*(exp(i*t+i*v-1/4*v2)+exp(-1/2*t-1/2-3/2*i+3/2*i*t+1/2*v+3/2*i*v-1/4*v2)+exp(2*i*t+2*i*v-1/4*v2)+exp(1/2*t-1/2+3/2*i+3/2*i*t-1/2*v+3/2*i*v-1/4*v2)程序中，我们使用 来代替手工求法中的 ， 用 来代替手工求法中的 。t v得到的结果图为：模糊函数幅度分布图 模糊函数信号项幅度分布图模糊函数交叉项幅度分布图当参数值 时，模糊函数表达式：12120,twAz=502186294043021/5629499534213120*exp(-5/2*t2)/pi(1/2)*10(1/2)*(exp(5*t-101/40+3/2*i+3/2*i*t-1/20*v+3/2*i*v-1/40*v2)+exp(-5*t-101/40-3/2*i+3/2*i*t+1/20*v+3/2*i*v-1/40*v2)+exp(i*t+i*v-1/40*v2)+exp(2*i*t+2*i*v-1/40*v2)程序中，我们使用 来代替t得到的结果图为：模糊函数幅度分布图 模糊函数信号项幅度分布图模糊函数交叉项幅度分布图三 令信号 是由 和 21-t4()jtjtzteX(t)=21-t4e()=相乘而得到的，Y(t)是线性调频 Chirp 信号，分别求出20jtjteX（ t）,Y(t)和 Z(t)的 WV 分布并画出三者的三维分布图。A 求解过程WV 分布的计算公式如下： 2(,)()2jfzWtfztted带入 X（t）,Y(t)和 Z(t)后得到 X（t）的 WV 分布如下式： 221212(,)()expexp2ejfx jfWfttdttedtjfdKtjf 得到 Y（t）的 WV 分布如下式： 2 212 20 0120(,)()expexpj2ejfy jfWftytdtjt tedjtfd 0Ktf得到 Z（t）的 WV 分布如下式：2212 00(,)()()expjfz jfWtfzttedyytedtjfKtf在本题中，Y（ t）频率 为 30HZ，调频斜率 ，有 Nyquist 采样定理得采样频0 率 。信号长度选为奇数点 901，衰减控制因子 ，信号持续时段为10sfHz 2。5B 程序实现运行程序得到 X（t）,Y(t)和 Z(t)的 WV 三维分布如下：X(t)的 WV 三维分布 Y(t)的 WV 三维分布Z(t)的 WV 三维分布附：实验代码1.LD 算法function a delte gama=LD_Calucation(R);p=length(R)-1;%阶数比自相关矩阵长度少一a=zeros(p+1,1);%初始化AR系数，该系数包含p+1个值，其中a(1)=1delte=R(1);%由0阶计算1阶的AR参数D=R(2);gama(1)=D/delte;delte=(1-gama(1)2)*delte;a(2)=-gama;a(1)=1;for k=1:p-1 %迭代知道计算出p阶的AR参数delte=R(1)+a(2:(k+1)*R(2:(k+1);D=a(1:(k+1)*R(k+2):-1:2);gama(k+1)=D/delte;delte=(1-gama(k+1)2)*delte;a(k+2)=-gama(k+1);a(2:(k+1)=a(2:(k+1)-gama(k+1)*a(k+1):-1:2);enda;delte;gama;2.非平稳信号 的 WV 分布12 2411()exp()exp()zttjttjt及模糊函数求解代码clear all;format long;syms x t aifa=input(Please input aifa:); %设置初始参数值t1=input(Please input t1:);t2=input(Please input t2(t1t2):);w1=input(Please input w1:);w2=input(Please input w2(w1w2):);z1=exp(-aifa/2*(t+x/2-t1)2+j*w1*(t+x/2)+exp(-aifa/2*(t+x/2-t2)2+j*w2*(t+x/2); %求取z(t+x/2)z2=exp(-aifa/2*(t-x/2-t1)2-j*w1*(t-x/2)+exp(-aifa/2*(t-x/2-t2)2-j*w2*(t-x/2); %求取z(t-x/2)的共轭z3=exp(-aifa/2*(t+x/2-t1)2+j*w1*(t+x/2)*exp(-aifa/2*(t-x/2-t1)2-j*w1*(t-x/2); %为求WV分布信号项作准备z4=exp(-aifa/2*(t+x/2-t2)2+j*w2*(t+x/2)*exp(-aifa/2*(t-x/2-t2)2-j*w2*(t-x/2);z=z1*z2;z_auto=z3+z4;Wz=simple(fourier(z)*(aifa/pi)(1/2) %求取WV分布Wz_auto=simple(fourier(z_auto)*(aifa/pi)(1/2); %求取WV分布的信号项Wz=inline(Wz); Wz_auto=inline(Wz_auto);t=0:0.01:5; %设置t,f的取值范围f=-5:0.01:5;T,F=meshgrid(t,f);W=double(Wz(T,2*pi*F); %求取在指定范围内的WV分布值W_auto=double(Wz_auto(T,2*pi*F); %求取在指定范围内WV分布信号项的值figure(name,WV分布); %画出WV分布图mesh(T,F,W);xlabel(time);ylabel(frequency);figure(name,WV分布的信号项); %画出WV分布信号项的三维图mesh(T,F,W_auto);xlabel(time);ylabel(frequency);W_cross=W-W_auto;figure(name,WV分布的交叉项); %画出WV分布交叉项的三维图mesh(T,F,W_cross);xlabel(time);ylabel(frequency);xlabel(time);ylabel(frequency);clear all;format long;syms w t vaifa=input(Please input aifa:); %设置初始参数值t1=input(Please input t1:);t2=input(Please input t2(t1t2):);w1=input(Please input w1:);w2=input(Please input w2(w1w2):);z1=exp(-aifa/2*(w+t/2-t1)2+j*w1*(w+t/2)+exp(-aifa/2*(w+t/2-t2)2+j*w2*(w+t/2); %求取z(w+t/2)z2=exp(-aifa/2*(w-t/2-t1)2-j*w1*(w-t/2)+exp(-aifa/2*(w-t/2-t2)2-j*w2*(w-t/2); %求取z(w-t/2)的共轭z3=exp(-aifa/2*(w+t/2-t1)2+j*w1*(w+t/2)*exp(-aifa/2*(w-t/2-t1)2-j*w1*(w-t/2); %为求模糊函数信号项作准备z4=exp(-aifa/2*(w+t/2-t2)2+j*w2*(w+t/2)*exp(-aifa/2*(w-t/2-t2)2-j*w2*(w-t/2);z=z1*z2;z_auto=z3+z4;Az=simple(ifourier(z,v)*(aifa/pi)(1/2)Az_auto=simple(ifourier(z_auto)*(aifa/pi)(1/2);Az=inline(Az); %求取模糊函数，在这里暂时还没有乘前面的系数项Az_auto=inline(Az_auto); %求取模糊函数的信号项，同样没有乘以系数t=-4:0.01:4; %设置t,v的取值范围v=-4:0.01:4;T,V=meshgrid(t,v);A=abs(double(Az(T,V); %求取在指定范围内的模糊函数值A_auto=abs(double(Az_auto(T,V); %求取在指定范围内模糊函数信号项的值figure(name,模糊函数); %画出模糊函数图mesh(T,V,A);xlabel(time interval);ylabel(variable v);figure(name,模糊函数中的信号项); %画出模糊函数信号项的三维图mesh(T,V,A_auto);xlabel(time interval);ylabel(variable v);A_cross=abs(double(Az(T,V)-double(Az_auto(T,V); %画出模糊函数交叉项的三维图figure(name,模糊函数中的交叉项);mesh(T,V,A_cross);xlabel(time interval);ylabel(variable v);ylabel(frequency);xlabel(time);ylabel(frequency);3. 信号 X(t),Y(t),Z(t)的 WV 分布求解代码clear all; alpha = 2*pi; beta = pi;f0 = 30; omiga0 = 2*pi*f0;fs = 100;T = 9.01; N = T*fs;L = (N-1)/2;n = -L:L; t = n/fs;x = (alpha/pi)0.25*exp( -alpha/2*t.2) ;y = exp(0.5i*beta*t.2 + 1i*omiga0*t);z = x.*y;m = -L:L;k = 0:L;Wv_x = zeros(N, L+1);Wv_y = zeros(N, L+1);Wv_z = zeros(N, L+1);for ii = 1:Nfor jj = 1:L+1nm1 = n(ii) + m; t_tau1 = nm1/fs;nm2 = n(ii) - m; t_tau2 = nm2/fs;x_tau1 = (alpha/pi)0.25*exp( -alpha/2*t_tau1.2);