function c = LS(a,omega,p,lr);

%LS             Regularized discrete cepstrum estimation
%               [c,e] = LS(a,omega,p[,lr]);
%               a               Squared harmonic magnitudes
%               omega           Fundamental frequency or list of frequencies
%               p               Cepstrum order
%               lr              Smoothing parameter
%
%This method is referred to as "LS" in the paper.
%
%_______________________________________________________________________
%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.

% Author: Olivier Cappé, 1994 - 2000
% ENST Dpt. TSI / LTCI (CNRS URA 820), Paris

% Input parameters
if (nargin < 4)
  lr = 3.5e-2;
end
% Number of squared amplitude values
L = length(a);
if (length(omega) == 1)
  % If frequencies are harmonic
  omega = omega*(1:L);
end;

% Construct matrix and vector for "amplitude matching part"
al = log(reshape(a,L,1));
M = 2*cos(reshape(omega,L,1)*(1:p));
M = [ones(L,1) M]; 

% Construct matrix for "regularization part"
B = diag(2*((0:p).^2));

% Least-square solution
c = (M'*M + lr*B)\M'*al;
