calc_noise_figure

PURPOSE ^

CALC_NOISE_FIGURE: calculate the noise amplification (NF) of an algorithm

SYNOPSIS ^

function [NF,SE] = calc_noise_figure( inv_model, hp, iterations)

DESCRIPTION ^

 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)
    if hp is [] => use the hp on the inv_model
  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

 [NF,SE] = calc_noise_figure( inv_model, vh, vi)
    interface provided for compatibility with the deprecated
    calc_noise_params

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 %    if hp is [] => use the hp on the inv_model
0010 %  NF = calculated NF. SE = standard error on NF
0011 %
0012 % hp is specified, it will be used for the hyperparameter.
0013 %    Otherwise the inv_model.hyperparameter will be used.
0014 %
0015 % Noise Figure must be defined for a specific measurement
0016 % In order to specify data, use
0017 %     inv_model.hyperparameter.tgt_data.meas_t1
0018 %     inv_model.hyperparameter.tgt_data.meas_t2
0019 %   to use a temporal solver (or the Kalman filter), the
0020 %   measurement to perform the NF calc must also be specified,
0021 %   using:
0022 %     inv_model.hyperparameter.tgt_data.meas_select
0023 %   otherwise, the middle measurement will be used
0024 %
0025 % In order to automatically simulate data, specify tgt_elems,
0026 %   containing a vector of elements to use
0027 %
0028 %     inv_model.hyperparameter.tgt_elems
0029 %
0030 % [NF,SE] = calc_noise_figure( inv_model, vh, vi)
0031 %    interface provided for compatibility with the deprecated
0032 %    calc_noise_params
0033 
0034 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0035 % $Id: calc_noise_figure.m 5112 2015-06-14 13:00:41Z aadler $
0036 
0037 % A normal definition of noise power is based on power:
0038 %      NF = SNR_z / SNR_x
0039 %    SNR_z = sumsq(z0) / var(z) = sum(z0.^2) / trace(Rn)
0040 %    SNR_x = sumsq(x0) / var(x) = sum(x0.^2) / trace(ABRnB'A')
0041 % but his doesn't work, because the signal spreads like
0042 % the amplitude, not like the power, thus
0043 %      NF = SNR_z / SNR_x
0044 %    SNR_z = sum(|z0|/len_z) / std(z/len_z)
0045 %    SNR_x = sum(|x0|/len_x) / std(x/len_x)
0046 
0047 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'), do_unit_test, return, end
0048 
0049 if nargin>=2 && numel(hp) == 1
0050    inv_model.hyperparameter.value= hp;
0051 % Remove function parameter because it will recurse
0052    try; inv_model.hyperparameter = rmfield(inv_model.hyperparameter,'func'); end
0053 end
0054 if nargin == 3 && numel(hp) > 1
0055     h_data = hp;
0056     c_data = iterations;
0057 else
0058     [inv_model, h_data, c_data] = process_parameters( inv_model );
0059 end
0060 
0061 % disable solution checking (slow)
0062 inv_model.inv_solve.calc_solution_error = false;
0063 
0064 %NF= nf_calc_use_matrix( inv_model, h_data, c_data);
0065 %NF= nf_calc_iterate( inv_model, h_data, c_data);
0066 if nargin<3; iterations= 10; end
0067 solver = inv_model.solve;
0068 if ischar(solver) && strcmp(solver, 'eidors_default')
0069     solver = eidors_default('get','inv_solve');
0070 end
0071 if isa(solver,'function_handle')
0072     solver = func2str(solver);
0073 end
0074 switch solver
0075     case {'inv_solve_backproj'
0076             'inv_solve_conj_grad'
0077             'inv_solve_diff_GN_one_step'
0078             'inv_solve_trunc_iterative'
0079             'inv_solve_TSVD'
0080             'solve_use_matrix'}
0081         NF = nf_calc_linear( inv_model, h_data, c_data);
0082         SE = 0;
0083     otherwise
0084         [NF,SE]= nf_calc_random( inv_model, h_data, c_data, iterations);
0085 end
0086 eidors_msg('@@ NF= %f', NF, 2);
0087 
0088 function [inv_model, h_data, c_data] = process_parameters( inv_model );
0089 
0090    if     isfield(inv_model.hyperparameter,'tgt_elems')
0091       [h_data, c_data]= simulate_targets( inv_model.fwd_model, ...
0092            inv_model.hyperparameter.tgt_elems);
0093    elseif isfield(inv_model.hyperparameter,'tgt_data')
0094       tgt_data= inv_model.hyperparameter.tgt_data;
0095       h_data= tgt_data.meas_t1;
0096       c_data= tgt_data.meas_t2;
0097    else
0098       error('unsure how to get data to measure signal');
0099    end
0100 
0101    % if hp is specified, then use that value
0102 
0103    if isfield(inv_model.hyperparameter,'func')
0104       funcname= inv_model.hyperparameter.func;
0105       if strcmp( class(funcname), 'function_handle')
0106          funcname= func2str(funcname);
0107       end
0108 
0109       if strcmp(funcname, 'choose_noise_figure')
0110          error('specifying inv_model.hp.func = choose_noise_figure will recurse');
0111       end
0112    end
0113 
0114 function params = nf_calc_linear(imdl, vh, vi )
0115 % params = GREIT_noise_params(imdl, homg_voltage, sig_voltage)
0116 %  params(1,:) = Noise Figure = SNR(image) / SNR(data)
0117 %
0118 %  see also: eval_GREIT_fig_merit or using test_performance
0119 
0120 % (C) 2008 Andy Adler. Licensed under GPL v2 or v3
0121 % $Id: calc_noise_figure.m 5112 2015-06-14 13:00:41Z aadler $
0122 
0123 % NOTE THAT WE ASSUME A LINEAR ALGORITHM FOR THIS MEASURE
0124 
0125 if 0 % old code with random noise
0126      % Keep this to validate while we test it
0127    Nnoise = 1000;
0128    noise = 0.01*std(vh)*randn(size(vh,1),Nnoise);
0129    vhn= vh*ones(1,Nnoise) + noise;
0130 else % use independent noise model on each channel
0131    noise = 0.01*std(vh)*speye(size(vh,1));
0132    vhn= vh*ones(1,size(vh,1)) + noise;
0133 end
0134 
0135 signal_y = calc_difference_data( vh, vi,  imdl.fwd_model);
0136 noise_y  = calc_difference_data( vh, vhn, imdl.fwd_model);
0137 
0138 signal_x = inv_solve(imdl, vh, vi);  
0139 signal_x = data_mapper(signal_x); signal_x = signal_x.elem_data;
0140 noise_x  = inv_solve(imdl, vh, vhn); 
0141 noise_x  = data_mapper(noise_x);  noise_x  = noise_x.elem_data;
0142 
0143 use_rec = 1;
0144 try 
0145    use_rec = ~imdl.prior_use_fwd_not_rec;
0146 end
0147 if use_rec
0148    try
0149       VOL = get_elem_volume(imdl.rec_model);
0150    catch
0151       VOL = get_elem_volume(imdl.fwd_model);
0152    end
0153 else
0154    VOL = get_elem_volume(imdl.fwd_model);
0155 end
0156 VOL = spdiags(VOL,0, length(VOL), length(VOL));
0157 
0158 signal_x = VOL*signal_x;
0159 noise_x = VOL*noise_x;
0160 
0161 signal_x = mean(abs(signal_x),1);
0162 noise_x  = mean(std(noise_x));
0163 snr_x = signal_x / noise_x;
0164 
0165 signal_y = mean(abs(signal_y),1);
0166 noise_y  = mean(std(noise_y)); 
0167 snr_y = signal_y / noise_y;
0168 
0169 params= [snr_y(:)./snr_x(:)]';
0170 
0171 try
0172 eidors_msg('@@ NF= %f (hp=%e)', NF, imdl.hyperparameter.value, 2);
0173 end
0174 
0175 
0176 % NOTES on the calculations: AA - Feb 20, 2012
0177 % SNR = mean(abs(x)); VAR =
0178 % ym= E[y]
0179 % Sy= E[(y-ym)*(y-ym)'] = E[y*y'] - ym*ym'
0180 % ny = sqrt(trace(Sy))
0181 % xm= E[x]  = E[R*y] = R*E[y] = R*ym
0182 % Sx= E[(x-xm)*(x-xm)'] = E[x*x'] - xm*xm'
0183 %   = E[R*ym*ym'*R'] = R*E[ym*ym']*R' = R*Sy*R'
0184 % nx = sqrt(trace(Sx))
0185 %
0186 % signal = mean(abs(x));
0187 %
0188 % In this case, these are exactly the same:
0189 %    noise_x  = mean(std(noise_x));
0190 %    noise_x  = sqrt(mean(noise_x.^2,2));
0191    
0192    
0193 function NF= nf_calc_use_matrix( inv_model, h_data, c_data)
0194 % To model std(z) we use z=z0+n
0195 % so that std(z) = sqrt(var(z)) = sqrt(1/L * E[n'*n])
0196 % we know a priori that the mean noise is zero, thus
0197 % we do not need to divide by L-1 in the variance
0198 % E[n'*n] = E[ trace(n*n') ] = trace( cov_N )
0199 % Thus, std(z) = sqrt( trace( cov_N )/L )
0200 %              = sqrt( mean( diag( cov_N )))
0201 % And,  std(x) = sqrt( mean( diag( cov_X )))
0202 %
0203 % To model cov_N, we consider independant noise
0204 %  on each channel, cov_N= N*N', N=diag( sigma_i )
0205 % And,              cov_X= X*X', X=reconst(N)
0206 % where X= reconst(mdl, z0,z0+N)
0207 %
0208 % To run efficiently mean(diag(cov_N))=mean(sum(N.^2,2))
0209 % The sum over X is actually weighted by the volume of
0210 %  each element, so VOL.^2*(sum(X.^2,2)
0211    try 
0212        VOL = get_elem_volume(inv_model.rec_model)';
0213    catch
0214        VOL = get_elem_volume(inv_model.fwd_model)';
0215    end
0216 
0217    % calculate signal
0218    d_len   = size(h_data,1);
0219    delta   = 1e-2* mean(h_data);
0220    c_noise = c_data*ones(1,d_len) + eye(d_len);
0221    h_full  = h_data*ones(1,d_len);
0222 
0223    sig_data = mean(abs( ...
0224          calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0225                        ));
0226    var_data = mean(sum( ...
0227          calc_difference_data( h_full, c_noise, inv_model.fwd_model ) ...
0228                        .^2, 2)); 
0229                       
0230 
0231    % calculate image
0232    % Note, this won't work if the algorithm output is not zero biased
0233    [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0234                                h_full, c_noise);
0235 
0236    i_len = length(img0);
0237    sig_img= VOL*abs(img0) / i_len;;
0238    var_img= VOL.^2*sum(img0n.^2 ,2) / i_len;
0239    
0240    NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img)  );
0241 
0242    % For the record, the expression for var_img is derived as:
0243    % Equiv expresssions for var_img % given: A= diag(pp.VOLUME);
0244    % var_img= trace(A*RM*diag(Rn.^2)*RM'*A');
0245    % vv=A*RM*diag(Rn);var_img=trace(vv*vv'); var_img= sum(sum(vv.^2));
0246    % var_img= VOL2* (RM*diag(Rn)).^2
0247    % var_img= VOL2* RM.^2 * Rn.^2
0248 
0249    
0250 % simulate homg data and a small target in centre
0251 function [h_data, c_data]= simulate_targets( fwd_model, ctr_elems)
0252 
0253    homg= 1; % homogeneous conductivity level is 1
0254 
0255    %Step 1: homogeneous image
0256    sigma= homg*ones( size(fwd_model.elems,1) ,1);
0257 
0258    img= eidors_obj('image', 'homogeneous image', ...
0259                    'elem_data', sigma, ...
0260                    'fwd_model', fwd_model );
0261    h_data=fwd_solve( img );
0262    h_data= h_data.meas;
0263 
0264    %Step 1: inhomogeneous image with contrast in centre
0265    delta = 1e-2;
0266    sigma(ctr_elems) = homg*(1 + delta);
0267    img.elem_data = sigma;
0268    c_data=fwd_solve( img );
0269    c_data= c_data.meas;
0270 
0271 function [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0272                                h_full, c_noise);
0273    if isa(inv_model.solve,'function_handle')
0274       solve= func2str(inv_model.solve);
0275    else
0276       solve= inv_model.solve;
0277    end
0278 
0279 % Test for special functions and solve them specially
0280    switch solve
0281    case 'ab_tv_diff_solve'
0282       error('Dont know how to calculate TV noise figure')
0283 
0284    case 'inv_solve_diff_kalman'
0285       inv_model.inv_solve_diff_kalman.keep_K_k1= 1;
0286       stablize = 6;
0287       img0 = inv_solve( inv_model, h_data, ...
0288                                    c_data*ones(1,stablize) );
0289       K= img0.inv_solve_diff_kalman.K_k1;
0290       img0.elem_data = K*calc_difference_data( h_data , c_data , inv_model.fwd_model);
0291       img0n.elem_data= K*calc_difference_data( h_full , c_noise, inv_model.fwd_model);
0292 
0293    otherwise
0294       img0 = inv_solve( inv_model, h_data, c_data);
0295       if nargin>4
0296       img0n= inv_solve( inv_model, h_full, c_noise);
0297       end
0298    end
0299 
0300    % Need elem or nodal data
0301    if isfield(img0,'node_data');
0302       img0 = img0.node_data;
0303    else
0304       img0 = img0.elem_data;
0305    end
0306 
0307    if isfield(img0n,'node_data');
0308       img0n = img0n.node_data;
0309    else
0310       img0n = img0n.elem_data;
0311    end
0312 
0313 % OLD CODE - iterate
0314 function NF= nf_calc_iterate( inv_model, h_data, c_data);
0315    try 
0316        VOL = get_elem_volume(inv_model.rec_model)';
0317    catch
0318        VOL = get_elem_volume(inv_model.fwd_model)';
0319    end
0320    % calculate signal
0321    d_len   = size(h_data,1);
0322    delta   = 1e-2* mean(h_data);
0323    sig_data = mean(abs( ...
0324          calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0325                        ));
0326    % calculate image
0327    % Note, this won't work if the algorithm output is not zero biased
0328 
0329    [img0] = get_images( inv_model, h_data, c_data);
0330 %  sig_img= mean(VOL'.*abs(img0.elem_data));
0331    sig_img= VOL*abs(img0) / length(img0);
0332 
0333    % Now do noise
0334    var_data= 0;
0335    var_img = 0;
0336    for i=1:d_len
0337       this_noise = -ones(d_len, size(c_data,2))/(d_len-1);
0338       this_noise(i,:) = 1;
0339       c_noise = c_data + this_noise;
0340       [imgn] = get_images( inv_model, h_data, c_noise);
0341       if 1
0342          var_data = var_data + mean(sum( ...
0343             calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0344                           .^2, 2)); 
0345 %        var_img= var_img +  mean( (VOL'.*imgn.elem_data).^2 );
0346          var_img= var_img +  (VOL.^2)*sum(imgn.elem_data.^2,2 ) / length(imgn.elem_data); 
0347       else
0348          % OLD APPROACH BASED ON variance, rather than matrix calcs
0349          var_data = var_data + var( ...
0350             calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0351                                  ); 
0352          var_img= var_img + var( VOL'.*imgn.elem_data ); 
0353       end
0354    end
0355    var_data = var_data / d_len;
0356    var_img  = var_img  / d_len;
0357    NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img)  );
0358 
0359 function [NF,SE]= nf_calc_random( rec, vh, vi, N_RUNS);
0360    eidors_cache('boost_priority',-2); % low priority values
0361 
0362    imgr= inv_solve(rec, vh, vi);
0363 
0364    if isfield(imgr,'node_data');
0365       img0 = imgr.node_data;
0366       try
0367           VOL = get_elem_volume(rec.rec_model, 1);
0368       catch
0369           VOL = get_elem_volume(rec.fwd_model, 1);
0370       end
0371    else
0372       n_els = length(rec.fwd_model);
0373       img0 = imgr.elem_data(1:n_els); % elem_data can also contain movement
0374       try
0375           VOL = get_elem_volume(rec.rec_model, 0);
0376       catch
0377           VOL = get_elem_volume(rec.fwd_model, 0);
0378       end
0379    end
0380 
0381    sig_ampl = mean( abs( VOL .* img0 )) / ...
0382               mean( abs( calc_difference_data( vh, vi, rec.fwd_model )));
0383 
0384 % Estimate Signal Amplitude
0385    for i=1:N_RUNS
0386       vn= addnoise(vh, vi, 1.0);
0387 
0388       imgr= inv_solve(rec, vh, vn);
0389 
0390       if isfield(imgr,'node_data'); img0 = imgr.node_data;
0391       else;                         img0 = imgr.elem_data(1:n_els);
0392       end
0393 
0394       noi_imag(i) = std( VOL .* img0 );
0395       noi_sgnl(i) = std( calc_difference_data( vh, vn, rec.fwd_model ));
0396    end
0397 
0398    noi_ampl = noi_imag./noi_sgnl;
0399    NF =  mean(noi_ampl/sig_ampl);
0400    SE =  std(noi_ampl/sig_ampl)/sqrt(N_RUNS);
0401    eidors_msg('NF= %f+/-%f', NF, SE, 3);
0402 
0403    eidors_cache('boost_priority',2);
0404 
0405 function noise= addnoise( vh, vi, SNR);
0406       if isstruct(vh); vh= vh.meas; end
0407       if isstruct(vi); vi= vi.meas; end
0408       noise = randn(size(vh));
0409       noise = noise*std(vh-vi)/std(noise);
0410       noise = vh + SNR*noise;
0411 
0412 function do_unit_test
0413     ll = eidors_msg('log_level',1);
0414     test1; % can we deal with c2f ?
0415     imdl = mk_common_model('a2t2',16); test2(imdl);
0416     imdl = mk_common_model('d2t2',16); test2(imdl);
0417     imdl = mk_common_model('a2c2',16); test2(imdl);
0418     imdl = mk_common_model('d2c2',16); test2(imdl);
0419     ll = eidors_msg('log_level',ll);
0420 
0421 function test1
0422     % big model with c2f and supposedly an NF of 0.5
0423     fmdl = mk_library_model('pig_23kg_16el');
0424     [fmdl.stimulation fmdl.meas_select] = mk_stim_patterns(16,1,'{ad}','{ad}');
0425     fmdl = mdl_normalize(fmdl, 1);  % Use normalized difference imaging
0426     opt.noise_figure = 0.5; opt.imgsz = [64 64];
0427     imdl = mk_GREIT_model(fmdl, 0.25, [], opt);
0428     % homogeneous measurement
0429     img = mk_image(fmdl,1);
0430     vh = fwd_solve(img);
0431     % inhomogeneous measurement
0432     select_fcn = inline('(x-0).^2+(y-0).^2+(z-0.5).^2<0.1^2','x','y','z');
0433     mfrac = elem_select(fmdl, select_fcn);
0434     img.elem_data = img.elem_data + mfrac*0.1;
0435     vi = fwd_solve(img);
0436     
0437     nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0438     
0439     % We use different measurements so naturally we don't get 0.5 here
0440     imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0441     imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0442     try
0443         % calc_noise_figure doens't support dual models
0444         nf2 = calc_noise_figure(imdl);
0445     catch
0446         nf2 = 0;
0447     end
0448     unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0449 
0450 function test2(imdl)
0451     fmdl = imdl.fwd_model;
0452     % homogeneous measurement
0453     img = mk_image(fmdl,1);
0454     vh = fwd_solve(img);
0455     % inhomogeneous measurement
0456     select_fcn = inline('(x-0).^2+(y-0).^2.^2<15^2','x','y','z');
0457     mfrac = elem_select(fmdl, select_fcn);
0458     img.elem_data = img.elem_data + mfrac*0.1;
0459     vi = fwd_solve(img);
0460 
0461     nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0462     eidors_msg(nf1,0);
0463 
0464 imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0465 imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0466 % calc_noise_figure doens't support dual models
0467 nf2 = calc_noise_figure(imdl,[],1000);
0468 unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0469

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005