Search code examples
numpyscipygaussiannoise

How to generate random samples of Gaussian distribution directly in the frequency domain through Python: NumPy/SciPy?

One can easily draw (pseudo-)random samples from a normal (Gaussian) distribution by using, say, NumPy:

import numpy as np
mu, sigma = 0, 0.1 # mean and standard deviation
s = np.random.normal(mu, sigma, 1000)

Now, consider the Fast Fourier transform of s:

from scipy.fftpack import fft
sHat = fft(s)

Considering, presumably, generating Gaussian random numbers directly in the frequency domain might be more clever(?) (thus, efficient?) for various applications, and "the Fourier transform of white noise is white noise", and reportedly sHat can be generated directly without the Fourier-transform of s as theoretically shown herein;

Could you please advice novices (like myself) on how to implement such a useful idea? Favourably, a reference to an available implementation that I could not find in the web?

The following is my attempt to code the above theoretical explanation:

import numpy as np 
from scipy.fftpack import ifft

N = 100
gaussComplex = np.full(shape=N, dtype=complex, fill_value=0.+0.j)

mu, sigma = 0, 1
s = np.random.normal(mu, sigma, N)

iters = np.arange(N) # 0..N-1

# 0..N/2-1
for i, item in enumerate(iters[:N/2]):
    gaussComplex[i] = complex(s[i], s[i+N/2])

conjugateGaussComplex = np.conjugate(gaussComplex)

# N/2..N-1
for i, item in enumerate(iters[N/2:]):
    gaussComplex[item] = conjugateGaussComplex[N-item]

sNew = ifft(gaussComplex)

A comparison of s and sNew reveals the following as I expect that they would be the same:

plt.plot(sHat.real, 'blue')
plt.plot(s, 'red')

enter image description here

Solution

  • the Fourier transform of white noise is white noise

    This is true, but it does not mean that they are axactly the same - otherwise there would not be much point in doing the FFT.

    If you plot s and the real part of fft(s) you will see that the transformed noise has much higher values.

    Conversely, if you fill gaussComplex with random values of standard deviation 1, the inversely transformed noise will have a much lower standard deviation. This is what you observe.

    To generate the FFT of Gaussian white noise in the frequency domain you need to scale it correctly. In your case, this should do the trick: 

    gaussComplex *= np.sqrt(N/2)

    You need to scale by sqrt(N) because this is the normalization factor of the FFT. Furthermore, you also need to scale by 1/sqrt(2) because you put noise with standard deviation of 1 in both, the real and imaginary part of the FFT. (However, the absolute value should have standard deviation of 1, hence the need to divide by sqrt(1+1).)

    Plotting s and sNew with the correct scaling applied results in something like that: enter image description here

    The result is not exactly the same (as expected), but the noise varies in the same range - both have a standard deviation of ~1.