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