TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
eeg2hist.m
Go to the documentation of this file.
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.
6 %
7 % [HDR]=EEG2HIST(FILENAME);
8 %
9 % input: FILENAME EEG-File
10 % CHAN Channel select
11 % output:
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
16 %
17 %
18 % see also: SOPEN, SLOAD, HIST2RES
19 %
20 % REFERENCES:
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.
24 % http://dx.doi.org/10.1016/S1388-2457(99)00172-8
25 %
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.
29 %
30 % [3] http://pub.ist.ac.at/~schloegl/lectures/Q/index.htm
31 %
32 % [4] A. Schlögl, Time Series Analysis toolbox for Matlab. 1996-2003
33 % http://pub.ist.ac.at/~schloegl/matlab/tsa/
34 
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://biosig.sf.net/
38 %
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.
43 %
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.
48 %
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/>.
51 
52 
53 MODE=0;
54 if nargin<2, CHAN=0; end;
55 if ischar(CHAN),
56  MODE = 1;
57  CHAN = 0;
58 end;
59 
60 
61 if 0,
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
64 H = histo2(s)
65 
66 
67 else
68 HDR = sopen(FILENAME,'r',CHAN,'UCAL'); % open EEG file in uncalibrated mode (no scaling of the data)
69 if HDR.FILE.FID<0,
70  fprintf(2,'EEG2HIST: couldnot open file %s.\n',FILENAME);
71  return;
72 end;
73 
74 HDR.FLAG.UCAL = 1; % do not scale
75 HDR.FLAG.OVERFLOWDETECTION = 0; % OFF
76 if CHAN<1, CHAN=1:HDR.NS; end;
77 
78 H.datatype='HISTOGRAM';
79 if all(HDR.GDFTYP==3)
80  [s,HDR]=sread(HDR);
81 
82  H.H = zeros(2^16,HDR.NS);
83  for l = 1: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;
86  else
87  H.H(:,l)=sparse(double(s(:,l)')+2^15+1,1,1,2^16,1);
88  end;
89  end;
90  H.X = [-2^15:2^15-1]'; %int16
91  %tmp = find(any(H.H,2));
92  %H.X = tmp-2^15-1; %int16
93  %H.H = H.H(tmp,:);
94 
95 elseif all(HDR.GDFTYP==4)
96  [s,HDR]=sread(HDR);
97 
98  H.H = zeros(2^16,HDR.NS);
99  for l = 1:HDR.NS,
100  if exist('OCTAVE_VERSION') > 2,
101  for k = s(:,l)'+1, H.H(k,l) = H.H(k,l)+1; end;
102  else
103  H.H(:,l)=sparse(s(:,l)'+1,1,1,2^16,1);
104  end;
105  end;
106  H.X = [0:2^16-1]';
107  %tmp = find(any(H.H,2));
108  %H.X = tmp-1; %uint16
109  %H.H = H.H(tmp,:);
110 
111 
112 elseif all(HDR.GDFTYP==(255+24)) %strcmp(HDR.TYPE,'BDF'),
113  [s,HDR] = sread(HDR);
114 
115  H.H = sparse(2^24,HDR.NS);
116  for l = 1: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;
119  end;
120  tmp = find(any(H.H,2));
121  H.X = tmp-2^23-1;
122  H.H = H.H(tmp,:);
123 
124 
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);
127 
128  if isfield(HDR.AS,'SPR')
129  bi=[0;cumsum(HDR.AS.SPR)];
130  else
131  bi=0:HDR.NS;
132  HDR.AS.spb = HDR.SPR;
133  end;
134  CHAN = HDR.InChanSelect;
135  ns=length(CHAN);
136 
137  H.H = zeros(2^16,ns);
138 
139  k=0;
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);
145  end;
146 
147  %%%%% HISTOGRAM
148  for l=1:ns,
149  h = zeros(2^16,1);
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;
152  else
153  h = sparse(S(bi(CHAN(l))+1:bi(CHAN(l)+1),:)+2^15+1,1,1,2^16,1);
154  end;
155  H.H(:,l) = H.H(:,l) + h;
156  end;
157 
158  k=k+NoBlks;
159  end; % WHILE
160  tmp = find(any(H.H,2));
161  H.X = [-2^15:2^15-1]';
162  % (tmp-2^15-1); %int16
163  %H.H = H.H(tmp,:);
164 
165 elseif ~strcmp(HDR.TYPE,'unknown')
166  [s,HDR]=sread(HDR);
167  H = histo2(s);
168  if any(HDR.GDFTYP>6)
169  fprintf(2,'WARNING EEG2HIST: data is not integer.\n');
170  end;
171 
172 else
173  fprintf(2,'EEG2HIST: format %s not implemented yet.\n',HDR.TYPE);
174 end;
175 HDR = sclose(HDR);
176 end;
177 
178 HDR.HIS = H;
179 
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
183 
184 if strcmp(HDR.TYPE,'GDF') || strcmp(HDR.TYPE,'alpha')
185  HDR.THRESHOLD = [HDR.DigMin(:),HDR.DigMax(:)];
186 end;
187 
188  if ~isfield(HDR,'THRESHOLD')
189  HDR.THRESHOLD = inf*ones(HDR.NS,1)./HDR.Cal(:) *[-1,1];
190  end;
191 
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);
195  h = H.H(:,K);
196 
197  mu = sumskipnan(t.*h)/H.N(K);
198  x = t-mu;
199  sd2= sumskipnan((x.^2).*h)./H.N(K);
200 
201  if 0,
202  elseif isfield(HDR,'THRESHOLD'),
203  MaxMin=HDR.THRESHOLD(K,[2,1])*HDR.Calib(K+1,K)+HDR.Calib(1,K);
204  else
205  %[max(t) min(t)],
206  MaxMin=[max(t) min(t)];
207  end;
208  xrange = [min(t(h>0)),max(t(h>0))];
209  xrange = xrange + [-1,1]*(diff(xrange)/2+eps);
210 
211  a(K)= subplot(ceil(size(H.H,2)/N),N,K);
212  tmp = diff(t);
213  dQ = min(tmp(abs(tmp)>eps));
214  tmp = sqrt(sum(h(h>0))/sqrt(2*pi*sd2)*dQ);
215 
216  h2 = h;
217  h2((t>min(MaxMin)) & (t<max(MaxMin)))=NaN;
218 
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);
222  title(HDR.Label{K});
223  end;
224 
225 if MODE,
226  figure(gcf); %set(gcf,'CurrentFigure');
227  drawnow;
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).');
231 
232  % modify Threshold
233  if ~isfield(HDR,'THRESHOLD'),
234  HDR.THRESHOLD = repmat(NaN,HDR.NS,2);
235  end;
236  b = 0;
237  K0 = 0;
238  while (b<2),
239  [x,y,b] = ginput(1);
240  if isempty(b) break; end;
241 
242  v=axis;
243  K = find(a==gca);
244  %min(K,size(H.X,2))
245  t = H.X(:,min(K,size(H.X,2)))*HDR.Calib(K+1,K)+HDR.Calib(1,K);
246  %HISTO=hist2pdf(HISTO);
247  h = H.H(:,K);
248 
249  if K~=K0, MaxMin = [NaN,NaN]; end;
250  MaxMin = [MaxMin,x];
251  tmp = h;
252  tmp((t>min(MaxMin)) & (t<max(MaxMin)))=NaN;
253  semilogy(t,[h+.01],'b-',t,tmp+.01,'-r');
254  ix = sort(MaxMin);
255  HDR.THRESHOLD(K,1:2) = (ix(1:2)-HDR.Calib(1,K))/HDR.Calib(K+1,K);
256  K0 = K;
257  v=[v(1:2) 1-eps max(h)]; axis(v);
258  title(HDR.Label{K});
259  end;
260 end;
261 fprintf(1,' <FINISHED>\n');
262 
263 
264 %%% complete return argument
265 HDR.HIS = H;
266 %HDR.RES = hist2res(H);
267 HDR.datatype = 'qc:histo';
268 
hist2res
function hist2res(in H, in fun)
sclose
function sclose(in HDR)
gdfdatatype
function gdfdatatype(in GDFTYP)
sload
function sload(in FILENAME, in varargin)
eeg2hist
function eeg2hist(in FILENAME, in CHAN)
sopen
function sopen(in arg1, in PERMISSION, in CHAN, in MODE, in arg5, in arg6)
sread
function sread(in HDR, in NoS, in StartPos)