numpy - 如何通过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" ，据说 sHat无需 s的傅立叶变换即可直接生成如理论所示 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)的实部您会看到转换后的噪声具有更高的值。

相反，如果您填写 gaussComplex对于标准偏差为 1 的随机值，逆变换的噪声将具有低得多的标准偏差。这就是你观察到的。

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

gaussComplex *= np.sqrt(N/2)

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

绘图 s和 sNew应用正确的缩放会导致类似的结果:


结果并不完全相同(如预期)，但噪声在相同范围内变化 - 两者的标准偏差都为 ~1。