NumPy中使用FFT时的频率单位

Question

我正在使用NumPy中的FFT函数进行一些信号处理。我有一个名为signal的数组，其中每个小时有一个数据点，总共有576个数据点。我在signal上使用以下代码来查看其傅里叶变换。

t = len(signal)
ft = fft(signal,n=t)
mgft=abs(ft)
plot(mgft[0:t/2+1])

我看到两个峰，但我不确定x轴的单位是什么，即它们如何映射成小时？任何帮助都将不胜感激。

- msohail

Answer 1

通过以下关系式，给定采样率 FSample 和变换块大小 N，可以计算出频率分辨率 deltaF、采样间隔 deltaT 和总捕获时间 capT：

deltaT = 1/FSample = capT/N
deltaF = 1/capT = FSample/N

需要注意的是，FFT返回值范围为0到FSample，等同于-FSample/2到FSample/2。在你的图表中，你已经把-FSample/2到0的部分去掉了。NumPy包含一个辅助函数fftfreq可以帮助你计算这些。

对于你的值deltaT = 1 小时和N = 576，你会得到deltaF = 0.001736 cycles/hour = 0.04167 cycles/day，从-0.5 cycles/hour到0.5 cycles/hour。因此，如果你在第48个（和第528个）bin处有幅度峰值，则相应的频率分量为48*deltaF = 0.0833 cycles/hour = 2 cycles/day.

一般来说，在计算FFT之前，你应该对时间域数据使用窗口函数来减少频谱泄漏。Hann窗口几乎从不是一个坏选择。你也可以使用rfft函数来跳过输出的-FSample/2, 0部分。那么，你的代码将是：

ft = np.fft.rfft(signal*np.hanning(len(signal)))
mgft = abs(ft)
xVals = np.fft.fftfreq(len(signal), d=1.0) # in hours, or d=1.0/24 in days
plot(xVals[:len(mgft)], mgft)

- mtrw

---

@JoeKington - 谢谢！今天我在文档中找不到它，我以为我是幻觉。 - mtrw

请注意，每个FFT结果bin都有一个宽度，因此在您的情况下可能存在约半小时左右的潜在误差范围。如果您使用窗口，则可能会更大，约为1小时或更多。 - hotpaw2

也许可以使用 xVals = np.fft.rfftfreq(len(signal), d=1.0) 然后 plot(xVals, mgft) 进行操作？ - LittleLittleQ

Answer 2

fft变换的结果并不映射到HOURS，而是映射到您的数据集中包含的频率。如果有您转换后的图形，我们可以看到峰值出现的位置，这将是有益的。

由于您没有进行任何窗口处理，可能在转换缓冲区的开头会出现峰值。

- Daniel Mošmondor

Answer 3

通常情况下，FFT 的频率维度单位与输入到 FFT 的数据的采样率的维度单位相同，例如：每米、每弧度、每秒，或者在您的情况下，每小时。

频率的缩放单位，每个 FFT 结果 bin 索引，为 N / theSampleRate，具有与上述相同的维度单位，其中 N 是完整 FFT 的长度（在严格实数数据的情况下，您可能只绘制了其中一半）。

请注意，每个 FFT 结果峰值 bin 都代表一个带有非零带宽的滤波器，因此您可能需要向映射到频率值的结果点添加一些不确定性或误差边界。如果需要并适用于源数据，甚至可以使用插值估计方法。

- hotpaw2