python - python 中的 FFT 无法绘制正确的频率

Question

import datetime
import shlex

import numpy as np
from scipy.fftpack import rfft, irfft, fftfreq

strWeatherFile = "data/weather.dat"

signal = []

with open(strWeatherFile, 'rt') as f:
    for line in f:
        arrDat = shlex.split(line)
        d = float(arrDat[3])
        if d < -30:  # if it's dummy data
            d = signal[len(signal) - 1]
        signal.append(d)

size = len(signal)

time   = np.linspace(1,size,size)

f_signal = rfft(signal)


import pylab as plt
plt.subplot(221)
plt.plot(time,signal)
plt.subplot(222)
plt.plot(time,f_signal)
plt.xlim(0,size)

plt.show()

这是天气数据，因此它应该显示频率为 365 的尖峰。但结果并不像预期的那样。 代码有问题吗？

数据来自以下链接: http://academic.udayton.edu/kissock/http/Weather/gsod95-current/CALOSANG.txt

Answer 1

数据文件包含 14 个值为 -99 的坏样本。我使用最近的好样本用线性插值替换了坏样本。文件中的样本总数为 7006。该数据针对未知的日常天气参数。

正如 Ben 所说，频率单位是 1/天而不是天，因此您不会在 365 个频率单位处看到峰值(假设数据是周期性的，周期为 365 天)。

下图显示了 FFT 的数据输入。已使用线性插值程序删除了不良样本。右端的零值由 FFT 自动添加，以便将输入样本的数量填充到下一个更高的 2 次幂。

下图显示了 FFT 的数据输入，去除了直流偏移 (62.3127)。

下图显示了满量程的 FFT 输出。 FFT 输入包括直流偏移。

下图放大了 FFT 输出的低频端。在图的左侧附近可以看到一个峰值。该峰值对应于您要查找的频率。 FFT 输入包括直流偏移。

下图显示了感兴趣的频率峰值。峰值位于 0.0026855 (22/8192) 个频率单位，对应于 372 天 (1/0.0026855) 的时间值。 FFT 输入包括直流偏移。

下图显示了使用 Blackman-Harris 92 dB 高分辨率窗口对输入数据进行窗口化后的 FFT。感兴趣的频率峰值比周围背景高 18 dB。在这种情况下开窗显着减少了频谱泄漏。窗口化后，峰值保持在 0.0026855 (22/8192) 个频率单位，对应于 372 天 (1/0.0026855) 的时间值。 FFT 输入包括直流偏移。

下图显示了从输入数据中去除直流偏移 (62.3127) 且未对输入数据应用窗口化后的 FFT。

下图显示了从输入数据中去除直流偏移 (62.3127) 并应用 Blackman-Harris 92 dB 高分辨率窗口后的 FFT。

FFT 和图表是用 Sooeet FFT calculator 完成的