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

Question

我的问题是关于快速傅立叶变换，因为这是我第一次使用它们。 因此，我有一组按年份(从 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 中的频谱泄漏误差。

下图显示了满量程的 FFT 输出。没有窗口功能。峰值位于 0.0917968 (47/512) 个频率单位，对应于 10.89 年 (1/0.0917968) 的时间值。

下图显示了满量程的 FFT 输出。应用 Kaiser-Bessel a=3.5 窗口。峰值保持在 0.0917968 (47/512) 个频率单位处的相同频率位置，对应于 10.89 年 (1/0.0917968) 的时间值。由于窗函数减少了光谱泄漏，峰值在背景上方更加清晰可见。

总之，我们可以高度肯定地说，原帖中提供的太阳黑子数据是周期性的，基本周期为 10.89 年。

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