现代信号处理大作业.docx

1、现代信号处理大作业现代信号处理大作业姓名：潘晓丹学号：0140349045 班级：A1403492作业1LD算法实现AR过程估计1.1AR模型p阶AR模型的差分方程为：px(n) y ai x( n i) =w(n),其中w(n)是均值为o的白噪声。i 二Yule-Walker方程递推AR过程的线性预测方法为：先求得观测数据的自相关函数，然后利用求得模型参数，再根据公式求得功率谱的估计。Yule-Walker方程可写成矩阵形式:J(0)rxx(1)rxx(2)rxx(-P)-1 1frxx(1)xx(O)rxx( 1)rxx(-P +1)ap(1)0rxx(2)rxx(1)rxx(0)xx(

2、-P + 2)ap(2)=0xx( p)rxx( P - 1)xx(P 2)rxx(0) 一ap( P)_i i0 一1.2LD算法介绍Levinson-Durbin算法可求解上述问题，其一般步骤为：1)计算观测值各自相关系数rxx(j ), j二0,1,，P ； 5二rxx(0) ； i=1 ；2)利用以下递推公式运算：irxx(i) + 迟 a(j )订-j )ai(i)二 kiai (j ) - ai J(j ) - ki ai J(i - j ), j - 1,2. i -12:i = : i 二(1 - ki )3) i=i+1 ，若ip，则算法结束；否则，返回 。1.3以AR模型:

3、matlab编程实现为例，Matlab程序代码如下:clear; clc;var = 1;noise = var*ra ndn (1,10000);p = 2;coefficient = 1 -0.5 0.5;x = filter(1,coefficient,noise);divide = linspace(-pi,pi,200);for ii = 1:200w = divide(ii);S1(ii) = var/(abs(1+coefficie nt(2:3)*exp(-j*w*(1:2)F2;enda_p var_p=Levinson_Durbin(x,p);for ii = 1:200w

4、 = divide(ii);Sxx(ii) = var_p/(abs(1+a_p(2:p+1)*exp(-j*w*(1:p)A2;endfigure;subplot(2,2,1); plot(divide,S1, b);grid onxlabel(w); ylabel(功率); title(AR 功率谱);subplot(2,2,2); plot(divide,Sxx, r-);grid onxlabel(w); ylabel(功率); title(L-D 算法估计);subplot(2,2,3); plot(divide,S1, b);hold onplot(divide,Sxx, r-);

5、hold offgrid onxlabel(w); ylabel(功率); title( AR 功率谱和算法比较 );子函数： Levinson_Durbin .mfunction a_p var_p = Levinson_Durbin(x,p)N = length(x);for ii=1:NRxx(ii) = x(1:N-ii+1)*(x(ii:N)/N;enda(1)=1;a(2)=-Rxx(2)/Rxx(1);for k=1:p-1 % Levinson-Durbin algorithmvar(k+1) = Rxx(0+1)+a(1+1:k+1)*Rxx(1+1:k+1); reflec

6、t_coefficient(k+1+1) = -a(0+1:k+1)*(fliplr(Rxx(2:k+1+1)/var(k+1); var(k+1+1) = (1-(reflect_coefficient(k+1+1)A2)*var(k+1);a_temp(1) = 1;for kk=1:ka_temp(kk+1) = a(kk+1)+reflect_coefficient(k+1+1)*a(k+1-kk+1); enda_temp(k+1+1) = reflect_coefficient(k+1+1);a = a_temp;enda_p = a;var_p = var(p+1);% pred

7、icti on coeffecie nts% predicti on error power1.4仿真结果1)阳功率i亂和LB登法比啟p=2时，仿真结果图如下预测系数：误差功率：var_p=1.01942)p=20时，仿真结果图如下预测系数:误差功率：var_p=0.99983)ARU谙和LD負:去比较Wp=50时，仿真结果图如下预测系数:误差功率：var_p=0.99551.5结果分析由不同阶数(P值)得到的仿真结果可得：当 P的阶数较低时，L-D算法估计AR模型对 功率谱估计的分辨率较低，有平滑的效果，从P=2的仿真结果可以看出估计得到的功率谱与 原始功率谱基本吻合，且曲线平滑没有毛刺；随

8、着阶数增大，采用 L-D算法进行估计后，得到的功率谱会产生振荡，从仿真可以看到，当阶数 P较高为50时，估计得到的功率谱与原始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。从LD算法得到的预测系数可得， 不论阶数P(2)取何值，通过该算法得到的预测参数与 原始目标函数中的一致，其余各个参数均接近 0。因此，L-D算法得到可计算得到较为精确的预测值或估计值。作业非平稳信号由两个高斯信号叠加而成，为12 :- 2z(t) =() exp( (t -ti) j -it) exp( (t -t2) j 2t)TL 2 2其中tl “21 J，2，分别求出z(t)的WV分布及其模糊函数，画

9、出二者的波形图，指出并分析其信号项和交叉项。2.1 WV分布由WV分布的定义知： Wt) z(t )z*(t )e j d .2 2i若记 Zi(t )=(】)4 exp( (t -ti)2 j .J),JT 21Z2(t )=(二)4 exp( (t -t2)2 j 2t)二 2则 W(t , ) = WU(t , ) WCross (t ,),2)Z2*(t-2 d)z/*(t )ed其中，W/uto(t, )二Zi(t )Zi (t )ej d . Z2(t2 2Wross(t, J 二Zi(t )Z2*(t )ej d. Z2(t2 2由此，我们计算得到：信号项：2 1 2 1W/ut

10、o(t, ) =2exp-:(t -t1) (；-； ) 2exp -:(t -12) (； ；，2)a a交叉项：2 1Wross (t , ) = 4 exp-： (t -t m) ( ；：， ；： m) COS( ；：， ；：，mX d (t - t m) d CLt 1 t 2 - 1其中，t m ,m ,t d - t 1 _ t 2 / _ 22 22.2模糊函数由模糊函数的定义知:*,可二：z(t )z*(t 0如 2 2i若记 Zi(t ) =()4exp( (t _t i)2 j it )兀 21Z2(t ) = ( 一)4 exp( (t t 2)2 j 2t )JT 2则

11、 Ae8) =Aauto(8) +Across(y，8)，其中，Aauto(,讨二Zi(t )Zi (t )ej 勺 dt ，Z2(t )Z2(t )ej dt2 2 2 2Across (,刃二 Zi(t 尹2 (t - ：)e心dt Z2(t ；)Zi (t - ；)e心dt2 2 2 2由此，我们计算得到：信号项：2 1 2i. - jt p) exp( j 2 - jt p)心o().exp(盲石)exp(交叉项:1 2 2 1 2 2Across( ” d) = expj ( F -t dt m)exp石u -d)盲(td) expw(d)-C td)2.3Matlab编程实现取11

12、 = 10,t 2 = 4, x = 8, ，24,=20，进行matlab编程如下clear all;format long;alpha仁20; t1=10;t2=4;w 1=8;w2=4;a=alpha1/2; td=t1-t2; omegad=w1-w2;tm=0.5*(t1+t2); omegam=0.5*(w1+w2);m=1; n=1;for t=0:0.1:8for omega=-6:0.1:12W_auto(m, n)=2*(exp(-a*(t-t1)A2-1/a*(omega-w1F2)+exp(-a*(t-t2)A2-1/a*(omega-w2)A2);W_cross( m

13、,n)=4*exp(-a*(t-tm)A2-1/a*(omega-omegam)A2)*cos(omega-omegam)*td+omegad*t)W(m,n)=W_auto(m,n)+W_cross(m,n);n=n+1;endm=m+1;n=1;endfigure;mesh(-6:0.1:12,0:0.1:8,W);xlabel(time );ylabel(frequency );title( WV 分布 );figure;mesh(-6:0.1:12,0:0.1:8,W_auto);xlabel(time );ylabel(frequency );title( WV 分布信号项 );fi

14、gure;mesh(-6:0.1:12,0:0.1:8,W_cross);xlabel(time );ylabel(frequency );title( WV 分布交叉项 );format long; %模糊函数a=10; t1=6;t2=2;w1=6;w2=2;td=t1-t2; wd=w1-w2;tm=0.5*(t1+t2); wm=0.5*(w1+w2);m=1;n=1;for t=-10:0.1:10for w=-10:0.1:10A_auto(m,n)=abs(exp(-a/4*tA2-1/(4*a)*wA2)*(exp(i*w1*t-i*t1*w)+exp(i*w2*t-i*t2*

15、w);A_cross(m,n)=abs(exp(i*wm*t+i*w*tm+i*wd*tm)*(exp(-1/(4*a)*(w+wd)A2-a/4*(t-td)A2)+exp(-1 /(4*a)*(w-wd)A2-a/4*(t+td)A2);A(m,n)=A_auto(m,n)+A_cross(m,n);n=n+1;endm=m+1;n=1;endfigure;mesh(-10:0.1:10,-10:0.1:10,A);xlabel(time );ylabel(frequency);title(模糊函数);figure;mesh(-10:0.1:10,-10:0.1:10,A_auto); x

16、label(time);ylabel(frequency);title(模糊函数信号项);figure;mesh(-10:0.1:10,-10:0.1:10,A_cross); xlabel(time);ylabel(frequency);title(模糊函数交叉项);2.4仿真结果及分析1)Z(t)函数的WV分布波形图如下2)WV分布信号项波形图如下卅山玮唁号曲3） WV分布交叉项波形图如下 .iciisncr根据结果可以看到， WV分布的两个信号项是分开的，分别以 （ti,1）和（t22）为中心;WV分布的交叉项则耦合在一起，以 （t mJ m）为中心。it4 ）模糊函数波形图5）模糊函数

17、信号项囈舉的劉百号顶frquncy 106）模糊函数交叉项frequency in -10由结果可以看出，模糊函数的两个信号项耦合在一起，以原点（ 0,0）为中心；而交叉项是作业三相乘得到。信号- 是由 一 与分别求出三者的 WV分布，并画出三维分布图。该计算过程和作业二类似，在此不再赘述。3.1 Matlab编程实现syms x ta=10; b=2; o=2;zt=(a/pi)A0.25*exp(-a*tA2/2+1i*b*tA2/2+1i*o*t);xt=(a/pi)A0.25*exp(-a*tA2/2);yt=exp(1i*b*tA2/2+1i*o*t);x1=(a/pi)A0.25*

18、exp(-a*(t+x/2)A2/2);x2=(a/pi)A0.25*exp(-a*(t-x/2)A2/2);x12=(a/pi)A0.5*exp(-a*(t-x/2)A2/2-a*(t-x/2)A2/2);y1=exp(1i*b*(t+x/2)A2/2+1i*o*(t+x/2);y2=exp(-1i*b*(t-x/2)A2/2-1i*o*(t-x/2);z1=(a/pi)A0.25*exp(-a*(t+x/2)A2/2+1i*b*(t+x/2)A2/2+1i*o*(t+x/2);z2=(a/pi)A0.25*exp(-a*(t-x/2)A2/2+-1i*b*(t-x/2)A2/2-1i*o*

19、(t-x/2);z12=(a/pi)A0.5*exp(-a*(t+x/2)A2/2+1i*b*(t+x/2)A2/2+1i*o*(t+x/2)-a*(t-x/2)A2/2+-1i*b*(t-x/2)A2/2-1i *o*(t-x/2);Wx=simple(fourier(x12)* (a/pi)A(1/2);Wz=simple(fourier(z12)* (a/pi)A(1/2);Wx=i nlin e(Wx)Wz=i nli ne(Wz)t=-6:0.01:6;f=-6:0.01:6;T,F=meshgrid(t,f);Wxwv= abs(Wx(T,2*pi*F);Wzwv= abs(Wz(T,2*pi*F);figure;mesh (T,F,Wxwv);axis(-6 6 -6 6 -1 3);xlabel(时间);ylabel(频率);title(X(t)WV 分布);figure;mesh (T,F,Wzwv);axis(-6 6 -6 6 -1 3);xlabel(时间);ylabel(频率);title( Z(t)WV 分布);3.2仿真结果及分析1)X(t)函数的 WV分布2)注意到y信号的WV分布是奇异函数，这里便并没有将其画出来。手算得到其WV分布为：( )3)Z(t)函数的 WV分布