如何通过Python:NumPy/SciPy在频域直接生成高斯分布的随机样本？

Question

可以使用 NumPy 轻松地从正态（高斯）分布中抽取（伪）随机样本：

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

现在，考虑s的快速傅里叶变换：

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

考虑到，据推测，直接在频域中生成高斯随机数对于各种应用可能更聪明（？）（因此，高效？），并且"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;

请问新手（比如我）如何实现这样一个有用的想法？有利的是，对我在网络上找不到的可用实现的引用？

以下是我对上述理论解释的代码尝试：

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)

s 和 sNew 的比较揭示了以下内容，因为我预计它们是相同的：

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

Answer 1

the Fourier transform of white noise is white noise

这是事实，但这并不意味着它们完全相同 - 否则进行 FFT 就没有多大意义。

如果绘制 s 和 fft(s) 的实部，您会发现转换后的噪声具有更高的值。

相反，如果用标准差为 1 的随机值填充 gaussComplex，则经过逆变换的噪声将具有更低的标准差。这就是你观察到的。

要在频域中生成高斯白噪声的 FFT，您需要对其进行正确缩放。在您的情况下，这应该可以解决问题：

gaussComplex *= np.sqrt(N/2)

您需要按 sqrt(N) 缩放，因为这是 FFT 的归一化因子。此外，您还需要按 1/sqrt(2) 缩放，因为您在 FFT 的实部和虚部中放置了标准偏差为 1 的噪声。（但是，绝对值应该有标准偏差 1，因此需要除以 sqrt(1+1)。）

使用正确的缩放比例绘制 s 和 sNew 会产生类似这样的结果：

结果并不完全相同（如预期），但噪音在同一范围内变化 - 两者的标准差均为 ~1。

如何通过Python:NumPy/SciPy在频域直接生成高斯分布的随机样本？

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

numpy

gaussian

noise

scipy