python - 使用 FFT 分析 Google 趋势时间序列的季节性

Question

我正在尝试使用快速傅里叶变换评估 Google 趋势时间序列的振幅谱。如果您在提供的数据中查看“饮食”数据 here它显示出非常强烈的季节性模式:

我想我可以使用 FFT 分析这种模式，它大概应该有一个为期 1 年的强峰值。

但是，当我像这样应用 FFT 时(a_gtrend_ham 是时间序列乘以汉明窗口):

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft, fftshift
import pandas as pd

gtrend = pd.read_csv('multiTimeline.csv',index_col=0)

gtrend.index = pd.to_datetime(gtrend.index, format='%Y-%m')

# Sampling rate
fs = 12 #Points per year

a_gtrend_orig = gtrend['diet: (Worldwide)']
N_gtrend_orig = len(a_gtrend_orig)
length_gtrend_orig = N_gtrend_orig / fs
t_gtrend_orig = np.linspace(0, length_gtrend_orig, num = N_gtrend_orig, endpoint = False)

a_gtrend_sel = a_gtrend_orig.loc['2005-01-01 00:00:00':'2017-12-01 00:00:00']
N_gtrend = len(a_gtrend_sel)
length_gtrend = N_gtrend / fs
t_gtrend = np.linspace(0, length_gtrend, num = N_gtrend, endpoint = False)

a_gtrend_zero_mean = a_gtrend_sel - np.mean(a_gtrend_sel)
ham = np.hamming(len(a_gtrend_zero_mean))
a_gtrend_ham = a_gtrend_zero_mean * ham

N_gtrend = len(a_gtrend_ham)
ampl_gtrend = 1/N_gtrend * abs(fft(a_gtrend_ham))
mag_gtrend = fftshift(ampl_gtrend)
freq_gtrend = np.linspace(-0.5, 0.5, len(ampl_gtrend))
response_gtrend = 20 * np.log10(mag_gtrend)
response_gtrend = np.clip(response_gtrend, -100, 100)

我得到的振幅谱没有显示任何主峰:

我对如何使用FFT获取数据序列的频谱的误解在哪里？

Answer 1

这是我认为您要完成的目标的干净实现。我包括图形输出和对其可能含义的简短讨论。

首先，我们使用 rfft() 因为数据是真实值。这节省了时间和精力(并降低了错误率)，否则会产生冗余的负频率。并且我们使用 rfftfreq() 来生成频率列表(同样，无需手动编码，并且使用 api 降低了错误率)。

对于您的数据，Tukey 窗口比 Hamming 和类似的基于 cos 或 sin 的窗口函数更合适。另请注意，我们在乘以窗函数之前减去中位数。 median() 是对基线的相当可靠的估计，肯定比 mean() 更可靠。

在图表中，您可以看到数据从其初始值快速下降，然后以低位结束。 Hamming 和类似窗口为此对中间采样太窄，不必要地衰减了大量有用的数据。

对于 FT 图，我们跳过零频率区间(第一个点)，因为它只包含基线，省略它可以更方便地缩放 y 轴。

您会注意到 FT 输出图中的一些高频分量。 我在下面包含了一个示例代码，它说明了这些高频分量的可能来源。

好的，这是代码:

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import rfft, rfftfreq
from scipy.signal import tukey

from numpy.fft import fft, fftshift
import pandas as pd

gtrend = pd.read_csv('multiTimeline.csv',index_col=0,skiprows=2)
#print(gtrend)

gtrend.index = pd.to_datetime(gtrend.index, format='%Y-%m')
#print(gtrend.index)

a_gtrend_orig = gtrend['diet: (Worldwide)']
t_gtrend_orig = np.linspace( 0, len(a_gtrend_orig)/12, len(a_gtrend_orig), endpoint=False )

a_gtrend_windowed = (a_gtrend_orig-np.median( a_gtrend_orig ))*tukey( len(a_gtrend_orig) )

plt.subplot( 2, 1, 1 )
plt.plot( t_gtrend_orig, a_gtrend_orig, label='raw data'  )
plt.plot( t_gtrend_orig, a_gtrend_windowed, label='windowed data' )
plt.xlabel( 'years' )
plt.legend()

a_gtrend_psd = abs(rfft( a_gtrend_orig ))
a_gtrend_psdtukey = abs(rfft( a_gtrend_windowed ) )

# Notice that we assert the delta-time here,
# It would be better to get it from the data.
a_gtrend_freqs = rfftfreq( len(a_gtrend_orig), d = 1./12. )

# For the PSD graph, we skip the first two points, this brings us more into a useful scale
# those points represent the baseline (or mean), and are usually not relevant to the analysis
plt.subplot( 2, 1, 2 )
plt.plot( a_gtrend_freqs[1:], a_gtrend_psd[1:], label='psd raw data' )
plt.plot( a_gtrend_freqs[1:], a_gtrend_psdtukey[1:], label='windowed psd' )
plt.xlabel( 'frequency ($yr^{-1}$)' )
plt.legend()

plt.tight_layout()
plt.show()

这是以图形方式显示的输出。在 1/year 和 0.14(恰好是 1/14 年的 1/2)处有很强的信号，并且有一组更高频率的信号，乍一看似乎很神秘。

我们看到窗口函数实际上在将数据带入基线方面非常有效，并且您看到应用窗口函数后 FT 中的相对信号强度没有太大改变。

如果你仔细观察数据，一年之内似乎有一些重复的变化。如果这些以某种规律性发生，则可以预期它们会在 FT 中显示为信号，而且实际上 FT 中是否存在信号通常用于区分信号和噪声。但正如将要展示的那样，对高频信号有更好的解释。

好的，现在这里有一个示例代码，说明了可以生成这些高频组件的一种方法。在这段代码中，我们创建了一个音调，然后我们创建了一组与音调频率相同的尖峰。然后我们对两个信号进行傅里叶变换，最后绘制原始数据和 FT 数据。

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import rfft, rfftfreq

t = np.linspace( 0, 1, 1000. )

y = np.cos( 50*3.14*t )

y2 = [ 1. if 1.-v < 0.01 else 0. for v in y ]

plt.subplot( 2, 1, 1 )
plt.plot( t, y, label='tone' )
plt.plot( t, y2, label='spikes' )
plt.xlabel('time')

plt.subplot( 2, 1, 2 )
plt.plot( rfftfreq(len(y),d=1/100.), abs( rfft(y) ), label='tone' )
plt.plot( rfftfreq(len(y2),d=1/100.), abs( rfft(y2) ), label='spikes' )
plt.xlabel('frequency')
plt.legend()

plt.tight_layout()
plt.show()

好的，这里是音调和尖峰的图表，然后是它们的傅立叶变换。请注意，尖峰产生的高频分量与我们数据中的高频分量非常相似。

换句话说，高频分量的起源很可能是在与原始数据中信号的尖峰特征相关的短时间尺度内。