資源描述：
《Levinson-Durbin算法和伯格(Burg)算法.doc》由會員上傳分享，免費在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫。
1、數(shù)字信號仿真實現(xiàn)題目:Levinson-Durbin算法和伯格(Burg)算法講課老師:學(xué)生姓名：所屬院系：信息科學(xué)與工程學(xué)院專業(yè)：信息與通信工程學(xué)號：完成日期：2015/5/12Levinson-Durbin算法%檢驗Levinson-Durbin算法clear;%清除內(nèi)存變量clc;%清屏close;%===========================================================%估計2階自回歸模型的功率譜%步驟1：建立信號模型，產(chǎn)生觀測數(shù)據(jù)xvar=1;u=var*randn(1,1000);%產(chǎn)生均值為0，方差為1的高斯
2、白噪聲u，數(shù)據(jù)長度為1000%figure(3);%plot(u);a0=[100.81];x=filter(1,a0,u);%信號模型，白噪聲通過線性系統(tǒng)H(z)=1/(1+a1*z^(-1)+a2*z^(-2))產(chǎn)生信號向量%畫出信號x的功率譜,Sxx(exp(j*w))=var/abs(1+sum(ak*exp(-j*w*k)))^2,求和范圍為從1到2，系數(shù)ak為給定%的模型參量0，0.81w=linspace(-pi,pi,2000);%將-pi到pi均分為2000等分formm=1:2000c=w(mm);S(mm)=var/(abs(1+a0(2:3)*
3、exp(-j*c*(1:2))'))^2;endsubplot(211);plot(w,S,'b');%根據(jù)已知參量畫出信號x的功率譜xlabel('角頻率/rad');ylabel('x的功率譜');title('2階自回歸AR模型的功率譜');%===========================================================%估計2階滑動平均模型的功率譜%步驟1：建立信號模型，產(chǎn)生觀測數(shù)據(jù)x%clear;var1=1;u1=var1*randn(1,1000);%產(chǎn)生均值為0，方差為1的高斯白噪聲u，數(shù)據(jù)長度為1000b0=[
4、111];x1=filter(b0,1,u1);%信號模型，白噪聲通過MA(2)階系統(tǒng)x(n)=u(n)+u(n-1)+u(n-2)%根據(jù)已知信號參量畫出信號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');%根據(jù)已知參量畫出信號x的功率譜xlabel('角頻率/rad');ylabel('x的功率譜');title('2階滑動平均MA模型的功率譜');運行結(jié)果
5、如下圖1所示：圖1Levinson-Durbin算法結(jié)果圖伯格(Burg)算法clear%取樣點%定義常數(shù)值N=32;a(1)=-0.;d2=0.;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個復(fù)數(shù)點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);
6、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:
7、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))^
