使用 FFT 分析谷歌趋势时间序列的季节性

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 图，我们跳过零频率 bin（第一个点），因为它只包含基线，省略它可以更方便地缩放 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 年和 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()

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

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