function f = lsf_saoudi(cor, iter);
% Syntax: f = lsf_saoudi(cor, iter);
% 
% computes Line Spectral Frequencies LSF from speech autocorrelation 
% coefficients using fast (not in Matlab) algorithm described in:
% Saoudi, Boucher, Guyader:
% "A new efficient algorithm to compute the LSP parameters for speech coding" 
% Signal Processing, no. 28, 1992, pp. 201-212.
% 
% cor is the vector of correlation coefficients c0...cP, where P is the
%   order of the prediction filter
% f is the vector of LSF parameters sorted in increasing order, size P
% iter (optional) is the precision of eigenvalue computation in eig_bisect
%   the max. error of eigenvalues is 4/2^iter, but for LSF it is modified
%   by the arccos. Default 9.
%
% - all vectors are COLUMN
% - The function works only for P even
%
% The function is calling:
%  levinson_anti.m   - calculation of alphap and alphaq parameters.
%  saoudi_diag.m - creation of diagonals using alphap and alphaq parameters.
%  eig_bisect.m   - finding eigenvalues of 3-diagonal matr.
%
% The implementation in Matlab is not "vectorized", so it may be very slow

if nargin==1,
  iter=9;	% default for eigenvalues
end
P=length(cor)-1;
if(P/2 ~= floor(P/2) )
 error ('P is not even');
end

%%% alpha parameters %%%
[alphap, alphaq] = levinson_anti (cor);

%%% diagonals creation %%%
[dp, ep, dq, eq] = saoudi_diag(alphap, alphaq);

%%% eigenvalues using bisection %%%
evp=eig_bisect (dp, ep, iter);
evq=eig_bisect (dq, eq, iter);

% for tests using func. eig. of Matlab.
%ep1=ep(2:P/2); eq1=eq(2:P/2);  % elimination of 1st dummy element
%JP=diag(dp) + diag(sqrt(ep1),+1) + diag(sqrt(ep1),-1); % putting dp to the main
%			% diag and sqrt of ep to secondary diags.
%JQ=diag(dq) + diag(sqrt(eq1),+1) + diag(sqrt(eq1),-1); % putting dp to the main
%			% diag and sqrt of ep to secondary diags.
%evp=eig(JP); evq=eig(JQ);

%%% conversion, sorting %%%
f=acos([evp; evq] * 0.5) / 2 / pi;
f=sort(f);