用傅里葉變換解碼時間序列：從頻域視角解析季節性模式
在眾多時間序列模型中，SARIMA（seasonal autoregressive integrated moving average，季節性自回歸積分滑動平均模型）能夠有效處理時間序列中的季節性成分。但是在實際應用中，如何準確識別和提取這些季節性模式一直是一個挑戰。
傳統上，識別季節性模式往往依賴于數據的可視化分析。但是我們可以使用傅里葉變換以及周期圖（Periodogram）這一強大工具，用一種更系統的方法來解決這個問題。
本文將詳細介紹一些基礎但重要的概念，這些方法對于每位研究時間序列的數據科學家都具有實用價值。
- 傅里葉變換的基本原理
- Python實現傅里葉變換
- 周期圖分析方法
問題概述
我們以AEP（American Electric Power）能源消耗數據集為例（數據集使用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()
從數據的初步可視化中可以明顯觀察到季節性模式的存在，但要精確捕獲所有這些模式并非易事。
傳統的手動分析方法通常需要多個時間尺度的可視化，如下所示：
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小時）
雖然對于能源消耗數據而言，這些季節性模式可以通過領域知識推斷，但僅依靠人工檢查存在以下局限性：
- 主觀性：容易忽略不明顯的模式
- 低效率：需要逐一檢查不同時間窗口
- 擴展性差：難以應用于大規模數據分析
作為數據科學家，我們需要一個能夠快速、準確識別時間序列中關鍵頻率的工具。這正是傅里葉變換發揮作用的地方。
1、傅里葉變換的基本原理
傅里葉變換是一種將信號從時域轉換到頻域的數學工具。在時域中，我們觀察數據隨時間的變化；而在頻域中，我們可以看到構成信號的各個頻率及其相對重要性。
根據傅里葉理論，任何滿足一定條件的函數f(x)都可以表示為不同頻率、振幅和相位的正弦函數的疊加。換言之，每個時間序列都可以分解為一系列基本波形的組合。
其中：
- F(f)表示頻域中的函數
- f(x)表示時域中的原始函數
- exp(-i2πf(x))是復指數函數，用作頻率分析器
函數F(f)的值表明頻率f在原始信號中的貢獻程度。
考慮一個由三個不同頻率（2Hz、3Hz和5Hz）的正弦波組成的復合信號：
應用傅里葉變換后，我們可以提取這些基本頻率：
頻域表示清晰地顯示出信號中存在的三個基本頻率（2Hz、3Hz和5Hz）。將信號分解為基本波形：
原始信號（藍色）是三個基本波形（紅色）的疊加。這種分解方法可以應用于任何時間序列，用于識別其主要頻率成分。
2、Python中的傅里葉變換實現
現在我們將這個理論應用到之前的AEP能源消耗數據中。
Python的numpy.fft模塊提供了計算離散信號傅里葉變換的工具。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
# 繪制傅里葉變換結果
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）的頻域表示。相比簡單的傅里葉變換，周期圖能更好地量化各頻率分量的強度，并且可以有效降低次要頻率的噪聲影響。
周期圖的數學定義為：
其中：
- P(f)是頻率f處的功率譜密度
- X(f)是信號的傅里葉變換
- N是樣本數量
Python實現如下：
power_spectrum = np.abs(X)**2 / N # 計算每個頻率的功率
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：對應24小時周期
- 4,380Hz：對應6個月周期
- 168Hz：對應每星期
除此之外，還發現了三個次要周期：
- 12Hz周期
- 84Hz周期（對應半周）
- 8,760Hz周期（對應年度周期）
我們也可以使用scipy.signal模塊中的periodogram函數獲得相同的結果：
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()總結
在時間序列分析中，準確識別季節性模式是一個關鍵任務。本文介紹的周期圖分析方法為我們提供了一個系統、高效的工具，可以準確識別時間序列中的各種周期性成分。這種方法不僅克服了傳統視覺分析的局限性，還能發現可能被忽略的次要周期模式，為時間序列分析提供了更全面的視角。
