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 4838 2015-03-30 07:41:23Z aadler $
0021 
0022 reqNF= inv_model.hyperparameter.noise_figure;
0023 try % remove the value field
0024    inv_model.hyperparameter = rmfield(inv_model.hyperparameter,'value');
0025 end
0026 
0027 copt.boost_priority = 3;
0028 HP = eidors_cache(@HP_for_NF_search,{reqNF,inv_model},copt);
0029 
0030 
0031 
0032 function [HP,NF,SE] = HP_for_NF_search(dNF,imdl)
0033    llv = eidors_msg('log_level');
0034    if llv==3; eidors_msg('log_level',2); end % HIDE A LOT OF MESSAGES
0035 
0036    hp= search1(dNF, imdl, 1);
0037 
0038    dx= hp-linspace(-0.7,0.7,5);
0039    hp= search2(dNF, imdl, dx);
0040 
0041    dx= hp-linspace(-0.2,0.2,5);
0042    hp= search2(dNF,imdl, dx);
0043 
0044    dx= hp-linspace(-0.1,0.1,5); 
0045    hp= search2(dNF,imdl, dx);
0046 
0047    dx= hp-linspace(-0.05,0.05,21);
0048    hp= search2(dNF,imdl, dx);
0049 
0050    HP= 10^-hp;
0051    eidors_msg('log_level',llv); %% Reset log_level
0052   
0053 function hp= search1(dNF, imdl, hp)
0054   [NF,SE]=calc_noise_figure( imdl, 10^(-hp));
0055    if     NF+3*SE < dNF; dir = 1;
0056    elseif NF-3*SE > dNF; dir = -1;
0057    else   dir = 0; end
0058    while  dir*NF+3*SE < dir*dNF %>
0059      hp= hp+0.5*dir;
0060      [NF,SE]=calc_noise_figure( imdl, 10^(-hp));
0061    end
0062 
0063 function hp=search2(dNF, imdl, dx)
0064    hp = [];
0065    it = 0;
0066    nf = zeros(1,length(dx));
0067    while isempty(hp) && it < 2
0068        it = it+1;
0069        %if it > 1 keyboard, end
0070        for k=1:length(dx)
0071            nf(k)=nf(k)+calc_noise_figure( imdl, 10^-dx(k), 10 );
0072        end
0073        log_nf = log10(nf/it);
0074        p= polyfit( dx, log_nf-log10(dNF), 3);
0075        hp = roots(p);
0076        %eliminate complex roots
0077        hp = hp(imag(hp)==0);
0078        %eliminate out of range roots
0079        hp = hp( hp<max(dx) & hp>min(dx) );  %USE if poly>1
0080        % pick the root with the smallest derivative
0081        if numel(hp) >1
0082            p2 = p.*(numel(p)-1:-1:0); p2(end) = [];
0083            [jnk, pos] = min(abs(polyval(p2,hp)));
0084            hp = hp(pos);
0085        end
0086    end
0087    if isempty(hp)
0088        %fallback
0089        [jnk,idx] = min(abs(log_nf-log10(dNF)));
0090        hp = dx(idx);
0091    end