% -------------------------------------------------------------------------
% levinson.m
% AUTHOR:   Markus Bantle
% DATE:     04-26-2010
%
% Levinson Algorithm for solving a linear system with a Toeplitz matrix
% Source: Golup, Van Loan: Matrix Computations 
%
% INPUT: r Vector defining the matrix T
%        b Right hand side
%
% OUTPUT: x = T^{-1}b
% -------------------------------------------------------------------------

function x = levinson(r,b)
r= r(:);
n = size(r,1);
b = b./ r(1);
r = r(2:n)./r(1);

                                                            
y(1) = -r(1);   x(1) = b(1);
beta = 1;       alpha = -r(1);

for k = 1:n-1
    beta = (1-alpha^2)*beta;
    mu = (b(k+1)-r(1:k)'*x(k:-1:1)')/beta;
    x(1:k) = x(1:k)+mu*y(k:-1:1);
    x(k+1) = mu;
    
    if k < n-1
        alpha = -(r(k+1)+r(1:k)'*y(k:-1:1)')/beta;
        y(1:k) = y(1:k)+alpha*y(k:-1:1);
        y(k+1) = alpha;
    end
end
x = x(:);
% whos
% pause