(numpy) FFT数组的幅度错误(?)

Question

我正在使用numpy和matplotlib来分析我模拟输出的数据。不过，我发现一个明显的不一致之处，搞不清楚原因。具体来说：

我有一个信号，它的能量大约是a^2~1。当我用rfft进行快速傅里叶变换（FFT）并计算傅里叶空间中的能量时，结果却明显更大。为了不详细说明我的数据，这里用一个简单的正弦波作为例子：

from pylab import *
xx=np.linspace(0.,2*pi,128)
a=np.zeros(128)
for i in range(0,128):
    a[i]=sin(xx[i])
aft=rfft(a)
print mean(abs(aft)**2),mean(a**2) 

原则上，这两个数字应该是相同的（至少在数值上是这样），但我从这段代码中得到的结果是：

62.523081632 0.49609375

我尝试查看numpy.fft的文档，但没有找到相关的信息。在这里搜索后找到了以下内容，但我无法理解那里的解释：

合成信号和过滤信号之间的FFT幅度差异很大

我到底漏掉了什么/误解了什么？如果能提供任何帮助或指引，我将非常感激。

谢谢！

Answer 1

在大多数快速傅里叶变换（FFT）库中，各种离散傅里叶变换（DFT）的形式并不是< a href="http://en.wikipedia.org/wiki/Orthogonal_matrix" rel="nofollow">正交的。numpy.fft库只在进行逆变换时应用必要的归一化。

看看< a href="http://en.wikipedia.org/wiki/Discrete_Fourier_transform" rel="nofollow">维基百科对DFT的描述；逆DFT有一个1/N的项，而DFT没有这个项（这里的N是变换的长度）。要制作一个正交版本的DFT，你需要把未归一化的DFT结果乘以1/√N。在这种情况下，变换就是正交的（也就是说，如果我们把正交DFT定义为F，那么逆DFT就是F的共轭转置，或者说厄米转置）。

在你的情况下，你可以通过简单地把aft乘以1.0/√(len(a))来得到正确的答案（注意N是通过变换的长度得出的；真实的FFT会丢弃大约一半的值，所以a的长度是重要的）。

我怀疑把归一化留到最后的原因是，在大多数情况下，这并不重要，因此可以节省进行两次归一化的计算成本。实际上，快速的< a href="http://www.fftw.org" rel="nofollow">FFTW库在任何方向上都不进行归一化，完全由用户来处理。

编辑：为了更清楚，以上的解释并不完全正确。通过那种简单的缩放是无法得到正确答案的，因为在你的情况下，直流分量会被加两次，尽管1.0/√(len(a))仍然是产生单位变换的正确缩放。

Answer 2

亨利在非归一化的部分说得对，但还有一点需要补充，因为你使用的是 rfft，而不是 fft。下面的内容和他的回答是一致的：

>>> x = np.linspace(0, 2 * np.pi, 128)
>>> y = 1 - np.sin(x)
>>> fft = np.fft.fft(y)
>>> np.mean((fft * fft.conj()).real)
191.49999999999991
>>> np.mean(y**2)
1.4960937500000004
>>> fft = fft / np.sqrt(len(fft))
>>> np.mean((fft * fft.conj()).real)
1.4960937499999991

不过，如果你现在尝试用 rfft，事情就不太一样了：

>>> rfft = np.fft.rfft(y)
>>> np.mean((rfft * rfft.conj()).real)
314.58462009358772
>>> rfft /= np.sqrt(len(rfft))
>>> np.mean((rfft * rfft.conj()).real)
4.8397633860551954
65
>>> np.mean((rfft * rfft.conj()).real) / len(rfft)
4.8397633860551954

不过，下面这个是可以正常工作的：

>>> (rfft[0] * rfft[0].conj() +
...  2 * np.sum(rfft[1:] * rfft[1:].conj())).real / len(y)
1.4960937873636722

当你使用 rfft 时，得到的结果并不是你数据的完整离散傅里叶变换（DFT），而只是正半部分，因为负半部分是对称的。要计算平均值，你需要把除了直流分量以外的每个值都考虑两次，这就是最后一行代码所做的事情。