function [f, G] = ceplo_eval(c, a, PHI, B, anl);

%OLC_eval       Computes the opposite of the penalized likelihood
%               criterion and its gradient
%               [f, G] = OLC_eval(c, a, PHI, B, anl)
%               c       Cepstrum vector
%               a       Squared harmonic magnitudes
%               PHI     Cepstral regression matrix
%               B       Diagonal of regularization matrix
%               anl     Apparent noise level at the harm. freq.
%
%This function is a subroutine of OLC and should not be called directly.
%
%_______________________________________________________________________
%For more information about this mfile please take a look at the paper
%"Estimation of the spectral envelope of voiced sounds using a penalized
%likelihood approach" by Marine Oudot, Olivier Cappé and Eric Moulines.

% Authors: Marine Oudot, Olivier Cappé, 1996 - 2000
% ENST Dpt. TSI / LTCI (CNRS URA 820), Paris

% Cepstrum order
p=length(c)-1;

% Regularisation part
reg = sum(B.*(c.^2));

% Envelope fit part
if max(PHI*c)>709
  fprintf(1, 'Warning: envelope may contain +Inf values \n')
end
s = exp(PHI*c);
u = 2*sqrt(s.*a./(anl.^2));

% Compute special functions for vector a
logI0 = zeros(size(u));
I1surI0 = zeros(size(u));
ind = (u < 3.75);
ind1 = find(ind);
ind2 = find(~ind);
% Directly compute I0(u) and I1(u) for small values of the argument
if (length(ind1)~=0)
  tmp = (u(ind1)/3.75).^2;
  bess0 = 1.0+tmp.*(3.5156229+tmp.*(3.0899424+tmp.*(1.2067492+tmp.*(0.2659732+tmp.*(0.360768e-1+tmp.*0.45813e-2)))));
  bess1 = u(ind1).*(0.5+tmp.*(0.87890594+tmp.*(0.51498869+tmp.*(0.15084934+tmp.*(0.2658733e-1+tmp.*(0.301532e-2+tmp.*0.32411e-3))))));
  logI0(ind1) = log(bess0);
  I1surI0(ind1) = bess1./bess0;
end
% Compute log(I0(u)) and I0(u) and I1(u) up to a common factor of
% exp(u)/sqrt(u) to prevent overflows
if (length(ind2)~=0)
  tmp = 3.75./u(ind2);
  bess0_part = 0.39894228+tmp.*(0.1328592e-1+tmp.*(0.225319e-2+tmp.*(-0.157565e-2+tmp.*(0.916281e-2+tmp.*(-0.2057706e-1+tmp.*(0.2635537e-1+tmp.*(-0.1647633e-1+tmp.*0.392377e-2)))))));
  logI0(ind2) = u(ind2)-0.5*log(u(ind2)) + log(bess0_part);
  bess1_part = 0.39894228+tmp.*(-0.3988024e-1+tmp.*(-0.362018e-2+tmp.*(0.163801e-2+tmp.*(-0.1031555e-1+tmp.*(0.2282967e-1+tmp.*(-0.2895312e-1+tmp.*(0.1787654e-1-tmp.*0.420059e-2)))))));
  I1surI0(ind2) = bess1_part./bess0_part;
end

% Value of - log-likeligood
ll = log(anl) + (s+a)./anl - logI0;
% Value of the criterion
f = sum(ll) + reg;

% Gradient of the likelihood part
gl = s./anl - u/2.*I1surI0;
% Add gradient for the regularization term
G = (PHI'*gl + 2*(B.*c));
