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