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

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 
0010 % Reference: Real-time management of faulty electrodes in
0011 %  electrical impedance tomography AE Hartinger, R Guardo,
0012 %  A Adler, H Gagnon. IEEE T BME 2008.
0013 
0014 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0015 % $Id: calc_reciproc_error.html 2819 2011-09-07 16:43:11Z aadler $
0016 
0017 try;    tau = inv_model.calc_reciproc_error.tau;
0018 catch;  tau = 2.5e-6;                            end
0019 
0020 fmdl = inv_model.fwd_model;
0021 
0022 data = calc_difference_data( 0, data, fmdl);
0023 
0024 idx = reciprocity_idx( fmdl );
0025 if any(isnan(idx));
0026    error('not all meas have reciprocity for this stim_pattern');
0027 end
0028 
0029 data = data/max(abs(data));
0030 recerr=  data - data(idx,:);
0031 recerr= mean( abs(recerr), 2);
0032 s2  = exp(-recerr.^2/tau); % eqn 8 from paper
0033 
0034 nmeas = length(idx);
0035 meas_icov= spdiags(s2(:),0,nmeas,nmeas);
0036 
0037 function oldcode
0038 % Use only with no_rotate
0039 Nel = length(inv_model.fwd_model.electrode);
0040 nframes= size(data,2);
0041 nmeas = size(data,1);
0042 if nmeas ~= Nel^2
0043    data2= zeros(Nel^2, nframes);
0044    mselect = inv_model.fwd_model.meas_select;
0045    data2(mselect,:)= data;
0046    data2= reshape( data2, Nel,Nel,nframes);
0047 else 
0048    data2= reshape( data, Nel,Nel,nframes);
0049 end
0050 
0051 ndata2= abs(data2)/max(abs(data2(:)));
0052 ndata2r= permute(ndata2, [2,1,3]); % Transpose
0053 recerr= mean( ndata2 - ndata2r, 3);
0054 e2    = recerr.^2; % eqn 8 from paper
0055 
0056 s2  = exp(-e2/tau);
0057 
0058 if isfield(inv_model.fwd_model,'meas_select')
0059    mselect = inv_model.fwd_model.meas_select;
0060    s2= s2(mselect);
0061    nmeas = length(find(mselect));
0062 end
0063 
0064 meas_icov= spdiags(s2(:),0,nmeas,nmeas);