1 function [HDR]=
eeg2hist(FILENAME,CHAN);
2 % EEG2HIST histogram analysis based on [1]
3 % It displays the histograms of the recorded data,
4 % and allows editing of the thresholds
for
5 % the automated overflow detection.
7 % [HDR]=EEG2HIST(FILENAME);
9 % input: FILENAME EEG-File
12 % HDR header information
13 % HDR.HIS histograms
for each channel
14 % HDR.RES summary statistics based on the histogram analysis [1]
15 % HDR.THRESHOLD Threshold values
for overflow detection
18 % see also: SOPEN, SLOAD, HIST2RES
21 % [1] A. Schlögl, B. Kemp, T. Penzel, D. Kunz, S.-L. Himanen,A. Värri, G. Dorffner, G. Pfurtscheller.
22 % Quality Control of polysomnographic Sleep Data by Histogram and EntropyAnalysis.
23 % Clin. Neurophysiol. 1999, Dec; 110(12): 2165 - 2170.
26 % [2] A. Schlögl, G. Klösch, W. Koele, J. Zeitlhofer, P.Rappelsberger, G. Pfurtscheller
27 % Qualitätskontrolle von Biosignalen,
28 % Jahrestagung der Österreichischen Gesellschaft für Klinische Neurophysiologie, 27. Nov. 1999, Vienna.
32 % [4] A. Schlögl, Time Series Analysis toolbox
for Matlab. 1996-2003
35 % $Id:
eeg2hist.m 2713 2011-06-16 08:55:51Z schloegl $
36 % Copyright (C) 2002,2003,2006,2007,2010,2011 by Alois Schloegl <alois.schloegl@gmail.com>
37 % This is part of the BIOSIG-toolbox http:
39 % BioSig is free software: you can redistribute it and/or modify
40 % it under the terms of the GNU General Public License as published by
41 % the Free Software Foundation, either version 3 of the License, or
42 % (at your option) any later version.
44 % BioSig is distributed in the hope that it will be useful,
45 % but WITHOUT ANY WARRANTY; without even the implied warranty of
46 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
47 % GNU General Public License
for more details.
49 % You should have received a copy of the GNU General Public License
50 % along with BioSig. If not, see <http://www.gnu.org/licenses/>.
54 if nargin<2, CHAN=0; end;
62 [s, HDR] =
sload(FILENAME,0,
'OVERFLOWDETECTION',
'OFF',
'UCAL',
'ON');
63 %%% FIXME:
if there is only a single value exceeding the threshold,
this is not visualized
68 HDR =
sopen(FILENAME,
'r',CHAN,
'UCAL'); % open EEG file in uncalibrated mode (no scaling of the data)
70 fprintf(2,
'EEG2HIST: couldnot open file %s.\n',FILENAME);
74 HDR.FLAG.UCAL = 1; %
do not scale
75 HDR.FLAG.OVERFLOWDETECTION = 0; % OFF
76 if CHAN<1, CHAN=1:HDR.NS; end;
78 H.datatype=
'HISTOGRAM';
82 H.H = zeros(2^16,HDR.NS);
84 if exist(
'OCTAVE_VERSION') > 2,
85 for k = double(s(:,l)
')+2^15+1, H.H(k,l) = H.H(k,l)+1; end;
87 H.H(:,l)=sparse(double(s(:,l)')+2^15+1,1,1,2^16,1);
90 H.X = [-2^15:2^15-1]
'; %int16
91 %tmp = find(any(H.H,2));
92 %H.X = tmp-2^15-1; %int16
95 elseif all(HDR.GDFTYP==4)
98 H.H = zeros(2^16,HDR.NS);
100 if exist('OCTAVE_VERSION
') > 2,
101 for k = s(:,l)'+1, H.H(k,l) = H.H(k,l)+1; end;
103 H.H(:,l)=sparse(s(:,l)
'+1,1,1,2^16,1);
107 %tmp = find(any(H.H,2));
108 %H.X = tmp-1; %uint16
112 elseif all(HDR.GDFTYP==(255+24)) %strcmp(HDR.TYPE,
'BDF'),
113 [s,HDR] =
sread(HDR);
115 H.H = sparse(2^24,HDR.NS);
117 H.H(:,l)=sparse(s(:,l)
'+2^23+1,1,1,2^24,1);
118 %for k = s(:,l)'+2^23+1, H.H(k,l) = H.H(k,l)+1; end;
120 tmp = find(any(H.H,2));
125 elseif (strcmp(HDR.TYPE,
'EDF') || strcmp(HDR.TYPE,
'GDF') || strcmp(HDR.TYPE,
'ACQ')) && all(HDR.GDFTYP(1)==HDR.GDFTYP) && (HDR.GDFTYP(1)<=3)
126 NoBlks=ceil(60*HDR.SampleRate/HDR.SPR);
128 if isfield(HDR.AS,
'SPR')
129 bi=[0;cumsum(HDR.AS.SPR)];
132 HDR.AS.spb = HDR.SPR;
134 CHAN = HDR.InChanSelect;
137 H.H = zeros(2^16,ns);
140 while (k<HDR.NRec) && ~feof(HDR.FILE.FID)
141 % READ BLOCKS of DATA
142 [S, count] = fread(HDR.FILE.FID,[HDR.AS.spb,NoBlks],
gdfdatatype(HDR.GDFTYP(1)));
143 if 0, count < HDR.AS.spb*NoBlks
144 fprintf(2,' Warning EEG2HIST: read error, only %i of %i read\n', count, HDR.AS.spb*NoBlks);
150 if exist('OCTAVE_VERSION','builtin'),
151 for k=reshape(S(bi(CHAN(l))+1:bi(CHAN(l)+1),:),1,HDR.SPR(l)*NoBlks)+2^15+1, h(k,l) = h(k,l)+1; end;
153 h = sparse(S(bi(CHAN(l))+1:bi(CHAN(l)+1),:)+2^15+1,1,1,2^16,1);
155 H.H(:,l) = H.H(:,l) + h;
160 tmp = find(any(H.H,2));
161 H.X = [-2^15:2^15-1]';
162 % (tmp-2^15-1); %int16
165 elseif ~strcmp(HDR.TYPE,'unknown')
169 fprintf(2,'WARNING EEG2HIST: data is not integer.\n');
173 fprintf(2,'EEG2HIST: format %s not implemented yet.\n',HDR.TYPE);
180 %%% complete histogram and display it
181 H.N = sumskipnan(H.H);
182 %H.X = [ones(size(H.X,1),1),repmat(H.X,1,length(CHAN))]*HDR.Calib; %int16
184 if strcmp(HDR.TYPE,'GDF') || strcmp(HDR.TYPE,'alpha')
185 HDR.THRESHOLD = [HDR.DigMin(:),HDR.DigMax(:)];
188 if ~isfield(HDR,'THRESHOLD')
189 HDR.THRESHOLD = inf*ones(HDR.NS,1)./HDR.Cal(:) *[-1,1];
192 N=ceil(sqrt(size(H.H,2)));
193 for K = 1:size(H.H,2);
194 t = H.X(:,min(K,size(H.X,2)))*HDR.Calib(K+1,K)+HDR.Calib(1,K);
197 mu = sumskipnan(t.*h)/H.N(K);
199 sd2= sumskipnan((x.^2).*h)./H.N(K);
202 elseif isfield(HDR,'THRESHOLD'),
203 MaxMin=HDR.THRESHOLD(K,[2,1])*HDR.Calib(K+1,K)+HDR.Calib(1,K);
206 MaxMin=[max(t) min(t)];
208 xrange = [min(t(h>0)),max(t(h>0))];
209 xrange = xrange + [-1,1]*(diff(xrange)/2+eps);
211 a(K)= subplot(ceil(size(H.H,2)/N),N,K);
213 dQ = min(tmp(abs(tmp)>eps));
214 tmp = sqrt(sum(h(h>0))/sqrt(2*pi*sd2)*dQ);
217 h2((t>min(MaxMin)) & (t<max(MaxMin)))=NaN;
219 semilogy(t,[h+.01,exp(-((t-mu).^2)/(sd2*2))/sqrt(2*pi*sd2)*sum(h(h>0))*dQ],'-',t,h2+.01,'r',mu+sqrt(sd2)*[-5 -3 -1 0 1 3 5]',tmp*ones(7,1),'+-',MaxMin,tmp,'rx');
220 %v=axis; v=[v(1:2) 1 max(h)]; axis(v);
221 v=axis; v=[xrange 1-eps max(h)]; axis(v);
226 figure(gcf); %set(gcf,
'CurrentFigure');
228 fprintf(1,
'\nEEG2HIST:>Now You can modify the thresholds with mouse clicks.');
229 fprintf(1,
'\nEEG2HIST:>When you are finished, PRESS ANY KEY on the keyboard.');
230 fprintf(1,
'\nEEG2HIST:> (make sure the focus is on the figure window).');
233 if ~isfield(HDR,
'THRESHOLD'),
234 HDR.THRESHOLD = repmat(NaN,HDR.NS,2);
240 if isempty(b)
break; end;
245 t = H.X(:,min(K,size(H.X,2)))*HDR.Calib(K+1,K)+HDR.Calib(1,K);
246 %HISTO=hist2pdf(HISTO);
249 if K~=K0, MaxMin = [NaN,NaN]; end;
252 tmp((t>min(MaxMin)) & (t<max(MaxMin)))=NaN;
253 semilogy(t,[h+.01],'b-',t,tmp+.01,'-r');
255 HDR.THRESHOLD(K,1:2) = (ix(1:2)-HDR.Calib(1,K))/HDR.Calib(K+1,K);
257 v=[v(1:2) 1-eps max(h)]; axis(v);
261 fprintf(1,
' <FINISHED>\n');
264 %%% complete
return argument
267 HDR.datatype =
'qc:histo';