## Copyright (C) 2016 Chris Adams ## ## 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 3 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, see . ## -*- texinfo -*- ## @deftypefn {} {@var{retval} =} statelevels (@var{input1}, @var{input2}) ## ## @seealso{} ## @end deftypefn #Uses histogram method. #1 Take a histogram #2 find ilow and ihigh, min and max non zero hist values #3 make two subhist #4 take mean of each for lower and upper function [S,H,B] = statelevels (A, varargin) #Setup arguments based on number of given arguments #arg1 is bins, arg2 is method #setup default values which can be later overwritten nBins = 100; method = 1; nargin = length(varargin); if(nargin>0) if(!isnumeric(varargin{1})) error('statelevels expects a number for argument 2'); end nBins = varargin{1}; end if(nargin>1) if(strcmp(varargin{2}, 'mode')) method = 1; elseif(strcmp(varargin{2}, 'mean')) method = 0; else error('Unexpected argument for method'); end end #Get size first for various uses sz = size(A); #Row column vector checking #Are any of the columns of size one? if(sz(1) != 1 && sz(2) != 1) error('statelevels expects a vector'); end #Switch to row if it's a clumn vector if(sz(2) == 1) X = A'; else X = A; end #Generate a histogram bin centres that are bounded by signal bounds upper = max(X); lower = min(X); #Generate hist vector based on the bounds in the histogram histVec = lower + ((1:nBins) - 0.5)' * (upper - lower) / nBins; #Get a histogram #Generate indexes for the histogram idx = ceil(nBins * (X-lower)/(upper-lower)); #Round zeros up to 1 idx(idx==0) = 1; #Pick only indexes in the range idx = idx(idx>=1 & idx<=nBins); #Initialise the histogram with zeros histogram = zeros(nBins, 1); #For every valid index increment the value if it hits for i=1:size(idx,2) histogram(idx(i)) = histogram(idx(i)) + 1; end #Find the bounds, lowest non zero and highest for i=1:nBins if(histogram(i) > 0) iLow = i; break; end end for i=nBins:-1:1 if(histogram(i) > 0) iHigh = i; break; end end #Define index pairs that define the lower and upper histogram slices. #We'll need these later lLow = iLow; lHigh = iLow + floor((iHigh - iLow)/2); uLow = iLow + floor((iHigh - iLow)/2); uHigh = iHigh; lHist = histogram(iLow:lHigh, 1); uHist = histogram(uLow:uHigh, 1); #What's the amplitude, size ratio (see ref) dy = (upper-lower) / size(histogram,1); #If mode (default) if (method==1) [~, iMax] = max(lHist); [~, iMin] = max(uHist); S1 = lower + dy * (lLow + iMax(1) - 1.5); S2 = lower + dy * (uLow + iMin(1) - 1.5); elseif (method==0) S1 = lower + dy * dot((lLow:lHigh)-0.5, lHist) / sum(lHist); S2 = lower + dy * dot((uLow:uHigh)-0.5, uHist) / sum(uHist); end S = [S1 S2]; H = histVec; B = histogram; endfunction