function [phi_kk] = spacf(r_k)
% SPACF - Computes the S(ample) P(artial) A(utocorrelation) F(unction) using the recursive method of Levinson (1946) and Durbin (1960)
% 
% See the book page 212
% 
% Inputs:
%   r_k : the sample autocorrelations for the given data set.  It is assumed that r_k(1) != 1.0.  
%         If this is true a warning is printed.         
% Outputs: 
%   phi_kk : the sample partial autocorrelation function 
% 
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

r_k = r_k(:).';

if( r_k(1)==1.0 )
  fprintf('WARNING r_k(1)==1.0 ... disgarding this term and continuing\n')
  r_k = r_k(2:end); 
end

phi_kk = []; 

% fill the element phi_11 by hand:
phi_kk(1) = r_k(1); 
phi_kj    = [ r_k(1) ]; 

for ii=2:length(r_k), 
  num        = r_k(ii) - dot( phi_kj, fliplr( r_k(1:(ii-1)) ) ); 
  den        =       1 - dot( phi_kj, r_k(1:(ii-1)) ); 
  phi_kk(ii) = num/den;
  phi_kjm1   = phi_kj - phi_kk(ii) * fliplr( phi_kj );
  phi_kj     = [ phi_kjm1, phi_kk(ii) ]; 
end