%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    lev.m     
%
%    Attila Vidacs
%
%    30/05/99
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  function [inverse,log_det] = lev(c)
%
%  Input :  c       :  first column of the Toeplitz matrix to be inverted
%
%  Output:  inverse :  inverse of toeplitz(c)
%           log_det :  logarithm of the determinant of toeplitz(c)

function [inverse,log_det] = lev(c)

% Levinson's algorithm

n = length(c);

log_detC = 0;
w(1) = 1;
alf(1) = c(1);

for i = 1:n-1

	beta(i) = flipud(c(2:i+1))'*w;
	
	K(i) = -beta(i)/alf(i);

	u = w./alf(i);
	v(2:i) = flipud(u(2:i));
	v(1) = 0;
	
	Lu = tril(toeplitz(u));
	Lv = tril(toeplitz(v));

	B = ( Lu*Lu' - Lv*Lv')/u(1);

	log_detC = log_detC + log(c(1) - c(2:i+1)'*B*c(2:i+1));
	
% --- i+1 ---
	
	w(i+1,1) = 0;
	w = w + K(i)*flipud(w);
	
	alf(i+1) = alf(i) + K(i)*beta(i);
end

%%% Inverse 

u = w./alf(n);
v(2:n) = flipud(u(2:n));
v(1) = 0;
	
Lu = tril(toeplitz(u));
Lv = tril(toeplitz(v));

inverse = ( Lu*Lu' - Lv*Lv')/u(1);
log_det = log_detC;

