function [estimate,nbias,sigma,descriptor]=entropy(x,descriptor,approach,base)
%ENTROPY   Estimates the entropy of stationary signals with
%          independent samples using various approaches.
%   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X) or
%   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR) or
%   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH) or
%   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH,BASE)
%
%   ESTIMATE    : The entropy estimate
%   NBIAS       : The N-bias of the estimate
%   SIGMA       : The standard error of the estimate
%   DESCRIPTOR  : The descriptor of the histogram, seel alse ENTROPY
%
%   X           : The time series to be analyzed, a row vector
%   DESCRIPTOR  : Where DESCRIPTOR=[LOWERBOUND,UPPERBOUND,NCELL]
%     LOWERBOUND: Lowerbound of the histogram
%     UPPERBOUND: Upperbound of the histogram
%     NCELL     : The number of cells of the histogram       
%   APPROACH    : The method used, one of the following ones:
%     'unbiased': The unbiased estimate (default)
%     'mmse'    : The minimum mean square error estimate
%     'biased'  : The biased estimate
%   BASE        : The base of the logarithm; default e
%
%   See also: http://www.cs.rug.nl/~rudy/matlab/

%   R. Moddemeijer 
%   Copyright (c) by R. Moddemeijer
%   $Revision: 1.1 $  $Date: 2001/02/05 08:59:36 $



if nargin <1
   disp('Usage: [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X)')
   disp('       [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR)')
   disp('       [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH)')
   disp('       [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH,BASE)')
   disp('Where: DESCRIPTOR = [LOWERBOUND,UPPERBOUND,NCELL]')
   return
end

% Some initial tests on the input arguments

[NRowX,NColX]=size(x);

if NRowX~=1
  error('Invalid dimension of X');
end;

if nargin>4
  error('Too many arguments');
end;

if nargin==1
  [h,descriptor]=histogram(x);
end;

if nargin>=2
  [h,descriptor]=histogram(x,descriptor);
end;

if nargin<3
  approach='unbiased';
end;

if nargin<4
  base=exp(1);
end;

lowerbound=descriptor(1);
upperbound=descriptor(2);
ncell=descriptor(3);

estimate=0;
sigma=0;
count=0;

for n=1:ncell
  if h(n)~=0 
    logf=log(h(n));
  else
    logf=0;
  end;
  count=count+h(n);
  estimate=estimate-h(n)*logf;
  sigma=sigma+h(n)*logf^2;
end;

% biased estimate

estimate=estimate/count;
sigma   =sqrt( (sigma/count-estimate^2)/(count-1) );
estimate=estimate+log(count)+log((upperbound-lowerbound)/ncell);
nbias   =-(ncell-1)/(2*count);

% conversion to unbiased estimate

if approach(1)=='u'
  estimate=estimate-nbias;
  nbias=0;
end;

% conversion to minimum mse estimate

if approach(1)=='m'
  estimate=estimate-nbias;
  nbias=0;
  lambda=estimate^2/(estimate^2+sigma^2);
  nbias   =(1-lambda)*estimate;
  estimate=lambda*estimate;
  sigma   =lambda*sigma;
end;

% base transformation

estimate=estimate/log(base);
nbias   =nbias   /log(base);
sigma   =sigma   /log(base);