从数据的FFT分析中提取含义

Question

我想问的是关于快速傅里叶变换（FFT）的事情，因为这是我第一次使用它。

我有一组按年份（从1700年到2009年）记录的数据，每一年都有一个对应的数值（读数）。当我把这些读数和年份画成图时，得到了下面的第一个图。现在，我的目标是用Python的FFT找出读数最高的主导周期（从图上看，大约在1940到1950年之间）。于是我进行了FFT，得到了它的幅度和功率谱（功率谱的图见第二个图）。功率谱显示主导频率在0.08到0.1之间（循环/年）。我的问题是，我该如何将这个频率与读数和年份联系起来？也就是说，我怎么知道这个主导频率对应的是哪一年（或者哪几年的时间段）？我该如何利用它来找到这个信息呢？

数据列表可以在这里找到：

http://www.physics.utoronto.ca/%7Ephy225h/web-pages/sunspot_yearly.txt

我写的代码是：

from pylab import *
from numpy import *
from scipy import *
from scipy.optimize import leastsq
import numpy.fft

#-------------------------------------------------------------------------------
# Defining the time array
tmin = 0
tmax = 100 * pi
delta = 0.1
t = arange(tmin, tmax, delta)

# Loading the data from the text file
year, N_sunspots = loadtxt('/Users/.../Desktop/sunspot_yearly.txt', unpack = True) # years and number of sunspots

# Ploting the data
figure(1)
plot(year, N_sunspots)
title('Number of Sunspots vs. Year')
xlabel('time(year)')
ylabel('N')

# Computing the FFT
N_w = fft(N_sunspots)

# Obtaining the frequencies 
n = len(N_sunspots)
freq = fftfreq(n) # dt is default to 1

# keeping only positive terms
N = N_w[1:len(N_w)/2.0]/float(len(N_w[1:len(N_w)/2.0]))
w = freq[1:len(freq)/2.0]

figure(2)
plot(w, real(N))
plot(w, imag(N))
title('The data function f(w) vs. frequency')
xlabel('frequency(cycles/year)')
ylabel('f(w)')
grid(True)

# Amplitude spectrum
Amp_spec = abs(N)

figure(3)
plot(w, Amp_spec)
title('Amplitude spectrum')
xlabel('frequency(cycles/year)')
ylabel('Amplitude')
grid(True)

# Power spectrum
Pwr_spec = abs(N)**2

figure(4)
plot(w, Pwr_spec 'o')
title('Power spectrum')
xlabel('frequency(cycles/year)')
ylabel('Power')
grid(True)

show()

Answer 1

下面的图展示了输入到快速傅里叶变换（FFT）中的数据。原始数据文件总共有309个样本。右边的零值是FFT自动添加的，用来把输入样本的数量填充到下一个更高的2的幂次（2^9 = 512）。

下面的图展示了应用了Kaiser-Bessel窗口函数（a=3.5）后的FFT输入数据。窗口函数可以减少FFT中的频谱泄漏错误，尤其是在输入信号在采样区间内不是周期性的时候，就像这个例子一样。

下面的图展示了没有应用窗口函数的FFT输出，显示在全尺度下。峰值出现在0.0917968（47/512）频率单位，这对应的时间值是10.89年（1/0.0917968）。

下面的图展示了应用了Kaiser-Bessel窗口（a=3.5）的FFT输出，显示在全尺度下。峰值仍然位于同样的频率位置0.0917968（47/512）频率单位，对应的时间值也是10.89年（1/0.0917968）。由于窗口函数减少了频谱泄漏，峰值在背景上更加明显。

总之，我们可以很有把握地说，原始帖子中提供的太阳黑子数据是周期性的，基本周期为10.89年。

FFT和图表是通过Sooeet FFT计算器完成的。