可以参考下面的函数
- 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);
复制代码
|