function [ax,lambda] = lspfit(a,wa,L,tresh)
% function [ax,lambda] = lspfit(a,wa,L,tresh)
%    Expand filter 'a' through interpolation of the symmetric and 
%    antisymmetric LSP polynomials with optimal weighting 'lambda',
%    that is, optimise 'lambda' (in function 'lspx') through golden
%    section -search, by iterative comparisons of absolute difference 
%    of autocorrelations of the original signal and the impulse 
%    response of the expanded filter.
%    Function returns the optimal LSP expanded filter 'ax' and the 
%    optimal 'lambda'.
%
%    This function is an appendix to the article "All-pole Modeling 
%    Technique Based on the Weighted Sum of the LSP Polynomials",
%    presented at ICASSP'02. However, the golden section -search
%    was not discussed in that paper, but more information is
%    available in
%       Bazaraa, Sherali and Shetty, "Nonlinear Programming - Theory
%       and Algorithms", Wiley, 1993.
%
%    See also: WLSP, LSP, LSPX

% Tom Bäckström / 11.10.2001

% check optinal input variables
if nargin < 4
  tresh = 0.000001;
  if nargin < 3
    L = 2*length(a);
  end;
end;

% golden section ratio

gold = sqrt(5/4)-0.5;

% autocorrelation vector

R = xcorr(wa,L,'coeff');
R = R(1+L+(0:L));

% construct the four points required for golden section
% search (left edge, left mid, right mid and right edge), 
% and set initial values

S = struct('value',0,'lambda',0);
ledge = S;
redge = S;
lmid = S;
rmid = S;
ledge.lambda = 0.001;
ledge.value = lspdiff(lspx(a,ledge.lambda),R,L);
redge.lambda = 0.999;
redge.value = lspdiff(lspx(a,redge.lambda),R,L);
rmid.lambda = gold;
rmid.value = lspdiff(lspx(a,rmid.lambda),R,L);
lmid.lambda = 1-gold;
lmid.value = lspdiff(lspx(a,lmid.lambda),R,L);


% perform golden section -search

while(1)
  
  % check if left triplet (left edge, left mid and right mid)
  % contains minimum
  
  if lmid.value < rmid.value
    
    % reset three of the four points to left triplet, and
    % calculate new value for the missing point
    
    redge.value = rmid.value; 
    redge.lambda = rmid.lambda;
    
    rmid.value = lmid.value; 
    rmid.lambda = lmid.lambda;
    
    lmid.lambda = redge.lambda - gold*(redge.lambda - ledge.lambda);
    lmid.value = lspdiff(lspx(a,lmid.lambda),R,L);

  else
    
    % reset three of the four points to right triplet, and
    % calculate new value for the missing point

    ledge.value = lmid.value; 
    ledge.lambda = lmid.lambda;
    
    lmid.value = rmid.value; 
    lmid.lambda = rmid.lambda;
    
    rmid.lambda = ledge.lambda + gold*(redge.lambda - ledge.lambda);
    rmid.value = lspdiff(lspx(a,rmid.lambda),R,L);
    
  end

  % check for end conditions
  
  if rmid.lambda - lmid.lambda < tresh
    
    % value of sufficient accuracy has been found -> exit
    
    disp('Accuracy threshold limit reached');
    break;
    
  elseif (redge.value < rmid.value) | (ledge.value < lmid.value)
    
    % objective function seems to be non-convex (this sometimes
    % happens due to noise in autocorrelation vector)
    
    disp('Noise level reached'); 
    break;
  end
end

% choose minimum of mid-points as global minimum
% and return associated lambda and coefficient vector

if lmid.value < rmid.value
  lambda = lmid.lambda;
  ax = lspx(a,lambda);
else  
  lambda = rmid.lambda;
  ax = lspx(a,lambda);
end;
