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

数字信号仿真实现题目:Levinson-Durbin算法和伯格(Burg)算法讲课老师:学生姓名：所属院系：信息科学与工程学院专业：信息与通信工程学号：完成日期：2015/5/12Levinson-Durbin算法%检验Levinson-Durbin算法clear;%清除内存变量clc;%清屏close;%===========================================================%估计2阶自回归模型的功率谱%步骤1：建立信号模型，产生观测数据xvar=1;u=var*randn(1,1000);%产生均值为0，方差为1的高斯白噪声u，数据长度为1000%figure(3);%plot(u);a0=[100.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等分formm=1:2000c=w(mm);S(mm)=var/(abs(1+a0(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=[111];x1=filter(b0,1,u1);%信号模型，白噪声通过MA(2)阶系统x(n)=u(n)+u(n-1)+u(n-2)%根据已知信号参量画出信号x的功率谱，S=var*abs(sum(bk*exp(-j*w*k)))^2formm=1:2000c=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所示：图1Levinson-Durbin算法结果图伯格(Burg)算法clear%取样点%定义常数值N=32;a(1)=-0.850848;d2=0.101043;f1=0.05;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);forn=2:Nz(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);%初值rxx=0;p0=zeros(1,11);ef=zeros(11,N);eb=zeros(11,N);a=zeros(10,10);forn=1:Nrxx=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);fork=1:pe1=0;e2=0;forn=(k+1):Ne1=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);endfori=1:(k-1)a(k,i)=a(k-1,i)+kk(k)*conj(a(k-1,k-i));enda(k,k)=kk(k);p0(k+1)=(1-abs(kk(k))^2)*p0(k);end%功率谱forj=1:nf+1pxx=0;fork=1:ppxx=pxx+a(10,k)*exp(-t*2*pi*f(j)*k);endpxxf(j)=d2/[(abs(1+pxx))^2];pxxf(j)=10*log10(pxxf(j));endplot(f,pxxf)title('伯格(Burg)算法');运行结果如下图2所示：图2Burg算法结果图