choose_noise_figure

 CHOOSE_NOISE_FIGURE: choose hyperparameter based on NF calculation
 HP= choose_noise_figure( inv_model );
 inv_model  => inverse model struct

 In order to use this function, it is necessary to specify
 inv_model.hyperparameter. has the following fields
 hpara.func         = @choose_noise_figure;
 hpara.noise_figure = NF Value requested
 hpara.tgt_elems    = vector of element numbers of contrast in centre

 The NF parameter is from the definition in Adler & Guardo (1996).

 SNR_z = sumsq(z0) / var(z) = sum(z0.^2) / trace(Rn)
 SNR_x = sumsq(A*x0) / var(A*x) = sum((A*x0).^2) / trace(ABRnB'A')
   where Rn = noise covariance and A_ii = area of element i
 NF = SNR_z / SNR_x

0001 function HP= choose_noise_figure( inv_model );
0002 % CHOOSE_NOISE_FIGURE: choose hyperparameter based on NF calculation
0003 % HP= choose_noise_figure( inv_model );
0004 % inv_model  => inverse model struct
0005 %
0006 % In order to use this function, it is necessary to specify
0007 % inv_model.hyperparameter. has the following fields
0008 % hpara.func         = @choose_noise_figure;
0009 % hpara.noise_figure = NF Value requested
0010 % hpara.tgt_elems    = vector of element numbers of contrast in centre
0011 %
0012 % The NF parameter is from the definition in Adler & Guardo (1996).
0013 %
0014 % SNR_z = sumsq(z0) / var(z) = sum(z0.^2) / trace(Rn)
0015 % SNR_x = sumsq(A*x0) / var(A*x) = sum((A*x0).^2) / trace(ABRnB'A')
0016 %   where Rn = noise covariance and A_ii = area of element i
0017 % NF = SNR_z / SNR_x
0018 
0019 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0020 % $Id: choose_noise_figure.m 6942 2024-06-18 14:53:21Z aadler $
0021 
0022 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0023 
0024 reqNF= inv_model.hyperparameter.noise_figure;
0025 try % remove the value field
0026    inv_model.hyperparameter = rmfield(inv_model.hyperparameter,'value');
0027 end
0028 
0029 copt.boost_priority = 3;
0030 HP = eidors_cache(@HP_for_NF_search,{reqNF,inv_model},copt);
0031 
0032 
0033 
0034 function [HP,NF,SE] = HP_for_NF_search(dNF,imdl)
0035    llv = eidors_msg('log_level');
0036    if llv==3; eidors_msg('log_level',2); end % HIDE A LOT OF MESSAGES
0037 
0038    hp= search1(dNF, imdl, 1);
0039 
0040    dx= hp-linspace(-0.7,0.7,5);
0041    hp= search2(dNF, imdl, dx);
0042 
0043    dx= hp-linspace(-0.2,0.2,5);
0044    hp= search2(dNF,imdl, dx);
0045 
0046    dx= hp-linspace(-0.1,0.1,5); 
0047    hp= search2(dNF,imdl, dx);
0048 
0049    dx= hp-linspace(-0.05,0.05,21);
0050    hp= search2(dNF,imdl, dx);
0051 
0052    HP= 10^-hp;
0053    eidors_msg('log_level',llv); %% Reset log_level
0054   
0055 function hp= search1(dNF, imdl, hp)
0056   [NF,SE]=calc_noise_figure( imdl, 10^(-hp));
0057    if     NF+3*SE < dNF; dir = 1;
0058    elseif NF-3*SE > dNF; dir = -1;
0059    else   dir = 0; end
0060    while  dir*NF+3*SE < dir*dNF %>
0061      hp= hp+0.5*dir;
0062      [NF,SE]=calc_noise_figure( imdl, 10^(-hp));
0063    end
0064 
0065 function hp=search2(dNF, imdl, dx)
0066    hp = [];
0067    it = 0;
0068    nf = zeros(1,length(dx));
0069    while isempty(hp) && it < 2
0070        it = it+1;
0071        %if it > 1 keyboard, end
0072        for k=1:length(dx)
0073            nf(k)=nf(k)+calc_noise_figure( imdl, 10^-dx(k), 10 );
0074        end
0075        log_nf = log10(nf/it);
0076        p= polyfit( dx, log_nf-log10(dNF), 3);
0077        hp = roots(p);
0078        %eliminate complex roots
0079        hp = hp(imag(hp)==0);
0080        %eliminate out of range roots
0081        hp = hp( hp<max(dx) & hp>min(dx) );  %USE if poly>1
0082        % pick the root with the smallest derivative
0083        if numel(hp) >1
0084            p2 = p.*(numel(p)-1:-1:0); p2(end) = [];
0085            [jnk, pos] = min(abs(polyval(p2,hp)));
0086            hp = hp(pos);
0087        end
0088    end
0089    if isempty(hp)
0090        %fallback
0091        [jnk,idx] = min(abs(log_nf-log10(dNF)));
0092        hp = dx(idx);
0093    end
0094 
0095 function do_unit_test
0096   do_unit_test1
0097   do_unit_test2
0098 
0099 function do_unit_test1
0100   inv2d=  mk_common_model('b2c',16);
0101   inv2d.hyperparameter.func = @choose_noise_figure;
0102   inv2d.hyperparameter.noise_figure= 0.5;
0103   inv2d.hyperparameter.tgt_elems= 1:4;
0104   HP = choose_noise_figure(inv2d);
0105   unit_test_cmp('test hp',HP,0.039338520257055,1e-12);
0106   inv2d.hyperparameter.noise_figure= 2.0;
0107   inv2d.hyperparameter.tgt_elems= 1:16;
0108   HP = choose_noise_figure(inv2d);
0109   unit_test_cmp('test hp',HP,0.008265463857691,1e-12);
0110 
0111 function do_unit_test2
0112   inv2d=  mk_common_model('b2c',16);
0113   inv2d.hyperparameter.func = @choose_noise_figure;
0114   inv2d.hyperparameter.noise_figure= 0.5;
0115   img = mk_image(inv2d); vh = fwd_solve(img);
0116   img.elem_data(1:4)=1.1;vi = fwd_solve(img);
0117   inv2d.hyperparameter.tgt_data.meas_t1 = vh;
0118   inv2d.hyperparameter.tgt_data.meas_t2 = vi;
0119   HP = choose_noise_figure(inv2d);
0120   unit_test_cmp('test hp',HP,0.039347649458332,1e-12);