levin_my.m

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

function [a, aa, rho]= levin_my(r,ip)
      
rho0=r(1);
aa(1,1)=-r(2)/r(1);
rho(1)=(1.-abs(aa(1,1))^2)*r(1);
a(1)=aa(1,1);
if ip == 1 
	return
end
      
for i=2:ip
	b=-r(i+1);
	for m=1:i-1
		b=b-aa(m,i-1)*r(i+1-m);
	end
      	aa(i,i)=b/rho(i-1);
      	for k=1:i-1
		aa(k,i)=aa(k,i-1)+aa(i,i)*conj(aa(i-k,i-1));
	end
      	rho(i)=(1.-abs(aa(i,i))^2)*rho(i-1);
end

for i=1:ip
	a(i)=aa(i,ip);
end