Levinson-Durbin-算法-和-伯格(Burg)算法(總6頁(yè)

1、數(shù)字信號(hào)仿真實(shí)現(xiàn)題 目: Levinson-Durbin 算法 和 伯格(Burg)算法講課老師: 學(xué)生姓名：所屬院系： 信息科學(xué)與工程學(xué)院專 業(yè)： 信息與通信工程學(xué) 號(hào)： 完成日期： 2015/5/12Levinson-Durbin 算法%檢驗(yàn)Levinson-Durbin算法clear;%清除內(nèi)存變量clc;%清屏close;%=%估計(jì)2階自回歸模型的功率譜%步驟1：建立信號(hào)模型，產(chǎn)生觀測(cè)數(shù)據(jù) xvar = 1;u = var*randn(1,1000);%產(chǎn)生均值為0，方差為1的高斯白噪聲u，數(shù)據(jù)長(zhǎng)度為1000% figure(3);% plot(u);a0 = 1 0 0.81;x =

2、 filter(1,a0,u);%信號(hào)模型，白噪聲通過線性系統(tǒng)H(z)=1/(1+a1*z(-1)+a2*z(-2)產(chǎn)生信號(hào)向量%畫出信號(hào)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等分for mm = 1:2000 c = w(mm); S(mm) = var/(abs(1+a0(2:3)*exp(-j*c*(1:2)2;endsubplot(211);plot(w,S,b);%根據(jù)已知參量畫出信號(hào)x

3、的功率譜xlabel(角頻率/rad);ylabel(x的功率譜);title(2階自回歸AR模型的功率譜);%=%估計(jì)2階滑動(dòng)平均模型的功率譜%步驟1：建立信號(hào)模型，產(chǎn)生觀測(cè)數(shù)據(jù) x% clear;var1 = 1;u1 = var1*randn(1,1000);%產(chǎn)生均值為0，方差為1的高斯白噪聲u，數(shù)據(jù)長(zhǎng)度為1000b0 = 1 1 1;x1 = filter(b0,1,u1);%信號(hào)模型，白噪聲通過MA(2)階系統(tǒng)x(n)=u(n)+u(n-1)+u(n-2)%根據(jù)已知信號(hào)參量畫出信號(hào)x的功率譜，S = var*abs(sum(bk*exp(-j*w*k)2for mm = 1:200

4、0 c = w(mm); S1(mm) = var1*(abs(b0*exp(-j*c*(0:2)2;endsubplot(212);plot(w,S1,b);%根據(jù)已知參量畫出信號(hào)x的功率譜xlabel(角頻率/rad);ylabel(x的功率譜);title(2階滑動(dòng)平均MA模型的功率譜);運(yùn)行結(jié)果如下圖1 所示 ：圖 1 Levinson-Durbin算法結(jié)果圖伯格(Burg)算法clear%取樣點(diǎn)%定義常數(shù)值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.*r

5、andn(1,N);u=ur+ui*i;%定義32個(gè)復(fù)數(shù)點(diǎn)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);%初值rxx=0;p0=zeros(1,11);ef=zeros(11,N);eb=zeros(11,N);

6、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 pxxf