function [alpha,gamma,edeltab,snr,core] = MLfit(infile,fitparameter,varargin)
 
% MLfit evaluate the parameters alpha (¦Á) and gamma (¦Ã) for the inverse power 
% background model using maximum likelihood background estimation. 
%
% Input:
%  infile         Input data, in form of input data matrix or input file name. The input data
%                   matrix is a 2xm matrix, in which the first row stands for the loss energy and 
%                   the second row are the recording background. The input file is a *.dm3 file. 
%                  
% fitparameter = [eprestart, eprewidth, signalstart, signalwidth]
% or fitparameter = [eprestart, eprewidth]
%
%  eprestart      starting energy for the background fitting window : eV
%  eprewidth    width for the backgound fitting window : eV
%  signalstart    starting energy for the signal window : eV. The signal windows is
%                       required to calculate the snr. Default: 467
%  signalwidth  width for the signal window : eV. Default: 12.
%
%  gamma0       intial guess of gamma0
%  min              termination tolerance, default: 10^(-14)
%
%
% Output:
%  alpha,gamma        parameters for the inverse power law background model. B=¦ÁE^(-¦Ã)
%  edeltab                  background fitting error at each energy channel
%  snr                         signal-to-noise ratio
%  core                       signals after removing background
%
% Example:
%  infile = 'S1-T2-2 EELS SI (dark ref crt) aligned sum.dm3' ;
%  eprestart = 454;
%  eprewidth = 10; 
%  fitparameter = [eprestart, eprewidth];
%  [alpha,gamma,snr] = MLfit(infile,fitparameter);
% 
% Ref: Unser et al, Optimal background estimation in EELS, J. Microsc.,
% 1987(145):245-256
% 
%  Note: DM3Import.m (by Robert A. McLeod) is required to input *.dm3 files. 
% 
% Zhu,Guozhen
% Oct 1st, 2011

% check inputs
msg = nargchk(1,4,nargin);
if ~isempty(msg)
    fprintf('Error! Invalid Input!')
    return
end
 
if isnumeric(infile)
    [m,n] = size(infile);
    if m==2
        if n<2||n==2
            fprintf('Input Error! Wrong Input EELS Signal!')
        end
    elseif m>2
        if n==2
            infile =infile';
        else
            fprintf('Input Error! Wrong Input EELS Signal!')
        end
    else
        fprintf('Input Error! Wrong Input EELS Signal!')
    end
    escale = infile(1,:);
    spectrum = infile(2,:); 
elseif isstr(infile)
    spectra_struct = DM3Import(infile);
    spectra =spectra_struct.spectra_data;  % unit:e-
    if size(spectra,2) > 1
        spectrum = spectra{1}.* spectra_struct.intensity.scale;
        edispersion = spectra_struct.xaxis{1}.scale;
        origin = spectra_struct.xaxis{1}.origin;
    else
        spectrum = spectra{1}.* spectra_struct.intensity.scale;
        edispersion = spectra_struct.xaxis.scale;
        origin = spectra_struct.xaxis.origin;
    end
    N = length( spectrum);
    escale = ( (0:N-1)-origin).* edispersion;   % unit:eV
else
    fprintf('Input Error! The first input should be either a 2xm or mx2 matrix or a *.dm3 file!')
end
 
gamma = 2;
min = 10^(-14);
eprewidth = [];
eprestart = [];
plot(escale,spectrum);
axisx = escale;   
spectrum_i = spectrum; 
signalstart = 467;        
signalwidth = 12;        
 
if isnumeric(fitparameter)
    [m,n] = size(fitparameter);
    if n == 1
        fitparameter = fitparameter';
    end
    [m,n] = size(fitparameter);
    if n== 2
        eprestart = fitparameter(1); 
        eprewidth = fitparameter(2);
    elseif n==4
        eprestart = fitparameter(1); 
        eprewidth = fitparameter(2);
        signalstart = fitparameter(3); 
        signalwidth = fitparameter(4);
    else
        fprintf('Input Error! No fitting parameter!')
    end
else
    fprintf('Input Error! No fitting parameter!')
end
        
if (nargin == 3)
    if isnumeric(varargin{nargin-1}) 
        gamma = varargin{nargin-1};
    else
        fprintf('Input Error! The third input is the initial value of gamma. ')
    end
elseif (nargin == 4)
    if isnumeric(varargin{nargin-1})&&isnumeric(varargin{nargin-2})
        min = varargin{nargin-1};
        gamma = varargin{nargin-2};
    else
        fprintf('Input Error! The third input is the initial value r and the fourth input is the termination tolerance.')
    end
elseif (nargin==1)
    temp_input = input('Please input the start pointing for fitting (unit: eV, e.g. 454):', 's');
    eprestart = str2num(temp_input);
    temp_input = input('Please input the wide of fitting window (unit: eV, e.g. 10):', 's');
    eprewidth = str2num(temp_input);
end
 
if ~isempty(eprestart)&&~isempty(eprewidth)
    efilter = ((escale >= eprestart) & (escale <=(eprestart + eprewidth)));
    escale = nonzeros(escale(efilter))';
    spectrum = nonzeros(spectrum(efilter))';
end

% maximum likelihood background estimation   
x = log(escale);
xw = sum(x.*spectrum)/sum(spectrum);
mr0 = sum(exp(-gamma.*x));
mr1 = sum(x.*exp(-gamma.*x));
fr = mr1/mr0-xw;
 
k = 0; 
while (abs(fr)>min)
    mr2 = sum(x.*x.*exp(-gamma.*x));
    dfdr = (mr1/mr0)^2-mr2/mr0;
    gamma = gamma-fr/dfdr;
    mr0 = sum(exp(-gamma.*x));
    mr1 = sum(x.*exp(-gamma.*x)); 
    fr = mr1/mr0-xw;
    k = k+1;     %  to avoid endless loop for bad experimental datasets. 
    if k >= 30
%         fprintf('%g\n', k);
        break;
    end
end
 
alpha = sum(spectrum)/mr0;  
 
back = alpha.*axisx.^(-gamma);    
hold on; plot(axisx, back,'k'); 

% error estimation
x = log(axisx);   
mlamda2 = sum(x(efilter).*x(efilter).*back(efilter))/sum(back(efilter));
mlamda1 = sum(x(efilter).*back(efilter))/sum(back(efilter));
b = back;
dbdr = -x.*back;
temp = back(efilter).*(x(efilter)-mlamda1).*(x(efilter)-mlamda1);
edeltab = sqrt((b.*b.*mlamda2+dbdr.*dbdr+2*b.*dbdr.*mlamda1)/sum(temp));

efilter2 = ((axisx >= signalstart) & (axisx <=(signalstart + signalwidth)));
mlamda2 = sum(x(efilter).*x(efilter).*back(efilter))/sum(back(efilter));
mlamda1 = sum(x(efilter).*back(efilter))/sum(back(efilter));
b = sum(back(efilter2));     
dbdr = sum(-x(efilter2).*back(efilter2)); 
temp = back(efilter2).*(x(efilter2)-mlamda1).*(x(efilter2)-mlamda1); 
edeltabsnr = (b*b*mlamda2+dbdr*dbdr+2*b*dbdr*mlamda1)/sum(temp);
 
t = sum(spectrum_i(efilter2));  
snr = (t-b)/sqrt(t+edeltabsnr);
 
core = spectrum_i - back;