1function [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$ 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 53MODE=0; 54if nargin<2, CHAN=0; end; 55if ischar(CHAN), 56 MODE = 1; 57 CHAN = 0; 58end; 59 60 61if 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 64H = histo2(s) 65 66 67else 68HDR = sopen(FILENAME,'r',CHAN,'UCAL'); % open EEG file in uncalibrated mode (no scaling of the data) 69if HDR.FILE.FID<0, 70 fprintf(2,'EEG2HIST: couldnot open file %s.\n',FILENAME); 71 return; 72end; 73 74HDR.FLAG.UCAL = 1; % do not scale 75HDR.FLAG.OVERFLOWDETECTION = 0; % OFF 76if CHAN<1, CHAN=1:HDR.NS; end; 77 78H.datatype='HISTOGRAM'; 79if 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 95elseif 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 112elseif 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 125elseif (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 165elseif ~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 172else 173 fprintf(2,'EEG2HIST: format %s not implemented yet.\n',HDR.TYPE); 174end; 175HDR = sclose(HDR); 176end; 177 178HDR.HIS = H; 179 180%%% complete histogram and display it 181H.N = sumskipnan(H.H); 182%H.X = [ones(size(H.X,1),1),repmat(H.X,1,length(CHAN))]*HDR.Calib; %int16 183 184if strcmp(HDR.TYPE,'GDF') || strcmp(HDR.TYPE,'alpha') 185 HDR.THRESHOLD = [HDR.DigMin(:),HDR.DigMax(:)]; 186end; 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 225if 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; 260end; 261fprintf(1,' <FINISHED>\n'); 262 263 264%%% complete return argument 265HDR.HIS = H; 266%HDR.RES = hist2res(H); 267HDR.datatype = 'qc:histo'; 268 269