ls_my.m

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

function [a,b,iflag] = ls_my(y,z,yz,n,ip,iq,epsilon) 
      
ipq=ip+iq;
l=1;
for j=1:ip+iq 
      	for i=1:j 
      		hh(l)=0.;
      		for k=1:n 
      			if i <= ip & k-i >= 1 
				hki=y(k-i);
			end
      			if i <= ip & k-i <= 0 
				hki=0.;
			end
      			if i > ip & k-(i-ip) >= 1 
				hki=z(k-(i-ip));
			end
      			if i > ip & k-(i-ip) <= 0 
				hki=0.;
			end
      			if j <= ip & k-j >= 1 
				hkj=y(k-j);
			end
      			if j <= ip & k-j <= 0 
				hkj=0.;
			end
      			if j > ip & k-(j-ip) >= 1 
				hkj=z(k-(j-ip));
			end
      			if j > ip & k-(j-ip) <= 0 
				hkj=0.;
			end
			hh(l)=hh(l)+conj(hki)*hkj;
		end
      		l=l+1;
	end
end
for i=1:ip+iq 
      	theta(i)=0.;
      	for k=1:n 
      		if i <= ip & k-i >= 1
			hki=y(k-i);
		end
      		if i <= ip & k-i <= 0
			hki=0.;
		end
      		if i > ip & k-(i-ip) >= 1
			hki=z(k-(i-ip));
		end
      		if i > ip & k-(i-ip) <= 0
			hki=0.;
		end
		theta(i)=theta(i)+conj(hki)*yz(k);
	end
end
[iflag theta] = cholesky(hh,theta,ipq,epsilon);
for i=1:ip+iq 
      	if i <= ip
		a(i)=-theta(i);
	end
	if i > ip
		b(i-ip)=theta(i);
	end
end