用傅里葉變換解碼時間序列：從頻域視角解析季節(jié)性模式
在眾多時間序列模型中，SARIMA（seasonal autoregressive integrated moving average，季節(jié)性自回歸積分滑動平均模型）能夠有效處理時間序列中的季節(jié)性成分。但是在實(shí)際應(yīng)用中，如何準(zhǔn)確識別和提取這些季節(jié)性模式一直是一個挑戰(zhàn)。
傳統(tǒng)上，識別季節(jié)性模式往往依賴于數(shù)據(jù)的可視化分析。但是我們可以使用傅里葉變換以及周期圖（Periodogram）這一強(qiáng)大工具，用一種更系統(tǒng)的方法來解決這個問題。
本文將詳細(xì)介紹一些基礎(chǔ)但重要的概念，這些方法對于每位研究時間序列的數(shù)據(jù)科學(xué)家都具有實(shí)用價值。
- 傅里葉變換的基本原理
- Python實(shí)現(xiàn)傅里葉變換
- 周期圖分析方法
問題概述
我們以AEP（American Electric Power）能源消耗數(shù)據(jù)集為例（數(shù)據(jù)集使用CC0許可）：
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv("data/AEP_hourly.csv", index_col=0)
df.index = pd.to_datetime(df.index)
df.sort_index(inplace=True)
fig, ax = plt.subplots(figsize=(20,4))
df.plot(ax=ax)
plt.tight_layout()
plt.show()
從數(shù)據(jù)的初步可視化中可以明顯觀察到季節(jié)性模式的存在，但要精確捕獲所有這些模式并非易事。
傳統(tǒng)的手動分析方法通常需要多個時間尺度的可視化，如下所示：
fig, ax = plt.subplots(3, 1, figsize=(20,9))
df_3y = df[(df.index >= '2006–01–01') & (df.index < '2010–01–01')]
df_3M = df[(df.index >= '2006–01–01') & (df.index < '2006–04–01')]
df_7d = df[(df.index >= '2006–01–01') & (df.index < '2006–01–08')]
ax[0].set_title('AEP energy consumption 3Y')
df_3y[['AEP_MW']].groupby(pd.Grouper(freq = 'D')).sum().plot(ax=ax[0])
for date in df_3y[[True if x % (24 * 365.25 / 2) == 0 else False for x in range(len(df_3y))]].index.tolist():
ax[0].axvline(date, color = 'r', alpha = 0.5)
ax[1].set_title('AEP energy consumption 3M')
df_3M[['AEP_MW']].plot(ax=ax[1])
for date in df_3M[[True if x % (24 * 7) == 0 else False for x in range(len(df_3M))]].index.tolist():
ax[1].axvline(date, color = 'r', alpha = 0.5)
ax[2].set_title('AEP energy consumption 7D')
df_7d[['AEP_MW']].plot(ax=ax[2])
for date in df_7d[[True if x % 24 == 0 else False for x in range(len(df_7d))]].index.tolist():
ax[2].axvline(date, color = 'r', alpha = 0.5)
plt.tight_layout()
plt.show()
不同時間尺度下的AEP能源消耗模式
通過多尺度分析，我們可以觀察到以下主要周期：
- 半年度周期（約180天）
- 周度周期（7天）
- 日度周期（24小時）
雖然對于能源消耗數(shù)據(jù)而言，這些季節(jié)性模式可以通過領(lǐng)域知識推斷，但僅依靠人工檢查存在以下局限性：
- 主觀性：容易忽略不明顯的模式
- 低效率：需要逐一檢查不同時間窗口
- 擴(kuò)展性差：難以應(yīng)用于大規(guī)模數(shù)據(jù)分析
作為數(shù)據(jù)科學(xué)家，我們需要一個能夠快速、準(zhǔn)確識別時間序列中關(guān)鍵頻率的工具。這正是傅里葉變換發(fā)揮作用的地方。
1、傅里葉變換的基本原理
傅里葉變換是一種將信號從時域轉(zhuǎn)換到頻域的數(shù)學(xué)工具。在時域中，我們觀察數(shù)據(jù)隨時間的變化；而在頻域中，我們可以看到構(gòu)成信號的各個頻率及其相對重要性。
根據(jù)傅里葉理論，任何滿足一定條件的函數(shù)f(x)都可以表示為不同頻率、振幅和相位的正弦函數(shù)的疊加。換言之，每個時間序列都可以分解為一系列基本波形的組合。
其中：
- F(f)表示頻域中的函數(shù)
- f(x)表示時域中的原始函數(shù)
- exp(-i2πf(x))是復(fù)指數(shù)函數(shù)，用作頻率分析器
函數(shù)F(f)的值表明頻率f在原始信號中的貢獻(xiàn)程度。
考慮一個由三個不同頻率（2Hz、3Hz和5Hz）的正弦波組成的復(fù)合信號：
應(yīng)用傅里葉變換后，我們可以提取這些基本頻率：
頻域表示清晰地顯示出信號中存在的三個基本頻率（2Hz、3Hz和5Hz）。將信號分解為基本波形：
原始信號（藍(lán)色）是三個基本波形（紅色）的疊加。這種分解方法可以應(yīng)用于任何時間序列，用于識別其主要頻率成分。
2、Python中的傅里葉變換實(shí)現(xiàn)
現(xiàn)在我們將這個理論應(yīng)用到之前的AEP能源消耗數(shù)據(jù)中。
Python的numpy.fft模塊提供了計(jì)算離散信號傅里葉變換的工具。FFT（快速傅里葉變換）是一種高效的算法，用于將離散信號分解為頻率成分：
from numpy import fft
X = fft.fft(df['AEP_MW'])
N = len(X)
frequencies = fft.fftfreq(N, 1)
periods = 1 / frequencies
fft_magnitude = np.abs(X) / N
mask = frequencies >= 0
# 繪制傅里葉變換結(jié)果
fig, ax = plt.subplots(figsize=(20, 3))
ax.step(periods[mask], fft_magnitude[mask]) # 僅繪制正頻率
ax.set_xscale('log')
ax.xaxis.set_major_formatter('{x:,.0f}')
ax.set_title('AEP energy consumption - Frequency-Domain')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Magnitude')
plt.show()
頻域可視化顯示了某些頻率具有較高的幅值，表明這些頻率在原始信號中具有重要作用。
3、周期圖分析
周期圖是信號功率譜密度（Power Spectral Density, PSD）的頻域表示。相比簡單的傅里葉變換，周期圖能更好地量化各頻率分量的強(qiáng)度，并且可以有效降低次要頻率的噪聲影響。
周期圖的數(shù)學(xué)定義為：
其中：
- P(f)是頻率f處的功率譜密度
- X(f)是信號的傅里葉變換
- N是樣本數(shù)量
Python實(shí)現(xiàn)如下：
power_spectrum = np.abs(X)**2 / N # 計(jì)算每個頻率的功率
fig, ax = plt.subplots(figsize=(20, 3))
ax.step(periods[mask], power_spectrum[mask])
ax.set_title('AEP energy consumption Periodogram')
ax.set_xscale('log')
ax.xaxis.set_major_formatter('{x:,.0f}')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.show()
從周期圖分析中，我們可以清晰地識別出以下主要頻率：
- 24Hz：對應(yīng)24小時周期
- 4,380Hz：對應(yīng)6個月周期
- 168Hz：對應(yīng)每星期
除此之外，還發(fā)現(xiàn)了三個次要周期：
- 12Hz周期
- 84Hz周期（對應(yīng)半周）
- 8,760Hz周期（對應(yīng)年度周期）
我們也可以使用scipy.signal模塊中的periodogram函數(shù)獲得相同的結(jié)果：
from scipy.signal import periodogram
frequencies, power_spectrum = periodogram(df['AEP_MW'], return_onesided=False)
periods = 1 / frequencies
fig, ax = plt.subplots(figsize=(20, 3))
ax.step(periods, power_spectrum)
ax.set_title('Periodogram')
ax.set_xscale('log')
ax.xaxis.set_major_formatter('{x:,.0f}')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.show()總結(jié)
在時間序列分析中，準(zhǔn)確識別季節(jié)性模式是一個關(guān)鍵任務(wù)。本文介紹的周期圖分析方法為我們提供了一個系統(tǒng)、高效的工具，可以準(zhǔn)確識別時間序列中的各種周期性成分。這種方法不僅克服了傳統(tǒng)視覺分析的局限性，還能發(fā)現(xiàn)可能被忽略的次要周期模式，為時間序列分析提供了更全面的視角。
