add_noise

 ADD_NOISE: Add a given SNR to EIDORS data
 v1_w_noise = add_noise( SNR, v1, v2, options)

 Usage:
  add_noise( SNR, v1 )           - add noise to v1 where signal = vv
  add_noise( SNR, v1, v2)        - add noise to v1 where signal = v1 - v2
  add_noise( SNR, v1, v2,'norm') - add noise to v1 where signal = (v1-v2)/v2

 SNR is defined in terms of power SNR =  || signal || / || noise ||

0001 function vv = add_noise( SNR, v1, v2, options)
0002 % ADD_NOISE: Add a given SNR to EIDORS data
0003 % v1_w_noise = add_noise( SNR, v1, v2, options)
0004 %
0005 % Usage:
0006 %  add_noise( SNR, v1 )           - add noise to v1 where signal = vv
0007 %  add_noise( SNR, v1, v2)        - add noise to v1 where signal = v1 - v2
0008 %  add_noise( SNR, v1, v2,'norm') - add noise to v1 where signal = (v1-v2)/v2
0009 %
0010 % SNR is defined in terms of power SNR =  || signal || / || noise ||
0011 
0012 % (C) 2010 Andy Adler. License: GPL v2 or v3. $Id: add_noise.html 2819 2011-09-07 16:43:11Z aadler $
0013 
0014 if isstr(SNR) && strcmp(SNR,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 if nargin>=2; try; v1 = v1.meas; end; end
0017 if nargin>=3; try; v2 = v2.meas; end; end
0018    
0019 if nargin == 2;
0020    signal = v1;
0021 elseif nargin==3;
0022    signal = v1 - v2;
0023 elseif strcmp( options, 'norm' )
0024    signal = (v1 - v2) ./ v1;
0025 else
0026    error('add_noise: input arguments not understood');
0027 end
0028    
0029 noise = randn(size(signal));
0030 
0031 % SNR = norm(signal)/norm(noise)
0032 % so  scale norm(noise) by norm(signal)/SNR
0033 
0034 noise = noise *norm(signal) / norm(noise) / SNR;
0035 
0036 vv = eidors_obj('data','from add_noise');
0037 vv.meas = v1 + noise;
0038 
0039 function do_unit_test
0040     v1 = 2.0*ones(10,1);
0041     v2 = 2.1*ones(10,1);
0042     
0043     v0 = add_noise( 2, v1);
0044     SNR_test(2, v0.meas - v1, v1);
0045 
0046     v0 = add_noise(.1, v1, v2);
0047     SNR_test(.1, v0.meas - v1, v1 - v2);
0048 
0049     v0 = add_noise(.3, v1, v2, 'norm' );
0050     SNR_test(.3, v0.meas - v1, (v1 - v2)./v1);
0051 
0052 function SNR_test(SNRspec, noi, sig)
0053     SNR = norm(sig)/norm(noi);
0054     if abs(SNR - SNRspec) < .001; disp('ok');
0055     else;                         disp('fail'); end