calc_meas_icov

 meas_icov = calc_meas_icov( inv_model )
 CALC_MEAS_ICOV: calculate inverse covariance of measurements
   The meas_icov is a matrix n_meas x n_meas of the
     inverse of measurement covariances. Normally measurements
     are assumed to be independant, so the meas_icov is 
     a diagonal matrix of 1/var(meas)

 calc_meas_icov can be called as
    meas_icov= calc_meas_icov( inv_model )

 meas_icov   is the calculated data prior
 inv_model    is an inv_model structure

 if:
    inv_model.meas_icov    is a function
          -> call it to calculate meas_icov
    inv_model.meas_icov    is a matrix
          -> return it as meas_icov
    inv_model.meas_icov    does not exist
          -> use I, or 1./homg (for normalized difference)

0001 function meas_icov = calc_meas_icov( inv_model )
0002 % meas_icov = calc_meas_icov( inv_model )
0003 % CALC_MEAS_ICOV: calculate inverse covariance of measurements
0004 %   The meas_icov is a matrix n_meas x n_meas of the
0005 %     inverse of measurement covariances. Normally measurements
0006 %     are assumed to be independant, so the meas_icov is
0007 %     a diagonal matrix of 1/var(meas)
0008 %
0009 % calc_meas_icov can be called as
0010 %    meas_icov= calc_meas_icov( inv_model )
0011 %
0012 % meas_icov   is the calculated data prior
0013 % inv_model    is an inv_model structure
0014 %
0015 % if:
0016 %    inv_model.meas_icov    is a function
0017 %          -> call it to calculate meas_icov
0018 %    inv_model.meas_icov    is a matrix
0019 %          -> return it as meas_icov
0020 %    inv_model.meas_icov    does not exist
0021 %          -> use I, or 1./homg (for normalized difference)
0022 
0023 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0024 % $Id: calc_meas_icov.m 4832 2015-03-29 21:13:53Z bgrychtol-ipa $
0025 
0026 
0027 if isfield(inv_model,'meas_icov')
0028    if isnumeric(inv_model.meas_icov);
0029       meas_icov = inv_model.meas_icov;
0030    else
0031       try inv_model.meas_icov = str2func(inv_model.meas_icov); end
0032       meas_icov= eidors_cache( inv_model.meas_icov, {inv_model});
0033    end
0034 else
0035    meas_icov= eidors_cache(@default_meas_icov,{inv_model} );
0036 end
0037 
0038  
0039 % Calculate a data prior for an assumption of uniform noise
0040 % on each channel
0041 %
0042 function meas_icov = default_meas_icov( inv_model )
0043 
0044    fwd_model= inv_model.fwd_model;
0045    fwd_model = mdl_normalize(fwd_model,mdl_normalize(fwd_model));
0046 
0047    n =  calc_n_meas( fwd_model );
0048 
0049    if ~mdl_normalize(fwd_model);
0050       meas_icov= speye( n );
0051    else
0052       homg_data=  solve_homg_image( fwd_model );
0053 % if we normalize, then small data get increased
0054 %   this means that noise on small data gets increased,
0055 %    so the covariance is large when data are small
0056 %   so the icov is small when data are small
0057 % sig = k/h -> std = k/h -> 1/std = kh
0058       meas_icov = sparse(1:n, 1:n, ( homg_data.meas ).^2 );
0059    end
0060 
0061 function n_meas = calc_n_meas( fwd_model )
0062 
0063    n_meas = 0;
0064    for i= 1:length(fwd_model.stimulation );
0065        n_meas = n_meas + size(fwd_model.stimulation(i).meas_pattern,1);
0066    end
0067 
0068 % create homogeneous image + simulate data
0069 function homg_data = solve_homg_image( fwd_mdl )
0070     n_elems= size( fwd_mdl.elems , 1);
0071     mat= ones( n_elems, 1);
0072     homg_img= eidors_obj('image', 'homogeneous image', ...
0073                          'elem_data', mat, 'fwd_model', fwd_mdl );
0074     homg_data=fwd_solve( homg_img);