Правильный способ добавления шума к сигналу

Ответ 1

Я добавляю еще один ответ, поскольку мне кажется, что Стивен не совсем прав, и предложение Horchler заглянуть внутрь функции awgn является хорошим.

Либо MATLAB, либо Octave (в панели инструментов связи) имеют функцию awgn, которая добавляет (белый гауссовский) шум для достижения желаемого уровня мощности сигнала к шуму; следующим является часть кода (из функции Octave):

  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);

Как вы можете видеть, кстати, p (мощность входных данных) вычисляется, ответ от Стивена не кажется совершенно правильным.

Вы можете попросить функцию вычислить общую мощность вашего массива данных и объединить это с желаемым значением s/n, которое вы предоставляете, для вычисления соответствующего уровня мощности добавленного шума. Вы делаете это, передавая строку "измеряется" среди дополнительных входов, например (см. здесь для документации Octave или здесь для документации MATLAB):

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

В конечном итоге это приводит к meas=1, и поэтому meas==1 является истинным в коде выше. Функция awgn затем использует переданный ей сигнал для вычисления мощности сигнала, и из этого и желаемого s/n затем вычисляет соответствующий уровень мощности для добавленного шума.

Как поясняется в документации,

По умолчанию snr и pwr считаются равными дБ и дБВт соответственно. Это поведение по умолчанию можно выбрать с типом, установленным в "ДБ". В случае, когда тип задан как "линейный", предполагается, что pwr в Ваттах и ​​snr - отношение.

Это означает, что вы можете передать отрицательное или 0 дБ значение snr. Результат также будет зависеть от других параметров, которые вы передаете, например, строка "измеренная".

В случае MATLAB я предлагаю прочитать документацию , в нем объясняется, как использовать функцию awgn в разных сценариях. Обратите внимание, что реализации в Octave и MATLAB не идентичны, вычисление мощности шума должно быть одинаковым, но могут быть разные варианты.

И вот соответствующая часть из wgn (вызывается выше 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

Если вы хотите проверить мощность вашего шума (np), функции awgn и awg предполагают следующие отношения:

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

где var(...,1) - дисперсия совокупности шумов y.

Ответ 2

Большинство ответов здесь забывают, что SNR указывается в децибелах. Поэтому вам не следует сталкиваться с ошибкой "деление на 0", потому что вы действительно должны делиться на 10^(targetSNR/10), которая никогда не является отрицательной или нулевой для реального targetSNR.

Ответ 3

Вы можете использовать randn() для генерации шумового вектора 'awgnNoise' требуемой длины. Затем, учитывая заданное значение SNR, вычислите мощность сигнала orignal и мощность вектора шума "awgnNoise". Получите правильный амплитудный масштабный коэффициент для вектора шума и просто масштабируйте его.

Следующий код - пример испорченного сигнала с белым шумом, если входной сигнал 1D и действительный.

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

Будьте осторожны с коэффициентом sqrt (2), когда вы имеете дело со сложным сигналом, если хотите генерировать реальную и часть изображения отдельно.

Ответ 4

Задача "не делить на 0" может быть легко решена, если вы добавите условие для проверки, является ли targetSNR равным 0, и делайте это только в том случае, если оно не равно 0. Когда целевое SNR равно 0, это означает, что это чистый шум.

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