使用 Python 进行 FFT - 意外的低频答案

【问题标题】：FFT using Python - unexpected low frequencies使用 Python 进行 FFT - 意外的低频
【发布时间】：2016-09-10 13:12:58
【问题描述】：

我仍在尝试使用 Python 中的 FFT 对此 data 进行频率分析。采样率为每分钟 1 个数据点。

我的代码是：

from scipy.fftpack import fft
df3 = pd.read_csv('Pressure - Dates by Minute.csv', sep=",", skiprows=0)
df3['Pressure FFT'] = df3['ATMOSPHERIC PRESSURE (hPa) mean'] - df3['ATMOSPHERIC PRESSURE (hPa) mean'].mean()
Pressure = df3['Pressure FFT']
Fs = 1/60
Ts = 1.0/Fs
n = len(Pressure)
k = np.arange(n)
T = n/Fs
t = np.arange(0,1,1/n) # time vector
frq = k/T # two sides frequency range
frq = frq[range(int(n/2))] # one side frequency range

Y = np.fft.fft(Pressure)/n # fft computing and normalization
Y = Y[range(int(n/2))]

fig, ax = plt.subplots(2, 1)
ax[0].plot(t,Pressure)
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Amplitude')
ax[1].plot(frq,abs(Y),'r') # plotting the spectrum
ax[1].set_xlabel('Freq (Hz)')
ax[1].set_ylabel('|Y(freq)|')

但结果给出：

所以我的问题是：

1) 为什么根本没有频率？数据明显是周期性的。

2) 为什么频谱这么低？ (0 - 0.009)

3) 也许我应该尝试不同的过滤技术？

有什么见解吗？

谢谢！！！

【问题讨论】：

FFT 函数返回的数组中的第一项具有 DC 分量，即原始数组中所有值的总和。这通常比周期性分量大几个数量级。尝试绘制 Y[1:] ，或在绘制之前执行 Y[0] = 0 ，你应该会看到你的频率出现了。
我在绘制之前尝试过 Y[0] = 0，但仍然没有频率。也许这是正常化？因为周期性行为是一天两次。
通常您希望使用对数刻度绘制 Y 轴。此外，您忘记在 FFT 之前应用合适的窗口函数。
窗口函数将有助于减少由样本窗口的起点和终点之间的不连续性引入的光谱泄漏和伪影 - 当您对 Y 轴使用对数刻度时，这些会更加明显.
看起来您有许多数据点比数据波动更接近。因此，您的所有结构都处于低频状态。我想，这就是你的答案。不过，总的来说，您要问几个关于 FFT 的基本问题。这与其说是一个编程问题，不如说是一个教程请求，在这里是题外话。

标签： python matplotlib scipy fft

【解决方案1】：

部分原因是那些频率似乎如此之低的是因为幅度绘图中的时间轴刻度均匀。如果您每60秒确实有一个样品，则X轴应在0到1690260秒之间（即20天！）。

通过眼睛，您似乎每50000秒（每天〜2秒）大约有一个小峰值，这将对应于约2x10⁻⁵Hz的频率。因此，根据X轴的尺度是多么大量的巨大，您的期限图对我来说看起来非常合理。

【讨论】：

【解决方案2】：

1) 为什么根本没有频率？数据明显是周期性的。

嗯，有频率内容，但由于其结构，它并不完全可见。尝试将绘制频谱的线从ax[1].plot(frq,abs(Y),'r') 更改为ax[1].semilogy(frq,abs(Y),'r')

这将导致：

我们现在应用了一个简单的转换来提升低值并限制高值。如需更多信息，请参阅this link。当然，移除 DC（就像您在代码的第 3 行所做的那样）也有帮助。

这看起来仍然有点模糊，但如果我们放大光谱的较低部分，我们会看到：

这显示了大约 2.3e-05 Hz 的峰值，对应于大约 12 小时。

2) 为什么频谱这么低？ (0 - 0.009)

因为您每 60 秒采样一次，因此您的采样频率为（大约）0.016 Hz。您的频谱包含 DC (0Hz) 和 0.0083Hz 之间的所有内容。更多信息请见this link

3) 也许我应该尝试不同的过滤技术？

如果您无法解决谐波，但看起来这里不需要它，您可以尝试加窗。

希望这会有所帮助。

【讨论】：