cholesky.m

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

function [iflag, b] = cholesky(a,b,n,epsilon)
      
iflag=0;
l=1;
d(1)=real(a(l));
for i=2:n
	for j=1:i-1
      		l=l+1;
      		xl(i,j)=conj(a(l))/d(j);
      		if j ~= 1
      			for k=1:j-1
				xl(i,j)=xl(i,j)-xl(i,k)*conj(xl(j,k))*d(k)/d(j);
			end
		end
	end
      	l=l+1;
      	d(i)=real(a(l));
      	for k=1:i-1
		d(i)=d(i)-d(k)*abs(xl(i,k))^2;
	end
      	if d(i) <= epsilon
      		iflag=-1;
      		return
	end
end
y(1)=b(1);
for k=2:n
	y(k)=b(k);
      	for j=1:k-1
    		y(k)=y(k)-xl(k,j)*y(j);
	end
end
b(n)=y(n)/d(n);
for k=n-1:-1:1
	b(k)=y(k)/d(k);
      	for j=k+1:n
		b(k)=b(k)-conj(xl(j,k))*b(j);
	end
end