Proper way to add noise to signal

Question

In many areas I have found that while adding noise, we mention some specification like zero mean and variance. I need to add AWGN, colored noise, uniform noise of varying SNR in Db. The following code shows the way how I generated and added noise. I am aware of the function awgn() but it is a kind of black box thing without knowing how the noise is getting added. So, can somebody please explain the correct way to generate and add noise. Thank you

SNR = [-10:5:30]; %in Db
snr = 10 .^ (0.1 .* SNR);

for I = 1:length(snr)
    noise = 1 / sqrt(2) * (randn(1, N) + 1i * randn(1, N));
    u = y + noise .* snr(I);
end

Answer 1

Most answers here forget that SNR is specified in decibels. Therefore, you shouldn't encounter 'division by 0' error, because you should really divide by 10^(targetSNR/10) which is never negative nor zero for real targetSNR.

Answer 2

This 'should not divide by 0' problem could be easily solved if you add a condition to check if targetSNR is 0 and do these only if it is not 0. When your target SNR is 0, it means it's pure noise.

function out_signal = addAWGN(signal, targetSNR)
sigLength = length(signal); % length
awgnNoise = randn(size(signal)); % orignal noise
pwrSig = sqrt(sum(signal.^2))/sigLength; % signal power
pwrNoise = sqrt(sum(awgnNoise.^2))/sigLength; % noise power
if targetSNR ~= 0
   scaleFactor = (pwrSig/pwrNoise)/targetSNR; %find scale factor
   awgnNoise = scaleFactor*awgnNoise; 
   out_signal = signal + awgnNoise; % add noise
else
   out_signal = awgnNoise; % noise only
end

Answer 3

I'm adding another answer since it strikes me that Steven's is not quite correct and Horchler's suggestion to look inside function awgn is a good one.

Either MATLAB or Octave (in the communications toolbox) have a function awgn that adds (white Gaussian) noise to attain a desired signal-to-noise power level; the following is the relevant portion of the code (from the Octave function):

  if (meas == 1)  % <-- if using signal power to determine appropriate noise power
    p = sum( abs( x(:)) .^ 2) / length(x(:));
    if (strcmp(type,"dB"))
      p = 10 * log10(p);
    endif
  endif

  if (strcmp(type,"linear"))
    np = p / snr;
  else   % <-- in dB
    np = p - snr;
  endif

  y = x + wgn (m, n, np, 1, seed, type, out);

As you can see by the way p (the power of the input data) is computed, the answer from Steven does not appear to be quite right.

You can ask the function to compute the total power of your data array and combine that with the desired s/n value you provide to compute the appropriate power level of the added noise. You do this by passing the string "measured" among the optional inputs, like this (see here for the Octave documentation or here for the MATLAB documentation):

     y = awgn (x, snr, 'measured')

This leads ultimately to meas=1 and so meas==1 being true in the code above. The function awgn then uses the signal passed to it to compute the signal power, and from this and the desired s/n it then computes the appropriate power level for the added noise.

As the documentation further explains

By default the snr and pwr are assumed to be in dB and dBW respectively. This default behavior can be chosen with type set to "dB". In the case where type is set to "linear", pwr is assumed to be in Watts and snr is a ratio.

This means you can pass a negative or 0 dB snr value. The result will also depend then on other options you pass, such as the string "measured".

For the MATLAB case I suggest reading the documentation, it explains how to use the function awgn in different scenarios. Note that implementations in Octave and MATLAB are not identical, the computation of noise power should be the same but there may be different options.

And here is the relevant part from wgn (called above by awgn):

  if (strcmp(type,"dBW"))
    np = 10 ^ (p/10);
  elseif (strcmp(type,"dBm"))
    np = 10 ^((p - 30)/10);
  elseif (strcmp(type,"linear"))
    np = p;
  endif

  if(!isempty(seed))
    randn("state",seed);
  endif

  if (strcmp(out,"complex"))
    y = (sqrt(imp*np/2))*(randn(m,n)+1i*randn(m,n)); % imp=1 assuming impedance is 1 Ohm
  else
    y = (sqrt(imp*np))*randn(m,n);
  endif

If you want to check the power of your noise (np), the awgn and awg functions assume the following relationships hold:

  np = var(y,1);        % linear scale
  np = 10*log10(np);    % in dB 

where var(...,1) is the population variance for the noise y.

Answer 4

You can use randn() to generate a noise vector 'awgnNoise' of the length you want. Then, given a specified SNR value, calculate the power of the orignal signal and the power of the noise vector 'awgnNoise'. Get the right amplitude scaling factor for the noise vector and just scale it.

The following code is an example to corrupt signal with white noise, assuming input signal is 1D and real valued.

function out_signal = addAWGN(signal, targetSNR)
sigLength = length(signal); % length
awgnNoise = randn(size(signal)); % orignal noise
pwrSig = sqrt(sum(signal.^2))/sigLength; % signal power
pwrNoise = sqrt(sum(awgnNoise.^2))/sigLength; % noise power

scaleFactor = (pwrSig/pwrNoise)/targetSNR; %find scale factor
awgnNoise = scaleFactor*awgnNoise; 
out_signal = signal + awgnNoise; % add noise

Be careful about the sqrt(2) factor when you deal with complex signal, if you want to generate the real and imag part separately.