burg.m

burg法估计AR（P)模型参数的算法。里面ef是前项误差bf是后项误差

clear;
close all;
[X ts bits]=wavread('c:\windows\media\ding.wav');
sound(X,ts,bits);
data=X(:,1);
len=length(data);
order=5;
ef=zeros(order+1,len);
eb=zeros(order+1,len);
%mse=zeros(order+1,1);
%r=zeros(order+1,1);
a=zeros(order+1,order+1);
a(:,1)=1;

ef(1,:)=data';
eb(1,:)=data';
mse(1)=sum(data.*data)/len;
for k=1:order
 ef_temp=ef(k,k+1:len);
 bf_temp=eb(k,k:len-1);

 r(k+1)=2*sum(ef_temp.*bf_temp)/sum(ef_temp.*ef_temp+bf_temp.*bf_temp);

 mse(k+1)=(1-r(k+1)*r(k+1))*mse(k);


 ef(k+1,:)=ef(k,:)-[0,r(k+1)*eb(k,1:len-1)]; 
 eb(k+1,:)=[0,eb(k,1:len-1)]-r(k+1)*ef(k,:); 

  for i=1:k
   a(k+1,i+1)=a(k,i+1)-r(k+1)*a(k,k-i+1);
  end
end

for m=1:order+1
    for n=1:m
        l(m,m-n+1)=a(m,n);
    end
end