高斯函数的傅里叶变换不是高斯函数？这是错误的！- Python

Question

我尝试使用Numpy的fft函数，但是当我给函数一个简单的高斯函数时，该高斯函数的fft不是高斯函数，它接近但是被分成两半，每半位于x轴的两端。

我计算的高斯函数是 y = exp(-x^2)

这是我的代码：

from cmath import *
from numpy import multiply
from numpy.fft import fft
from pylab import plot, show

""" Basically the standard range() function but with float support """
def frange (min_value, max_value, step):
    value = float(min_value)
    array = []
    while value < float(max_value):
        array.append(value)
        value += float(step)
    return array


N = 256.0 # number of steps
y = []
x = frange(-5, 5, 10/N)

# fill array y with values of the Gaussian function   
cache = -multiply(x, x)
for i in cache: y.append(exp(i))

Y = fft(y)

# plot the fft of the gausian function
plot(x, abs(Y))
show()

结果并不完全正确，因为高斯函数的快速傅里叶变换应该是一个高斯函数本身...

- chutsu

---

你可能想要查看numpy.arange()这个函数。 - Sven Marnach

Answer 1

np.fft.fft 返回的结果遵循所谓的“标准顺序”(来自文档)：

如果 A = fft(a, n)，那么A[0] 包含零频率项（信号的平均值），对于实数输入它始终是纯实数。然后，A[1:n/2] 包含正频率项，A[n/2+1:] 包含负频率项，按照逐渐变小的负频率排序。

函数np.fft.fftshift将结果重新排列为大多数人期望的顺序（也适合绘图）:

例程np.fft.fftshift（A）转换和它们的频率，以便将零频率分量放在中间...

因此，使用np.fft.fftshift:

import matplotlib.pyplot as plt
import numpy as np

N = 128
x = np.arange(-5, 5, 10./(2 * N))
y = np.exp(-x * x)
y_fft = np.fft.fftshift(np.abs(np.fft.fft(y))) / np.sqrt(len(y))
plt.plot(x,y)
plt.plot(x,y_fft)
plt.show()

- unutbu

---

谢谢，我在寻找fftshift，但是我找不到它。 - Steve Tjoa

为什么要有sqrt(2N)？ - rainman

@omehoque：定义DFT有许多不同的约定。它们在常数上都是相同的。不同的群体（例如工程师、物理学家、数学家）倾向于不同的约定。上面使用的约定受到数学家的青睐，因为它使DFT和逆DFT公式对称。这使得Parseval定理特别美丽：(y**2).sum()等于(y_fft**2).sum()。它也是绘图的一个好选择，因为它使y和y_fft大致相同大小。 - unutbu

NumPy 定义了 DFT 这样。将该公式除以 1/sqrt(N) 就是 这个答案 所称的单位缩放约定。之所以在上面使用 1/sqrt(2*N) 是因为 len(y) 是 2*N，而不是 N。 - unutbu

我找到了一个关于幅度的相关问题：https://scicomp.stackexchange.com/q/11233/32010。有一个答案指出了采样频率（在这种情况下是10./2/128）。但我不明白它从哪里来。 - Jason

更新自己（我希望他们允许编辑超过5分钟的评论）：unutbu评论中的第二个链接谈论了采样频率。 - Jason

Answer 2

你的结果与高斯分布完全不同，甚至不能分成两半。

要得到你期望的结果，你需要将自己的高斯分布放置在索引0处，结果也将以此为中心位置。尝试以下代码：

from pylab import *
N = 128
x = r_[arange(0, 5, 5./N), arange(-5, 0, 5./N)]
y = exp(-x*x)
y_fft = fft(y) / sqrt(2 * N)
plot(r_[y[N:], y[:N]])
plot(r_[y_fft[N:], y_fft[:N]])
show()

绘图命令将数组分成两半并交换它们以获得更好的图片。

- Sven Marnach

---

我是否需要将数组分成两半，然后对每个信号进行FFT？或者这只适用于高斯函数，如果是这样，为什么？ - chutsu

@chutsu：我不知道你想做什么。只要你不在意高斯傅里叶变换的图像看起来如何，那么放置原点的位置可能并不重要。 - Sven Marnach

Answer 3

该内容显示的中心（即平均值）在系数索引为零位置。因此，右半部分出现在左侧，左半部分出现在右侧。

编辑：请查看以下代码：

import scipy
import scipy.signal as sig
import pylab
x = sig.gaussian(2048, 10)
X = scipy.absolute(scipy.fft(x))
pylab.plot(x)
pylab.plot(X)
pylab.plot(X[range(1024, 2048)+range(0, 1024)])

最后一行将从向量的中心开始绘制X，然后环绕到开头。

- Steve Tjoa

---

那么我该如何更改上述内容以解决这个问题？这只影响高斯函数吗？其他波形/信号呢？如果我的数学不太好，对不起... - chutsu

Answer 4

继Sven Marnach的回答之后，一个更简单的版本如下：

from pylab import *
N = 128

x = ifftshift(arange(-5,5,5./N))

y = exp(-x*x)
y_fft = fft(y) / sqrt(2 * N)

plot(fftshift(y))
plot(fftshift(y_fft))

show()

这将产生与上面相同的图形。

关键是（这对我来说似乎很奇怪），NumPy假定的数据排序方式---在频率和时间域中---是将“零”值放在第一位。这不是我从其他FFT实现（如C中的FFTW3库）中所期望的。

这在unutbu和Steve Tjoa的答案中稍微有些欺骗性，因为他们在绘制之前取FFT的绝对值，从而抹去了由于未使用时间上的“标准顺序”而导致的相位问题。

- Bernard

---

我知道这篇文章已经8年了，但这个答案是我找到的唯一可接受的答案。其他所有的答案都通过取绝对值来掩盖相移问题，甚至在被接受的答案中也是如此（比如在这篇文章中）。还可以参见这里：https://stackoverflow.com/questions/56154992/fourier-transform-of-a-gaussian-function-in-python、https://scicomp.stackexchange.com/questions/11233/amplitude-of-discrete-fourier-transform-of-gaussian-is-incorrect，以及许多无法放在评论中的内容。感谢您不回避这个问题。取绝对值就是完全错误的。 - snar

Answer 5

傅里叶变换隐含地无限重复，因为它是一个信号的变换，该信号隐含地无限重复。请注意，当您传递要变换的y时，实际上并没有提供x值，因此被转换的高斯函数是以0和256之间的中位数128为中心的一个。

还要记住，f(x)的翻译是F(x)的相位变化。

- Phil H