1function [Y0,EVENT] = OAHE(S,Fs,varargin)
2% OAHE detectes obstructive Apnea/Hypopnea event
3%
4%    [Y,EVENT] = OAHE(X,Fs)
5%    [Y,EVENT] = OAHE(filename,CHAN)
6%   ... = OAHE(... ,'-o',outputFilename)
7%   ... = OAHE(... ,'-e',eventFilename)
8%
9% INPUT:
10%       X   respiratory channel
11%       Fs  sampleing rate
12%       filename        source filename
13%       CHAN            respiratory channels for calculating OAHE
14%	outputFilename
15%		name of file for storing the resulting data with the
16%		detected spikes and bursts in GDF format.
17%	eventFilename
18%		filename to store event inforamation in GDF format. this is similar to
19%		the outputFile, except that the signal data is not included and is, therefore,
20%		much smaller than the outputFile
21%
22% OUTPUT:
23%       Y       detection trace
24%       EVENT   event structure as used in BIOSIG
25%
26% see also: SVIEWER, SLOAD
27%
28%
29% REFERENCES:
30% [1] AASM Task Force. Sleep-Related Breathing Disorders in Adults:
31%       Recommendations for Syndrome Definition, and Measurement Techniques in Clinical Research. Sleep, 22(5), 1999.
32% [2] Meoli AL, Casey KR, Clark RW, Coleman JA Jr, Fayle RW, Troell RJ, Iber C; Clinical Practice Review Committee.
33%       Hypopnea in Sleep-Disordered Breathing in Adults. Sleep. 2001 Jun 15;24(4):469-70.
34% [3] Penzel T.,  Brandenburg U., Fischer J., Jobert M., Kurella B., Mayer G., Nioewerth H.J., Peter J.H., Pollmächer T., Schäfer T., Steinberg R., Trowitzsch E., Warmuth R., Weeß H.-G., Wölk C., Zulley J.,
35%       Empfehlungen zur computerunterstützen Aufzeichnung und Auswertung von Polygraphien. Somnologie, 2, 42-48, 1998.
36% [4] Ross SD, Sheinhait IA, Harrison KJ, Kvasz M, Connelly JE, Shea SA, Allen IE.
37%       Systematic review and meta-analysis of the literature regarding the diagnosis of sleep apnea. Sleep. 2000 Jun 15;23(4):519-32.
38% [5] Schafer T. Methodik der Atmungsmessung im Schlaf: Kapneographie zur Beurteilung der Ventilation.
39%       Biomed Tech (Berl). 2003 Jun; 48(6):170-5.
40% [6] Thurnheer R, Xie X, Bloch KE. Accuracy of nasal cannula pressure recordings for assessment of ventilation during sleep.
41%       Am J Respir Crit Care Med. 2001 Nov 15;164(10 Pt 1):1914-9.
42% [7] Whitney CW, Gottlieb DJ, Redline S, Norman RG, Dodge RR, Shahar E, Surovec S, Nieto FJ.
43%       Reliability of scoring respiratory disturbance indices and sleep staging. Sleep. 1998 Nov 1;21(7):749-57.
44% [8] Schlögl A, Kemp B, Penzel T, Kunz D, Himanen SL, Varri A, Dorffner G, Pfurtscheller
45%       G. Quality control of polysomnographic sleep data by histogram and entropy analysis.
46%       Clin Neurophysiol. 1999 Dec;110(12):2165-70.
47
48%	$Revision: 1.3 $
49%	$Id$
50%	Copyright (C) 2003,2005 by Alois Schloegl <alois.schloegl@gmail.com>
51
52% This library is free software; you can redistribute it and/or
53% modify it under the terms of the GNU Library General Public
54% License as published by the Free Software Foundation; either
55% Version 3 of the License, or (at your option) any later version.
56%
57% This library is distributed in the hope that it will be useful,
58% but WITHOUT ANY WARRANTY; without even the implied warranty of
59% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
60% Library General Public License for more details.
61%
62% You should have received a copy of the GNU Library General Public
63% License along with this library; if not, write to the
64% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
65% Boston, MA  02111-1307, USA.
66
67%%%%% default settings %%%%%
68outFile = [];
69evtFile = [];
70
71%%%%% analyze input arguments %%%%%
72k = 1;
73while k <= length(varargin)
74	if ischar(varargin{k})
75		if (strcmp(varargin{k},'-o'))
76			k = k + 1;
77			outFile = varargin{k};
78		elseif (strcmp(varargin{k},'-e'))
79			k = k + 1;
80			evtFile = varargin{k};
81		else
82			warning(sprintf('unknown input argument <%s>- ignored',varargin{k}));
83		end;
84	end;
85	k = k+1;
86end;
87
88if ischar(S)
89        [S,HDR]=sload(S,Fs);
90        Fs = HDR.SampleRate;
91else
92        HDR.InChanSelect = 1:size(X,1);
93end;
94
95[nr,nc]=size(S);
96
97nseg  = ceil (nr/(120*Fs));
98nseg2 = floor(nr/(  5*Fs));
99
100S = [repmat(nan,125*Fs,nc); S; repmat(nan, nseg*120*Fs-nr, nc)];
101for k = 1:nc,
102        X = reshape(S(:,k),5*Fs,24*nseg+25);
103
104        hi_baseline = max(X);
105        lo_baseline = min(X);
106
107        ix = hankel(1:nseg2,nseg2:nseg2+25)';
108
109        blh = median(hi_baseline(ix(1:24,:)));
110        bll = median(lo_baseline(ix(1:24,:)));
111
112        mh  = max(hi_baseline(ix(25:26,:)));
113        ml  = min(lo_baseline(ix(25:26,:)));
114
115        Y   = real((blh-bll)' > 2*(mh-ml)');
116        Y(any(isnan([mh;ml;blh;bll]))) = NaN;
117
118        Y0(:,k) = Y;
119end;
120
121%%%%% generate event encoding
122Y = Y0;
123ix = isnan(Y(1,:));
124Y(1,ix) = 0;
125for k=2:length(Y),
126	ix = isnan(Y(k,:));
127	Y(k,ix) = Y(k-1,ix);
128end;
129Y(end+1,:)=0;
130
131EVENT.SampleRate = Fs;
132N = 0;
133for k=1:nc;
134	tmp1 = diff(Y(:,k));
135        tmp2 = find(tmp1>0);
136	L    = length(tmp2);
137	EVENT.TYP(N+1:N+L,1) = hex2dec('0401');
138	EVENT.POS(N+1:N+L,1) = tmp2*Fs;
139	EVENT.DUR(N+1:N+L,1) = (find(tmp1<0)-tmp2)*Fs;
140	EVENT.CHN(N+1:N+L,1) = HDR.InChanSelect(k);
141end;
142
143
144%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145%	Output
146%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147
148if ~isempty(outFile)
149		%%% write data to output
150		HDR.TYPE  = 'GDF';
151		HDR.VERSION = 3;
152		%[p,f,e]=fileparts(fn);
153		HDR.FILE = [];
154		HDR.FileName  = outFile;
155		HDR.FILE.Path = '';
156		HDR.PhysMax   = max(s);
157		HDR.PhysMin   = min(s);
158		HDR.DigMax(:) =  2^15-1;
159		HDR.DigMin(:) = -2^15;
160		HDR.GDFTYP(:) = 3;
161		HDR.FLAG.UCAL = 0;
162		HDR.NRec = size(s,1);
163		HDR.SPR = 1;
164		HDR.NS  = size(s,2);
165		HDR.Dur = 1/HDR.SampleRate;
166		HDR = rmfield(HDR,'AS');
167		HDR.EVENT = EVENT;
168		HDR = sopen(HDR,'w');
169		if (HDR.FILE.FID < 0)
170			fprintf(2,'Warning can not open file <%s> - GDF file can not be written\n',HDR.FileName);
171		else
172			HDR = swrite(HDR,s);
173			HDR = sclose(HDR);
174		end;
175end;
176
177if ~isempty(evtFile)
178		%%% write data to output
179		HDR.TYPE  = 'EVENT';
180		HDR.VERSION = 3;
181		%[p,f,e]=fileparts(fn);
182		HDR.FILE = [];
183		HDR.FileName  = evtFile;
184		HDR.FILE.Path = '';
185		HDR.NRec = 0;
186		HDR.SPR = 0;
187		HDR.Dur = 1/HDR.SampleRate;
188		HDR = rmfield(HDR,'AS');
189		HDR = sopen(HDR, 'w');
190		if (HDR.FILE.FID<0)
191			fprintf(2,'Warning can not open file <%s> - EVT file can not be written\n', HDR.FileName);
192		else
193			HDR = sclose(HDR);
194 		end;
195end;
196
197
198