1 function [R]=hist2res(H,fun)
2 % Evaluates Histogram data
6 % estimates fun-statistic
9 % 'std' standard deviation
11 % 'sem' standard error of the mean
12 % 'rms' root mean square
13 % 'meansq' mean of squares
15 % 'sumsq' sum of squares
16 % 'CM#' central moment of order #
18 % 'kurtosis' excess coefficient (Fisher kurtosis)
20 % see also: NaN/statistic
23 % [1] C.L. Nikias and A.P. Petropulu "Higher-Order Spectra Analysis" Prentice Hall, 1993.
24 % [2] C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
25 % [3] http://www.itl.nist.gov/
26 % [4] http://mathworld.wolfram.com/
28 % This program is free software; you can redistribute it and/or
29 % modify it under the terms of the GNU General Public License
30 % as published by the Free Software Foundation; either version 2
31 % of the License, or (at your option) any later version.
33 % This program is distributed in the hope that it will be useful,
34 % but WITHOUT ANY WARRANTY; without even the implied warranty of
35 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 % GNU General Public License for more details.
38 % You should have received a copy of the GNU General Public License
39 % along with this program; if not, write to the Free Software
40 % Foundation, Inc., 51 Franklin Street - Fifth Floor, Boston, MA 02110-1301, USA.
42 % $Id: hist2res.m 9387 2011-12-15 10:42:14Z schloegl $
43 % Copyright (c) 1996-2002,2006 by Alois Schloegl <alois.schloegl@gmail.com>
44 % This function is part of the NaN-toolbox
45 % http://pub.ist.ac.at/~schloegl/matlab/NaN/
48 if strcmp(H.datatype,'HISTOGRAM'),
50 elseif strcmp(H.datatype,'qc:histo')
52 if isfield(H,'THRESHOLD'),
55 TH = repmat([-inf,inf],HDR.NS,1);
59 % remove overflowing samples
60 HIS.N = sumskipnan(HIS.H);
61 for k = 1:size(HIS.H,2);
62 t = HIS.X(:,min(k,size(HIS.X,2)));
63 HIS.H(xor(t<=min(TH(k,:)), t>=max(TH(k,:))),k) = 0;
65 Nnew = sumskipnan(HIS.H);
66 R.ratio_lost = 1-Nnew./HIS.N;
69 % scale into physical values
72 %for k=1:length(HDR.InChanSelect),
73 % HIS.X(:,k) = t(:,min(size(t,2),k))*HDR.Calib(k+1,k)+HDR.Calib(1,k);
75 HIS.X = [ones(size(HIS.X,1),1),repmat(HIS.X,1,size(HIS.H,2)./size(HIS.X,2))]*H.Calib;
79 fprintf(2,'ERROR: arg1 is not a histogram\n');
82 if nargin<2, fun=[]; end;
84 global FLAG_implicit_unbiased_estimation;
85 %%% check whether FLAG was already defined
86 if ~exist('FLAG_implicit_unbiased_estimation','var'),
87 FLAG_implicit_unbiased_estimation=[];
89 %%% set DEFAULT value of FLAG
90 if isempty(FLAG_implicit_unbiased_estimation),
91 FLAG_implicit_unbiased_estimation=logical(1);
94 sz = size(H.H)./size(H.X);
95 R.N = sumskipnan(H.H,1);
96 R.SUM = sumskipnan(H.H.*repmat(H.X,sz),1);
97 R.SSQ = sumskipnan(H.H.*repmat(H.X.*H.X,sz),1);
98 %R.S3P = sumskipnan(H.H.*repmat(H.X.^3,sz),1); % sum of 3rd power
99 R.S4P = sumskipnan(H.H.*repmat(H.X.^4,sz),1); % sum of 4th power
100 %R.S5P = sumskipnan(H.H.*repmat(H.X.^5,sz),1); % sum of 5th power
105 R.SSQ0 = R.SSQ-R.SUM.*R.MEAN; % sum square of mean removed
107 if FLAG_implicit_unbiased_estimation,
108 n1 = max(R.N-1,0); % in case of n=0 and n=1, the (biased) variance, STD and STE are INF
113 R.VAR = R.SSQ0./n1; % variance (unbiased)
114 R.STD = sqrt(R.VAR); % standard deviation
115 R.SEM = sqrt(R.SSQ0./(R.N.*n1)); % standard error of the mean
116 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
117 R.Coefficient_of_variation = R.STD./R.MEAN;
120 x = repmat(H.X,sz) - repmat(R.MEAN,size(H.X,1),1);
121 R.CM3 = sumskipnan(H.H.*(x.^3),1)./n1;
122 R.CM4 = sumskipnan(H.H.*(x.^4),1)./n1;
123 %R.CM5 = sumskipnan(H.H.*(x.^5),1)./n1;
125 R.SKEWNESS = R.CM3./(R.STD.^3);
126 R.KURTOSIS = R.CM4./(R.VAR.^2)-3;
127 R.MAD = sumskipnan(H.H.*abs(x),1)./R.N; % mean absolute deviation
129 H.PDF = H.H./H.N(ones(size(H.H,1),1),:);
130 status=warning('off');
131 R.ENTROPY = -sumskipnan(H.PDF.*log2(H.PDF),1);
133 R.QUANT = repmat(min(diff(H.X,[],1)),1,size(H.H,2)/size(H.X,2));
136 R.RANGE = R.MAX-R.MIN;
140 if strncmp(fun,'CM',2)
141 oo = str2double(fun(3:length(fun)));
142 R = sumskipnan(H.PDF.*(x.^oo),1);