Description of calc_noise

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
