Gaussian noise generation for a given SNR ?

Question

I am trying to add a Gaussian noise, normal distributed to a signal I have simulated (sig_noiseFree), to get a noisy signal (sig_noisy). lease see the code.

I wrote the function my self according to what I understood. Can please anyone confirm that I am doing this right ? especially the noise addition and SNR calculation ?

I have a second question : I understand that the variable "k" controls the noise level, how can I set the SNR and generate the suitable k so I get the noisy signal with the SNR wanted ? Thank you in advance.

def SNR(sig, noise,dt):
    Signal = np.sum(np.abs(np.fft.fft(sig)*dt)**2)/len(np.fft.fft(sig))
    Noise = np.sum(np.abs(np.fft.fft(noise)*dt)**2)/len(np.fft.fft(noise))
    return (10 * np.log10(Signal/Noise))


def Gauss_noise(k,sig):
     return np.random.normal(scale=k*np.max(sig), size=len(sig))

def noisy_sig(sig,k, dt):
    return np.fft.ifft(np.fft.fft(sig)*dt+np.fft.fft(np.random.normal(scale= k*np.max(sig), size=len(sig)))*dt)/dt   

def main(argv):
    k = 2e-2                
    noise_g = Gauss_noise(k,sig_noiseFree)
    sig_noisy = noisy_sig(sig_noiseFree,k, dt)

if __name__ == "__main__":
main(sys.argv)

Answer 1

Since nobody answered this question yet, i'm going to do my best to do it. Be aware i'm no expert on DSP.

I believe your first function is calculating the signal and noise average power based on each one's Fourier Transforms. So you're returning the SNR in decibels at the end.

\begin{equation} SNR_{db} = 10 * \log_{10}{\frac{P_{signal}}{P_{noise}}} \end{equation}

Your gaussian noise function generates the noise based on a scaling factor k of the signal max amplitude. Since you want to scale the amplitude of the noise based on your signal, i believe you want a relationship of:

\begin{equation} k=\frac{A_{noise}}{A_{signal}} \end{equation}

With each A meaning RMS amplitude. That can be generated used the following equation:

\begin{equation} k = \sqrt{\frac{1}{SNR_{linear}}} \end{equation}

or, to be more clear:

\begin{equation} k = \sqrt{\frac{P_{noise}}{P_{signal}}} =\sqrt{\frac{A_{noise}^2}{A_{signal}^2}} \end{equation}

If you were given the SNR in decibels and was asked to generate a noise based on it, you can use the following equation:

\begin{equation} k = \frac{1}{10^{\frac{SNR_{db}}{10}}} \end{equation}

In the third function you're generating the output signal by adding the frequency components of each signal, but if it's just an additive gaussian noise, you could just add the noise to the signal. (I'm not exactly sure on this).

You could also generate the linear SNR from your SNR in decibels, I've used this function in one of my projects once:

def linear_snr_from_db(snrdb):
    return 10.0 ** (snrdb/10.0)

I got all of the equations based on the wikipedia page for SNR and my little bit of experience.

Answer 2

It looks like one change to note and make in the code based on the theory Eduardo presented above, would be to change max (amplitude) value of signal to rms (amplitude) in your code, yes? on the line where k*np.max(sig), you will need to create python code to calculate rms amplitude and replace that with np.max(sig) in the line:

return np.random.normal(scale=k*np.max(sig), size=len(sig))