calc_reciproc_error

 CALC_RECIPROC_ERROR: CALCULATE RECIPROCITY ERROR MATRIX
 meas_icov = calc_reciproc_error(inv_model, data )

 Calculate meas_icov from data
 
 specify tau as 
    inv_model.calc_reciproc_error.tau
 as tau increases, effect is less. There is almost
 no effect when tau > 1e-4

 also accepts a fwd_model parameter

 Reference: Real-time management of faulty electrodes in
  electrical impedance tomography AE Hartinger, R Guardo,
  A Adler, H Gagnon. IEEE T BME 2008.

0001 function meas_icov = calc_reciproc_error(inv_model, data )
0002 % CALC_RECIPROC_ERROR: CALCULATE RECIPROCITY ERROR MATRIX
0003 % meas_icov = calc_reciproc_error(inv_model, data )
0004 %
0005 % Calculate meas_icov from data
0006 %
0007 % specify tau as
0008 %    inv_model.calc_reciproc_error.tau
0009 % as tau increases, effect is less. There is almost
0010 % no effect when tau > 1e-4
0011 %
0012 % also accepts a fwd_model parameter
0013 %
0014 % Reference: Real-time management of faulty electrodes in
0015 %  electrical impedance tomography AE Hartinger, R Guardo,
0016 %  A Adler, H Gagnon. IEEE T BME 2008.
0017 
0018 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0019 % $Id: calc_reciproc_error.m 3040 2012-06-06 15:20:40Z aadler $
0020 
0021 try;    tau = inv_model.calc_reciproc_error.tau;
0022 catch;  tau = 2.5e-6;                            end
0023 
0024 switch inv_model.type
0025   case 'inv_model'; fmdl = inv_model.fwd_model;
0026   case 'fwd_model'; fmdl = inv_model;
0027   otherwise;        error('calc_reciproc_error: require model input');
0028 end
0029 
0030 data = calc_difference_data( 0, data, fmdl);
0031 
0032 idx = reciprocity_idx( fmdl );
0033 if any(isnan(idx));
0034    error('not all meas have reciprocity for this stim_pattern');
0035 end
0036 
0037 data = data/max(abs(data));
0038 recerr=  data - data(idx,:);
0039 recerr= mean( abs(recerr), 2);
0040 s2  = exp(-recerr.^2/tau); % eqn 8 from paper
0041 
0042 nmeas = length(idx);
0043 meas_icov= spdiags(s2(:),0,nmeas,nmeas);
0044 
0045 function oldcode
0046 % Use only with no_rotate
0047 Nel = length(inv_model.fwd_model.electrode);
0048 nframes= size(data,2);
0049 nmeas = size(data,1);
0050 if nmeas ~= Nel^2
0051    data2= zeros(Nel^2, nframes);
0052    mselect = inv_model.fwd_model.meas_select;
0053    data2(mselect,:)= data;
0054    data2= reshape( data2, Nel,Nel,nframes);
0055 else 
0056    data2= reshape( data, Nel,Nel,nframes);
0057 end
0058 
0059 ndata2= abs(data2)/max(abs(data2(:)));
0060 ndata2r= permute(ndata2, [2,1,3]); % Transpose
0061 recerr= mean( ndata2 - ndata2r, 3);
0062 e2    = recerr.^2; % eqn 8 from paper
0063 
0064 s2  = exp(-e2/tau);
0065 
0066 if isfield(inv_model.fwd_model,'meas_select')
0067    mselect = inv_model.fwd_model.meas_select;
0068    s2= s2(mselect);
0069    nmeas = length(find(mselect));
0070 end
0071 
0072 meas_icov= spdiags(s2(:),0,nmeas,nmeas);