TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
sread.m
Go to the documentation of this file.
1 function [S,HDR,time] = sread(HDR,NoS,StartPos)
2 % SREAD loads selected segments of signal file
3 %
4 % [S,HDR] = sread(HDR [,NoS [,StartPos]] )
5 % NoS Number of seconds, default = 1 (second)
6 % StartPos Starting position, if not provided the following data is read continously from the file.
7 % no reposition of file pointer is performed
8 %
9 % HDR=sopen(Filename,'r',CHAN);
10 % [S,HDR] = sread(HDR, NoS, StartPos)
11 % reads NoS seconds beginning at StartPos
12 %
13 % [S,HDR] = sread(HDR, inf)
14 % reads til the end starting at the current position
15 %
16 % [S,HDR] = sread(HDR, N*HDR.Dur)
17 % reads N trials of an BKR file
18 %
19 %
20 % See also: fread, SREAD, SWRITE, SCLOSE, SSEEK, SREWIND, STELL, SEOF
21 
22 % This program is free software; you can redistribute it and/or
23 % modify it under the terms of the GNU General Public License
24 % as published by the Free Software Foundation; either version 3
25 % of the License, or (at your option) any later version.
26 
27 % $Id: sread.m 2997 2012-06-18 15:15:49Z schloegl $
28 % (C) 1997-2005,2007,2008,2011 by Alois Schloegl <alois.schloegl@gmail.com>
29 % This is part of the BIOSIG-toolbox http://biosig.sf.net/
30 
31 S = [];
32 time = [];
33 
34 if nargin<2,
35  NoS = inf;
36 end;
37 
38 if ~isnumeric(NoS) || (NoS<0),
39  fprintf(HDR.FILE.stderr,'Error SREAD: NoS must be non-negative number\n');
40  return;
41 end;
42 if (nargin>=3)
43  if (StartPos<0),
44  fprintf(HDR.FILE.stderr,'Error SREAD: StartPos must be non-negative\n');
45  return;
46  end;
47  tmp = HDR.SampleRate*StartPos;
48  if tmp ~= round(tmp),
49  % fprintf(HDR.FILE.stderr,'Warning SREAD: StartPos yields non-integer position\n');
50  StartPos = round(tmp)/HDR.SampleRate;
51  end;
52 else
53  StartPos = HDR.FILE.POS/HDR.SampleRate;
54 end;
55 
56 tmp = HDR.SampleRate*NoS;
57 if tmp ~= round(tmp),
58  fprintf(HDR.FILE.stderr,'Warning SREAD: NoS yields non-integer position [%f, %f]\n',NoS,HDR.SampleRate);
59  NoS = round(tmp)/HDR.SampleRate;
60 end;
61 
62 % define HDR.out.EVENT. This is used by EEGLAB.
63 ix = (HDR.EVENT.POS >= StartPos*HDR.SampleRate) & (HDR.EVENT.POS <= (StartPos+NoS)*HDR.SampleRate);
64 HDR.out.EVENT.POS = HDR.EVENT.POS(ix)-StartPos;
65 HDR.out.EVENT.TYP = HDR.EVENT.TYP(ix);
66 if isfield(HDR.EVENT,'CHN')
67  if ~isempty(HDR.EVENT.CHN)
68  HDR.out.EVENT.CHN = HDR.EVENT.CHN(ix);
69  end;
70 end;
71 if isfield(HDR.EVENT,'DUR')
72  if ~isempty(HDR.EVENT.DUR)
73  HDR.out.EVENT.DUR = HDR.EVENT.DUR(ix);
74  end;
75 end;
76 
77 STATUS = 0;
78 if isfield(HDR,'THRESHOLD') % save THRESHOLD status (will be modified in BKR);
79  THRESHOLD = HDR.THRESHOLD;
80 end
81 FLAG_CALIB_DONE = 0;
82 
83 if 0,
84 elseif strcmp(HDR.TYPE,'BDF'),
85  if nargin==3,
86  HDR.FILE.POS = round(HDR.SampleRate*StartPos);
87  end;
88 
89  nr = min(HDR.NRec*HDR.SPR-HDR.FILE.POS, NoS*HDR.SampleRate);
90  block1 = floor(HDR.FILE.POS/HDR.SPR);
91  ix1 = HDR.FILE.POS - block1*HDR.SPR; % starting sample (minus one) within 1st block
92  nb = ceil((HDR.FILE.POS+nr)/HDR.SPR)-block1;
93  fp = HDR.HeadLen + block1*HDR.AS.bpb;
94  status = fseek(HDR.FILE.FID, fp, 'bof');
95 
96  count = 0;
97  if HDR.NS==0,
98 
99  else
100  if (HDR.AS.spb*nb<=2^22), % faster access
101  S = repmat(NaN,HDR.SPR*nb,length(HDR.InChanSelect));
102  [s,c] = fread(HDR.FILE.FID,[3*HDR.AS.spb, nb],'uint8');
103  s = reshape(2.^[0,8,16]*reshape(s(:),3,c/3),[HDR.AS.spb, nb]);
104  c = c/3;
105  for k = 1:length(HDR.InChanSelect),
106  K = HDR.InChanSelect(k);
107  if (HDR.AS.SPR(K)>0)
108  S(:,k) = rs(reshape(s(HDR.AS.bi(K)+1:HDR.AS.bi(K+1),:),HDR.AS.SPR(K)*nb,1),HDR.AS.SPR(K),HDR.SPR);
109  else
110  S(:,k) = NaN;
111  end;
112  end;
113  S = S(ix1+1:ix1+nr,:);
114  count = nr;
115 
116  if HDR.FLAG.OVERFLOWDETECTION && isfield(HDR,'BDF') && isfield(HDR.BDF,'Status') && isfield(HDR.BDF.Status,'Channel'),
117  % BDF overflow detection is based on Status bit20
118  K = HDR.BDF.Status.Channel;
119  OVERFLOW = ~bitand(reshape(s(HDR.AS.bi(K)+1:HDR.AS.bi(K+1),:),HDR.AS.SPR(K)*nb,1),2^19);
120  OVERFLOW = rs(OVERFLOW,HDR.AS.SPR(K),HDR.SPR);
121  OVERFLOW = OVERFLOW(ix1+1:ix1+nr,:);
122  S(OVERFLOW>0,:)=NaN;
123  end;
124  else
125  S = repmat(NaN,nr,length(HDR.InChanSelect));
126  while (count<nr);
127  len = ceil(min([(nr-count)/HDR.SPR,2^22/HDR.AS.spb]));
128  [s,c] = fread(HDR.FILE.FID,[3*HDR.AS.spb, len],'uint8=>uint8');
129  s1 = zeros(HDR.SPR*c/(3*HDR.AS.spb),length(HDR.InChanSelect));
130  for k = 1:length(HDR.InChanSelect),
131  K = HDR.InChanSelect(k);
132  tmp = 2.^[0,8,16]*double(reshape(s(HDR.AS.bi(K)*3+1:HDR.AS.bi(K+1)*3,:),3,HDR.AS.SPR(K)*c/HDR.AS.bpb));
133  if (HDR.AS.SPR(K)>0)
134  s1(:,k) = rs(tmp',HDR.AS.SPR(K),HDR.SPR);
135  else
136  s1(:,k) = NaN;
137  end;
138  end;
139  if HDR.FLAG.OVERFLOWDETECTION && isfield(HDR,'BDF') && isfield(HDR.BDF,'Status') && isfield(HDR.BDF.Status,'Channel'),
140  % BDF overflow detection is based on Status bit20
141  K = HDR.BDF.Status.Channel;
142  tmp = 2.^[0,8,16]*double(reshape(s(HDR.AS.bi(K)*3+1:HDR.AS.bi(K+1)*3,:),3,HDR.AS.SPR(K)*c/HDR.AS.bpb));
143  OVERFLOW = rs(~bitand(tmp',2^19),HDR.AS.SPR(K),HDR.SPR);
144  s1(OVERFLOW>0,:)=NaN;
145  end;
146  ix2 = min(nr-count, size(s1,1)-ix1);
147  S(count+1:count+ix2,:) = s1(ix1+1:ix1+ix2,:);
148  count = count+ix2;
149  ix1 = 0;
150  end;
151  end;
152  HDR.FILE.POS = HDR.FILE.POS + count;
153  S = S - 2^24*(S>=2^23);
154  end
155 
156 elseif strcmp(HDR.TYPE,'EDF') || strcmp(HDR.TYPE,'GDF') || strcmp(HDR.TYPE,'BDF') || strcmp(HDR.TYPE,'ACQ'),
157  % experimental, might replace SDFREAD.M
158  if nargin==3,
159  HDR.FILE.POS = round(HDR.SampleRate*StartPos);
160  end;
161 
162  nr = min(HDR.NRec*HDR.SPR-HDR.FILE.POS, NoS*HDR.SampleRate);
163  S = repmat(NaN,nr,length(HDR.InChanSelect));
164 
165  block1 = floor(HDR.FILE.POS/HDR.SPR);
166  ix1 = HDR.FILE.POS- block1*HDR.SPR; % starting sample (minus one) within 1st block
167  nb = ceil((HDR.FILE.POS+nr)/HDR.SPR)-block1;
168  fp = HDR.HeadLen + block1*HDR.AS.bpb;
169  STATUS = fseek(HDR.FILE.FID, fp, 'bof');
170  count = 0;
171  if HDR.NS==0,
172  elseif all(HDR.GDFTYP==HDR.GDFTYP(1)),
173  if (HDR.AS.spb*nb<=2^24), % faster access
174  S = [];
175  %[HDR.AS.spb, nb,block1,ix1,HDR.FILE.POS],
176  [s,c] = fread(HDR.FILE.FID,[HDR.AS.spb, nb],gdfdatatype(HDR.GDFTYP(1)));
177  for k = 1:length(HDR.InChanSelect),
178  K = HDR.InChanSelect(k);
179  if (HDR.AS.SPR(K)>0)
180  S(:,k) = rs(reshape(s(HDR.AS.bi(K)+1:HDR.AS.bi(K+1),:),HDR.AS.SPR(K)*nb,1),HDR.AS.SPR(K),HDR.SPR);
181  else
182  S(:,k) = NaN;
183  end;
184  end;
185  S = S(ix1+1:ix1+nr,:);
186  count = nr;
187  else
188  S = repmat(NaN,[nr,length(HDR.InChanSelect)]);
189  while (count<nr);
190  len = ceil(min([(nr-count)/HDR.SPR,2^22/HDR.AS.spb]));
191  [s,c] = fread(HDR.FILE.FID,[HDR.AS.spb, len],gdfdatatype(HDR.GDFTYP(1)));
192  s1 = zeros(HDR.SPR*c/HDR.AS.spb,length(HDR.InChanSelect));
193  for k = 1:length(HDR.InChanSelect),
194  K = HDR.InChanSelect(k);
195  if (HDR.AS.SPR(K)>0)
196  tmp = reshape(s(HDR.AS.bi(K)+1:HDR.AS.bi(K+1),:),HDR.AS.SPR(K)*c/HDR.AS.spb,1);
197  s1(:,k) = rs(tmp,HDR.AS.SPR(K),HDR.SPR);
198  end;
199  end;
200  ix2 = min(nr-count, size(s1,1)-ix1);
201  S(count+1:count+ix2,:) = s1(ix1+1:ix1+ix2,:);
202  count = count+ix2;
203  ix1 = 0;
204  end;
205  end;
206  else
207  fprintf(2,'SREAD (GDF): different datatypes - this might take some time.\n');
208 
209  S = repmat(NaN,[nr,length(HDR.InChanSelect)]);
210  while (count<nr);
211  s = [];
212  for k=1:length(HDR.AS.TYP),
213  [s0,tmp] = fread(HDR.FILE.FID,[HDR.AS.c(k), 1],gdfdatatype(HDR.AS.TYP(k)));
214  s = [s;s0];
215  end;
216 
217  s1 = repmat(NaN,[HDR.SPR,length(HDR.InChanSelect)]);
218  for k = 1:length(HDR.InChanSelect),
219  K = HDR.InChanSelect(k);
220  if (HDR.AS.SPR(K)>0)
221  s1(:,k) = rs(s(HDR.AS.bi(K)+1:HDR.AS.bi(K+1),:),HDR.AS.SPR(K),HDR.SPR);
222  end;
223  end;
224  ix2 = min(nr-count, size(s1,1)-ix1);
225  S(count+1:count+ix2,:) = s1(ix1+1:ix1+ix2,:);
226  count = count+HDR.SPR;
227  ix1 = 0;
228  end;
229  end;
230  if strcmp(HDR.TYPE,'GDF') % read non-equidistant sampling channels of GDF2.0 format
231  if (HDR.VERSION>1.94) %& isfield(HDR.EVENT,'VAL'),
232  for k = 1:length(HDR.InChanSelect),
233  ch = HDR.InChanSelect(k);
234  if (HDR.AS.SPR(ch)==0),
235  ix = find((HDR.EVENT.TYP==hex2dec('7fff')) & (HDR.EVENT.CHN==ch));
236  pix= HDR.EVENT.POS(ix)-HDR.FILE.POS;
237  ix1= find((pix > 0) & (pix <= count));
238  S(pix(ix1),k)=HDR.EVENT.DUR(ix(ix1));
239  end;
240  end;
241  end;
242  end
243  HDR.FILE.POS = HDR.FILE.POS + count;
244 
245 
246 elseif strcmp(HDR.TYPE,'AINF'),
247  if nargin==3,
248  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.SampleRate*StartPos*2*(HDR.NS+2),'bof');
249  HDR.FILE.POS = HDR.SampleRate*StartPos;
250  end;
251 
252  %[S,count] = fread(HDR.FILE.FID,[HDR.NS+2,HDR.SampleRate*NoS],'int16');
253  nr = min(HDR.SampleRate*NoS, HDR.SPR*HDR.NRec-HDR.FILE.POS);
254  S = [];
255  time = [];
256  count = 0;
257  while (count<nr),
258  [s,c] = fread(HDR.FILE.FID, [HDR.NS+2, min(nr-count,floor(2^24/HDR.NS))], 'int16');
259  if nargout>2,
260  time = [time; [s(1:2,:)'+2^16*(s(1:2,:)'<0)]*(2.^[16;0])];
261  end;
262  S = [S; s(2+HDR.InChanSelect,:)'];
263  count = count + c/(HDR.NS+2);
264  end;
265  HDR.FILE.POS = HDR.FILE.POS + count;
266 
267 
268 elseif strmatch(HDR.TYPE,{'BKR'}),
269  if nargin==3,
270  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.SampleRate*HDR.NS*StartPos*2,'bof');
271  HDR.FILE.POS = HDR.SampleRate*StartPos;
272  end;
273  [S,count] = fread(HDR.FILE.FID,[HDR.NS,HDR.SampleRate*NoS],'int16');
274  if ~isempty(S),
275  S = S(HDR.InChanSelect,:)';
276  HDR.FILE.POS = HDR.FILE.POS + count/HDR.NS;
277  else
278  S = zeros(0,length(HDR.InChanSelect));
279  end;
280  THRESHOLD(HDR.AS.TRIGCHAN,:)=NaN; % do not apply overflow detection for Trigger channel
281 
282 
283 elseif strmatch(HDR.TYPE,{'AIF','CNT','SND','WAV','Sigma'})
284  if nargin==3,
285  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.SampleRate*HDR.AS.bpb*StartPos,'bof');
286  HDR.FILE.POS = HDR.SampleRate*StartPos;
287  end;
288 
289  maxsamples = min(HDR.SPR*HDR.NRec - HDR.FILE.POS, HDR.SampleRate*NoS);
290  [S,count] = fread(HDR.FILE.FID,[HDR.NS,maxsamples], gdfdatatype(HDR.GDFTYP));
291 
292  S = S(HDR.InChanSelect,:)';
293  HDR.FILE.POS = HDR.FILE.POS + count/HDR.NS;
294 
295  if ~HDR.FLAG.UCAL,
296  if isfield(HDR.FILE,'TYPE')
297  if HDR.FILE.TYPE==1,
298  S = mu2lin(S);
299  end;
300  end;
301  end;
302 
303 
304 elseif strmatch(HDR.TYPE,{'BLSC2','CFWB','CNT','DEMG','DDT','ET-MEG','ISHNE','Nicolet','RG64'}),
305 
306  tc = strcmp(HDR.TYPE,'CFWB') && isfield(HDR,'FLAG') && isfield(HDR.FLAG,'TimeChannel') && HDR.FLAG.TimeChannel;
307 
308  if nargin==3,
309  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.AS.bpb*round(HDR.SampleRate*StartPos),'bof');
310  HDR.FILE.POS = HDR.SampleRate*StartPos;
311  end;
312  if strcmpi(HDR.FLAG.OUTPUT,'single'),
313  DT = [gdfdatatype(HDR.GDFTYP),'=>single'];
314  elseif any(HDR.GDFTYP==[1:6,16]) && ~exist('OCTAVE_VERSION','builtin'),
315  % preserve data type
316  DT = ['*',gdfdatatype(HDR.GDFTYP)];
317  else
318  % convert to double
319  DT = [gdfdatatype(HDR.GDFTYP)];
320  end;
321  maxsamples = min(HDR.SampleRate*NoS, HDR.NRec*HDR.SPR-HDR.FILE.POS);
322  S = []; count = 0;
323  while maxsamples>0,
324  % the maximum block size of 2^23 is a heuristical value
325  [s,c] = fread(HDR.FILE.FID, [HDR.NS+tc,maxsamples], DT);
326  c = c/(HDR.NS+tc);
327  count = count + c;
328  maxsamples = maxsamples - c;
329  if c>0,
330  S = [S; s(HDR.InChanSelect+tc,:)'];
331  else
332  fprintf(HDR.FILE.stderr,'Warning SREAD(%s): could not read %i samples, only %i samples read\n',HDR.TYPE,maxsamples+count,count);
333  break;
334  end;
335  end;
336  HDR.FILE.POS = HDR.FILE.POS + count;
337 
338 
339 elseif strcmp(HDR.TYPE,'EPL'),
340  if nargin==3,
341  HDR.FILE.POS = HDR.SampleRate*StartPos;
342  end;
343  startblock = floor(HDR.FILE.POS/HDR.SPR);
344  STATUS = fseek(HDR.FILE.FID, HDR.HeadLen + startblock*HDR.AS.bpb, 'bof'); % fseek needed because HDR.FILE.POS can be changed by SSEEK
345  curblock = startblock;
346  endpos = min(HDR.FILE.POS+NoS*HDR.SampleRate, HDR.NRec*HDR.SPR);
347 
348  [datablock,count] = fread(HDR.FILE.FID, [(1+HDR.NS)*HDR.SPR,ceil(endpos/HDR.SPR)-startblock], 'int16');
349  datablock = reshape(datablock(256+1:end,:),HDR.NS,HDR.SPR*size(datablock,2))'; % remove mark track, and reshape data
350 
351  S = datablock(HDR.FILE.POS-startblock*HDR.SPR+1:endpos-startblock*HDR.SPR,HDR.InChanSelect);
352  HDR.FILE.POS = HDR.FILE.POS + size(S,1);
353 
354 
355 elseif strcmp(HDR.TYPE,'SMA'),
356  if nargin==3,
357  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.SampleRate*HDR.AS.bpb*StartPos,'bof');
358  HDR.FILE.POS = HDR.SampleRate*StartPos;
359  end;
360  tmp = min(NoS*HDR.SampleRate,(HDR.NRec*HDR.SPR-HDR.FILE.POS));
361  [S,count] = fread(HDR.FILE.FID,[HDR.NS,tmp],'float32'); % read data frame
362  tmp = HDR.NS*tmp;
363  if count < tmp,
364  fprintf(HDR.FILE.stderr,'Warning SREAD SMA: only %i out of %i samples read\n',count/HDR.NS,tmp/HDR.NS);
365  end;
366  S = S(HDR.InChanSelect,:)';
367  HDR.FILE.POS = HDR.FILE.POS + count/HDR.NS;
368 
369  HDR.SMA.events = diff(sign([HDR.Filter.T0',S(HDR.SMA.EVENT_CHANNEL,:)]-HDR.SMA.EVENT_THRESH))>0;
370  HDR.EVENT.POS = find(HDR.SMA.events);
371  HDR.EVENT.TYP = HDR.SMA.events(HDR.EVENT.POS);
372 
373  if size(S,2) > 0,
374  HDR.Filter.T0 = S(HDR.SMA.EVENT_CHANNEL,size(S,2))';
375  end;
376 
377 
378 elseif strcmp(HDR.TYPE,'RDF'),
379  S = [];
380  if nargin>2,
381  HDR.FILE.POS = StartPos;
382  end;
383  POS = HDR.FILE.POS;
384 
385  NoSeg = min(NoS,length(HDR.Block.Pos)-HDR.FILE.POS);
386  count = 0;
387  S = zeros(NoSeg*HDR.SPR, length(HDR.InChanSelect));
388 
389  for k = 1:NoSeg,
390  STATUS = fseek(HDR.FILE.FID,HDR.Block.Pos(POS+k),-1);
391 
392  % Read nchans and block length
393  tmp = fread(HDR.FILE.FID,34+220,'uint16');
394 
395  %STATUS = fseek(HDR.FILE.FID,2,0);
396  nchans = tmp(2); %fread(HDR.FILE.FID,1,'uint16');
397  %fread(HDR.FILE.FID,1,'uint16');
398  block_size = tmp(4); %fread(HDR.FILE.FID,1,'uint16');
399  %ndupsamp = fread(HDR.FILE.FID,1,'uint16');
400  %nrun = fread(HDR.FILE.FID,1,'uint16');
401  %err_detect = fread(HDR.FILE.FID,1,'uint16');
402  %nlost = fread(HDR.FILE.FID,1,'uint16');
403  nevents = tmp(9); %fread(HDR.FILE.FID,1,'uint16');
404  %STATUS = fseek(HDR.FILE.FID,50,0);
405 
406  [data,c] = fread(HDR.FILE.FID,[nchans,block_size],'int16');
407  %S = [S; data(HDR.InChanSelect,:)']; % concatenate data blocks
408  S((k-1)*HDR.SPR+(1:c/nchans),:) = data(HDR.InChanSelect,:)';
409  count = count + c;
410  end;
411  HDR.FILE.POS = HDR.FILE.POS + NoSeg;
412 
413 
414 elseif strcmp(HDR.TYPE,'LABVIEW'),
415  if nargin==3,
416  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.SampleRate*HDR.AS.bpb*StartPos,'bof');
417  HDR.FILE.POS = HDR.SampleRate*StartPos;
418  end;
419  [S,count] = fread(HDR.FILE.FID,[HDR.NS,HDR.SampleRate*NoS],'int32');
420  if count,
421  S = S(HDR.InChanSelect,:)';
422  HDR.FILE.POS = HDR.FILE.POS + count/HDR.NS;
423  end;
424 
425 
426 elseif strcmp(HDR.TYPE,'alpha'),
427  if nargin==3,
428  POS = HDR.SampleRate*StartPos*HDR.AS.bpb/HDR.SPR;
429  if POS~=ceil(POS),
430  fprintf(HDR.FILE.stderr,'warning SREAD (alpha): starting position is non-integer (%f)\n',POS);
431  end
432  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen + POS,'bof');
433  HDR.FILE.POS = HDR.SampleRate*StartPos;
434  end;
435 
436  nr = min(HDR.SampleRate*NoS, HDR.NRec*HDR.SPR-HDR.FILE.POS)*HDR.AS.bpb;
437  if (nr - round(nr)) < .01,
438  nr = round(nr);
439  else
440  fprintf(HDR.FILE.stderr,'Error SREAD (alpha): can not deal with odd number of samples \n');
441  return;
442  end
443 
444  if HDR.Bits==12,
445  [s,count] = fread(HDR.FILE.FID,[3,nr/3],'uint8');
446  s(1,:) = s(1,:)*16 + floor(s(2,:)/16);
447  s(3,:) = s(3,:)+ mod(s(2,:),16)*256;
448  s = reshape(s([1,3],:),2*size(s,2),1);
449  s = s - (s>=2^11)*2^12;
450  nr = floor(length(s)/HDR.NS);
451  S = reshape(s(1:nr*HDR.NS),HDR.NS,nr);
452  count = count*2/3;
453 
454  elseif HDR.Bits==16,
455  [S,count] = fread(HDR.FILE.FID,[HDR.NS,nr],'int16');
456 
457  elseif HDR.Bits==32,
458  [S,count] = fread(HDR.FILE.FID,[HDR.NS,nr],'int32');
459  end;
460 
461  if count,
462  S = S(HDR.InChanSelect,:)';
463  HDR.FILE.POS = HDR.FILE.POS + count/HDR.NS;
464  end;
465 
466 
467 elseif strcmp(HDR.TYPE,'MIT'),
468  if nargin==3,
469  STATUS = fseek(HDR.FILE.FID,HDR.SampleRate*HDR.AS.bpb*StartPos,'bof');
470  tmp = HDR.SampleRate*StartPos;
471  if HDR.FILE.POS~=tmp,
472  HDR.mode8.accu = zeros(1,length(HDR.InChanSelect));
473  HDR.mode8.valid= 0;
474  end;
475  end;
476 
477  DataLen = NoS*HDR.SampleRate/HDR.SPR;
478  if HDR.VERSION == 212,
479  [A,count] = fread(HDR.FILE.FID, [1,DataLen*HDR.AS.bpb], 'uint8'); % matrix with 3 rows, each 8 bits long, = 2*12bit
480  DataLen = floor(count/HDR.AS.bpb);
481  if (count~= DataLen*HDR.AS.bpb) && isfinite(DataLen),
482  fprintf(HDR.FILE.stderr,'Warning SREAD (MIT): non-integer block length %i,%f\n',count, DataLen*HDR.AS.bpb);
483  %HDR = sseek(HDR,HDR.FILE.POS,'bof');
484  %return;
485  A = A(1:DataLen*HDR.AS.bpb);
486  end;
487  S = [A(1:3:end) + mod(A(2:3:end),16)*256; A(3:3:end) + floor(A(2:3:end)/16)*256];
488  clear A;
489  S = S - 2^12*(S>=2^11); % 2-th complement
490  S = reshape(S,HDR.AS.spb,prod(size(S))/HDR.AS.spb)';
491 
492  elseif HDR.VERSION == 310,
493  [A,count] = fread(HDR.FILE.FID, [HDR.AS.bpb/2, DataLen], 'uint16');
494  A = A'; DataLen = count/HDR.AS.bpb*2;
495  for k = 1:ceil(HDR.AS.spb/3),
496  k1=3*k-2; k2=3*k-1; k3=3*k;
497  S(:,3*k-2) = floor(mod(A(:,k*2-1),2^12)/2);
498  S(:,3*k-1) = floor(mod(A(:,k*2),2^12)/2);
499  S(:,3*k ) = floor(A(:,k*2-1)*(2^-11)) + floor(A(:,k*2)*(2^-11))*2^5;
500  S = mod(S(:,1:HDR.AS.spb),2^10);
501  S = S - 2^10*(S>=2^9); % 2-th complement
502  end;
503  S = S;
504 
505  elseif HDR.VERSION == 311,
506  [A,count] = fread(HDR.FILE.FID, [HDR.AS.bpb/4, DataLen], 'uint32');
507  A = A'; DataLen = count/HDR.AS.bpb*4;
508  for k = 1:ceil(HDR.AS.spb/3),
509  S(:,3*k-2) = mod(A(:,k),2^10);
510  S(:,3*k-1) = mod(floor(A(:,k)*2^(-11)),2^10);
511  S(:,3*k) = mod(floor(A(:,k)*2^(-22)),2^10);
512  S = S(:,1:HDR.AS.spb);
513  S = S - 2^10*(S>=2^9); % 2-th complement
514  end;
515  S = S';
516 
517  elseif HDR.VERSION == 8,
518  [S,count] = fread(HDR.FILE.FID, [HDR.AS.spb,DataLen], 'int8');
519  S = S'; DataLen = count/HDR.AS.spb
520 
521  elseif HDR.VERSION == 80,
522  [S,count] = fread(HDR.FILE.FID, [HDR.AS.spb,DataLen], 'uint8');
523  S = S'-128; DataLen = count/HDR.AS.spb;
524 
525  elseif HDR.VERSION == 160,
526  [S,count] = fread(HDR.FILE.FID, [HDR.AS.spb,DataLen], 'uint16');
527  S = S'-2^15; DataLen = count/HDR.AS.spb;
528 
529  elseif HDR.VERSION == 16,
530  [S,count] = fread(HDR.FILE.FID, [HDR.AS.spb,DataLen], 'int16');
531  S = S'; DataLen = count/HDR.AS.spb;
532 
533  elseif HDR.VERSION == 61,
534  [S,count] = fread(HDR.FILE.FID, [HDR.AS.spb,DataLen], 'int16');
535  S = S'; DataLen = count/HDR.AS.spb;
536 
537  else
538  fprintf(2, 'ERROR MIT-ECG: format %i not supported.\n',HDR.VERSION);
539 
540  end;
541  if any(HDR.AS.SPR>1),
542  A = S;
543  S = zeros(size(A,1)*HDR.SPR,length(HDR.InChanSelect));
544  for k = 1:length(HDR.InChanSelect),
545  ch = HDR.InChanSelect(k);
546  ix = HDR.AS.bi(ch)+1:HDR.AS.bi(ch+1);
547  S(:,k)=rs(reshape(A(:,ix)',size(A,1)*HDR.AS.SPR(ch),1),HDR.AS.SPR(ch),HDR.SPR);
548  end
549  else
550  S = S(:,HDR.InChanSelect);
551  end
552  if HDR.VERSION == 8,
553  if HDR.FILE.POS==0,
554  HDR.mode8.accu = zeros(1,length(HDR.InChanSelect));
555  HDR.mode8.valid= 1;
556  end;
557  if ~HDR.mode8.valid;
558  fprintf(2,'Warning SREAD: unknown offset (TYPE=MIT, mode=8) \n');
559  else
560  S(1,:) = S(1,:) + HDR.mode8.accu;
561  end;
562  S = cumsum(S);
563  HDR.mode8.accu = S(size(S,1),:);
564  end;
565  HDR.FILE.POS = HDR.FILE.POS + DataLen;
566 
567 
568 elseif strcmp(HDR.TYPE,'TMS32'),
569  if nargin==3,
570  HDR.FILE.POS = round(HDR.SampleRate*StartPos);
571  end;
572 
573  blockN = floor(HDR.FILE.POS/HDR.SPR);
574  ix1 = HDR.FILE.POS- blockN*HDR.SPR; % starting sample (minus one) within 1st block
575  fp = HDR.HeadLen + blockN*HDR.AS.bpb;
576  status = fseek(HDR.FILE.FID, fp, 'bof');
577 
578  nr = min(HDR.AS.endpos-HDR.FILE.POS, NoS*HDR.SampleRate);
579  S = repmat(NaN,nr,length(HDR.InChanSelect));
580  count = 0;
581 
582  fread(HDR.FILE.FID,86,'uint8');
583  while (count<nr) && ~feof(HDR.FILE.FID),
584  if all(HDR.GDFTYP==HDR.GDFTYP(1))
585  [s,c] = fread(HDR.FILE.FID,[HDR.NS,HDR.SPR],gdfdatatype(HDR.GDFTYP(1)));
586  else
587  s = repmat(NaN,HDR.NS,HDR.SPR);
588  c = 0;
589  for k1 = 1:HDR.SPR,
590  for k2 = 1:HDR.NS,
591  [s(k2,k1),c2] = fread(HDR.FILE.FID,1,gdfdatatype(HDR.GDFTYP(k2)));
592  end;
593  end;
594  end;
595  ix2 = min(nr-count, size(s,2)-ix1);
596  S(count+1:count+ix2,:) = s(HDR.InChanSelect, ix1+1:ix1+ix2)';
597  count = count + ix2;
598  ix1 = 0; % reset starting index,
599  fread(HDR.FILE.FID,86,'uint8');
600  end;
601  HDR.FILE.POS = HDR.FILE.POS + count;
602 
603 
604 elseif strcmp(HDR.TYPE,'EGI'),
605  if nargin==3,
606  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.AS.bpb*StartPos,'bof');
607  HDR.FILE.POS = HDR.SampleRate*StartPos;
608  end;
609 
610  if HDR.FLAG.TRIGGERED,
611  NoS = min(NoS,(HDR.NRec-HDR.FILE.POS));
612  S = zeros(NoS*HDR.SPR,length(HDR.InChanSelect))+NaN;
613  for i = (1:NoS),
614  SegmentCatIndex(HDR.FILE.POS+i) = fread(HDR.FILE.FID,1,'uint16');
615  SegmentStartTime(HDR.FILE.POS+i) = fread(HDR.FILE.FID,1,'uint32');
616 
617  [s,count] = fread(HDR.FILE.FID, [HDR.NS + HDR.EGI.N, HDR.SPR*HDR.NRec], gdfdatatype(HDR.GDFTYP));
618  tmp = (HDR.NS + HDR.EGI.N) * HDR.SPR;
619  if isfinite(tmp) && (count < tmp),
620  fprintf(HDR.FILE.stderr,'Warning SREAD EGI: only %i out of %i samples read\n',count,tmp);
621  end;
622  HDR.FILE.POS = HDR.FILE.POS + count/tmp;
623 
624  if (HDR.EGI.N > 0),
625  [HDR.EVENT.POS,HDR.EVENT.CHN,HDR.EVENT.TYP] = find(s(HDR.NS+1:size(s,1),:)');
626  HDR.EVENT.DUR = ones(size(HDR.EVENT.POS));
627  end
628  S((i-1)*HDR.SPR + (1:size(s,2)),:) = s(HDR.InChanSelect,:)';
629  end;
630  else
631  [S,count] = fread(HDR.FILE.FID,[HDR.NS + HDR.EGI.N, HDR.SampleRate*NoS],gdfdatatype(HDR.GDFTYP));
632  tmp = HDR.SampleRate * NoS;
633  if isfinite(tmp) && (count < tmp),
634  fprintf(HDR.FILE.stderr,'Warning SREAD EGI: only %i out of %i samples read\n',count,tmp);
635  end;
636  HDR.FILE.POS = HDR.FILE.POS + round(count/(HDR.NS + HDR.EGI.N));
637  S = S(HDR.InChanSelect,:)';
638  end;
639 
640 
641 elseif strcmp(HDR.TYPE,'AVG'),
642  S = repmat(nan,HDR.SPR,HDR.NS);
643  count = 0;
644  for i = 1:HDR.NS,
645  [tmp,c] = fread(HDR.FILE.FID,5,'uint8'); % no longer used
646  count = count + c;
647  [S(:,i), c] = fread(HDR.FILE.FID,HDR.SPR,'float');
648  count = count + c*4;
649  end
650  S = S(:,HDR.InChanSelect);
651  HDR.FILE.POS = HDR.FILE.POS + count/HDR.AS.bpb;
652 
653 
654 elseif strcmp(HDR.TYPE,'COH'),
655  warning('.COH data not tested yet')
656  if (prod(size(NoS))==1) && (nargin>2),
657  rows = NoS; cols = StartPos;
658  elseif prod(size(NoS))==2
659  rows = NoS(1); cols = NoS(2);
660  else
661  fprintf(HDR.FILE.stderr,'Error SREAD mode=COH: invalid arguments.\n');
662  end;
663 
664  STATUS = fseek(HDR.FILE.FID,HDR.COH.directory(rows,cols)+8,'bof'); % skip over a small unused header of 8 bytes
665  sr = fread(HDR.FILE.FID, HDR.SPR, 'float32'); % read real part of coherence
666  si = fread(HDR.FILE.FID, HDR.SPR, 'float32'); % read imag part of coherence
667  S = sr + i * si;
668 
669 
670 elseif strcmp(HDR.TYPE,'CSA'),
671  warning('.CSA data not tested yet')
672  S = fread(HDR.FILE.FID, [HDR.NRec*(HDR.SPR+6)*HDR.NS], 'float32');
673 
674 
675 elseif strcmp(HDR.TYPE,'EEG'),
676  if nargin>2,
677  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.AS.bpb*StartPos,'bof');
678  end;
679 
680  NoS = min(NoS, HDR.NRec-HDR.FILE.POS);
681  S = zeros(NoS*HDR.SPR, length(HDR.InChanSelect));
682  count = 0;
683  for i = 1:NoS, %h.compsweeps,
684  h.sweep(i).accept = fread(HDR.FILE.FID,1,'uchar');
685  tmp = fread(HDR.FILE.FID,2,'ushort');
686  h.sweep(i).ttype = tmp(1);
687  h.sweep(i).correct = tmp(2);
688  h.sweep(i).rt = fread(HDR.FILE.FID,1,'float32');
689  tmp = fread(HDR.FILE.FID,2,'ushort');
690  h.sweep(i).response = tmp(1);
691  h.sweep(i).reserved = tmp(2);
692 
693  [signal,c] = fread(HDR.FILE.FID, [HDR.NS,HDR.SPR], gdfdatatype(HDR.GDFTYP));
694 
695  S(i*HDR.SPR+(1-HDR.SPR:0),:) = signal(HDR.InChanSelect,:)';
696  count = count + c;
697  end;
698  HDR.FILE.POS = HDR.FILE.POS + count/HDR.AS.spb;
699 
700 
701 elseif strcmp(HDR.TYPE,'MFER'),
702  if (HDR.FRAME.N ~= 1),
703  fprintf(2,'Warning MWFOPEN: files with more than one frame not implemented, yet.\n');
704  return;
705  end
706 
707  N = 1;
708  if ~isfield(HDR,'data'),
709  STATUS = fseek(HDR.FILE.FID,HDR.FRAME.POS(N),'bof');
710  [tmp,count] = fread(HDR.FILE.FID,HDR.FRAME.sz(N,1:2),gdfdatatype(HDR.FRAME.TYP(N)));
711  if isnan(HDR.NRec),
712  HDR.NRec = count/(HDR.SPR*HDR.NS);
713  end;
714 
715  if count==(HDR.SPR*HDR.NS), %% alternate mode format
716  tmp = reshape(tmp,[HDR.SPR,HDR.NS]);
717  else
718  tmp = reshape(tmp,[HDR.SPR,HDR.NS,HDR.NRec]); % convert into 3-Dim
719  tmp = permute(tmp,[1,3,2]); % re-order dimensions
720  tmp = reshape(tmp,[HDR.SPR*HDR.NRec,HDR.NS]); % make 2-Dim
721  end;
722  HDR.data = tmp;
723  end;
724 
725  if nargin>2,
726  HDR.FILE.POS = HDR.SampleRate*StartPos;
727  end;
728 
729  nr = min(HDR.SampleRate*NoS,size(HDR.data,1)-HDR.FILE.POS);
730  S = HDR.data(HDR.FILE.POS + (1:nr), HDR.InChanSelect);
731  HDR.FILE.POS = HDR.FILE.POS + nr;
732 
733 
734 elseif strcmp(HDR.TYPE,'BCI2000'),
735  if nargin==3,
736  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.SampleRate*HDR.AS.bpb*StartPos,'bof');
737  HDR.FILE.POS = HDR.SampleRate*StartPos;
738  end;
739  [S,count] = fread(HDR.FILE.FID,[HDR.NS,HDR.SampleRate*NoS],HDR.BCI2000.GDFTYP,HDR.BCI2000.StateVectorLength);
740  if count,
741  S = S(HDR.InChanSelect,:)';
742  HDR.FILE.POS = HDR.FILE.POS + count/HDR.NS;
743  end;
744 
745 
746 elseif strcmp(HDR.TYPE,'native') || strcmp(HDR.TYPE,'SCP'),
747  if nargin>2,
748  HDR.FILE.POS = HDR.SampleRate*StartPos;
749  end;
750 
751  nr = min(round(HDR.SampleRate * NoS), size(HDR.data,1) - HDR.FILE.POS);
752  S = HDR.data(HDR.FILE.POS+1:HDR.FILE.POS+nr,:);
753  HDR.FILE.POS = HDR.FILE.POS + nr;
754 % FLAG_CALIB_DONE = 1;
755 
756 elseif strcmp(HDR.TYPE,'NEX'),
757  %% hack: read NEX data once, and transform into "native"
758  for k = 1:HDR.NEX.NS,
759  fseek(HDR.FILE.FID, HDR.NEX.offset(k), 'bof');
760  if 0, % elseif HDR.NEX.type(k)==1,
761 
762  elseif HDR.NEX.type(k)==2, % interval
763  tmp = fread(HDR.FILE.FID, [HDR.NEX.nf(k),2], 'int32');
764  HDR.EVENT.POS = [HDR.EVENT.POS; tmp(:,1)];
765  HDR.EVENT.DUR = [HDR.EVENT.DUR; tmp(:,2)];
766  HDR.EVENT.CHN = [HDR.EVENT.CHN; repmat(k,size(tmp,1),1)];
767 
768  elseif HDR.NEX.type(k)==3, % waveform
769  HDR.HeadLen = HDR.NEX.offset(k);
770  HDR.NEX6.ts{k} = fread(HDR.FILE.FID, [HDR.NEX.nf(k),1], 'int32');
771  HDR.NEX6.data{k} = fread(HDR.FILE.FID, [HDR.NEX.SPR(k), HDR.NEX.nf(k)], 'int16');
772 
773  elseif HDR.NEX.type(k)==5, % continous variable
774  HDR.NEX5.ts{k} = fread(HDR.FILE.FID, [HDR.NEX.nf(k), 2], 'int32')';
775  HDR.data(:,HDR.AS.chanreduce(k)) = fread(HDR.FILE.FID, [HDR.NEX.SPR(k), 1], 'int16');
776 
777  elseif HDR.NEX.type(k)==6, % marker
778  ts = fread(HDR.FILE.FID, [1,HDR.NEX.nf(k)], 'int32');
779  names = zeros(1,64);
780  m = zeros(HDR.NEX.SPR(k), nl, nm);
781  for j=1:nm
782  names(j, :) = fread(HDR.FILE.FID, [1 64], 'uint8');
783  for p = 1:HDR.NEX.SPR(k)
784  m(p, :, j) = fread(HDR.FILE.FID, [1 nl], 'uint8');
785  end
786  end
787  HDR.NEX.names = names;
788  HDR.NEX.m = m;
789  else
790  ts = fread(HDR.FILE.FID, [1,HDR.NEX.nf(k)], 'int32');
791  HDR.NEX0.ts{k}=ts;
792  end;
793  end
794  fclose(HDR.FILE.FID);
795  HDR.FILE.OPEN = 0;
796  HDR.TYPE = 'native';
797  HDR.data = HDR.data(:,HDR.InChanSelect);
798 
799  % sequence for reading "native" format
800  if nargin>2,
801  HDR.FILE.POS = HDR.SampleRate*StartPos;
802  end;
803 
804  nr = min([round(HDR.SampleRate * NoS), size(HDR.data,1) - HDR.FILE.POS]);
805  S = HDR.data(HDR.FILE.POS + (1:nr), :);
806  HDR.FILE.POS = HDR.FILE.POS + nr;
807 
808 
809 elseif strcmp(HDR.TYPE,'PLEXON'),
810  %% hack: read PLEXON data once, and transform into "native"
811  HDR.data = repmat(NaN,max(HDR.PLX.adcount),HDR.NS);
812  NRec = 0;
813  nET = 0;
814  ET = repmat(NaN,1024,6);
815  wav = zeros(1024,HDR.PLX.wavlen);
816 
817  tscount=zeros(5,130);
818  wfcount=zeros(5,130);
819  evcount=zeros(1,300);
820  adcount=zeros(1,212);
821 
822  adpos = zeros(1,HDR.NS);
823  typ_ubyte = fread(HDR.FILE.FID,2,'int16');
824  while ~feof(HDR.FILE.FID),
825  TYP = typ_ubyte(1);
826  ubyte= typ_ubyte(2);
827  POS = fread(HDR.FILE.FID,1,'int32');
828  tmp = fread(HDR.FILE.FID,4,'int16');
829  CHN = tmp(1)+1;
830  unit = tmp(2);
831  nwf = tmp(3);
832  DUR = tmp(4);
833  if nwf>0,
834  wf = fread(HDR.FILE.FID,[1,DUR],'int16');
835  end;
836  if TYP<5,
837  nET = nET + 1;
838  if size(ET,1) < nET; % memory allocation
839  ET = [ET ; repmat(NaN,size(ET))];
840  end;
841  ET(nET,:) = [TYP,ubyte*2^32+POS,CHN,DUR,unit,nwf];
842  end;
843  if TYP==1, % spike
844  tscount(unit+1,CHN) = tscount(unit+1,CHN)+1;
845  if DUR>0,
846  wfcount(unit+1,CHN) = wfcount(unit+1,CHN)+1;
847  end;
848  elseif TYP==4, % events,
849  evcount(CHN) = evcount(CHN) + 1;
850  elseif TYP==5, % continous
851  t1 = adpos(CHN) + 1;
852  t2 = adpos(CHN) + DUR;
853  HDR.data(t1:t2,CHN) = wf';
854  adpos(CHN) = t2;
855  end;
856  NRec = NRec+1;
857  typ_ubyte = fread(HDR.FILE.FID,2,'int16');
858  end
859  fclose(HDR.FILE.FID);
860  HDR.FILE.OPEN = 0;
861 
862  HDR.PLX2.tscount = tscount;
863  HDR.PLX2.wfcount = wfcount;
864  HDR.PLX2.evcount = evcount;
865  HDR.PLX2.adcount = adcount;
866 
867  ET = ET(1:nET,:);
868  HDR.EVENT.ET = ET;
869  HDR.EVENT.TYP = ET(:,1);
870  HDR.EVENT.POS = ET(:,2);
871  HDR.EVENT.CHN = ET(:,3);
872  HDR.EVENT.DUR = ET(:,4);
873  HDR.EVENT.unit= ET(:,5);
874 
875  HDR.TYPE = 'native';
876  HDR.data = HDR.data(:,HDR.InChanSelect);
877 
878  % sequence for reading "native" format
879  if nargin>2,
880  HDR.FILE.POS = HDR.SampleRate*StartPos;
881  end;
882 
883  nr = min(round(HDR.SampleRate * NoS), size(HDR.data,1) - HDR.FILE.POS);
884  S = HDR.data(HDR.FILE.POS + (1:nr), :);
885  HDR.FILE.POS = HDR.FILE.POS + nr;
886 
887 
888 elseif strcmp(HDR.TYPE,'SCP'),
889  % this branch is not in use yet.
890  % its experiemental to improve performance
891 
892  % decompress data in first call
893  HT = HDR.Huffman.HT;
894 
895  dd = [0:255]';
896  ACC = zeros(size(dd));
897  c = 0;
898  for k2 = 1:8,
899  ACC = ACC + (dd>127).*(2^c);
900  dd = mod(dd*2, 256);
901  c = c + 1;
902  end;
903 
904  S2 = [];
905  for k3 = 5:6,
906  if (k3==5) && isfield(HDR,'SCP5');
907  SCP = HDR.SCP5;
908  elseif (k3==6) && isfield(HDR,'SCP6');
909  SCP = HDR.SCP6;
910  else
911  SCP = [];
912  end;
913 
914  if ~isempty(SCP),
915  if ~isfield(HDR,'SCP2'),
916  S2 = SCP.data(:,HDR.InChanSelect);
917 
918  elseif HDR.SCP2.NHT==19999,
919  HuffTab = HDR.Huffman.DHT;
920  S2 = zeros(0,length(HDR.InChanSelect));
921  for k0 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k0); %HDR.NS,
922  s2 = SCP.data{k};
923  s2 = [s2; repmat(0,ceil(max(HDR.SCP2.HT(:,4))/8),1)];
924  k1 = 0;
925  l2 = 0;
926  accu = 0;
927  c = 0;
928  x = [];
929  HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==1),3:7);
930  while (l2 < HDR.LeadPos(k,2)),
931  while ((c < max(HT(:,2))) && (k1<length(s2)-1));
932  k1 = k1 + 1;
933  dd = s2(k1);
934  accu = accu + ACC(dd+1)*(2^c);
935  c = c + 8;
936  end;
937 
938  ixx = 1;
939  acc = accu - 2^32*floor(accu*(2^(-32))); % bitand returns NaN if accu >= 2^32
940  while (bitand(acc,2^HT(ixx,1)-1) ~= HT(ixx,5)),
941  ixx = ixx + 1;
942  end;
943 
944  dd = HT(ixx,2) - HT(ixx,1);
945  if HT(ixx,3)==0,
946  HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==HT(ixx,5)),3:7);
947  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: Switching Huffman Tables is not tested yet.\n');
948  elseif (dd==0),
949  l2 = l2 + 1;
950  x(l2) = HT(ixx,4);
951  else %if (HT(ixx,3)>0),
952  l2 = l2 + 1;
953 
954  tmp = floor(accu*(2^(-HT(ixx,1)))); %
955  %tmp = bitshift(accu,-HT(ixx,1));
956  tmp = tmp - (2^dd)*floor(tmp*(2^(-dd))); %
957  %tmp = bitand(tmp,2^dd)
958 
959  % reverse bit-pattern
960  if dd==8,
961  tmp = ACC(tmp+1);
962  else
963  tmp = dec2bin(tmp);
964  tmp = [char(repmat('0',1,dd-length(tmp))),tmp];
965  tmp = bin2dec(tmp(length(tmp):-1:1));
966  end
967  x(l2) = tmp-(tmp>=(2^(dd-1)))*(2^dd);
968  end;
969  accu = floor(accu*2^(-HT(ixx,2)));
970  c = c - HT(ixx,2);
971  end;
972 
973  if k0==1,
974  S2 = x';
975  elseif size(x,2)==size(S2,1),
976  S2(:,k0) = x';
977  else
978  fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
979  return;
980  end;
981  end;
982 
983 
984  elseif (HDR.SCP2.NHT==19999), % alternative decoding algorithm.
985  HuffTab = HDR.Huffman.DHT;
986 
987  for k0 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k0); %HDR.NS,
988  tmp = SCP.data{k};
989  accu = [tmp(4)+256*tmp(3)+65536*tmp(2)+2^24*tmp(1)];
990  %accu = bitshift(accu,HDR.SCP2.prefix,32);
991  c = 0; %HDR.SCP2.prefix;
992  l = 4;
993  l2 = 0;
994  clear x;
995  Ntmp = length(tmp);
996  tmp = [tmp; zeros(4,1)];
997  while c <= 32, %1:HDR.SPR(k),
998  ixx = 1;
999  while (bitand(accu,HDR.Huffman.mask(ixx)) ~= HDR.Huffman.PREFIX(ixx)),
1000  ixx = ixx + 1;
1001  end;
1002 
1003  if ixx < 18,
1004  c = c + HDR.Huffman.prefix(ixx);
1005  %accu = bitshift(accu, HDR.Huffman.prefix(ixx),32);
1006  accu = mod(accu.*(2^HDR.Huffman.prefix(ixx)),2^32);
1007  l2 = l2 + 1;
1008  x(l2) = HuffTab(ixx,1);
1009 
1010  elseif ixx == 18,
1011  c = c + HDR.Huffman.prefix(ixx) + 8;
1012  %accu = bitshift(accu, HDR.Huffman.prefix(ixx),32);
1013  accu = mod(accu.*(2^HDR.Huffman.prefix(ixx)),2^32);
1014  l2 = l2 + 1;
1015 
1016  acc1 = mod(floor(accu*2^(-24)),256);
1017  %accu = bitshift(accu, 8, 32);
1018  accu = mod(accu*256, 2^32);
1019 
1020  x(l2) = acc1-(acc1>=2^7)*2^8;
1021  acc2 = 0;
1022  for kk = 1:8,
1023  acc2 = acc2*2 + mod(acc1,2);
1024  acc1 = floor(acc1/2);
1025  end;
1026 
1027  elseif ixx == 19,
1028  c = c + HDR.Huffman.prefix(ixx);
1029  %accu = bitshift(accu, HDR.Huffman.prefix(ixx),32);
1030  accu = mod(accu.*(2^HDR.Huffman.prefix(ixx)),2^32);
1031  l2 = l2 + 1;
1032  while (c > 7) && (l < Ntmp),
1033  l = l+1;
1034  c = c-8;
1035  accu = accu + tmp(l)*2^c;
1036  end;
1037 
1038  acc1 = mod(floor(accu*2^(-16)),2^16);
1039  %accu = bitshift(accu, 16, 32);
1040  accu = mod(accu.*(2^16), 2^32);
1041 
1042  x(l2) = acc1-(acc1>=2^15)*2^16;
1043  acc2 = 0;
1044  for kk= 1:16,
1045  acc2 = acc2*2+mod(acc1,2);
1046  acc1 = floor(acc1/2);
1047  end;
1048  %x(l2) = acc2;
1049  c = c + 16;
1050  end;
1051 
1052  while (c > 7) && (l < Ntmp),
1053  l = l+1;
1054  c = c-8;
1055  accu = accu + tmp(l)*(2^c);
1056  end;
1057  end;
1058 
1059  x = x(1:end-1)';
1060  if k==1,
1061  S2=x;
1062  elseif size(x,1)==size(S2,1),
1063  S2(:,k0) = x;
1064  else
1065  fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
1066  return;
1067  end;
1068  end;
1069  elseif (HDR.SCP2.NHT==1) && (HDR.SCP2.NCT==1) && (HDR.SCP2.prefix==0),
1070  S2 = SCP.data(:,HDR.InChanSelect);
1071 
1072  elseif HDR.SCP2.NHT~=19999,
1073  fprintf(HDR.FILE.stderr,'Warning SOPEN SCP-ECG: user specified Huffman Table not supported\n');
1074  return;
1075  else
1076  HDR.SCP2,
1077  end;
1078 
1079  % Decoding of Difference encoding
1080  if SCP.FLAG.DIFF==2,
1081  for k1 = 3:size(S2,1);
1082  S2(k1,:) = S2(k1,:) + [2,-1] * S2(k1-(1:2),:);
1083  end;
1084  elseif SCP.FLAG.DIFF==1,
1085  S2 = cumsum(S2);
1086  end;
1087  S2 = S2 * SCP.Cal;
1088 
1089  if (k3==5) && isfield(HDR,'SCP5');
1090  % HDR.SCP5.data = S2;
1091 
1092  elseif (k3==6) && isfield(HDR,'SCP6');
1093  if HDR.SCP6.FLAG.bimodal_compression,
1094  F = HDR.SCP5.SampleRate/HDR.SCP6.SampleRate;
1095  HDR.SampleRate = HDR.SCP5.SampleRate;
1096  HDR.FLAG.F = F;
1097 
1098  tmp=[HDR.SCP4.PA(:,1);HDR.LeadPos(1,2)]-[1;HDR.SCP4.PA(:,2)+1];
1099  if ~all(tmp==floor(tmp))
1100  tmp,
1101  end;
1102  t = (1:HDR.N) / HDR.SampleRate;
1103  S1 = zeros(HDR.N, HDR.NS);
1104 
1105 
1106  p = 1;
1107  k2 = 1;
1108  pa = [HDR.SCP4.PA;NaN,NaN];
1109  flag = 1;
1110  for k1 = 1:HDR.N,
1111  if k1 == pa(p,2)+1,
1112  flag = 1;
1113  p = p+1;
1114  accu = S2(k2,:);
1115  elseif k1 == pa(p,1),
1116  flag = 0;
1117  k2 = ceil(k2);
1118  end;
1119 
1120  if flag,
1121  S1(k1,:) = ((F-1)*accu + S2(floor(k2),:)) / F;
1122  k2 = k2 + 1/F;
1123  else
1124  S1(k1,:) = S2(k2,:);
1125  k2 = k2 + 1;
1126  end;
1127  end;
1128 
1129  HDR.SCP.S2 = S2;
1130  HDR.SCP.S1 = S1;
1131  S2 = S1;
1132  end;
1133 
1134  if HDR.FLAG.ReferenceBeat,
1135  for k = find(~HDR.SCP4.type(:,1)'),
1136  t1 = (HDR.SCP4.type(k,2):HDR.SCP4.type(k,4));
1137  t0 = t1 - HDR.SCP4.type(k,3) + HDR.SCP4.fc0;
1138  S2(t1,:) = S2(t1,:) + HDR.SCP5.data(t0,:);
1139  end;
1140  end;
1141  end;
1142  end;
1143  end;
1144  HDR.data = S2;
1145  HDR.TYPE = 'native'; % decompression is already applied
1146 
1147  if nargin>2,
1148  HDR.FILE.POS = HDR.SampleRate*StartPos;
1149  end;
1150 
1151  nr = min(round(HDR.SampleRate * NoS), size(HDR.data,1) - HDR.FILE.POS);
1152  S = HDR.data(HDR.FILE.POS + (1:nr), :);
1153  HDR.FILE.POS = HDR.FILE.POS + nr;
1154 
1155 
1156 elseif strcmp(HDR.TYPE,'SIGIF'),
1157  if nargin==3,
1158  HDR.FILE.POS = StartPos;
1159  end;
1160 
1161  S = [];
1162  for k = 1:min(NoS,HDR.NRec-HDR.FILE.POS),
1163  HDR.FILE.POS = HDR.FILE.POS + 1;
1164  STATUS = fseek(HDR.FILE.FID, HDR.Block.Pos(HDR.FILE.POS), 'bof');
1165  if HDR.FLAG.TimeStamp,
1166  HDR.Frame(k).TimeStamp = fread(HDR.FILE.FID,[1,9],'uint8');
1167  end;
1168 
1169  if HDR.FLAG.SegmentLength,
1170  HDR.Block.Length(k) = fread(HDR.FILE.FID,1,'uint16'); %#26
1171  STATUS = fseek(HDR.FILE.FID,HDR.Block.Length(k)*H1.Bytes_per_Sample,'cof');
1172  else
1173  tmp = HDR.Segment_separator-1;
1174  [dat,c] = fread(HDR.FILE.FID,[HDR.NS,HDR.Block.Length/HDR.NS],gdfdatatype(HDR.GDFTYP));
1175  [tmpsep,c] = fread(HDR.FILE.FID,1,gdfdatatype(HDR.GDFTYP));
1176 
1177  if (tmpsep~=HDR.Segment_separator);
1178  fprintf(HDR.FILE.stderr,'Error SREAD Type=SIGIF: blockseparator not found\n');
1179  end;
1180  end;
1181  S = [S; dat(HDR.InChanSelect,:)'];
1182  end;
1183 
1184 
1185 elseif strcmp(HDR.TYPE,'CTF'),
1186  if nargin>2,
1187  STATUS = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.NS*HDR.SPR*4*StartPos,'bof');
1188  HDR.FILE.POS = StartPos;
1189  end;
1190 
1191  nr = min(NoS, HDR.NRec - HDR.FILE.POS);
1192 
1193  S = []; count = 0;
1194  for k = 1:nr,
1195  %[tmp,c] = fread(HDR.FILE.FID, 1, 'int32')
1196  [s,c] = fread(HDR.FILE.FID, [HDR.SPR, HDR.NS], 'int32');
1197  S = [S; s(:,HDR.InChanSelect)];
1198  count = count + c;
1199  end;
1200 
1201  HDR.FILE.POS = HDR.FILE.POS + count/(HDR.SPR*HDR.NS);
1202 
1203 
1204 elseif strcmp(HDR.TYPE,'EEProbe-CNT'),
1205  if nargin>2,
1206  HDR.FILE.POS = HDR.SampleRate*StartPos;
1207  end;
1208 
1209  nr = min(HDR.SampleRate*NoS, HDR.SPR*HDR.NRec-HDR.FILE.POS);
1210  if exist('read_eep_cnt','file')==3,
1211  tmp = read_eep_cnt(HDR.FileName, HDR.FILE.POS+1, HDR.FILE.POS+nr);
1212  sz = size(tmp.data);
1213  S = tmp.data(HDR.InChanSelect,:)';
1214  clear tmp;
1215  if HDR.FLAG.UCAL,
1216  S = S*diag(1./HDR.Cal(HDR.InChanSelect));
1217  end;
1218  HDR.FILE.POS = HDR.FILE.POS + sz(2);
1219 
1220  elseif exist('OCTAVE_VERSION','builtin')
1221  fprintf(HDR.FILE.stderr,'ERROR SREAD (EEProbe): Reading EEProbe-file format is not supported.\n');
1222  return;
1223  else
1224  fprintf(HDR.FILE.stderr,'ERROR SREAD (EEProbe): Cannot open EEProbe-file, because read_eep_cnt.mex not installed. \n');
1225  fprintf(HDR.FILE.stderr,'ERROR SREAD (EEProbe): You can downlad it from http://www.smi.auc.dk/~roberto/eeprobe/\n');
1226  %% ftp://ftp.fcdonders.nl/pub/fieldtrip/external/eeprobe.zip
1227  return;
1228  end
1229 
1230 elseif strcmp(HDR.TYPE,'EEProbe-AVR'),
1231  if nargin>2,
1232  HDR.FILE.POS = HDR.SampleRate*StartPos;
1233  end;
1234 
1235  nr = min(HDR.SPR-HDR.FILE.POS,NoS*HDR.SampleRate);
1236 
1237  S = HDR.EEP.data(HDR.FILE.POS+(1:nr),:);
1238  HDR.FILE.POS = HDR.FILE.POS + nr;
1239 
1240 
1241 elseif strncmp(HDR.TYPE,'BrainVision',11), %Brainvision
1242  if strncmpi(HDR.BV.DataFormat, 'binary',5)
1243  tc = strcmp(HDR.TYPE,'BrainVisionVAmp');
1244  NS = HDR.NS+tc;
1245  if strncmpi(HDR.BV.DataOrientation, 'multiplexed',6),
1246  if nargin>2,
1247  STATUS = fseek(HDR.FILE.FID,StartPos*HDR.SampleRate*HDR.AS.bpb,'bof');
1248  HDR.FILE.POS = HDR.SampleRate*StartPos;
1249  end;
1250 
1251  nr = min(HDR.SampleRate*NoS, HDR.SPR*HDR.NRec - HDR.FILE.POS);
1252  if (length(HDR.InChanSelect)*2>HDR.NS)
1253  [s,c] = fread(HDR.FILE.FID, [NS, nr], ['*',gdfdatatype(HDR.GDFTYP)]);
1254  count = c/NS;
1255  S = s(HDR.InChanSelect,:)';
1256  else
1257  S = [];
1258  count = 0;
1259  while (count<nr),
1260  [s,c] = fread(HDR.FILE.FID, [NS, min(nr-count,floor(2^24/NS))], gdfdatatype(HDR.GDFTYP));
1261  if ~c, break; end;
1262  S = [S; s(HDR.InChanSelect,:)'];
1263  count = count + c/NS;
1264  end;
1265  end;
1266  HDR.FILE.POS = HDR.FILE.POS + count;
1267  elseif strncmpi(HDR.BV.DataOrientation, 'vectorized',6),
1268  S = [];
1269  nr = min(HDR.SampleRate*NoS, HDR.AS.endpos-HDR.FILE.POS);
1270 
1271  count = 0;
1272  for chan = 1:length(HDR.InChanSelect);
1273  STATUS = fseek(HDR.FILE.FID, HDR.HeadLen + HDR.FILE.POS + HDR.AS.bpb*HDR.SPR*(chan-1)/NS, 'bof');
1274  [s,count] = fread(HDR.FILE.FID, [nr,1], gdfdatatype(HDR.GDFTYP));
1275  if count ~= nr,
1276  fprintf(2,'ERROR READ BV-bin-vec: \n');
1277  return;
1278  end;
1279  S(:,chan) = s;
1280  end
1281  HDR.FILE.POS = HDR.FILE.POS + count;
1282  end;
1283 
1284  elseif strncmpi(HDR.BV.DataFormat, 'ascii',5)
1285  %%%% OBSOLETE: supported by 'native' %%%%
1286  if nargin>2,
1287  HDR.FILE.POS = HDR.SampleRate*StartPos;
1288  end;
1289  nr = min(HDR.SampleRate*NoS, HDR.AS.endpos-HDR.FILE.POS);
1290  S = HDR.BV.data(HDR.FILE.POS+(1:nr),HDR.InChanSelect);
1291 
1292  end
1293 
1294 
1295 elseif strcmp(HDR.TYPE,'SierraECG'), %% SierraECG 1.03 *.open.xml from PHILIPS
1296  if ~isfield(HDR,'data');
1297  [HDR.data,status] = str2double(HDR.XML.waveforms.parsedwaveforms);
1298  if any(status)
1299  error('SREAD: compressed SierraECG (Philips) format not supported')
1300  end;
1301  HDR.data = reshape(HDR.data,length(HDR.data)/HDR.NS,HDR.NS);
1302  HDR.SPR = size(HDR.data,1);
1303  else
1304  % base64 - decoding
1305  base64 = ['A':'Z','a':'z','0':'9','+','/'];
1306  decode64 = repmat(nan,256,1);
1307  decode64(abs(base64)) = 0:63;
1308  tmp = decode64(HDR.XML.waveforms.parsedwaveforms);
1309  tmp(isnan(tmp)) = [];
1310  n = length(tmp);
1311  tmp = reshape([tmp;zeros(mod(n,4),1)], 4, ceil(n/4));
1312  t1 = tmp(1,:)*4 + floor(tmp(2,:)/16);
1313  t2 = mod(tmp(2,:),16)*16 + floor(tmp(3,:)/4);
1314  t3 = mod(tmp(3,:),4)*64 + tmp(4,:);
1315  tmp = reshape([t1,t2,t3], ceil(n/4)*3, 1);
1316  end;
1317  if nargin>2,
1318  HDR.FILE.POS = HDR.SampleRate*StartPos;
1319  end;
1320  nr = min(HDR.SampleRate*NoS, HDR.SPR-HDR.FILE.POS);
1321  S = HDR.data(HDR.FILE.POS+(1:nr),HDR.InChanSelect);
1322  HDR.FILE.POS = HDR.FILE.POS + nr;
1323 
1324 
1325 elseif strcmp(HDR.TYPE,'ATF');
1326  if HDR.FILE.OPEN,
1327  fseek(HDR.FILE.FID,HDR.HeadLen,-1);
1328  t = fread(HDR.FILE.FID,[1,inf],'uint8');
1329  fclose(HDR.FILE.FID);
1330  HDR.FILE.OPEN=0;
1331  [HDR.ATF.NUM,status,HDR.ATF.STR] = str2double(char(t));
1332  end;
1333  S = HDR.ATF.NUM;
1334 
1335 
1336 elseif strcmp(HDR.TYPE,'FEPI3');
1337  if nargin==3,
1338  HDR.FILE.POS = HDR.SampleRate*StartPos;
1339  end;
1340  Duration = min(NoS*HDR.SampleRate,(HDR.AS.endpos-HDR.FILE.POS));
1341 
1342  S = repmat(NaN,Duration,length(HDR.InChanSelect));
1343  ix = find([HDR.FEPI.SEG(:,2) > HDR.FILE.POS] & [HDR.FEPI.SEG(:,1) <= HDR.FILE.POS+Duration]);
1344  for k = 1:length(ix),
1345  fid = fopen(fullfile(HDR.FILE.Path,[HDR.FEPI.ListOfDataFiles{k},'.bin']));
1346  if k==1,
1347  fseek(fid,(HDR.FILE.POS-HDR.FEPI.SEG(ix(k),1))*2*HDR.NS,-1);
1348  end;
1349  data= fread(fid,[HDR.NS,inf],'int16')';
1350  ix1 = HDR.FEPI.SEG(ix(k),1)-HDR.FEPI.SEG(ix(1),1);
1351  ix2 = min(size(data,1),Duration+HDR.FILE.POS-HDR.FEPI.SEG(ix(k),1)+1);
1352  S(ix1+1:ix1+ix2,:) = data(1:ix2,HDR.InChanSelect);
1353  fclose(fid);
1354  end;
1355  HDR.FILE.POS = HDR.FILE.POS+Duration;
1356 
1357 
1358 elseif strcmp(HDR.TYPE,'WG1'), %walter-graphtek
1359  % code from Robert Reijntjes, Amsterdam, NL
1360  % modified by Alois Schloegl 19. Feb 2005
1361  if nargin==3,
1362  HDR.FILE.POS = round(HDR.SampleRate*StartPos);
1363  end;
1364 
1365  ix1 = mod(HDR.FILE.POS, HDR.SPR); % starting sample (minus one) within 1st block
1366  fp = HDR.HeadLen + floor(HDR.FILE.POS/HDR.SPR)*HDR.AS.bpb;
1367  status = fseek(HDR.FILE.FID, fp, 'bof');
1368 
1369  nr = min(HDR.AS.endpos-HDR.FILE.POS, NoS*HDR.SampleRate);
1370  S = repmat(NaN,nr,length(HDR.InChanSelect));
1371  count = 0;
1372  endloop= 0;
1373  c = 1;
1374  offset = 0;
1375  while ~endloop & (c>0) && (offset(1)~=(hex2dec('AEAE5555')-2^32)) && (count<nr);
1376  [offset,c] = fread(HDR.FILE.FID, HDR.WG1.szOffset, 'int32');
1377  [databuf,c] = fread(HDR.FILE.FID,[HDR.WG1.szBlock,HDR.NS+HDR.WG1.szExtra],'uint8');
1378  dt = HDR.WG1.conv(databuf(:,1:HDR.NS)+1);
1379  if any(dt(:)==HDR.WG1.unknownNr),
1380  %dt(dt==HDR.WG1.unknownNr) = NaN;
1381  dt(:) = NaN;
1382  %fprintf(HDR.FILE.stderr,'Warning SREAD (WG1): error in reading datastream');
1383  end;
1384  dt(1,:) = dt(1,:) + offset(1:HDR.NS)';
1385  dt = cumsum(dt,1);
1386 
1387  ix2 = min(nr-count, size(dt,1)-ix1);
1388  S(count+1:count+ix2,:) = dt(ix1+1:ix1+ix2, HDR.InChanSelect);
1389  count = count + ix2;
1390  ix1 = 0; % reset starting index,
1391 
1392  k = 0;
1393  while (k<HDR.WG1.szExtra) && ~endloop,
1394  endloop = ~isempty(strfind(databuf(:,HDR.NS+k)',[85,85,174,174]));
1395  k = k+1;
1396  end;
1397  end;
1398  %S = S(1:count,:);
1399  HDR.FILE.POS = HDR.FILE.POS + count;
1400 
1401 
1402 elseif strcmp(HDR.TYPE,'XML-FDA'), % FDA-XML Format
1403  if ~isfield(HDR,'data');
1404  tmp = HDR.XML.component.series.derivation;
1405  if isfield(tmp,'Series');
1406  tmp = tmp.Series.component.sequenceSet.component;
1407  else % Dovermed.CO.IL version of format
1408  tmp = tmp.derivedSeries.component.sequenceSet.component;
1409  end;
1410  for k = 1:length(HDR.InChanSelect);
1411  HDR.data(:,k) = str2double(tmp{HDR.InChanSelect(k)+1}.sequence.value.digits)';
1412  end;
1413  HDR.SPR = size(HDR.data,1);
1414  end;
1415  if nargin>2,
1416  HDR.FILE.POS = HDR.SampleRate*StartPos;
1417  end;
1418  nr = min(HDR.SampleRate*NoS, HDR.SPR-HDR.FILE.POS);
1419  S = HDR.data(HDR.FILE.POS+(1:nr),:);
1420  HDR.FILE.POS = HDR.FILE.POS + nr;
1421 
1422 % using XML4MAT instead of XMLTREE
1423 % str2double(HDR.XML0{end}.component{1}.series{end}.derivation{:}.derivedSeries{5}.component{1}.sequenceSet{5}.component{1}.sequence{2}.value{3}.digits)'
1424 
1425 
1426 elseif strcmp(HDR.TYPE,'FIF'),
1427  % some parts of this code are from Robert Oostenveld,
1428  if ~(exist('rawdata')==3 & exist('channames')==3)
1429  error('cannot find Neuromag import routines on your Matlab path (see http://boojum.hut.fi/~kuutela/meg-pd)');
1430  end
1431  if nargin<3,
1432  StartPos = HDR.FILE.POS/HDR.SampleRate;
1433  end
1434  if nargin>2,
1435  HDR.FILE.POS = HDR.SampleRate*StartPos;
1436  end;
1437 
1438  t1 = rawdata('goto', HDR.FILE.POS/HDR.SPR);
1439  t2 = t1;
1440  dat = [];
1441  count = 0;
1442  status = 'ok';
1443 
1444  while (t2<(StartPos + NoS)) && ~strcmp(status,'eof'),
1445  [buf, status] = rawdata('next');
1446  if 0
1447  elseif strcmp(status, 'ok')
1448  count = count + size(buf,2);
1449  dat = [dat; buf(HDR.InChanSelect,:)'];
1450  elseif strcmp(status, 'eof')
1451  elseif strcmp(status, 'skip')
1452  elseif strcmp(status, 'error')
1453  error('error reading selected data from fif-file');
1454  else
1455  error('undefined status code return from RAWDATA(FIF-file)');
1456  end
1457  t2 = rawdata('t');
1458  end
1459  t = t1*HDR.SampleRate+1:t2*HDR.SampleRate;
1460  ix = (t>StartPos*HDR.SampleRate) & (t<=(StartPos+NoS)*HDR.SampleRate);
1461  S = dat(ix,HDR.InChanSelect);
1462  HDR.FILE.POS = t2*HDR.SampleRate;
1463 
1464 elseif strcmp(HDR.TYPE,'EVENT'),
1465  s = [];
1466 
1467 elseif strncmp(HDR.TYPE,'IMAGE:',6),
1468  % forward call to IREAD
1469  [S,HDR] = iread(HDR);
1470  return;
1471 
1472 else
1473  fprintf(2,'Error SREAD: %s-format not supported yet.\n', HDR.TYPE);
1474  return;
1475 end;
1476 
1477 
1478 %%% TOGGLE CHECK - checks whether HDR is kept consist %%%
1479 if 0,
1480 global SREAD_TOGGLE_CHECK
1481 if isfield(HDR.FLAG,'TOGGLE');
1482  if HDR.FLAG.TOGGLE~=SREAD_TOGGLE_CHECK,
1483  fprintf(HDR.FILE.stderr,'Warning SREAD: [s,HDR]=sread(HDR, ...) \nYou forgot to pass HDR in %i call(s) of SREAD\n',SREAD_TOGGLE_CHECK-HDR.FLAG.TOGGLE);
1484  end;
1485 else
1486  HDR.FLAG.TOGGLE=0;
1487  SREAD_TOGGLE_CHECK=0;
1488 end;
1489 SREAD_TOGGLE_CHECK = SREAD_TOGGLE_CHECK+1;
1490 HDR.FLAG.TOGGLE = HDR.FLAG.TOGGLE+1;
1491 end;
1492 
1493 if STATUS,
1494  fprintf(HDR.FILE.stderr,'WARNING SREAD: something went wrong. Please send the files %s and BIOSIGCORE to <a.schloegl@ieee.org>',HDR.FileName);
1495  save biosigcore.mat
1496 end;
1497 
1498 if isempty(S),
1499 
1500 elseif isfield(HDR,'THRESHOLD') && HDR.FLAG.OVERFLOWDETECTION,
1501  ix = (S~=S);
1502  for k=1:length(HDR.InChanSelect),
1503  TH = THRESHOLD(HDR.InChanSelect(k),:);
1504  %ix(:,k) = (S(:,k)<=TH(1)) | (S(:,k)>=TH(2));
1505  ix = (S(:,k)<=TH(1)) | (S(:,k)>=TH(2));
1506  S(ix,k)=NaN;
1507  end
1508  if exist('double','builtin')
1509  S = double(S);
1510  end;
1511 % S(ix>0) = NaN;
1512 elseif HDR.FLAG.OVERFLOWDETECTION,
1513  % no HDR.THRESHOLD defined
1514  warning('no Threshold defined');
1515 elseif isfield(HDR,'THRESHOLD'),
1516  % automated overflow detection has been turned off
1517 end;
1518 
1519 if ~HDR.FLAG.UCAL,
1520  % S = [ones(size(S,1),1),S]*HDR.Calib;
1521  % perform the previous function more efficiently and
1522  % taking into account some specialities related to Octave sparse
1523  % data.
1524  if isempty(S), % otherwise Octave 2.1.64 could break below,
1525  if size(S,2)~=length(HDR.InChanSelect),
1526  fprintf(HDR.FILE.stderr,'Warning SREAD (%s): number of columns (%i) incorrect!\n',HDR.TYPE,size(S,2));
1527  end;
1528  S = zeros(0,length(HDR.InChanSelect));
1529  end;
1530 
1531  %if ~issparse(HDR.Calib); %
1532  if FLAG_CALIB_DONE,
1533 
1534  elseif strcmpi(HDR.FLAG.OUTPUT,'single')
1535  tmp = single(zeros(size(S,1),size(HDR.Calib,2)));
1536  for k = 1:size(S,1),
1537  tmp(k,:) = [1,S(k,:)] * HDR.Calib;
1538  end;
1539  S = tmp;
1540  clear tmp;
1541 
1542  elseif 1, % exist('OCTAVE_VERSION','builtin')
1543  % force octave to do a sparse multiplication
1544  % the difference is NaN*sparse(0) = 0 instead of NaN
1545  % this is important for the automatic overflow detection
1546  Calib = HDR.Calib;
1547  tmp = S;
1548  S = zeros(size(S,1),size(Calib,2)); % memory allocation
1549 
1550  for k = 1:size(Calib,2),
1551  chan = find(Calib(2:end,k));
1552  S(:,k) = double(tmp(:,chan)) * full(Calib(1+chan,k)) + Calib(1,k);
1553  end;
1554  else
1555  % S = [ones(size(S,1),1),S]*HDR.Calib;
1556  % the following is the same as above but needs less memory.
1557  S = full(double(S) * HDR.Calib(2:end,:));
1558  for k = 1:size(HDR.Calib,2),
1559  S(:,k) = S(:,k) + full(HDR.Calib(1,k));
1560  end;
1561  end;
1562 end;
1563 
1564 
1565 
gdfdatatype
function gdfdatatype(in GDFTYP)
function
function()
iread
function iread(in HDR, in CHAN, in StartPos)
sopen
function sopen(in arg1, in PERMISSION, in CHAN, in MODE, in arg5, in arg6)
str2double
function str2double(in s, in cdelim, in rdelim, in ddelim)
rs
function rs(in y1, in T, in f2)
sread
function sread(in HDR, in NoS, in StartPos)