toepeqn.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 67 行

function [a,iflag] = toepeqn(r,ip,iq,epsilon)
      
iflag=0;
aa(1,1)=-r(iq+2)/r(iq+1);
if ip == 1
	for i=1:ip
		a(i)=aa(i,ip);
	end
	return
end

bb(1,1)=-r(iq)/r(iq+1);
rho=(1.-aa(1,1)*bb(1,1))*r(iq+1);
if abs(rho)<= epsilon
      iflag=-1;
      return
end

for k=2:ip
	c=-r(iq+1+k);
	for l=1:k-1
		lag=iq+k-l;
		if lag >= 0 
			corr=r(lag+1);
		end
		if lag < 0 
			corr=conj(r(-lag+1));
		end
		c=c-aa(l,k-1)*corr;
	end
	aa(k,k)=c/rho;
	for i=1:k-1
		aa(i,k)=aa(i,k-1)+aa(k,k)*bb(k-i,k-1);
	end
	if k == ip
		break;
	end
	lag=iq-k;
	if lag >= 0 
		d=-r(lag+1);
	else
		d=-conj(r(-lag+1));
	end
	for l=1:k-1
		lag=iq-k+l;
		if lag >= 0
			corr=r(lag+1);
		end
		if lag < 0
			corr=conj(r(-lag+1));
		end
		d=d-bb(l,k-1)*corr;
	end
	bb(k,k)=d/rho;
	for i=1:k-1
		bb(i,k)=bb(i,k-1)+bb(k,k)*aa(k-i,k-1);
	end
	rho=(1.-aa(k,k)*bb(k,k))*rho;
	if abs(rho) <= epsilon
		iflag=-1;
		return
	end
end
for i=1:ip
	a(i)=aa(i,ip);
end