eidors_readdata

 EIDORS readdata - read data files from various EIT equipment
    manufacturers

 [vv, auxdata, stim ]= eidors_readdata( fname, format, frame_range, extra )
    frame_range is the list of frames to load from the file (ie. 1:100).
    if the frames are beyond the number in the file, no data is returned
    to get all frames, use frame_range= []

 Currently the list of supported file formats is:
    - "EIT" Files, may be one of 
          - Draeger Pulmovista (2008+)
          - GoeIIMF/Carefusion (2010+)
          - Swisstom BB2/Pioneer (2010+)
       The function will attempt to autodetect the format
       
    - MCEIT (GoeIIMF / Viasys) "get" file format 
        format = "GET" or "MCEIT"
        Note that the data is "untwisted" to correspond to a "no_rotate_meas" stim pattern
    - MCEIT (GoeIIMF) "get" file format
        format = "GET-RAW"
        Data in original order, corresponds to "rotate_meas" stim pattern

    - Draeger (Pulmovista) "EIT" file format (2008+)
        format = "DRAEGER-EIT"

    - Draeger "get" file format (- 2007 - format for Draeger equipment)
        format = "DRAEGER-GET"

    - Swisstom EIT equipment "EIT" (Pioneer set and BB2):
           'LQ1' (2010 - 2011)
           'LQ2' (2013 - 2014)
           'LQ4' (2015+)

    - Dixtal file format, from Dixtal inc, Brazil
        format = 'DIXTAL_encode' extract encoder from provided Dll
           - Output: [encodepage] = eidors_readdata( path,'DIXTAL_encode');
              where path= '/path/to/Criptografa_New.dll' (provided with system)
        format = 'DIXTAL'
           - output: [vv] = eidors_readdata( fname, 'DIXTAL', [], encodepage );
    - New Carefusion "EIT" file format
        format = "GOEIIMF-EIT" or "carefusion"

0001 function [vv, auxdata, stim ]= eidors_readdata( fname, format, frame_range, extra )
0002 % EIDORS readdata - read data files from various EIT equipment
0003 %    manufacturers
0004 %
0005 % [vv, auxdata, stim ]= eidors_readdata( fname, format, frame_range, extra )
0006 %    frame_range is the list of frames to load from the file (ie. 1:100).
0007 %    if the frames are beyond the number in the file, no data is returned
0008 %    to get all frames, use frame_range= []
0009 %
0010 % Currently the list of supported file formats is:
0011 %    - "EIT" Files, may be one of
0012 %          - Draeger Pulmovista (2008+)
0013 %          - GoeIIMF/Carefusion (2010+)
0014 %          - Swisstom BB2/Pioneer (2010+)
0015 %       The function will attempt to autodetect the format
0016 %
0017 %    - MCEIT (GoeIIMF / Viasys) "get" file format
0018 %        format = "GET" or "MCEIT"
0019 %        Note that the data is "untwisted" to correspond to a "no_rotate_meas" stim pattern
0020 %    - MCEIT (GoeIIMF) "get" file format
0021 %        format = "GET-RAW"
0022 %        Data in original order, corresponds to "rotate_meas" stim pattern
0023 %
0024 %    - Draeger (Pulmovista) "EIT" file format (2008+)
0025 %        format = "DRAEGER-EIT"
0026 %
0027 %    - Draeger "get" file format (- 2007 - format for Draeger equipment)
0028 %        format = "DRAEGER-GET"
0029 %
0030 %    - Swisstom EIT equipment "EIT" (Pioneer set and BB2):
0031 %           'LQ1' (2010 - 2011)
0032 %           'LQ2' (2013 - 2014)
0033 %           'LQ4' (2015+)
0034 %
0035 %    - Dixtal file format, from Dixtal inc, Brazil
0036 %        format = 'DIXTAL_encode' extract encoder from provided Dll
0037 %           - Output: [encodepage] = eidors_readdata( path,'DIXTAL_encode');
0038 %              where path= '/path/to/Criptografa_New.dll' (provided with system)
0039 %        format = 'DIXTAL'
0040 %           - output: [vv] = eidors_readdata( fname, 'DIXTAL', [], encodepage );
0041 %    - New Carefusion "EIT" file format
0042 %        format = "GOEIIMF-EIT" or "carefusion"
0043 
0044 %    - Sheffield MK I "RAW" file format
0045 %        format = "RAW" or "sheffield"
0046 %    - ITS (International Tomography Systems)
0047 %        format = "ITS" or "p2k"
0048 %    - IIRC (Impedance Imaging Research Center, Korea)
0049 %        format = "txt" or "IIRC"
0050 %    - University of Cape Town formats
0051 %        format = "UCT_SEQ"  UCT sequence file
0052 %           - Output: [ stimulation, meas_select]= eidors_readdata(fname, 'UCT_SEQ')
0053 %        format = "UCT_CAL"  UCT calibration file
0054 %           - Output: [vv, no_cur_caldata_raw ]= eidors_readdata( fname, 'UCT_CAL' )
0055 %                 where no_cur_caldata_raw is data captured with no current
0056 %        format = "UCT_DATA"  UCT data frame file
0057 %           - Output: [vv]= eidors_readdata( fname, 'UCT_DATA' )
0058 %
0059 % Usage
0060 % [vv, auxdata, stim ]= eidors_readdata( fname, format )
0061 %     vv      = measurements - data frames in each column
0062 %     auxdata = auxillary data - if provided by system
0063 %     stim    = stimulation structure, to be used with
0064 %                fwd_model.stimulation.
0065 %     fname = file name
0066 %     stim.framerate = acquisition rate (frames/sec) if available
0067 %
0068 %  if format is unspecified, an attempt to autodetect is made
0069 
0070 % (C) 2005-09 Andy Adler. License: GPL version 2 or version 3
0071 % $Id: eidors_readdata.m 5945 2019-05-24 20:16:47Z aadler $
0072 
0073 % TODO:
0074 %   - output an eidors data object
0075 %   - test whether file format matches specified stimulation pattern
0076 %   - todo MCEIT provides curr and volt on driven electrodes.
0077 %       how can this be provided to system?
0078 
0079 if ~exist(fname,'file')
0080    error([fname,' does not exist']);
0081 end
0082 
0083 if nargin < 2
0084 % unspecified file format, autodetect
0085    dotpos = find(fname == '.');
0086    if isempty( dotpos ) 
0087       error('file format unspecified, can`t autodetect');
0088    else
0089       dotpos= dotpos(end);
0090       format= fname( dotpos+1:end );
0091    end
0092 end
0093 
0094 
0095 auxdata = []; % default, can be overriden if the format has it
0096 fmt = pre_proc_spec_fmt( format, fname );
0097 switch fmt
0098    case 'mceit';
0099       [vv,curr,volt,auxdata_out,auxtime,rawdata] = mceit_readdata( fname );
0100       auxdata.auxdata = auxdata_out;
0101       auxdata.auxtime = auxtime;
0102       auxdata.curr    = curr;
0103       auxdata.volt    = volt;
0104 
0105 
0106       if strcmp(lower(format),'get-raw')
0107           vv= rawdata(1:208,:);
0108           stim = mk_stim_patterns(16,1,[0,1],[0,1], {'no_meas_current','rotate_meas'}, 1);
0109       else
0110           stim = basic_stim(16);
0111       end
0112    case 'draeger-get'
0113       vv = draeger_get_readdata( fname );
0114 
0115       stim = basic_stim(16);
0116 
0117    case 'draeger-eit'
0118      [fr] = read_draeger_header( fname );
0119      % Currently Draeger equipment uses this pattern with 5mA injection
0120      stim = mk_stim_patterns(16,1,[0,1],[0,1],{'rotate_meas','no_meas_current'},.005);
0121      [stim(:).framerate] = deal(fr);
0122      [vv] = read_draeger_file( fname );
0123      auxdata = vv; % the actual voltages
0124      % Based on the draeger *get file converter
0125      ft = [.00098242, .00019607];% estimated AA: 2016-04-07
0126      vv = ft(1)*vv(1:208,:) - ft(2)*vv(322+(1:208),:);
0127 
0128    case {'raw', 'sheffield'}
0129       vv = sheffield_readdata( fname );
0130 
0131       stim = basic_stim(16);
0132    case {'p2k', 'its'}
0133       vv = its_readdata( fname );
0134 
0135       stim = 'UNKNOWN';
0136    case {'txt','iirc'}
0137       vv = iirc_readdata( fname );
0138 
0139       stim = basic_stim(16);
0140    case 'uct_seq'
0141       [vv,auxdata] = UCT_sequence_file( fname );
0142 
0143       stim = 'UNKNOWN';
0144    case 'uct_cal'
0145       [vv,auxdata] = UCT_calibration_file( fname );
0146 
0147       stim = 'UNKNOWN';
0148    case 'uct_data'
0149       [vv] = UCT_LoadDataFrame( fname );
0150 
0151       stim = 'UNKNOWN';
0152    case {'carefusion','goeiimf-eit'}
0153       [vv, auxdata_out, auxtime] = carefusion_eit_readdata( fname );
0154       auxdata.auxdata = auxdata_out;
0155       auxdata.auxtime = auxtime;
0156   
0157       stim = mk_stim_patterns(16,1,[0,1],[0,1], {'no_meas_current','rotate_meas'}, .005);
0158 
0159    case 'lq1'
0160       [vv] = landquart1_readdata( fname );
0161 
0162       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0163       
0164    case {'lq2','lq3'}
0165       [vv, elecImps, tStampsAbs, tStampsRel] = landquart2_readdata( fname );
0166       auxdata.elec_impedance = elecImps;
0167       auxdata.t_abs = tStampsAbs;
0168       auxdata.t_rel = tStampsRel;
0169 
0170       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0171      
0172    case 'lq4pre'
0173       [vv] = landquart4pre_readdata( fname );
0174 
0175       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0176 
0177    case 'lq4'
0178       [vv, evtlist, elecImps, tStampsAbs, tStampsRel] = landquart4_readdata( fname );
0179       auxdata.event = evtlist;
0180       auxdata.elec_impedance = elecImps;
0181       auxdata.t_abs = tStampsAbs;
0182       auxdata.t_rel = tStampsRel;
0183 
0184       stim = mk_stim_patterns(32,1,[0,5],[0,5],{'no_rotate_meas','no_meas_current'},.005);
0185       
0186    case 'dixtal_encode'
0187       [vv] = dixtal_read_codepage( fname );
0188       stim = 'N/A';
0189 
0190    case 'dixtal'
0191       [vv] = dixtal_read_data( fname, frame_range, extra );
0192       auxdata = vv(1025:end,:);
0193       vv      = vv([1:1024],:);
0194       stim= mk_stim_patterns(32,1,[0,4],[0,4], ... 
0195               {'meas_current','no_rotate_meas'}, 1);
0196 
0197    otherwise
0198       error('eidors_readdata: file "%s" format unknown', format);
0199 end
0200 
0201 function stim = basic_stim(N_el);
0202    stim= mk_stim_patterns(16,1,[0,1],[0,1], ... \
0203           {'no_meas_current','no_rotate_meas'}, 1);
0204 
0205 function fmt = pre_proc_spec_fmt( format, fname );
0206    fmt= lower(format);
0207    if strcmp(fmt,'get')
0208       if is_get_file_a_draeger_file( fname)
0209          fmt= 'draeger-get';
0210       else
0211          fmt= 'mceit';
0212       end
0213    end
0214 
0215    if strcmp(fmt,'get-raw')
0216       fmt= 'mceit';
0217    end
0218    
0219 
0220    if strcmp(fmt,'eit')
0221       draeger =   is_eit_file_a_draeger_file( fname );
0222       swisstom=   is_eit_file_a_swisstom_file( fname );
0223       carefusion= is_eit_file_a_carefusion_file( fname );
0224       if carefusion
0225          eidors_msg('"%s" appears to be in GOEIIMF/Carefusion format',fname,3);
0226          fmt= 'carefusion';
0227       elseif draeger
0228          eidors_msg('"%s" appears to be in Draeger format',fname,3);
0229          fmt= 'draeger-eit';
0230       elseif swisstom
0231          fmt= sprintf('lq%d',swisstom);
0232          if swisstom == 3.5
0233             fmt= 'lq4pre';
0234          end
0235          eidors_msg('"%s" appears to be in %s format',fname,upper(fmt),3);
0236       else
0237          error('EIT file specified, but it doesn''t seem to be a Carefusion file')
0238       end
0239    end
0240 
0241 function df= is_get_file_a_draeger_file( fname)
0242    fid= fopen(fname,'rb');
0243    d= fread(fid,[1 26],'uchar');
0244    fclose(fid);
0245    df = all(d == '---Draeger EIT-Software---');
0246 
0247 function df= is_eit_file_a_draeger_file( fname );
0248    fid= fopen(fname,'rb');
0249    d= fread(fid,[1 80],'uchar');
0250    fclose(fid);
0251    ff = findstr(d, '---Draeger EIT-Software');
0252    if ff;
0253       df = 1;
0254       eidors_msg('Draeger format: %s', d(ff(1)+(0:30)),4);
0255    else 
0256       df = 0;
0257    end
0258 
0259 function df= is_eit_file_a_swisstom_file( fname );
0260    fid= fopen(fname,'rb');
0261    d= fread(fid,[1 4],'uchar');
0262    fclose(fid);
0263    
0264 % Note that there are two possible LQ4 formats, either with
0265 %     d=[0,0,0,4] or with d = [4,0,0,0]
0266 %  we call the first one version 3.5
0267    if any(d(4)==[2,3,4]) && all(d([1,2,3]) == 0);
0268       df = d(4);
0269       if d(4)==4; df = 3.5; end
0270 %     eidors_msg('Swisstom format: LQ%d', df, 4);
0271    elseif any(d(1:4) == [4,0,0,0])
0272       df = 4;
0273    else 
0274       df = 0;
0275    end
0276 
0277 function df = is_eit_file_a_carefusion_file( fname )
0278    fid= fopen(fname,'rb');
0279    d= fread(fid,[1 80],'uchar');
0280    fclose(fid);
0281    df = 0;
0282    if d(1:2) ~= [255,254]; return; end
0283    if d(4:2:end) ~= 0; return; end
0284    str= setstr(d(3:2:end));
0285    tests1= '<?xml version="1.0" ?>';
0286    tests2= '<EIT_File>';
0287    if ~any(findstr(str,tests1)); return; end
0288    if ~any(findstr(str,tests2)); return; end
0289    
0290    df= 1;
0291    return
0292 
0293 function [vv, auxdata, auxtime] = carefusion_eit_readdata( fname );
0294    fid= fopen(fname,'rb');
0295    d= fread(fid,[1 180],'uchar');
0296    str= setstr(d(3:2:end));
0297    outv = regexp(str,'<Header>(\d+) bytes</Header>','tokens');
0298    if length(outv) ~=1; 
0299       error('format problem reading carefusion eit files');
0300    end
0301    header = str2num(outv{1}{1});
0302 
0303 % NOTE: we're throwing away the header completely. This isn't
0304 % the 'right' way to read a file.
0305 
0306    fseek(fid, header, -1); % From the beginning
0307    d_EIT= fread(fid,[inf],'float32');
0308 
0309    fseek(fid, header, -1); % From the beginning
0310    d_struct= fread(fid,[inf],'int32');  % In the struct, meaning of every int: 1 Type, 2 Waveform channel No., 3 sequence No., 4 size in byte, 5 count
0311    fclose(fid);
0312    pt_EIT=1;               % pointer of the EIT data
0313    vv=[];
0314    auxdata=[];
0315    
0316    while (pt_EIT<=length(d_struct))
0317       switch d_struct(pt_EIT)
0318          case 3, % type=3,  EIT transimpedance measurement
0319            % d_struct(pt_EIT+2), Sequence No., start from 0
0320            vv(1:256,1+d_struct(pt_EIT+2))= d_EIT(pt_EIT+6:pt_EIT+6+255);
0321            pt_EIT=pt_EIT+6+256;
0322          case 7, %type=7, EIT image vector
0323            pt_EIT=pt_EIT+912+6;        % drop the image
0324          case 8, %type=8, Waveform data
0325             aux_seg_len = d_struct(pt_EIT+3)/4*d_struct(pt_EIT+4);
0326            if (d_struct(pt_EIT+1) == 0)     % waveform channel 0 is the analog input
0327                aux_segment = d_EIT(pt_EIT+6 : pt_EIT+6+aux_seg_len-1);
0328                auxdata = [auxdata; aux_segment];
0329            else
0330                % drop all other waveform channels (see Waves/Waveform in header file for what they are)
0331            end  
0332            pt_EIT=pt_EIT+6+aux_seg_len;           
0333          case 10,% type = 10, unknown type
0334            %drop
0335            pt_EIT=pt_EIT+6+d_struct(pt_EIT+3)/4*d_struct(pt_EIT+4);
0336          otherwise,
0337            eidors_msg('WARNING: unknown type in carefusion file type');
0338            pt_EIT=pt_EIT+6+d_struct(pt_EIT+3)/4*d_struct(pt_EIT+4);
0339        end
0340    end
0341    vv=vv(47+[1:208],:);
0342    auxtime = (cumsum([1 13*ones(1,15)])-1)/208;             % relative time as fractions of the EIT frame: assuming uniform sampling - not sure this is 100% correct(!)
0343    auxtime = reshape(repmat(1:size(vv,2), length(auxtime), 1),[],1) + repmat(auxtime, 1, size(vv,2))';
0344    
0345 
0346 % Read the old <2006 Draeger "get" file format
0347 function [vv,curr,volt,auxdata] = draeger_get_readdata( fname );
0348    fid= fopen(fname,'rb');
0349    emptyctr=0;
0350    while emptyctr<2
0351      line= fgetl(fid);
0352      if isempty(line)
0353         emptyctr= emptyctr+1;
0354      else
0355         eidors_msg('data not understood',0);
0356         emptyctr=0;
0357      end
0358    end
0359    d= fread(fid,inf,'float',0,'ieee-le');
0360    fclose(fid);
0361 pause
0362 
0363    if rem( length(d), 256) ~=0
0364       eidors_msg('File length strange - cropping file',0);
0365       d=d(1:floor(length(d)/256)*256);
0366    end
0367 
0368    dd= reshape( d, 256, length(d)/256);
0369    vv= untwist(dd);
0370 
0371    curr=0.00512*dd(209:224,:);  % Amps
0372    volt=12*dd(225:240,:); %Vrms
0373    auxdata= dd(241:255,:);
0374    auxdata= auxdata(:);
0375     
0376 function jnk
0377 %[adler@adler01 sept07]$ xxd Sch_Pneumoperitoneum_01_001.get | head -30
0378 %0000000: 2d2d 2d44 7261 6567 6572 2045 4954 2d53  ---Draeger EIT-S
0379 %0000010: 6f66 7477 6172 652d 2d2d 0d0a 2d2d 2d50  oftware---..---P
0380 %0000020: 726f 746f 636f 6c20 4461 7461 2d2d 2d2d  rotocol Data----
0381 %0000030: 2d0d 0a0d 0a44 6174 653a 2020 3138 2d30  -....Date:  18-0
0382 %0000040: 322d 3230 3034 0d0a 5469 6d65 3a20 2031  2-2004..Time:  1
0383 %0000050: 333a 3138 2050 4d0d 0a0d 0a46 696c 656e  3:18 PM....Filen
0384 %0000060: 616d 653a 2020 2020 2020 2020 2053 6368  ame:         Sch
0385 %0000070: 5f50 6e65 756d 6f70 6572 6974 6f6e 6575  _Pneumoperitoneu
0386 %0000080: 6d5f 3031 5f30 3031 2e67 6574 0d0a 4453  m_01_001.get..DS
0387 %0000090: 5020 5365 7269 616c 204e 722e 3a20 2020  P Serial Nr.:
0388 %00000a0: 4549 5430 322f 3035 2f30 3030 360d 0a0d  EIT02/05/0006...
0389 %00000b0: 0a46 7265 7175 656e 6379 205b 487a 5d3a  .Frequency [Hz]:
0390 %00000c0: 2020 2020 3937 3635 362e 330d 0a47 6169      97656.3..Gai
0391 %00000d0: 6e3a 2020 2020 2020 2020 2020 2020 2020  n:
0392 %00000e0: 2020 2031 3134 350d 0a41 6d70 6c69 7475     1145..Amplitu
0393 %00000f0: 6465 3a20 2020 2020 2020 2020 2020 2031  de:            1
0394 %0000100: 3030 300d 0a53 616d 706c 6520 5261 7465  000..Sample Rate
0395 %0000110: 205b 6b48 7a5d 3a20 2020 2035 3030 300d   [kHz]:    5000.
0396 %0000120: 0a50 6572 696f 6473 3a20 2020 2020 2020  .Periods:
0397 %0000130: 2020 2020 2020 2020 2032 300d 0a46 7261           20..Fra
0398 %0000140: 6d65 733a 2020 2020 2020 2020 2020 2020  mes:
0399 %0000150: 2020 2020 2031 300d 0a0d 0a0d 0a
0400 %                                       8bc33e       10........>
0401 %0000160: 3f 0548633e bf20933d e192393d a568ea  ?.Hc>. .=..9=.h.
0402 %0000170: 3c 2530f63c 27e6043d 2043ad3c 25ce93  <%0.<'..= C.<%..
0403 %0000180: 3c aebcce3c 8714643d e3533d3e 65b6e1  <...<..d=.S=>e..
0404 %0000190: 3e 7a62103f 81c4143e 8c99813d 35921d  >zb.?...>...=5..
0405 %00001a0: 3d 8e6b0d3d 690cf93c 9910713c 3c9289  =.k.=i..<..q<<..
0406 %00001b0: 3c f6736f3c 2291453d ad1cab3d 386f15  <.so<".E=...=8o.
0407 %00001c0: 3e e82a143f 2a952e3f f568493e f08a8c  >.*.?*..?.hI>...
0408 %00001d0: 3d e43e0e3d 7040253d 19f4af3c 67fd93  =.>.=p@%=...<g..
0409 
0410 function vv = sheffield_readdata( fname );
0411    fid=fopen(fname,'rb','ieee-le');
0412    draw=fread(fid,[104,inf],'float32');
0413    fclose(fid);
0414    ldat = size(draw,2);
0415 
0416    [x,y]= meshgrid( 1:16, 1:16);
0417    idxm= y-x;
0418  % HW gain table
0419    gtbl = [0,849,213,87,45,28,21,19,21,28,45,87,213,849,0];
0420    idxm(idxm<=0)= 1;
0421    idxm= gtbl(idxm);
0422 
0423  % Multiply by gains
0424    draw = draw .* (idxm(idxm>0) * ones(1,ldat)); 
0425 
0426    vv= zeros(16*16, size(draw,2));
0427    vv(find(idxm),:) = draw;
0428 
0429  % Add reciprocal measurements
0430    vv= reshape(vv,[16,16,ldat]);
0431    vv= vv + permute( vv, [2,1,3]);
0432    vv= reshape(vv,[16*16,ldat]);
0433 
0434    
0435 
0436 function [vv,curr,volt,auxdata,auxtime,rawdata] = mceit_readdata( fname );
0437 
0438    fid= fopen(fname,'rb');
0439    d= fread(fid,inf,'float');
0440    fclose(fid);
0441 
0442    if rem( length(d), 256) ~=0
0443       eidors_msg('File length strange - cropping file',0);
0444       d=d(1:floor(length(d)/256)*256);
0445    end
0446 
0447    dd= reshape( d, 256, length(d)/256);
0448    rawdata = dd;
0449    no_reciprocity = (dd(39,1)==0);      %104 measurements per frame
0450    if no_reciprocity
0451        dd=transfer104_208(dd);
0452    end
0453    vv= untwist(dd);
0454 
0455    curr=0.00512*dd(209:224,:);  % Amps
0456    volt=12*dd(225:240,:); %Vrms
0457    auxdata= dd(241:256,:);
0458    auxdata= auxdata(:);
0459    %input impedance=voltage./current-440;        Ohm
0460    
0461    if no_reciprocity
0462      % remove every 14th, 15th and 16th sample as they are always zero
0463      auxdata(14:16:end) = []; auxdata(14:15:end) = []; auxdata(14:14:end) = [];
0464      % only 104 measurements were performed with NON-UNIFORM sampling of AUX in between EIT samples
0465      % Sampling procedure was thought to be: 13 AUX1 13 AUX2 12 AUX3 11 AUX4 10 AUX5 9 ... 2 AUX13 1 [END OF FRAME]
0466 %      auxtime = (cumsum([1 13 12:-1:2]) - 1) / 104;
0467      % However, this seems to be a little different, finally the sampling points were determined
0468      % empirically by measuring a sawtooth signal (linear ramp) and fitting the timings
0469      auxtime = [0.0000,0.1390,0.2535,0.3617,0.4642,0.5486,0.6367,0.7110,0.7780,0.8354,0.8865,0.9304,0.9711];  % mean fit at 25 Hz
0470 %      auxtime = [0.0000,0.1460,0.2789,0.3895,0.4875,0.5788,0.6608,0.7356,0.7986,0.8569,0.9045,0.9454,0.9824,]; % mean fit at 44 Hz
0471    else
0472      % all 208 measurements were performed
0473      auxtime = (cumsum([1 13*ones(1,15)])-1)/208;             % relative time as fractions of the EIT frame
0474    end
0475    % time of AUX signal relative to frame number (starting with 1 - Matlab style)
0476    auxtime = reshape(repmat(1:size(dd,2), length(auxtime), 1),[],1) + repmat(auxtime, 1, size(dd,2))';
0477    
0478 
0479 function array208=transfer104_208(array104),
0480 % The order in the 208-vector is
0481 % Index   Inject    Measure Pair
0482 % 1         1-2        3-4        (alternate pair is 3-4, 1-2, at location 39)
0483 % 2         1-2        4-5
0484 % 3         1-2        5-6
0485 % ¡­
0486 % 13        1-2        15-16
0487 % 14        2-3        4-5
0488 % 15        2-3        5-6
0489 % ¡­
0490 % 26        2-3        16-1
0491 % 27        3-4        5-6
0492 % ¡­
0493 % 39        3-4        1-2
0494 % 40        4-5        6-7
0495 % ¡­
0496 
0497 ind=[39,51,63,75,87,99,111,123,135,147,159,171,183,52,64,76, ...
0498      88,100,112,124,136,148,160,172,184,196,65,77,89,101,113, ...
0499      125,137,149,161,173,185,197,78,90,102,114,126,138,150,162, ...
0500      174,186,198,91,103,115,127,139,151,163,175,187,199,104,116, ...
0501      128,140,152,164,176,188,200,117,129,141,153,165,177,189, ...
0502      201,130,142,154,166,178,190,202,143,155,167,179,191,203, ...
0503      156,168,180,192,204,169,181,193,205,182,194,206,195,207,208];
0504 ro=[1:13, 14:26, 27:38, 40:50, 53:62, 66:74, 79:86, 92:98, ...
0505     105:110,118:122,131:134,144:146,157:158,170:170];
0506 
0507 [x,y]=size(array104);
0508 if x~=256 && y~=256,
0509     eidors_msg(['eidors_readdata: expectingin an input array ', ...
0510                 'of size 208*n']);
0511     return;
0512 elseif y==256,
0513     array104=array104';
0514     y=x;
0515 end
0516 array208=array104;
0517 for i=1:y,
0518     array208(ind,i)=array104(ro,i);    
0519 end
0520    
0521 
0522 % measurements rotate with stimulation, we untwist them
0523 function vv= untwist(dd);
0524    elec=16;
0525    pos_i= [0,1];
0526    ELS= rem(rem(0:elec^2-1,elec) - ...
0527         floor((0:elec^2-1)/elec)+elec,elec)';
0528    ELS=~any(rem( elec+[-1 0 [-1 0]+pos_i*[-1;1] ] ,elec)' ...
0529         *ones(1,elec^2)==ones(4,1)*ELS')';
0530    twist= [               0+(1:13), ...
0531                          13+(1:13), ...
0532            39-(0:-1:0),  26+(1:12), ...
0533            52-(1:-1:0),  39+(1:11), ...
0534            65-(2:-1:0),  52+(1:10), ...
0535            78-(3:-1:0),  65+(1: 9), ...
0536            91-(4:-1:0),  78+(1: 8), ...
0537           104-(5:-1:0),  91+(1: 7), ...
0538           117-(6:-1:0), 104+(1: 6), ...
0539           130-(7:-1:0), 117+(1: 5), ...
0540           143-(8:-1:0), 130+(1: 4), ...
0541           156-(9:-1:0), 143+(1: 3), ...
0542           169-(10:-1:0),156+(1: 2), ...
0543           182-(11:-1:0),169+(1: 1), ...
0544           195-(12:-1:0), ...
0545           208-(12:-1:0) ];
0546     vv= zeros(256,size(dd,2));
0547     vv(ELS,:)= dd(twist,:);
0548    %vv= dd(1:208,:);
0549 
0550 % Read data from p2k files (I T S system)
0551 % FIXME: this code is very rough, it works for
0552 %   only eight ring data records
0553 function  vv = its_readdata( fname ) 
0554    fid= fopen( fname, 'rb', 'ieee-le');
0555    vv=[];
0556 
0557    % don't know how to interpret header
0558    header= fread(fid, 880, 'uchar');
0559    frameno= 0;
0560    rings= 8;
0561    while( ~feof(fid) )
0562        frameno= frameno+1;
0563        % don't know how to interpret frame header
0564        framehdr= fread(fid, 40);
0565        data= fread(fid, 104*rings, 'double');
0566        vv= [vv, data];
0567    end
0568 
0569    if 0 % convert a ring to 208
0570       ringno= 1;
0571       ld= size(vv,2);
0572       vx= [zeros(1,ld);vv( ringno*104 + (-103:0) ,: )];
0573       idx= ptr104_208;
0574       vv= vx(idx+1,:);
0575    end
0576 
0577 % pointer to convert from 104 to 208 meas patterns
0578 function idx= ptr104_208;
0579     idx= zeros(16);
0580     idx(1,:)= [0,0,1:13,0];
0581     ofs= 13;
0582     for i= 2:14
0583         mm= 15-i;
0584         idx(i,:) = [zeros(1,i+1), (1:mm)+ofs ];
0585         ofs= ofs+ mm;
0586     end
0587 
0588     idx= idx + idx';
0589     
0590 function vv = iirc_readdata( fname );
0591     fid= fopen( fname, 'r');
0592     while ~feof(fid)
0593        line = fgetl(fid);
0594        if isempty(line)
0595            continue;
0596        end
0597        
0598        num= regexp(line,'Channel : (\d+)');
0599        if ~isempty(num)
0600            channels= str2num( line(num(1):end ) );
0601            continue;
0602        end
0603        
0604        num= regexp(line,'Frequency : (\d+)kHz');
0605        if ~isempty(num)
0606            freqency= str2num( line(num(1):end ) );
0607            continue;
0608        end
0609 
0610        num= regexp(line,'Scan Method : (\w+)');
0611        if ~isempty(num)
0612            scan_method=  line(num(1):end );
0613            continue;
0614        end
0615 
0616        num= regexp(line,'Study : (\w+)');
0617        if ~isempty(num)
0618            study=  line(num(1):end);
0619            continue;
0620        end
0621            
0622        if strcmp(line,'Data');
0623            data= fscanf(fid,'%f',[4,inf])';
0624            continue;
0625        end
0626     end
0627     vv= data(:,1) + 1i*data(:,2);
0628     if length(vv) ~= channels^2
0629         error('eidors_readdata: data length wrong')
0630     end
0631 
0632 % stimulation is the fwd_model stimulation data structure
0633 % meas_select indicates if data is NOT measures on current electrode
0634 function [stimulations,meas_select] = UCT_sequence_file( fname );
0635    % (c) Tim Long
0636    % 21 January 2005
0637    % University of Cape Town
0638 
0639 
0640    % open the file
0641    fid = fopen(fname, 'rt');
0642 
0643    % check to see if file opened ok
0644    if fid == -1
0645          errordlg('File not found','ERROR')  
0646          return;
0647    end
0648 
0649 
0650    tline = fgetl(fid);             % get the spacer at top of text file
0651 
0652    % the measurement and injection pattern is stored as follows:
0653    % I1V1:  db #$00,#$0F,#$00,#$00,#$10,#$00,#$00,#$21,#$00 ...
0654    % I2V2:  db #$11,#$0F,#$11,#$11,#$10,#$11,#$11,#$21,#$11 ...
0655    % etc
0656    % need to put all the bytes in a vector
0657 
0658    % tokenlist will store the list of bytes as strings
0659    tokenlist = [];
0660 
0661 
0662    tline = fgetl(fid);             % get first line of data
0663 
0664    while length(tline) ~= 0
0665        
0666        % the first few characters in the line are junk
0667        rem = tline(11:end);            % extract only useful data
0668        
0669        % extract each byte
0670        while length(rem) ~=0
0671            [token, rem] = strtok(rem, ',');
0672            tokenlist = [tokenlist; token];
0673        end
0674        
0675        % get the next line in sequence file
0676        tline = fgetl(fid);
0677    end
0678 
0679    fclose(fid);
0680 
0681    % got everything in string form... need to covert to number format
0682 
0683    drive_lay = [];
0684    drive_elec = [];
0685    sense_lay = [];
0686 
0687    injection_no = 1;
0688    % for each injection
0689    for i=1:3:length(tokenlist)
0690        
0691        % get injection layer
0692        tsource_layer = tokenlist(i,3);
0693        tsink_layer = tokenlist(i,4);
0694        source_layer = sscanf(tsource_layer, '%x');
0695        sink_layer = sscanf(tsink_layer, '%x');
0696        
0697        drive_lay = [drive_lay; [source_layer sink_layer]];
0698          
0699        
0700        % get drive pair
0701        tsource_elec = tokenlist(i+1,3);
0702        tsink_elec = tokenlist(i+1,4);
0703        source_elec = sscanf(tsource_elec, '%x');
0704        sink_elec = sscanf(tsink_elec, '%x');
0705        
0706        drive_elec = [drive_elec; [source_elec sink_elec]];
0707        
0708        
0709        % get sense layer pair
0710        tpos_layer = tokenlist(i+2,3);
0711        tneg_layer = tokenlist(i+2,4);
0712        pos_sense_layer = sscanf(tpos_layer, '%x');
0713        neg_sense_layer = sscanf(tneg_layer, '%x');
0714        
0715        sense_lay = [sense_lay; [pos_sense_layer neg_sense_layer]];
0716    end
0717 
0718    n_elec = size(sense_lay,1);
0719    n_inj  = size(drive_lay,1);       % every injection
0720    elecs_per_plane = 16; % FIXED FOR THE UCT DEVICE
0721    raw_index = 0;
0722    meas_select = [];
0723 
0724    for i=1:n_inj      % for every injection
0725        stimulations(i).stimulation= 'mA';
0726        
0727        % find the injection electrodes
0728        e_inj_p = drive_lay(i, 1) * elecs_per_plane + drive_elec(i,1) + 1;
0729        e_inj_n = drive_lay(i, 2) * elecs_per_plane + drive_elec(i,2) + 1;
0730 
0731        % create the stimulation pattern for this injection
0732        inj = zeros(n_elec, 1);
0733        inj(e_inj_p) = 1;
0734        inj(e_inj_n) = -1;
0735        stimulations(i).stim_pattern = inj;
0736      
0737        % the UCT instrument always makes 16 measurements per injection.
0738        % the +ve and -ve electrodes are always adjacent, but might be on
0739        % different planes
0740        meas_pat = [];
0741        for e = 0:15
0742            raw_index = raw_index + 1;
0743            meas = zeros(1, n_elec);   % the measurement electrodes for this sample
0744 
0745            % find the measurement electrodes for this measurement (+ve elec is
0746            % next to -ve electrode)
0747            e_meas_p = sense_lay(i,1) * elecs_per_plane + mod(e+1,elecs_per_plane) + 1;
0748            e_meas_n = sense_lay(i,2) * elecs_per_plane + e + 1;
0749 
0750            % if either of the drive electrodes are equal to any of the sense
0751            % electrodes, we must not include this sample
0752            if any( e_meas_p == [e_inj_p, e_inj_n] ) | ...
0753               any( e_meas_n == [e_inj_p, e_inj_n] ) 
0754                continue;
0755            end
0756 
0757            meas(e_meas_p) = -1;
0758            meas(e_meas_n) = 1;
0759            meas_select = [meas_select;raw_index];
0760            
0761            % add this measurement to the measurement pattern
0762            meas_pat = [meas_pat; meas];
0763 
0764        end     % for each injection there are actually only 13-16 measurements
0765        stimulations(i).meas_pattern = sparse(meas_pat);
0766    end
0767 
0768 function [cur_data,no_cur_data] = UCT_calibration_file( fname );
0769    fid = fopen(fname, 'rb');
0770    mag_num = fread(fid, 1, 'int32');
0771    version = fread(fid, 1, 'int32');
0772 %  comments = fread(fid, 2048, 'char'); MATLAB CHANGED CHAR to 2 bytes
0773    comments = fread(fid, 2048, 'uint8');
0774    comments = setstr(comments');
0775    no_of_layers = fread(fid, 1, 'int32');
0776    uppa_dac = fread(fid, 1, 'float64');
0777    lowa_dac = fread(fid, 1, 'float64');
0778    raw_frame_size = fread(fid, 1, 'int32');
0779 
0780    no_cur_data = [];
0781    cur_data = [];
0782 
0783    for i=1:no_of_layers
0784        no_cur_data = [no_cur_data;fread(fid, raw_frame_size, 'float64')];
0785        cur_data = [cur_data;fread(fid, raw_frame_size, 'float64')];
0786    end
0787    fclose(fid);
0788 
0789 %  no_cur_data = UCT_ShuffleData( no_cur_data);
0790 %  cur_data =    UCT_ShuffleData( cur_data);
0791 
0792 function [v_raw] = UCT_LoadDataFrame(infilename)
0793 % [v_raw] = LoadDataFrame(infilename)
0794 %
0795 % Loads the data from the tomography file
0796 %
0797 % (c) Tim Long
0798 % 21 January 2005
0799 % University of Cape Town
0800 
0801 
0802 % Open the file
0803 fid = fopen(infilename, 'rb');
0804 
0805 if fid == -1
0806       errordlg('File not found','ERROR')  
0807       return;
0808 end
0809 
0810 %%% new file changes
0811 magic_number = fread(fid, 1, 'uint32');
0812 
0813 if magic_number ~= 2290649224
0814    disp('UCT File: wrong file type'); 
0815 end
0816 version = fread(fid, 1, 'uint32');
0817 foffset = fread(fid, 1, 'uint32');
0818 no_of_layers = fread(fid, 1, 'uint32');
0819 frame_size = fread(fid, 1, 'uint32');
0820 fseek(fid, foffset-8, 'cof');
0821 
0822 %%% end of new file changes
0823 %%% old file stuff
0824 % no_of_layers = fread(fid, 1, 'uint32');
0825 % frame_size = fread(fid, 1, 'uint32');
0826 
0827 frame_no = fread(fid, 1, 'uint32');
0828 
0829 v_raw = [];
0830 while feof(fid)==0
0831     v_raw = [v_raw, fread(fid, frame_size*no_of_layers, 'float64')];
0832     frame_no = fread(fid, 1, 'uint32');
0833 end
0834 
0835 fclose(fid);
0836 % UCT_ShuffleData??
0837 
0838 % Read data from the file format develped by Pascal
0839 % Gaggero at CSEM Landquart in Switzerland.
0840 function [vv] = landquart1_readdata( fname );
0841    FrameSize = 1024;
0842 
0843    enableCounter=1;
0844    counter=0;
0845 
0846    fid=fopen(fname);
0847 
0848    codec=fread(fid,1,'int'); %read codec
0849    if codec~=1; error('Codec unexpected value'); end
0850 
0851    nbFileInHeader=fread(fid,1,'int'); %read codec
0852    
0853    for i=1:nbFileInHeader
0854        lenghtFile = fread(fid,1,'int64'); 
0855        jnk= fread(fid,lenghtFile,'int8');% discard the file
0856    end
0857    
0858    
0859    vv= []; 
0860    while 1; 
0861        type=fread(fid,1,'int'); %check if not End of file
0862        if type == -1; break ; end
0863        
0864        nel= fread(fid,1,'int');
0865        sz= fread(fid,1,'int');
0866        iq=fread(fid,[2,sz/8],'int')';
0867 
0868        vv = [vv, iq*[1;1j]]; %I+ 1j*Q vectors
0869        
0870        for j=1:nel-1
0871            sz= fread(fid,1,'int');
0872            jnk = fread(fid,sz/4,'int');
0873        end
0874        
0875    end
0876    
0877 % Read data from the file format develped by Swisstom, Landquart, Switzerland.
0878 function [vv, elecImps, tStampsAbs, tStampsRel] = landquart2_readdata( fname )
0879    if ~exist('OCTAVE_VERSION')
0880    [fid msg]= fopen(fname,'r','ieee-be','UTF-8');
0881    else
0882    [fid msg]= fopen(fname,'r','ieee-be');
0883    end
0884    try
0885       format_version = fread(fid,1,'int32','ieee-be');
0886       if format_version ~= 3
0887          error('unsupported file format version');
0888       else
0889          % get expected number of frames and preallocate variables
0890          fseek(fid,16,'cof');
0891          nFrames = fread(fid, 1, 'int32', 'ieee-be');
0892          tStampsAbs = nan(1, nFrames);
0893          tStampsRel = nan(1, nFrames);
0894          viPayload = nan(64, nFrames);
0895          iqPayload = nan(2048, nFrames);
0896          
0897          header_size = 2264; 
0898          fseek(fid,header_size + 8,'bof');
0899          
0900          %%% get frame size and length of payload at end of frame
0901          frame_length = fread(fid, 1, 'int32', 'ieee-be') + 12;
0902          %%% move back to start of 1. frame
0903          fseek(fid,header_size,'bof');
0904          
0905          %%% Read frames
0906          i = 1;
0907          while fseek(fid, 1,'cof') ~= -1
0908             fseek(fid, -1,'cof');
0909             % read absolute timestamp (milliseconds resolution)
0910             tStampsAbs(i) = fread(fid,1,'int64','ieee-be');
0911             pl = fread(fid,1,'int32','ieee-le');    % payload length (i.e. frame length)
0912             frame_header = fread(fid,15,'int32','ieee-le'); 
0913             % read relative timestamp (microseconds resolution)
0914             tStampsRel(i) = frame_header(5);    
0915             
0916             % drop some unintersting parts
0917             fseek(fid,12, 'cof');
0918             
0919             % get electrode impedance measurement
0920             viPayload(:,i) = fread(fid,64,'int32','ieee-le');
0921             % get effective payload (voltage measurements)
0922             iqPayload(:,i) = fread(fid,2048,'int32','ieee-le');
0923             fseek(fid,header_size + i*frame_length,'bof');
0924             i = i +1;
0925          end
0926 
0927       end
0928    catch err
0929       fclose(fid);
0930       rethrow(err);
0931    end
0932    fclose(fid);
0933    
0934    i = i-1;
0935    if i ~= nFrames       
0936       % remove data which were preallocated but not read
0937       if i < nFrames       
0938          tStampsAbs(i+1:end) = [];
0939          tStampsRel(i+1:end) = [];
0940          viPayload(:,i+1:end) = [];
0941          iqPayload(:,i+1:end) = [];
0942       end
0943       eidors_msg('"%s": expected %.0f frames but read %.0f',fname, nFrames, i ,3);
0944    end
0945    
0946    % convert into matlab date format, e.g. read via datestr(tStampAbs(1))
0947    tStampsAbs = tStampsAbs/(24*3600*1E3) + datenum(1970, 1, 1, 1, 0, 0);
0948    % strictly speaking we would need to correct for DST, but as this is only
0949    % supported by ML R2014b or higher we avoid it to be backwards compatible
0950 %    tStampsAbs = tStampsAbs + isdst(datetime(datevec(tStampsAbs(1)), 'TimeZone', 'local'))/24;   % add one hour if DST
0951    
0952    % this is just a simple guess
0953    amplitudeFactor = 2.048 / (2^20 * 360 * 1000);
0954    vv = amplitudeFactor * (iqPayload(1:2:end,:) + 1i*iqPayload(2:2:end,:));
0955    
0956    elecImps = viPayload(1:2:end,:) + 1i*viPayload(2:2:end,:);
0957 
0958 % Read data from the file format develped by Swisstom, Landquart, Switzerland.
0959 function [vv, evtlist, elecImps, tStampsAbs, tStampsRel] = landquart4_readdata( fname )
0960    evtlist = [];
0961    if ~exist('OCTAVE_VERSION')
0962    [fid msg]= fopen(fname,'r','ieee-le','UTF-8');
0963    else
0964    [fid msg]= fopen(fname,'r','ieee-le');
0965    end
0966    try
0967       format_version = fread(fid,1,'int32','ieee-le');
0968       if format_version ~= 4
0969          error('unsupported file format version');
0970       else          
0971          header_size = fread(fid,1,'int32', 'ieee-le');   
0972          eit_frame_offset = 328; % 60 + 12 + 256 = 328
0973          iq_payload = 2048;
0974          vi_payload = 64;
0975          
0976          % get expected number of frames and preallocate variables
0977          fseek(fid,16,'cof');
0978          nFrames = fread(fid, 1, 'int32', 'ieee-le');
0979          tStampsAbs = nan(1, nFrames);
0980          tStampsRel = nan(1, nFrames);
0981          viPayload = nan(vi_payload, nFrames);
0982          iqPayload = nan(iq_payload, nFrames);
0983 
0984          fseek(fid,header_size,'bof');
0985          
0986          %%% Read frames
0987          i = 1;
0988          evti = 1;
0989          while fseek(fid, 1,'cof') ~= -1
0990             fseek(fid, -1,'cof');
0991             % drop frame header and some payload:
0992             % read absolute timestamp (milliseconds resolution)
0993             tStampsAbs(i) = fread(fid,1,'int64','ieee-le'); 
0994             ft = fread(fid,1,'int32','ieee-le'); 
0995             pl = fread(fid,1,'int32','ieee-le');
0996             if ft == 1
0997                % event
0998                evtlist(evti).timestamp = tStampsAbs(i) ;
0999                evtlist(evti).frame = i ;
1000                evti = evti + 1;
1001                if pl > 0
1002                   evtlist(evti).eventId = fread(fid, 1, 'int32', 'ieee-le');
1003                 end
1004             elseif ft == 0
1005                 frame_header = fread(fid,15,'int32','ieee-le'); 
1006                 % read relative timestamp (microseconds resolution)
1007                 tStampsRel(i) = frame_header(5);   
1008                 
1009                 % drop some unintersting parts
1010                 fseek(fid,12, 'cof');
1011                
1012                 % get electrode impedance measurement
1013                 viPayload(:,i) = fread(fid,vi_payload,'int32','ieee-le');
1014                 % get effective payload (voltage measurements)
1015                 iqPayload(:,i) = fread(fid,iq_payload,'int32','ieee-le');
1016                 fseek(fid,pl-4*iq_payload-eit_frame_offset,'cof');
1017                 i = i+1;
1018             elseif pl > 0
1019                fseek(fid,pl,'cof'); 
1020             else
1021                % nothing to do
1022             end
1023          end         
1024       end
1025    catch err
1026       fclose(fid);
1027       rethrow(err);
1028    end
1029    fclose(fid);
1030    
1031    i = i-1;
1032    if i ~= nFrames       
1033       % remove data which were preallocated but not read
1034       if i < nFrames       
1035          tStampsAbs(i+1:end) = [];
1036          tStampsRel(i+1:end) = [];
1037          viPayload(:,i+1:end) = [];
1038          iqPayload(:,i+1:end) = [];
1039       end
1040       eidors_msg('"%s": expected %.0f frames but read %.0f',fname, nFrames, i ,3);
1041    end
1042    
1043    % convert into matlab date format, e.g. read via datestr(tStampAbs(1))
1044    tStampsAbs = tStampsAbs/(24*3600*1E3) + datenum(1, 1, 1, 0, 0, 0);
1045    % strictly speaking we would need to correct for DST, but as this is only
1046    % supported by ML R2014b or higher we avoid it to be backwards compatible
1047 %    tStampsAbs = tStampsAbs + isdst(datetime(datevec(tStampsAbs(1)), 'TimeZone', 'local'))/24;   % add one hour if DST
1048 
1049    % this is just a simple guess
1050    amplitudeFactor = 2.048 / (2^20 * 360 * 1000);
1051    vv = amplitudeFactor * (iqPayload(1:2:end,:) + 1i*iqPayload(2:2:end,:));
1052    
1053    %% elecImps depends on a factor which varies with the SBC version,
1054    %%   however the following value is close for most versions
1055    impedanceFactor =  2.048 / (2^12 * 0.173 * 0.003) / 2^15; % = 0.9633 / 2^15;
1056    elecImps = impedanceFactor * (viPayload(1:2:end,:) + 1i*viPayload(2:2:end,:));
1057 
1058 function [vv] = landquart4pre_readdata( fname )
1059    if ~exist('OCTAVE_VERSION')
1060    [fid msg]= fopen(fname,'r','ieee-le','UTF-8');
1061    else
1062    [fid msg]= fopen(fname,'r','ieee-le');
1063    end
1064    try
1065       format_version = fread(fid,1,'int32','ieee-le');
1066       if format_version ~= 4
1067          error('unsupported file format version');
1068       else
1069          header_size = fread(fid,1,'int32', 'ieee-le'); 
1070          
1071          fseek(fid,header_size + 12,'bof');
1072          %%% get frame size and length of payload at end of frame
1073          frame_length = fread(fid, 1, 'int32', 'ieee-le') + 16;
1074          %%% move back to start of 1. frame
1075          fseek(fid,header_size,'bof');
1076          
1077          %%% Read frames
1078          i = 1;
1079          while fseek(fid, 1,'cof') ~= -1
1080             fseek(fid, -1,'cof');
1081             % drop frame header and some payload:
1082             % (16 + 60) + 12 + 256 = 344
1083             fseek(fid,344, 'cof');
1084             iqPayload(:,i) = fread(fid,2048,'int32','ieee-le');
1085             fseek(fid,header_size + i*frame_length,'bof');
1086             i = i +1;
1087          end
1088 
1089       end
1090    catch err
1091       fclose(fid);
1092       rethrow(err);
1093    end
1094    fclose(fid);
1095 
1096    % this is just a simple guess
1097    amplitudeFactor = 2.048 / (2^20 * 360 * 1000);
1098    vv = amplitudeFactor * (iqPayload(1:2:end,:) + 1i*iqPayload(2:2:end,:));
1099    
1100 %  Output: [encodepage] = eidors_readdata( path,'DIX_encode');
1101 %   where path= '/path/to/Criptografa_New.dll' (provided with system)
1102 function [vv] = dixtal_read_codepage( fname );
1103    %Fname must be the Criptografa_New dll
1104    fid = fopen(fname,'rb');
1105    b1234= fread(fid, [1,4], 'uint8');
1106    if b1234~= [77, 90, 144, 0];
1107       error('This does not appear to be the correct Dll');
1108    end
1109    fseek(fid, hex2dec('e00'), 'bof');
1110    encodepage1 = fread(fid, hex2dec('1000'), 'uint8');
1111    fclose(fid);
1112    encodepage1 = flipud(reshape(encodepage1,4,[]));
1113    vv = encodepage1(:);
1114 
1115 %        format = 'DIXTAL'
1116 %           - output: [vv] = eidors_readdata( fname, 'DIXTAL', encodepage );
1117 function [vv] = dixtal_read_data( file_name, frame_range, encodepage1 );
1118 
1119    % Get encodepage from the file
1120    fid=fopen(file_name,'rb');
1121    n1 = fgetl(fid);
1122    n2 = fgetl(fid);
1123    n3 = fgetl(fid); 
1124    nframes = str2num( n3 );
1125    
1126    fseek(fid, hex2dec('800'), 'bof');
1127    encodepage2 = fread(fid, hex2dec('1000'), 'uint8');
1128    fclose(fid);
1129 
1130    encodepage = bitxor(encodepage1, encodepage2);
1131 
1132    % Read the file
1133    dplen = hex2dec('1090');
1134    start = hex2dec('2800');
1135    fid=fopen(file_name,'rb');
1136    if isempty(frame_range)
1137       frame_range = 0:nframes;
1138    else
1139       frame_range = frame_range - 1;
1140       frame_range(frame_range>nframes) = [];
1141    end      
1142    vv= zeros(dplen/4, length(frame_range));
1143    k= 1;
1144    for i = frame_range
1145       status= fseek(fid, start + i*dplen, 'bof');
1146       if status~=0; break; end
1147       datapage = fread(fid, dplen, 'uint8');
1148       if length(datapage)== 0; 
1149          vv= vv(:,1:k-1); break; % frame number inside file is wrong
1150       end
1151       vv(:,k) = proc_dixtal_data( datapage, encodepage );
1152       k=k+1;
1153    end
1154    fclose(fid);
1155 
1156 
1157 function  data = proc_dixtal_data( b, encodepage );
1158 
1159    cryptolen = 1024*4;
1160    b(1:cryptolen) = bitxor(b(1:cryptolen), encodepage);
1161 
1162    b1 = b(1:4:end,:); %b1 = bitand(b1,127);
1163    b2 = b(2:4:end,:);
1164    b3 = b(3:4:end,:);
1165    b4 = b(4:4:end,:);
1166    bb  = [b1,b2,b3,b4]*(256.^[3;2;1;0]);
1167 
1168    sgnbit = bitand( bb,  2^31);
1169    sgnbit = +1*(sgnbit == 0) - 1*(sgnbit == 2^31);
1170 
1171    expbit = bitand( bb, 255*2^23 ) /2^23; 
1172    expbit = 2.^(expbit - 127);
1173 
1174    fracbt = bitand( bb, 2^23-1)/2^23 + 1;
1175    data = sgnbit .* expbit .* fracbt;
1176 
1177 function [fr] = read_draeger_header( filename );
1178 %READ_DRAEGER_HEADER   Opens and reads the header portion of the Draeger .eit
1179 %file. Current parameter returned is frame rate.
1180 %
1181 % function [fr,s] = ReadDraegerHeader(filename)
1182 % fr        double      scalar      frame rate in hertz
1183 
1184 % Determine file version
1185 K0 = 2048;   % bytes to read in char class
1186 K1='Framerate [Hz]:'; % Draeger frame rate line in header
1187 K2=15;                % Length of K1 string
1188 K3=13;                % Following F1 Field for delimiter (look for ^M)
1189 
1190 % Open file for reading in little-endian form
1191 fid = fopen(filename,'r','l');
1192 if fid == -1
1193     error('Error read_draeger_header: can''t open file');
1194 end
1195 header = fread(fid,K0,'*char')';
1196 index = strfind(header,K1);
1197 if ~isempty(index)
1198     [tok,rem]= strtok(header(index+K2:end),K3);
1199     fr = str2num(tok); 
1200 else
1201     error('Error read_draeger_header: frame rate unspecified');
1202 end
1203 
1204 if isempty(fr)
1205     error('Error read_draeger_header: Frame rate could not be read');
1206 end
1207 
1208 fclose(fid);
1209 
1210 
1211 function vd = read_draeger_file(filename)
1212 %READDRAEGERFILE   Opens and reads a Draeger .eit file. Returns an array
1213 %containing the voltage data arranged per frame.
1214 %
1215 % Input:
1216 % filename  char        1xN
1217 %
1218 % Output:
1219 % vd        double      MxN         M = volt measurements per frame (mV)
1220 %                                   N = number of frames
1221 
1222 
1223    ff = dir(filename);
1224    filelen = ff.bytes;
1225 
1226    % Open file for reading in little-endian form
1227    fid = fopen(filename,'r','l');
1228    if fid == -1
1229        error('Error read_draeger_file: file could not be opened');
1230    end
1231 
1232    % Determine file version
1233    K0 = 128;   % bytes to read in char class
1234 
1235    header = fread(fid,K0,'*char')';
1236 if 0 % THis analysis doesn't seem useful
1237    if ~isempty(strfind(header,'V3.2'))
1238        version = '3.2';
1239    elseif ~isempty(strfind(header,'V4.01'))
1240        version = '4.01'; 
1241    elseif ~isempty(strfind(header,'V1.10')) % 2014!!
1242        version = '1.10'; 
1243    else
1244        error('Error read_draeger_file: unknown file version');
1245    end
1246 end
1247 
1248    % Read Format version
1249    fseek(fid,0,-1); % beginning of file
1250    hdr = fread(fid, 8, 'uint8');
1251    offset1 =  (256.^(0:3))*hdr(5:8) + 16; % "good" data starts at offset #2
1252    switch sprintf('%02X-',header(1:4));
1253       case '1F-00-00-00-';
1254          type='Draeger v31';  SPC = 4112;
1255       case '20-00-00-00-';
1256          type='Draeger v32';  SPC = 5200;
1257       case '33-00-00-00-';
1258          type='Draeger v51';  SPC = 5495;
1259       otherwise;
1260          error('File "%s" is format version %d, which we don''t know how to read', ...
1261                hdr(1))
1262    end
1263 
1264    len = floor( (filelen-offset1)/SPC );
1265    vd= zeros(600,len);
1266    ss= offset1 + (0:len)*SPC;
1267 
1268    for k= 1:length(ss);
1269       if fseek(fid,ss(k),-1)<0; break; end
1270       vd(:,k)=fread(fid,600,'double', 0, 'ieee-le');
1271    end
1272 
1273    % End of function
1274    fclose(fid);
1275