从数据的 FFT 分析中提取意义答案

【问题标题】：Extracting meaning from an FFT analysis of data从数据的 FFT 分析中提取意义
【发布时间】：2014-04-28 09:21:19
【问题描述】：

我的问题是关于快速傅里叶变换，因为这是我第一次使用它们。因此，我有一组按年份（从 1700 年到 2009 年）和每年对应于某个值（读数）的数据。当我根据年份绘制读数时，它给了我下面的第一个图。现在，我的目标是使用 FFT 和 python 找到读数最高的主导时期（从图中看来，它似乎在 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()

【问题讨论】：

标签： python fft frequency

【解决方案1】：

下图显示了输入到 FFT 的数据。原始数据文件共包含 309 个样本。右端的零值由 FFT 自动添加，以将输入样本数填充到下一个更高的 2 次幂 (2^9 = 512)。

下图显示了输入到 FFT 的数据，其中应用了 Kaiser-Bessel a=3.5 窗函数。当 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 完成的

【讨论】：