Python 数学建模——傅里叶变换时间序列分析

前言

  在数学建模过程中，得到一个序列 $x_1,\cdots ,x_n$，我们首先要进行数据分析，其中就包括分析数据的周期性。这里的周期性不是数学上严格的周期性，即不必 $x_i$ 严格等于 $x_{i+T}$，大致相等就可以。
  比如下面是 $1700-1988$ 年间，每一年太阳黑子数量的记录图。太阳黑子数量随着时间是否有一定的周期性？

  光看图片的话，太阳黑子几乎都是在 $5$ 年内持续增长，随后 $5$ 年内持续下降，周而复始，大概有一个 $10$ 年的周期。但是光肉眼看还是太主观了，有没有系统一点的办法？
  傅里叶变换是将时域数据 $f(x)$ 变换到频域 $F(\mathrm{j}\omega )$ 的变换法，根据数据在频域的图象 $F(\mathrm{j}\omega )$ 的极大值，我们可以轻松看出 $f(x)$ 有哪些主要的频率分量。而这些分量就对应 $f(x)$ 主要的周期。
  本文主要介绍傅里叶变换时间序列分析的 Python 实现，对其原理不过多讲述，可以参考其他文献。

原理

傅里叶变换时间序列分析的用途：可以用来找时间序列的周期，或者主观地观察到了周期后，可以用这个理论验证一下。找到周期后，可以联系问题情境解释原因。如果需要用到 Prophet 的话，还可添加自定义周期性（add_seasonality，详见Python 数学建模——Prophet 时间序列预测_多变量prophet-CSDN博客）。

  傅里叶变换分为连续傅里叶变换（CFT）和离散傅里叶变换（DFT）。由于计算机中只能存储离散信号，计算 CFT 非常不方便，实际都是用的离散傅里叶变换，也称快速傅里叶变换（FFT）：$$y_k=\sum _{n=0}^{N-1}{e}^{-2\pi \mathrm{j}kn/N}{x}_n$$  其中 $y_1,\cdots ,y_n$ 是 $x_1,\cdots ,x_n$ 的离散傅里叶变换。

Python 库函数实现

单周期函数

  调用 Python 库函数scipy.fft.fft, scipy.fft.fftfreq可以很好地进行傅里叶变换：

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

N = 400
T = 1.0 / 1000

# 生成一个周期为 1/114 的正弦函数
x = np.linspace(0, N*T, N, endpoint = False)
y = np.sin(114 * 2 * np.pi * x)

# z 是对序列进行离散傅里叶变换的结果，abs 对复数取模
z = np.abs(fft(y))

# fftfreq 根据我们的采样周期和采样点个数，返回对应的频率值序列
# 其返回值按照 0 → 最大正数 → 最小负数 → 0 排列， [:N//2] 取前半部分（正频率）
fx = fftfreq(N, T)[:N//2]

plt.rcParams['font.family'] = 'Euclid'
plt.plot(fx,2.0 / N * z[:N//2])
plt.grid()
plt.show()

  在上面的代码中，我们生成了一个正弦函数 $y=\sin (114\times 2\pi x)$，这个函数是一个周期为 $T=1/114$ 的周期函数，因此它具有频率 $f=114$。代码中， $z$ 是 $y$ 的离散序列傅里叶变换，最后作图的结果如下所示：

  可以明显地看出，在傅里叶变换之后，频谱图象有一个峰值。理论上，这个峰值对应的频率应该就是 $f=114$。利用代码print(fx[np.argmax(z)])求出这个极值点是 $f_0=115$。

关于为什么 $f_0=115$ 而不是 $114$，这是因为程序求出来的 $z$ 并不是一个函数，而是函数上的一系列函数值。这些函数值对应的横坐标间隔比较大，影响了 $f$ 的精度。

多周期函数

  下面我们尝试 $y=\sin (114\times 2\pi x)+\sin (514\times 2\pi x)$，也就是两个周期函数的叠加，看看能不能同时找出这两个频率分量。

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
import pandas as pd
from scipy.signal import argrelextrema
N = 500
T = 1.0 / 2500

plt.rcParams['font.family'] = 'Euclid'

# 生成一个周期 1/114 的正弦函数和 1/514 的正弦函数的叠加
x = np.linspace(0, N*T, N, endpoint = False)
y = np.sin(114 * 2 * np.pi * x) + np.sin(514 * 2 * np.pi * x)

# z 是对序列进行离散傅里叶变换的结果，abs 对复数取模
z = np.abs(fft(y))
fx = fftfreq(N, T)
lst = list(zip(fx, 2.0 / N * z))
l = len(fx)
plt.plot(fx[:l//2],2.0 / N * z[:l//2])
plt.plot(fx[l//2:],2.0 / N * z[l//2:])
plt.show()
print(fx[argrelextrema(z,np.greater,order=1)])
# [ 115.  515. -515. -115.]

  作图结果如下所示，通过print(fx[argrelextrema(z,np.greater,order=1)])得出来的极大值点也是 $f_{1,2,3,4}=115,515,-515,-115$。负频率的出现是因为我们采用傅里叶变换的指数形式，频率值不完全与 $114,514$ 重合的原因也和上面相同。

真实数据挑战

  上面两个函数都是已知周期，用傅里叶变换去验证，看看 Python 库好不好用，效果真不真。现在给出一个未知周期的函数，让傅里叶变换发挥作用。当然，我们就用之前提到的太阳黑子数据。

点击下载路径，Ctrl+S下载太阳黑子数据sunspots.csv表单。

  使用太阳黑子数据sunspots.csv如上图（左）所示。认为它具有一定周期性，编写以下代码分析太阳黑子的周期性，绘制频谱图如上图（右）所示。根据图象可以得出结论，太阳黑子数量在短期内具有大概 $10$ 天的周期性 $f\approx 0.1$，在长期内具有大概 $100$ 天的周期性 $f\approx 0.01$。

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
import pandas as pd
from scipy.signal import argrelextrema


df = pd.read_csv("sunspots.csv")
y = df['counts']
y -= y.mean() # 减去均值以消去直流分量

# z 是对序列进行离散傅里叶变换的结果，abs 对复数取模
z = np.abs(fft(y.values))[:N//2]
fx = fftfreq(len(y), 1)[:N//2]
plt.plot(fx,2.0 / N * z)

plt.rcParams['font.family'] = 'Euclid'
plt.show()
# 根据图像找出最尖的那几个峰对应频率
print(fx[argrelextrema(z,lambda x,y: 2.0 / N * x > 13,order=1)])
# [0.01038062 0.08304498 0.0899654  0.10034602]

有关库函数argrelextrema的使用，参见Python 基本库用法：数学建模_数模代码库python-CSDN博客。上面最后一行代码是找出频谱图中纵坐标大于 $13$ 的点对应的横坐标（频率）。

参考文献：
Fourier Transforms (scipy.fft) — SciPy v1.14.1 Manual
时间序列分析之：傅里叶变换找周期_傅里叶变换去数据周期-CSDN博客