hist2res.m

Go to the documentation of this file.

function [R]=hist2res(H,fun)
% Evaluates Histogram data
% [R]=hist2res(H)
%
% [y]=hist2res(H,fun)
%   estimates fun-statistic
%
% fun   'mean'  mean
%   'std'   standard deviation
%   'var'   variance
%   'sem'   standard error of the mean
%   'rms'   root mean square
%   'meansq' mean of squares
%   'sum'   sum
%   'sumsq' sum of squares
%   'CM#'   central moment of order #
%   'skewness' skewness 
%   'kurtosis' excess coefficient (Fisher kurtosis)
%
% see also: NaN/statistic
%
% REFERENCES:
% [1] C.L. Nikias and A.P. Petropulu "Higher-Order Spectra Analysis" Prentice Hall, 1993.
% [2] C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
% [3] http://www.itl.nist.gov/
% [4] http://mathworld.wolfram.com/
 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 
%   $Id: hist2res.m 2202 2009-10-27 12:06:45Z schloegl $
%   Copyright (c) 1996-2002,2006 by Alois Schloegl <a.schloegl@ieee.org>
%       This is part of the BIOSIG-toolbox http://biosig.sf.net/
 
 
if strcmp(H.datatype,'HISTOGRAM')
 
elseif strcmp(H.datatype,'qc:histo')
    HDR = H; 
    if isfield(H,'THRESHOLD'),
        TH  = H.THRESHOLD;
    else
        TH = repmat([-inf,inf],HDR.NS,1); 
    end;
    HIS = H.HIS; 
 
    % remove overflowing samples
    HIS.N = sumskipnan(HIS.H); 
    for k = 1:size(HIS.H,2);
        t = HIS.X(:,min(k,size(HIS.X,2))); 
        HIS.H(xor(t<=min(TH(k,:)), t>=max(TH(k,:))),k) = 0; 
    end; 
    Nnew = sumskipnan(HIS.H); 
    R.ratio_lost = 1-Nnew./HIS.N;
    HIS.N = Nnew; 
      
    % scale into physical values
    if H.FLAG.UCAL,
        %t = HIS.X;
        %for k=1:length(HDR.InChanSelect),
        %   HIS.X(:,k) = t(:,min(size(t,2),k))*HDR.Calib(k+1,k)+HDR.Calib(1,k);
        %end;
        HIS.X = [ones(size(HIS.X,1),1),repmat(HIS.X,1,size(HIS.H,2)./size(HIS.X,2))]*H.Calib;
    end;    
    H = HIS; 
else
        fprintf(2,'ERROR: arg1 is not a histogram\n');
        return;
end;
if nargin<2, fun=[]; end;
 
global FLAG_implicit_unbiased_estimation; 
%%% check whether FLAG was already defined 
if exist('FLAG_implicit_unbiased_estimation')~=1,
    FLAG_implicit_unbiased_estimation=[];
end;
%%% set DEFAULT value of FLAG
if isempty(FLAG_implicit_unbiased_estimation),
    FLAG_implicit_unbiased_estimation=logical(1);
end;
 
sz  = size(H.H)./size(H.X);
R.N     = sumskipnan(H.H,1);
R.SUM   = sumskipnan(H.H.*repmat(H.X,sz),1);
R.SSQ   = sumskipnan(H.H.*repmat(H.X.*H.X,sz),1);
%R.S3P  = sumskipnan(H.H.*repmat(H.X.^3,sz),1); % sum of 3rd power
R.S4P   = sumskipnan(H.H.*repmat(H.X.^4,sz),1); % sum of 4th power
%R.S5P  = sumskipnan(H.H.*repmat(H.X.^5,sz),1); % sum of 5th power
 
R.MEAN  = R.SUM./R.N;
R.MSQ   = R.SSQ./R.N;
R.RMS   = sqrt(R.MSQ);
R.SSQ0  = R.SSQ-R.SUM.*R.MEAN;      % sum square of mean removed
 
if FLAG_implicit_unbiased_estimation,
    n1  = max(R.N-1,0);         % in case of n=0 and n=1, the (biased) variance, STD and STE are INF
else
    n1  = R.N;
end;
 
R.VAR   = R.SSQ0./n1;               % variance (unbiased) 
R.STD   = sqrt(R.VAR);              % standard deviation
R.SEM   = sqrt(R.SSQ0./(R.N.*n1));  % standard error of the mean
R.SEV   = sqrt(n1.*(n1.*R.S4P./R.N+(R.N.^2-2*R.N+3).*(R.SSQ./R.N).^2)./(R.N.^3)); % standard error of the variance
R.Coefficient_of_variation = R.STD./R.MEAN;
 
R.CM2   = R.SSQ0./n1;
x       = repmat(H.X,sz) - repmat(R.MEAN,size(H.X,1),1);
R.CM3   = sumskipnan(H.H.*(x.^3),1)./n1;
R.CM4   = sumskipnan(H.H.*(x.^4),1)./n1;
%R.CM5  = sumskipnan(H.H.*(x.^5),1)./n1;
 
R.SKEWNESS = R.CM3./(R.STD.^3);
R.KURTOSIS = R.CM4./(R.VAR.^2)-3;
R.MAD = sumskipnan(H.H.*abs(x),1)./R.N; % mean absolute deviation
 
H.PDF = H.H./H.N(ones(size(H.H,1),1),:);
status=warning('off'); 
R.ENTROPY = -sumskipnan(H.PDF.*log2(H.PDF),1);
warning(status); 
R.QUANT = repmat(min(diff(H.X,[],1)),1,size(H.H,2)/size(H.X,2));
R.MAX = max(H.X); 
R.MIN = min(H.X); 
R.RANGE = R.MAX-R.MIN;
 
if ~isempty(fun),
        fun=upper(fun);
        if strncmp(fun,'CM',2) 
                oo = str2double(fun(3:length(fun)));
                R = sumskipnan(H.PDF.*(x.^oo),1);
        else                
        R = getfield(R,fun);
    end;
end;