1function [LIM1,LIM2,LIM3,h0] = hist2limits(H,TH) 2% HIST2LIMITS returns the threshold for detecting saturation artifacts. 3% 4% Saturation thresholds can be obtained from the histogram [1]. This 5% routine tries to obtain the saturation threshold in an automated way. 6% 7% The routine was tested with the histograms of 528 recordings with 8% three respiratory channels each. 9% 10% 11% Reference(s): 12% [1] A. Schlögl, B. Kemp, T. Penzel, D. Kunz, S.-L. Himanen,A. Värri, G. Dorffner, G. Pfurtscheller. 13% Quality Control of polysomnographic Sleep Data by Histogram and Entropy Analysis. 14% Clin. Neurophysiol. 1999, Dec; 110(12): 2165 - 2170. 15 16 17% Copyright (C) 1999-2003,2021 by Alois Schloegl <alois.schloegl@gmail.com> 18 19% This program is free software; you can redistribute it and/or 20% modify it under the terms of the GNU General Public License 21% as published by the Free Software Foundation; either version 2 22% of the License, or (at your option) any later version. 23% 24% This program is distributed in the hope that it will be useful, 25% but WITHOUT ANY WARRANTY; without even the implied warranty of 26% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27% GNU General Public License for more details. 28% 29% You should have received a copy of the GNU General Public License 30% along with this program; if not, write to the Free Software 31% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 32 33 34 35 36if isnumeric(H) 37 H = histo2(H); 38end; 39 40for k = 1:size(H.H,2), 41 h = H.H(:,k); 42 x = H.X(:,min(k,size(H.X,2))); 43 44 tmp = find(h>0); 45 h(tmp([1,length(tmp)])) = 0; % remove max and min values; (necessary for H recordings and some B) 46 N = sum(h); 47 48% Lim = H.X(find(h>0),1); % calculate limit values of remaining Histogram 49 Lim = x(find(h>0)); % calculate limit values of remaining Histogram 50 Lim = [max(Lim),min(Lim)]; 51 LIM1(:,k) = sort(mean(Lim)+([1;-1]*TH/2)*abs(diff(Lim))); % take range between 10% and 90% of total range. 52 53 mu = x'*h/N; 54 sd = sqrt(((x-mu)'.^2*h)/N); 55 56 LIM2(:,k) = sort(sd*[1;-1]*norminv(2/N)+mu); 57 58 h0 = sum(h)*normpdf(x,mu,sd); 59 60 r0 = (h./max(h0,1e-2)); 61 [s,ix] = sort(-r0); 62 tmp = ix(x(ix)>mu); 63 LIM3(:,k) = [NaN;NaN]; 64 if 1;r0(tmp(1))>1e2; 65 LIM3(1,k) = x(tmp(1)); 66 end; 67 tmp = ix(x(ix)<mu); 68 if 1;r0(tmp(1))>1e3; 69 LIM3(2,k) = x(tmp(1)); 70 end; 71 72 LIM3(:,k) = sort(mean(LIM3(:,k))+([1;-1]*TH/2)*abs(diff(LIM3(:,k)))); % take range between 10% and 90% of total range. 73 74 75 LIM0(:,k) = [max(LIM1(1,k),LIM2(1,k));min(LIM1(2,k),LIM2(2,k))]; 76 LIM = LIM1; 77end; 78 79 80