lsmywe.m

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

function [a,b,sig2,iflag] = lsmywe(x,n,ip,iq,m,epsilon)
      
lag=m+1;
r = correlat(n,lag,0,x,x);
for i=1:m-iq
	for j=1:ip
		if iq+i-j >= 0
			rr(i,j)=r(iq+1+i-j);
		end
		if iq+i-j < 0
			rr(i,j)=conj(r(j-i-iq+1));
		end
	end
end
l=1;
for j=1:ip
	for i=1:j
		aa(l)=0;
		for k=1:m-iq
			aa(l)=aa(l)+conj(rr(k,i))*rr(k,j);
		end
		l=l+1;
	end
end
for i=1:ip
	a(i)=0;
	for k=1:m-iq
		a(i)=a(i)-conj(rr(k,i))*r(iq+1+k);
	end
end
[iflag a] = cholesky(aa,a,ip,epsilon);
if iflag == -1
	return
end
y = czeroflt(x,n,ip,a);
ll=round(n/5);
[b sig2] = durbin(y,n,iq,ll);