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

我这里可能没有一个很好的答案，但我想让你的代码更容易让别人复制，这样他们就不需要手动下载和处理文件了：

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

url = 'http://academic.udayton.edu/kissock/http/Weather/gsod95-current/CALOSANG.txt'
response = urllib2.urlopen(url)
html = response.read()
rows = [[float(c) if '.' in c else int(c)
         for c in row.split()] 
         for row in html.splitlines()]

signal = []
for line in rows:
    d = float(line[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()

Answer 2

你的快速傅里叶变换（FFT）应该是根据频率来绘制，而不是时间。

import shlex
import pylab as plt
import numpy as np
from scipy.fftpack import rfft, rfftfreq

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.arange(1,size+1)

f_signal = rfft(signal)
freq = rfftfreq(size, d=1.)


plt.ion()
plt.close("all")
plt.figure()
plt.plot(time,signal)
plt.figure()
plt.plot(freq,abs(f_signal))
plt.show()

你不会期待在365这个时间点看到一个尖峰，因为那是时间值。对于一个周期为365的信号，它对应的频率是1/365，或者说是0.00274。

Answer 3

这个数据文件里有14个坏样本，它们的值是-99。我用最近的好样本通过线性插值的方法替换了这些坏样本。文件里总共有7006个样本。这些数据是关于一个未知的每日天气参数。

正如Ben所说，频率单位是1/天，而不是天，所以你不会在365个频率单位的位置看到一个峰值（假设数据是周期性的，周期为365天）。

下面的图展示了输入到FFT（快速傅里叶变换）中的数据。坏样本已经通过线性插值的方式去掉了。右边的零值是FFT自动添加的，目的是为了把输入样本的数量填充到下一个更高的2的幂次。

下面的图展示了去掉直流偏移（62.3127）后的FFT输入数据。

下面的图展示了完整比例的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计算器完成的。