calc_noise_figure

 CALC_NOISE_FIGURE: calculate the noise amplification (NF) of an algorithm
 [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
    inv_model  => inverse model object
    hp         => value of hyperparameter to use (if not spec
         then use the value of inv_model.hyperparameter.value)
    iterations => number of iterations (default 10)
       (for calculation of noise figure using random noise)
  NF = calculated NF. SE = standard error on NF

 hp is specified, it will be used for the hyperparameter.
    Otherwise the inv_model.hyperparameter will be used.

 Noise Figure must be defined for a specific measurement
 In order to specify data, use
     inv_model.hyperparameter.tgt_data.meas_t1
     inv_model.hyperparameter.tgt_data.meas_t2
   to use a temporal solver (or the Kalman filter), the
   measurement to perform the NF calc must also be specified,
   using:
     inv_model.hyperparameter.tgt_data.meas_select
   otherwise, the middle measurement will be used

 In order to automatically simulate data, specify tgt_elems,
   containing a vector of elements to use
 
     inv_model.hyperparameter.tgt_elems

0001 function [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
0002 % CALC_NOISE_FIGURE: calculate the noise amplification (NF) of an algorithm
0003 % [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
0004 %    inv_model  => inverse model object
0005 %    hp         => value of hyperparameter to use (if not spec
0006 %         then use the value of inv_model.hyperparameter.value)
0007 %    iterations => number of iterations (default 10)
0008 %       (for calculation of noise figure using random noise)
0009 %  NF = calculated NF. SE = standard error on NF
0010 %
0011 % hp is specified, it will be used for the hyperparameter.
0012 %    Otherwise the inv_model.hyperparameter will be used.
0013 %
0014 % Noise Figure must be defined for a specific measurement
0015 % In order to specify data, use
0016 %     inv_model.hyperparameter.tgt_data.meas_t1
0017 %     inv_model.hyperparameter.tgt_data.meas_t2
0018 %   to use a temporal solver (or the Kalman filter), the
0019 %   measurement to perform the NF calc must also be specified,
0020 %   using:
0021 %     inv_model.hyperparameter.tgt_data.meas_select
0022 %   otherwise, the middle measurement will be used
0023 %
0024 % In order to automatically simulate data, specify tgt_elems,
0025 %   containing a vector of elements to use
0026 %
0027 %     inv_model.hyperparameter.tgt_elems
0028 %
0029 
0030 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0031 % $Id: calc_noise_figure.html 2819 2011-09-07 16:43:11Z aadler $
0032 
0033 % A normal definition of noise power is based on power:
0034 %      NF = SNR_z / SNR_x
0035 %    SNR_z = sumsq(z0) / var(z) = sum(z0.^2) / trace(Rn)
0036 %    SNR_x = sumsq(x0) / var(x) = sum(x0.^2) / trace(ABRnB'A')
0037 % but his doesn't work, because the signal spreads like
0038 % the amplitude, not like the power, thus
0039 %      NF = SNR_z / SNR_x
0040 %    SNR_z = sum(|z0|/len_z) / std(z/len_z)
0041 %    SNR_x = sum(|x0|/len_x) / std(x/len_x)
0042 
0043 
0044 if nargin>=2
0045    inv_model.hyperparameter.value= hp;
0046 % Remove function parameter because it will recurse
0047    try; inv_model.hyperparameter = rmfield(inv_model.hyperparameter,'func'); end
0048 end
0049 [inv_model, h_data, c_data] = process_parameters( inv_model );
0050 
0051 %NF= nf_calc_use_matrix( inv_model, h_data, c_data);
0052 %NF= nf_calc_iterate( inv_model, h_data, c_data);
0053 if nargin<3; iterations= 10; end
0054 [NF,SE]= nf_calc_random( inv_model, h_data, c_data, iterations);
0055 eidors_msg('calculating NF=%f', NF, 2);
0056 
0057 function [inv_model, h_data, c_data] = process_parameters( inv_model );
0058 
0059    if     isfield(inv_model.hyperparameter,'tgt_elems')
0060       [h_data, c_data]= simulate_targets( inv_model.fwd_model, ...
0061            inv_model.hyperparameter.tgt_elems);
0062    elseif isfield(inv_model.hyperparameter,'tgt_data')
0063       tgt_data= inv_model.hyperparameter.tgt_data;
0064       h_data= tgt_data.meas_t1;
0065       c_data= tgt_data.meas_t2;
0066    else
0067       error('unsure how to get data to measure signal');
0068    end
0069 
0070    % if hp is specified, then use that value
0071 
0072    if isfield(inv_model.hyperparameter,'func')
0073       funcname= inv_model.hyperparameter.func;
0074       if strcmp( class(funcname), 'function_handle')
0075          funcname= func2str(funcname);
0076       end
0077 
0078       if strcmp(funcname, 'choose_noise_figure')
0079          error('specifying inv_model.hp.func = choose_noise_figure will recurse');
0080       end
0081    end
0082 
0083 
0084 function NF= nf_calc_use_matrix( inv_model, h_data, c_data)
0085 % To model std(z) we use z=z0+n
0086 % so that std(z) = sqrt(var(z)) = sqrt(1/L * E[n'*n])
0087 % we know a priori that the mean noise is zero, thus
0088 % we do not need to divide by L-1 in the variance
0089 % E[n'*n] = E[ trace(n*n') ] = trace( cov_N )
0090 % Thus, std(z) = sqrt( trace( cov_N )/L )
0091 %              = sqrt( mean( diag( cov_N )))
0092 % And,  std(x) = sqrt( mean( diag( cov_X )))
0093 %
0094 % To model cov_N, we consider independant noise
0095 %  on each channel, cov_N= N*N', N=diag( sigma_i )
0096 % And,              cov_X= X*X', X=reconst(N)
0097 % where X= reconst(mdl, z0,z0+N)
0098 %
0099 % To run efficiently mean(diag(cov_N))=mean(sum(N.^2,2))
0100 % The sum over X is actually weighted by the volume of
0101 %  each element, so VOL.^2*(sum(X.^2,2)
0102    VOL = get_elem_volume(inv_model.fwd_model)';
0103 
0104    % calculate signal
0105    d_len   = size(h_data,1);
0106    delta   = 1e-2* mean(h_data);
0107    c_noise = c_data*ones(1,d_len) + eye(d_len);
0108    h_full  = h_data*ones(1,d_len);
0109 
0110    sig_data = mean(abs( ...
0111          calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0112                        ));
0113    var_data = mean(sum( ...
0114          calc_difference_data( h_full, c_noise, inv_model.fwd_model ) ...
0115                        .^2, 2)); 
0116                       
0117 
0118    % calculate image
0119    % Note, this won't work if the algorithm output is not zero biased
0120    [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0121                                h_full, c_noise);
0122 
0123    i_len = length(img0);
0124    sig_img= VOL*abs(img0) / i_len;;
0125    var_img= VOL.^2*sum(img0n.^2 ,2) / i_len;
0126    
0127    NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img)  );
0128 
0129    % For the record, the expression for var_img is derived as:
0130    % Equiv expresssions for var_img % given: A= diag(pp.VOLUME);
0131    % var_img= trace(A*RM*diag(Rn.^2)*RM'*A');
0132    % vv=A*RM*diag(Rn);var_img=trace(vv*vv'); var_img= sum(sum(vv.^2));
0133    % var_img= VOL2* (RM*diag(Rn)).^2
0134    % var_img= VOL2* RM.^2 * Rn.^2
0135 
0136    
0137 % simulate homg data and a small target in centre
0138 function [h_data, c_data]= simulate_targets( fwd_model, ctr_elems)
0139 
0140    homg= 1; % homogeneous conductivity level is 1
0141 
0142    %Step 1: homogeneous image
0143    sigma= homg*ones( size(fwd_model.elems,1) ,1);
0144 
0145    img= eidors_obj('image', 'homogeneous image', ...
0146                    'elem_data', sigma, ...
0147                    'fwd_model', fwd_model );
0148    h_data=fwd_solve( img );
0149    h_data= h_data.meas;
0150 
0151    %Step 1: inhomogeneous image with contrast in centre
0152    delta = 1e-2;
0153    sigma(ctr_elems) = homg*(1 + delta);
0154    img.elem_data = sigma;
0155    c_data=fwd_solve( img );
0156    c_data= c_data.meas;
0157 
0158 function [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0159                                h_full, c_noise);
0160    if isa(inv_model.solve,'function_handle')
0161       solve= func2str(inv_model.solve);
0162    else
0163       solve= inv_model.solve;
0164    end
0165 
0166 % Test for special functions and solve them specially
0167    switch solve
0168    case 'ab_tv_diff_solve'
0169       error('Dont know how to calculate TV noise figure')
0170 
0171    case 'inv_kalman_diff'
0172       inv_model.inv_kalman_diff.keep_K_k1= 1;
0173       stablize = 6;
0174       img0 = inv_solve( inv_model, h_data, ...
0175                                    c_data*ones(1,stablize) );
0176       K= img0.inv_kalman_diff.K_k1;
0177       img0.elem_data = K*calc_difference_data( h_data , c_data , inv_model.fwd_model);
0178       img0n.elem_data= K*calc_difference_data( h_full , c_noise, inv_model.fwd_model);
0179 
0180    otherwise
0181       img0 = inv_solve( inv_model, h_data, c_data);
0182       if nargin>4
0183       img0n= inv_solve( inv_model, h_full, c_noise);
0184       end
0185    end
0186 
0187    % Need elem or nodal data
0188    if isfield(img0,'node_data');
0189       img0 = img0.node_data;
0190    else
0191       img0 = img0.elem_data;
0192    end
0193 
0194    if isfield(img0n,'node_data');
0195       img0n = img0n.node_data;
0196    else
0197       img0n = img0n.elem_data;
0198    end
0199 
0200 % OLD CODE - iterate
0201 function NF= nf_calc_iterate( inv_model, h_data, c_data);
0202    VOL = get_elem_volume(inv_model.fwd_model)';
0203 
0204    % calculate signal
0205    d_len   = size(h_data,1);
0206    delta   = 1e-2* mean(h_data);
0207    sig_data = mean(abs( ...
0208          calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0209                        ));
0210    % calculate image
0211    % Note, this won't work if the algorithm output is not zero biased
0212 
0213    [img0] = get_images( inv_model, h_data, c_data);
0214 %  sig_img= mean(VOL'.*abs(img0.elem_data));
0215    sig_img= VOL*abs(img0) / length(img0);
0216 
0217    % Now do noise
0218    var_data= 0;
0219    var_img = 0;
0220    for i=1:d_len
0221       this_noise = -ones(d_len, size(c_data,2))/(d_len-1);
0222       this_noise(i,:) = 1;
0223       c_noise = c_data + this_noise;
0224       [imgn] = get_images( inv_model, h_data, c_noise);
0225       if 1
0226          var_data = var_data + mean(sum( ...
0227             calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0228                           .^2, 2)); 
0229 %        var_img= var_img +  mean( (VOL'.*imgn.elem_data).^2 );
0230          var_img= var_img +  (VOL.^2)*sum(imgn.elem_data.^2,2 ) / length(imgn.elem_data); 
0231       else
0232          % OLD APPROACH BASED ON variance, rather than matrix calcs
0233          var_data = var_data + var( ...
0234             calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0235                                  ); 
0236          var_img= var_img + var( VOL'.*imgn.elem_data ); 
0237       end
0238    end
0239    var_data = var_data / d_len;
0240    var_img  = var_img  / d_len;
0241    NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img)  );
0242 
0243 function [NF,SE]= nf_calc_random( rec, vh, vi, N_RUNS);
0244    eidors_cache('boost_priority',-2); % low priority values
0245 
0246    imgr= inv_solve(rec, vh, vi);
0247 
0248    if isfield(imgr,'node_data');
0249       img0 = imgr.node_data;
0250       VOL = get_elem_volume(rec.fwd_model, 1);
0251    else
0252       img0 = imgr.elem_data;
0253       VOL = get_elem_volume(rec.fwd_model, 0);
0254    end
0255 
0256    sig_ampl = mean( abs( VOL .* img0 )) / ...
0257               mean( abs( calc_difference_data( vh, vi, rec.fwd_model )));
0258 
0259 % Estimate Signal Amplitude
0260    for i=1:N_RUNS
0261       vn= addnoise(vh, vi, 1.0);
0262 
0263       imgr= inv_solve(rec, vh, vn);
0264 
0265       if isfield(imgr,'node_data'); img0 = imgr.node_data;
0266       else;                         img0 = imgr.elem_data;
0267       end
0268 
0269       noi_imag(i) = std( VOL .* img0 );
0270       noi_sgnl(i) = std( calc_difference_data( vh, vn, rec.fwd_model ));
0271    end
0272 
0273    noi_ampl = noi_imag./noi_sgnl;
0274    NF =  mean(noi_ampl/sig_ampl);
0275    SE =  std(noi_ampl/sig_ampl)/sqrt(N_RUNS);
0276    eidors_msg('NF= %f+/-%f', NF, SE, 1);
0277 
0278    eidors_cache('boost_priority',2);
0279 
0280    function noise= addnoise( vh, vi, SNR);
0281       if isstruct(vh); vh= vh.meas; end
0282       if isstruct(vi); vi= vi.meas; end
0283       noise = randn(size(vh));
0284       noise = noise*std(vh-vi)/std(noise);
0285       noise = vh + SNR*noise;
0286 
0287