现代信号处理大作业姓名：潘晓丹学号：0140349045班级：A1403492作业1LD 算法实现AR 过程估计1.1 AR 模型p 阶AR 模型的差分方程为：)()()(1n w i n x a n x pii =-+∑=，其中)(n w 是均值为0的白噪声。

AR 过程的线性预测方法为：先求得观测数据的自相关函数，然后利用Yule-Walker 方程递推求得模型参数，再根据公式求得功率谱的估计。

Yule-Walker 方程可写成矩阵形式：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--+-+--000)()2()1(1)0()2()1()()2()0()1()2()1()1()0()1()()2()1()0(2σp a a a r p r p r p r p r r r r p r r r r p r r r r p p p xx xx xx xxxx xx xx xx xx xx xx xx xx xx xx xx 1.2 LD 算法介绍Levinson-Durbin 算法可求解上述问题，其一般步骤为：1) 计算观测值各自相关系数p j j r xx ,,1,0),( =；)0(0xx r =ρ；i=1； 2) 利用以下递推公式运算：)1(1,...,2,1),()()()()()()(21111111i i i i i i i ii i i jxx i xx i k i j j i a k j a j a k i a j i r j a i r k -=-=--==-⋅+-=-----=-∑ρρρ3) i=i+1，若i>p ，则算法结束；否则，返回(2)。

1.3 matlab 编程实现以AR 模型：xn=12xn-1-12xn-2+w(n)为例，Matlab 程序代码如下：clear; clc;var = 1;noise = var*randn(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+coefficient(2:3)*exp(-j*w*(1:2))'))^2;end[a_p var_p]=Levinson_Durbin(x,p);for ii = 1:200w = divide(ii);Sxx(ii) = var_p/(abs(1+a_p(2:p+1)*exp(-j*w*(1:p))'))^2;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--');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)';reflect_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))^2)*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; % prediction coeffecientsvar_p = var(p+1); % prediction error power1.4 仿真结果1）p=2时，仿真结果图如下预测系数：a20,a21,a22=[1,-0.5068,0.5031]误差功率：var_p=1.01942）p=20时，仿真结果图如下预测系数：a20,a21,a22,a23,a24,……=[1,-0.5098,0.4999,-0.0066,0.0060,-0.0179,0.0193,……]误差功率：var_p=0.99983）p=50时，仿真结果图如下预测系数：a20,a21,a22,a23,a24,……=[1,-0.4951,0.5178,-0.0145,0.0117,-0.0169,0.0141,……]误差功率：var_p=0.99551.5 结果分析由不同阶数（P值）得到的仿真结果可得：当P的阶数较低时，L-D算法估计AR模型对功率谱估计的分辨率较低，有平滑的效果，从P=2的仿真结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；随着阶数增大，采用L-D算法进行估计后，得到的功率谱会产生振荡，从仿真可以看到，当阶数P较高为50时，估计得到的功率谱与原始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。

从LD算法得到的预测系数可得，不论阶数P(>2)取何值，通过该算法得到的预测参数与原始目标函数中的一致，其余各个参数均接近0。

因此，L-D算法得到可计算得到较为精确的预测值或估计值。

作业二非平稳信号由两个高斯信号叠加而成，为12241122()()[exp(())exp(())]22z t t t j t t t j t αααωωπ=--++--+其中2121,ωω>>t t ，分别求出()z t 的WV 分布及其模糊函数，画出二者的波形图，指出并分析其信号项和交叉项。

2.1 WV 分布由WV 分布的定义知：τττωωτd e t z t z t W j -∞∞--+=⎰)2()2(),(*若记))(2exp()()(121411t j t t t z ωαπα+--=,))(2exp()()(222412t j t t t z ωαπα+--=则 ),(),(),(ωωωt W t W t W cross auto +=, 其中,ττττττωωτωτd e t z t z d e t z t z t W j j auto -∞∞--∞∞--++-+=⎰⎰)2()2()2()2(),(*22*11ττττττωωτωτd e t z t z d e t z t z t W j j cross -∞∞--∞∞--++-+=⎰⎰)2()2()2()2(),(*12*21由此，我们计算得到： 信号项：)](1)(exp[2)](1)(exp[2),(222121ωωααωωααω----+----=t t t t t W auto 交叉项：])()cos[()](1)(exp[4),(2d m d m m m cross t t t t t t W ωωωωωααω-+-----=其中，21212121,,2,2ωωωωωω-=-=+=+=d d m m t t t t t t2.2 模糊函数由模糊函数的定义知：dt e t z t z A t j θττθτ)2()2(),(*-+=⎰∞∞-若记))(2exp()()(121411t j t t t z ωαπα+--=,))(2exp()()(222412t j t t t z ωαπα+--=则 ),(),(),(θτθτθτcross auto A A A +=, 其中，dt e t z t z dt e t z t z A tj t j auto θθττττθτ)2()2()2()2(),(*22*11-++-+=⎰⎰∞∞-∞∞-dt e t z t z dt et z t z A t j tj cross θθττττθτ)2()2()2()2(),(*12*21-++-+=⎰⎰∞∞-∞∞-由此，我们计算得到： 信号项：)]exp())[exp(414exp(),(221122θτωθτωθαταθτjt j jt j A auto -+---= 交叉项：]})(4)(41exp[])(4)(41)]{exp[(exp[),(2222d d d d m d m m cross t t t t j A +---+--+-++=ταωθαταωθαωθτωθτ其中，21212121,,2,2ωωωωωω-=-=+=+=d d m m t t t t t t2.3 Matlab 编程实现取20,4,8,4,102121=====αωωt t ，进行matlab 编程如下 clear all ;format long ;alpha1=20; t1=10;t2=4;w1=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)^2-1/a*(omega-w1)^2)+exp(-a*(t-t2)^2-1/a*(omega-w2)^2));W_cross(m,n)=4*exp(-a*(t-tm)^2-1/a*(omega-omegam)^2)*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分布信号项');figure;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:10 A_auto(m,n)=abs(exp(-a/4*t^2-1/(4*a)*w^2)*(exp(i*w1*t-i*t1*w)+exp(i*w2*t-i*t2*w)));A_cross(m,n)=abs(exp(i*wm*t+i*w*tm+i*wd*tm)*(exp(-1/(4*a)*(w+wd)^2-a/4*(t-td)^2)+exp(-1/(4 *a)*(w-wd)^2-a/4*(t+td)^2)));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); xlabel('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 分布交叉项波形图如下根据结果可以看到，WV 分布的两个信号项是分开的，分别以),(和),(2211ωωt t 为中心； WV 分布的交叉项则耦合在一起，以),(m m t ω为中心。