理解FFT输出

Question

我需要帮助理解DFT / FFT计算的输出结果。

我是一名有经验的软件工程师，需要解释一些智能手机加速度计读数，例如查找主频率。不幸的是，我在十五年前大学时睡了大部分的电子工程课程，但我最近几天一直在阅读DFT和FFT（显然收效甚微）。

请不要回答“去上一门电子工程课程”。如果我的雇主愿意付钱，我实际上正在计划这样做。 :)

所以这就是我的问题：

我以32 Hz的采样率捕获了一个信号。 这是32个点的1秒样本，我在Excel中将其绘制成图表。

然后我从哥伦比亚大学获取了一些用Java编写的FFT代码（在遵循“Java中可靠且快速的FFT”的建议后）。

此程序的输出如下。 我相信它正在运行就地FFT，因此它在输入和输出中重复使用相同的缓冲区。

Before: 

Re: [0.887  1.645  2.005  1.069  1.069  0.69  1.046  1.847  0.808  0.617  0.792  1.384  1.782  0.925  0.751  0.858  0.915  1.006  0.985  0.97  1.075  1.183  1.408  1.575  1.556  1.282  1.06  1.061  1.283  1.701  1.101  0.702  ]

Im: [0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ]

After: 

Re: [37.054  1.774  -1.075  1.451  -0.653  -0.253  -1.686  -3.602  0.226  0.374  -0.194  -0.312  -1.432  0.429  0.709  -0.085  0.0090  -0.085  0.709  0.429  -1.432  -0.312  -0.194  0.374  0.226  -3.602  -1.686  -0.253  -0.653  1.451  -1.075  1.774  ]

Im: [0.0  1.474  -0.238  -2.026  -0.22  -0.24  -5.009  -1.398  0.416  -1.251  -0.708  -0.713  0.851  1.882  0.379  0.021  0.0  -0.021  -0.379  -1.882  -0.851  0.713  0.708  1.251  -0.416  1.398  5.009  0.24  0.22  2.026  0.238  -1.474  ]

所以，此时我无法理解输出结果。我理解DFT的概念，例如实部是组成余弦波的振幅，虚部是组成正弦波的振幅。我也可以按照这本伟大的书《科学家与工程师数字信号处理指南》中的图表进行操作：

我的具体问题如下：

1. 从FFT的输出中，我如何找到“最常出现的频率”？这是分析加速度计数据的一部分。我应该阅读实数（余弦）数组还是虚数（正弦）数组？

2. 我有一个32点的时间域输入。FFT的输出不应该是一个16个元素的实数组和一个16个元素的虚数组吗？为什么程序会给我两个大小都为32的实数和虚数数组输出？

3. 关于上一个问题，我如何解析输出数组中的索引？考虑到我的32个采样值采样频率为32 Hz，我的理解是16个元素的数组输出应该平均分布在采样率的1/2（即32 Hz / 2 = 16 Hz），那么每个数组元素是否表示（32 Hz * 1/2）/16 = 1Hz？

4. 为什么FFT输出具有负值？我认为这些值代表正弦波的振幅。例如，Real [3]的输出为-1.075应该意味着频率为3的余弦波振幅为-1.075。对吗？振幅怎么可能是负数？

- stackoverflowuser2010

---

你想从加速度计读数中计算什么：速度、距离？加速度计读数的噪声遵循高斯分布，我无法看出拟合正弦波如何能解决这个问题。 - Ali

应该移除Java标签，因为它比特定的编程语言更为通用。 - user3791372

看哥伦比亚大学的源代码，它根本不高效。这是一个简单、未优化的 Cooley-Tucky 实现，并且比特反转是手动完成的，而不是使用现有的库函数。 - Mark Jeronimus

@MarkJeronimus：你能推荐一个高效的Java FFT实现吗？如果我没记错的话，我选择哥伦比亚大学的代码是因为FFTW库在Android智能手机上运行过于复杂。 - stackoverflowuser2010

我发现了一些零散的“优化”实现，但它们基本上是每个N大小一个算法，因此如果您需要一系列大小，则需要所有这些例程。在实践中，我主要使用英特尔集成性能原语（是的，从Java通过JNA），但那是非免费的。在家里，我使用基本上与您链接的相同算法，但是在2005年使用教科书从头编写。它只是FFT（快速傅里叶变换），没有什么“快速”的东西可以证明名为“快速FFT”。 - Mark Jeronimus

Answer 1

1. 你不应该寻找复数的实部或虚部（也就是你的实数和虚数数组），而应该寻找频率的幅度，其定义为sqrt（real * real + imag * imag）。这个数字总是正的。现在你只需要寻找最大值即可（忽略数组中的第一个条目。那是你的直流偏移量，并不携带任何频率相关信息）。

2. 你获得32个实数和32个虚数输出，因为你使用的是复合复合FFT。请记住，你已经通过用零虚部扩展它将32个样本转换为64个值（或32个复杂值）。这导致了对称的FFT输出，其中频率结果出现两次。一次准备在输出0到N / 2中使用，一次镜像在输出N / 2到N中。在你的情况下，最容易的方法是简单地忽略输出N / 2到N。你不需要它们，它们只是你计算FFT时的产物。

3. 频率到fft-bin方程式为(bin_id * freq/2) / (N/2)，其中freq是你的采样频率（即32 Hz），N是你的FFT大小。在你的情况下，这简化为每个bin 1 Hz。bins N/2到N代表负频率（奇怪的概念，我知道）。对于你的情况，它们不包含任何重要信息，因为它们只是前N / 2个频率的镜像。

4. 每个bin的实部和虚部组成一个复数。如果实部和虚部是负数，而频率本身的幅度是正的（请参见我对问题1的答案），那么这是可以接受的。我建议您阅读有关复数的资料。解释它们工作原理（以及它们为什么有用）超出了在单个stackoverflow问题中可以解释的范围。

注意：您可能还想了解自相关是什么，以及如何使用它来找到信号的基本频率。我有一种感觉，这正是你想要的。

- Nils Pipenbrinck

---

关于1：我看到了这个Matlab页面，显示了一个频谱(http://www.mathworks.com/help/techdoc/ref/fft.html)。在那个页面上有一个标题为“y(t)的单边幅度谱”的图表。它是否绘制了您建议的频率幅度，即sqrt(real^2 + img^2)？ 关于3：我仍然不明白每个bin 2Hz的结果。在我的情况下，N=32且freq=32，对吗？因此有N/2=32/2=16个bin，最高频率（奈奎斯特）为freq/2=32/2=16 Hz，导致每16个bin 16 Hz，每个bin 1 Hz？ - stackoverflowuser2010

是的，这个图表展示了频谱的幅度 - |Y(f)|。绝对值条表示幅度。箱子宽度 = 采样率 / FFT 大小。您的采样率为32赫兹，FFT大小为32。是的，您关于箱子宽度的想法是正确的！ - Matt Montag

Answer 2

你已经得到了一些不错的答案，但我想补充一点，你在进行FFT之前需要对时间域数据应用一个窗函数，否则你将会得到一些不好的谱线，这是由于频谱泄漏所导致的。请参考窗口函数和频谱泄漏。

- Paul R

Answer 3

1) 在实数数组中寻找除第一个（即DC分量）以外具有最高值的指数。为了获得有意义的结果，您可能需要比32 Hz更高的采样率和更大的窗口大小。

2) 两个数组的后一半是前一半的镜像。例如，请注意实数数组的最后一个元素（1.774）与第二个元素（1.774）相同，并且虚数数组的最后一个元素（1.474）是第二个元素的负数。

3) 在32 Hz的采样率下，您可以拾取的最大频率为16 Hz（奈奎斯特极限），因此每个步骤为2 Hz。如前所述，请记住第一个元素是0 Hz（即DC偏移）。

4) 当然，负振幅是有道理的。这只是意味着信号被“翻转”了——标准FFT基于余弦函数，通常在t = 0时的值为1，因此，在时间t = 0时值为-1的信号将具有负振幅。

- user149341

---

谢谢回复。(1)您的意思是我可以忽略虚数（正弦）数组，如果是这样，为什么？毕竟正弦分量必须很重要吧？(2)这种镜像现象是为什么？这只是FFT算法的结果吗？大多数人都会忽略镜像的一半吗？(3)您是如何计算2Hz的步长的？我了解16Hz的奈奎斯特极限，所以如果有16个（非镜像）数组元素，则每个元素必须是16 Hz / 16 = 1 Hz吗？(4)要找到主频率，我只需要在输出数组中取振幅值的绝对值吗？ - stackoverflowuser2010

你不应该在真实的数组中寻找最高值，也不能忽略正弦/I数组。相反，你需要考虑复合复杂向量的幅度。镜像出现是因为一半的输入（即I数组）都是零，所以结果自由度减半。如果数据严格为实数，则可以忽略它。 - hotpaw2

(3) 对我来说，步长为2Hz的值目前仍然是隐含的。 我们有16个箱子，由长度为16的数组表示。我们需要描述从0Hz到16Hz的所有频率。我假设每个箱子描述了该范围的一部分，不是吗？ - krafter

@krafter 我认为它被减半是因为你无法从单个值中推断出频率（因为没有重复）。 - JVE999

Answer 4

请注意，“最常出现的频率”可能会分散到多个FFT bin中，即使使用了窗口函数。因此，您可能需要使用更长的窗口、多个窗口或插值来更好地估计任何频谱峰值的频率。

- hotpaw2