现代信号处理大作业_1217

1、现代信号处理大作业 姚苏洋现代信号处理大作业姚苏洋 11303491711.使用matlab实现LD迭代算法1.1 Levinson-Durbin算法功率谱估计大致可以分为经典谱估计和现代功率谱估计，经典谱估计方法存在着以下几点缺陷：A). 数据加窗或自相关加窗，都隐含着假定在窗外未观测到的数据或自相关系数为零，该假设不切实际。B). 要性能好往往需要较长的数据，但实际数据长度有限C). 窗函数容易造成谱的模糊。采用AR模型的现代谱估计方法可以克服这些不足。其中LD递推算法可以在计算机上方便实现。LD递推算法具体计算步骤如下：（1） Yule-Walker方程的矩阵形式(1)所示： (1)系数

2、矩阵,为Hermitian矩阵，对角线上元素相同，即为Topliez矩阵。（2） P-1阶Yule-Walker方程为： (2)其中，为误差功率。写成联立方程： (3)取共轭得： (4)变量替换，并利用得： (5)表示成矩阵： (6)求解得： (7) (8) (9) (10) (11)（3） 当k=1时，即一阶递推为： 求解可得： (12)（4） 对于时，递推为: , , (13) (14)矩阵Rx已知，可得到各阶AR模型系数为： (15) (16) (17) (18)1.2实验结果 假设p=5,使用matlab求得递推结果矩阵A如图1所示。画出每次迭代矩阵A的对应值如图2所示。图1 LD算法

3、得到递推结果A图2 多次迭代输出2.令信号由三个不同频率的复正弦信号首尾相连而成。其中，。（1）试求x(t)的WV分布，并画出三维WV分布图（2）指出并分析其WV分布的信号项和交叉项。2.1 原理分析已知WV分布公式如(19)所示： (19) (20)求得的wv分布为：(21)其中： (22)2.2 实验结果使用Matlab时频分析工具箱，可以很方便地帮助我们画出WV分布以及信号项，交叉项。将信号看作三个分段信号的叠加，分别用信号1,2,3表示，其结果分别如图3至图9所示。图3 信号1信号项 图4 信号2信号项图5 信号3信号项 图6 x(t) WV分布图7 信号1 信号2 交叉项 图8 信号

4、1 信号3 交叉项图9 信号2 信号3 交叉项3.一非平稳信号由两个高斯信号叠加而成：分别求出z(t)的WV分布及模糊函数，画出二者的波形图，指出并分析其信号项和交叉项。3.1原理分析根据WV分布的定义式，可以得到的WV分布为; (23)所以的WV分布的信号项为: (24)交叉项为： (25) (26)其中,所以z(t)的WV分布为： (27)（1） 由模糊函数的定义： (28) 可以计算得出z1(t)的模糊函数为： (29)模糊函数信号项为： (30)交叉项为： (31) 其中，所以可得z(t)的模糊函数为： (32)3.2 实验结果图10 信号1wv分布信号项 图11 信号2 wv分布交叉

5、项如图12至图17所示，分别为wv分布信号项，交叉项，模糊函数，信号项，模糊函数交叉项。 图12 WV分布交叉项 图14 模糊函数图15 模糊函数交叉项图16 信号1 信号项图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); for i=2:N newregister(1)=mod(sum(fbconnection.*reg

6、ister),2); for j=2:n newregister(j)=register(j-1); 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);for i =1:p for j =1:p if i = j Rxx(i,j) = Rxx_0; elseif i<j Rxx(i,j) = Rxx_ii(

7、1+j-i); elseif i>j Rxx(i,j) = Rxx_jj(1+i-j); end endend% LD algorithmp = p-1;A = zeros(p,p+1);for 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(

8、loop-1,k)*Rxx(1,loop-k); end delta_p = Rxx(loop,loop)+pre_sum; K = -delta_p/sigma; sigma = sigma*(1-abs(K).2); 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,:);end2 WV分布w1=400; w2=200;w3=100;t = 1:1024;x1 = exp(1j*2*pi*w1*t(1:256)/1024);

9、x2 = exp(1j*2*pi*w2*t(257:512)/1024); x3 = exp(1j*2*pi*w3*t(513:1024)/1024); sig_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

10、auto% sig_temp = sig_3.'% tfrwv(sig_temp);% total sig autosig = sig_1+sig_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.'% tfrwv(sig_temp);% sig2 sig 3 crosssig_temp = sig_2,sig_3.'tfrwv(

11、sig_temp);3 模糊函数t = 1:1024;a = 0.5;t0 = 64;wn = 100;sig_1 = (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 si