function [Beta,r,ierr,L1_norm]=spefac(rhs,limit,epsilon) 
%function [Beta,r,ierr,L1_norm]=spefac(rhs,limit,epsilon) 
%
% The polynomial Beta(z)=z^n + b_1*z^(n-1) + ... + b_n
% and the real-valued parameter r are computed such that 
%                
%               r*Beta(z)*Beta(1/z) = rhs 
% 
% An iterative Newton algorithm by V. Kucera, 
% Discrete Linear Control (1979) is used. 
% Convergence problems may occur for root clusters close 
%  to the unit circle. 
% 
% Inputs: 
%  1.     rhs - "the right hand side" 
%  2.   limit - max number of iterations (default 300) 
%  3. epsilon - accuracy limit for termination  (default 1e-10) 
% 
% Outputs: 
%  1.    Beta - spectral factor, with zeros inside the unit circle 
%  2.       r - real-valued scalar
%  3.    ierr - Error codes: 
%               a) ierr=0, convergence has occured. 
%               b) ierr=1, no convergence within limited iterations. 
%               c) ierr=2, right hand side (rhs) not positive definite
%
%  4. L1_norm - the L1-norm of the error (r*Beta(z)*Beta(1/z) - rhs) 
%
%
%               [Beta,r,ierr,L1_norm]=spefac(rhs,limit,epsilon) 
% 
%
% M. Sternad  (1995-08) 
 
if nargin < 2 ,limit=300;end 
if nargin < 3 ,epsilon=1E-10;end 
 
stop=0; 
ierr=0; 
 
n=(length(rhs)+1)/2; 
 
rr = rhs(n:length(rhs)); 
if rr(1) < max(abs(rr)),ierr=2;stop=1;end 
 
if rhs(n) > 0 
  A=rhs(n:2*n-1)/sqrt(rhs(n)); 
  stop=0; 
else 
  ierr=1; 
  stop=1; 
end 
 
% start iteration 
 
i=0; 
while ~stop 
    
  i=i+1; 
  B=rhs(n:2*n-1); 
  W=A; 
  k=0; 
  while k~=n-1 
    for j=1:n-k,rhs1(j)=A(n-k+1-j);end 
    Al(k+1)=A(n-k)/rhs1(n-k); 
    if abs(Al(k+1)) < 0 
      ierr=2; 
      stop=1; 
    else 
      A(1:n-k-1)=A(1:n-k-1)-Al(k+1)*rhs1(1:n-k-1); 
      E(n-k)=2*B(n-k)/A(1); 
      if k<n-2 
        for j=1:n-k-2,B(j+1)=B(j+1)-E(n-k)*A(n-k-j)/2;end 
      end 
      k=k+1; 
    end 
  end 
  E(1)=B(1)/A(1); 
  while k>0 
    k=k-1; 
    for j=1:n-k,P(j)=E(n-k+1-j);end 
    E(1:n-k)=E(1:n-k)-Al(k+1)*P(1:n-k); 
  end 
   
  A=(E+W)/2; 
  if sum(A.^2)-B(1) < epsilon*B(1) 
    ierr=0; 
    stop=1; 
  elseif i==limit 
    ierr=1; 
    stop=1; 
 end 
  
end 
 
Beta = A/A(1); 
r = A(1)^2; 
L1_norm = norm(r*xcorr(Beta)-rhs,1);