calc_noise_params

 params = GREIT_noise_params(imdl, homg_voltage, sig_voltage)
  params(1,:) = Noise Figure = SNR(image) / SNR(data)

  see also: eval_GREIT_fig_merit or using test_performance

0001 function params = calc_noise_params(imdl, vh, vi )
0002 % params = GREIT_noise_params(imdl, homg_voltage, sig_voltage)
0003 %  params(1,:) = Noise Figure = SNR(image) / SNR(data)
0004 %
0005 %  see also: eval_GREIT_fig_merit or using test_performance
0006 
0007 % (C) 2008 Andy Adler. Licensed under GPL v2 or v3
0008 % $Id: calc_noise_params.m 6036 2019-12-30 19:17:42Z aadler $
0009 
0010 % NOTE THAT WE ASSUME A LINEAR ALGORITHM FOR THIS MEASURE
0011 
0012 warning('EIDORS:deprecated','CALC_NOISE_PARAMS is deprecated as of 02-Jul-2012. Use CALC_NOISE_FIGURE instead.');
0013 
0014 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'), do_unit_test, return, end
0015 
0016 if 0 % old code with random noise
0017      % Keep this to validate while we test it
0018    Nnoise = 1000;
0019    noise = 0.01*std(vh)*randn(size(vh,1),Nnoise);
0020    vhn= vh*ones(1,Nnoise) + noise;
0021 else % use independent noise model on each channel
0022    noise = 0.01*std(vh)*speye(size(vh,1));
0023    vhn= vh*ones(1,size(vh,1)) + noise;
0024 end
0025 
0026 signal_y = calc_difference_data( vh, vi,  imdl.fwd_model);
0027 noise_y  = calc_difference_data( vh, vhn, imdl.fwd_model);
0028 
0029 signal_x = inv_solve(imdl, vh, vi);  signal_x = signal_x.elem_data;
0030 noise_x  = inv_solve(imdl, vh, vhn); noise_x  = noise_x.elem_data;
0031 
0032 try 
0033     VOL = get_elem_volume(imdl.rec_model);
0034 catch
0035     VOL = get_elem_volume(imdl.fwd_model);
0036 end
0037 VOL = spdiags(VOL,0, length(VOL), length(VOL));
0038 
0039 signal_x = VOL*signal_x;
0040 noise_x = VOL*noise_x;
0041 
0042 signal_x = mean(abs(signal_x),1);
0043 noise_x  = mean(std(noise_x));
0044 snr_x = signal_x / noise_x;
0045 
0046 signal_y = mean(abs(signal_y),1);
0047 noise_y  = mean(std(noise_y)); 
0048 snr_y = signal_y / noise_y;
0049 
0050 params= [snr_y(:)./snr_x(:)]';
0051 
0052 % NOTES on the calculations: AA - Feb 20, 2012
0053 % SNR = mean(abs(x)); VAR =
0054 % ym= E[y]
0055 % Sy= E[(y-ym)*(y-ym)'] = E[y*y'] - ym*ym'
0056 % ny = sqrt(trace(Sy))
0057 % xm= E[x]  = E[R*y] = R*E[y] = R*ym
0058 % Sx= E[(x-xm)*(x-xm)'] = E[x*x'] - xm*xm'
0059 %   = E[R*ym*ym'*R'] = R*E[ym*ym']*R' = R*Sy*R'
0060 % nx = sqrt(trace(Sx))
0061 %
0062 % signal = mean(abs(x));
0063 %
0064 % In this case, these are exactly the same:
0065 %    noise_x  = mean(std(noise_x));
0066 %    noise_x  = sqrt(mean(noise_x.^2,2));
0067 
0068 function do_unit_test
0069 ll = eidors_msg('log_level',1);
0070 imdl = mk_common_model('a2t2',16); test2(imdl);
0071 imdl = mk_common_model('d2t2',16); test2(imdl);
0072 imdl = mk_common_model('a2c2',16); test2(imdl);
0073 imdl = mk_common_model('d2c2',16); test2(imdl);
0074 test1; % can we deal with c2f ?
0075 ll = eidors_msg('log_level',ll);
0076 
0077 function test1
0078 % big model with c2f and supposedly an NF of 0.5
0079 fmdl = mk_library_model('pig_23kg_16el');
0080 [fmdl.stimulation fmdl.meas_select] = mk_stim_patterns(16,1,'{ad}','{ad}');
0081 fmdl = mdl_normalize(fmdl, 1);  % Use normalized difference imaging
0082 opt.noise_figure = 0.5; opt.imgsz = [64 64];
0083 imdl = mk_GREIT_model(fmdl, 0.25, [], opt);
0084 % homogeneous measurement
0085 img = mk_image(fmdl,1);
0086 vh = fwd_solve(img);
0087 % inhomogeneous measurement
0088 select_fcn = inline('(x-0).^2+(y-0).^2+(z-0.5).^2<0.1^2','x','y','z');
0089 mfrac = elem_select(fmdl, select_fcn);
0090 img.elem_data = img.elem_data + mfrac*0.1;
0091 vi = fwd_solve(img);
0092 
0093 nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0094 
0095 % We use different measurements so naturally we don't get 0.5 here
0096 imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0097 imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0098 try
0099     % calc_noise_figure doens't support dual models
0100     nf2 = calc_noise_figure(imdl);
0101 catch
0102     nf2 = 0;
0103 end
0104 unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0105 
0106 function test2(imdl)
0107 fmdl = imdl.fwd_model;
0108 % homogeneous measurement
0109 img = mk_image(fmdl,1);
0110 vh = fwd_solve(img);
0111 % inhomogeneous measurement
0112 select_fcn = inline('(x-0).^2+(y-0).^2.^2<15^2','x','y','z');
0113 mfrac = elem_select(fmdl, select_fcn);
0114 img.elem_data = img.elem_data + mfrac*0.1;
0115 vi = fwd_solve(img);
0116 
0117 nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0118 eidors_msg(nf1,0);
0119 
0120 imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0121 imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0122 % calc_noise_figure doens't support dual models
0123 nf2 = calc_noise_figure(imdl,[],1000);
0124 unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);