Levinson-Durbin-算法-和-伯格(Burg)算法

1、精选优质文档-倾情为你奉上精选优质文档-倾情为你奉上专心-专注-专业专心-专注-专业精选优质文档-倾情为你奉上专心-专注-专业数字信号仿真实现题 目: Levinson-Durbin 算法 和 伯格(Burg)算法讲课老师: 学生姓名：所属院系： 信息科学与工程学院专 业： 信息与通信工程学 号： 完成日期： 2015/5/12Levinson-Durbin 算法%检验Levinson-Durbin算法clear;%清除内存变量clc;%清屏close;%=%估计2阶自回归模型的功率谱%步骤1：建立信号模型，产生观测数据 xvar = 1;u = var*randn(1,1000);%产生均值

2、为0，方差为1的高斯白噪声u，数据长度为1000% figure(3);% plot(u);a0 = 1 0 0.81;x = filter(1,a0,u);%信号模型，白噪声通过线性系统H(z)=1/(1+a1*z(-1)+a2*z(-2)产生信号向量%画出信号x的功率谱,Sxx(exp(j*w)=var/abs(1+sum(ak*exp(-j*w*k)2,求和范围为从1到2，系数ak为给定%的模型参量0，0.81w = linspace(-pi,pi,2000);%将-pi到pi均分为2000等分for mm = 1:2000 c = w(mm); S(mm) = var/(abs(1+a

3、0(2:3)*exp(-j*c*(1:2)2;endsubplot(211);plot(w,S,b);%根据已知参量画出信号x的功率谱xlabel(角频率/rad);ylabel(x的功率谱);title(2阶自回归AR模型的功率谱);%=%估计2阶滑动平均模型的功率谱%步骤1：建立信号模型，产生观测数据 x% clear;var1 = 1;u1 = var1*randn(1,1000);%产生均值为0，方差为1的高斯白噪声u，数据长度为1000b0 = 1 1 1;x1 = filter(b0,1,u1);%信号模型，白噪声通过MA(2)阶系统x(n)=u(n)+u(n-1)+u(n-2)%

4、根据已知信号参量画出信号x的功率谱，S = var*abs(sum(bk*exp(-j*w*k)2for mm = 1:2000 c = w(mm); S1(mm) = var1*(abs(b0*exp(-j*c*(0:2)2;endsubplot(212);plot(w,S1,b);%根据已知参量画出信号x的功率谱xlabel(角频率/rad);ylabel(x的功率谱);title(2阶滑动平均MA模型的功率谱);运行结果如下图1 所示 ：图 1 Levinson-Durbin算法结果图伯格(Burg)算法clear%取样点%定义常数值N=32;a(1)=-0.;d2=0.;f1=0.05

5、;f2=0.40;f3=0.42;ur=0.5*d2.*randn(1,N);ui=0.5*d2.*randn(1,N);u=ur+ui*i;%定义32个复数点z(1)=u(1);x(1)=6+z(1);for n=2:N z(n)=-a(1)*z(n-1)+u(n); x(n)=2*cos(2*pi*f1*(n-1)+2*cos(2*pi*f2*(n-1)+2*cos(2*pi*f3*(n-1)+z(n);end%定义f范围fmin=-0.5;fstep=0.001;fmax=0.5;f=fmin:fstep:fmax;nf=(fmax-fmin)/fstep; t=sqrt(-1);%初值

6、rxx=0;p0=zeros(1,11);ef=zeros(11,N);eb=zeros(11,N);a=zeros(10,10);for n=1:N rxx=rxx+(abs(x(n)2;endrxx=(1/N)*rxx;p0(1)=rxx;ef(1,:)=x;eb(1,:)=x;ef(1,1)=0;eb(1,32)=0;%算法p=10;kk=zeros(1,10);for k=1:p e1=0;e2=0; for n=(k+1):N e1=e1+ef(k,n)*(conj(eb(k,n-1); e2=e2+(abs(ef(k,n)2+abs(eb(k,n-1)2); kk(k)=(-2)*e1/e2; ef(k+1,n)=ef(k,n)+kk(k)*eb(k,n-1); eb(k+1,n)=eb(k,n-1)+conj(kk(k)*ef(k,n); end for i=1:(k-1) a(k,i)=a(k-1,i)+kk(k)*conj(a(k-1,k-i); end a(k,k)=kk(k); p0(k+1)=(1-abs(kk(k)2)*p0(k);end%功率谱for j=1:nf+1 pxx=0; for k=1:p pxx=pxx+a(10,k)*exp(-t*2*pi*f(j)*k); end px