burg.m

这是一个用于语音信号处理的工具箱

function [a, sig2] = burg(x,n,ip)
      
rho0=0.;
for i=1:n
	rho0=rho0+abs(x(i))^2/n;
end
for i=2:n
	efk1(i)=x(i);
	ebk1(i-1)=x(i-1);
end
for k=1:ip
	sumn=0.;
	sumd=0.;
	for i=k+1:n
		sumn=sumn+efk1(i)*conj(ebk1(i-1));
	    	sumd=sumd+abs(efk1(i))^2+abs(ebk1(i-1))^2;
	end
	aa(k,k)=-2.*sumn/sumd;
	if k == 1
		rho(k)=(1.-abs(aa(k,k))^2)*rho0;
	end
	if k > 1
		rho(k)=(1.-abs(aa(k,k))^2)*rho(k-1);
	end
	if ip == 1 
		break;
	end
	if k ~= 1
		for j=1:k-1
			aa(j,k)=aa(j,k-1)+aa(k,k)*conj(aa(k-j,k-1));
		end
	end    
	for i=k+2:n
		efk(i)=efk1(i)+aa(k,k)*ebk1(i-1);
	    	ebk(i-1)=ebk1(i-2)+conj(aa(k,k))*efk1(i-1);
	end
	for i=k+2:n
		efk1(i)=efk(i);
		ebk1(i-1)=ebk(i-1);
	end
end
sig2=rho(ip);
for i=1:ip
	a(i)=aa(i,ip);
end