【気象データで時系列分析③】フーリエ変換

概要

シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、フーリエ変換について扱います。

フーリエ変換の詳しい数式などの話は他記事に譲り，本記事では使い方メインの解説をしています．

フーリエ変換

時系列解析においてフーリエ変換（Fourier Transform） は，データの周期性を分析したい場合や，ある特定の周期成分を除去したい場合などに用いられる手法です．

簡単に言うと（離散）フーリエ変換とは，時系列データ$\{a_m\}_{m=0}^{n}$をいろいろな周期成分をもつ$\sin$や$\cos$の線形和で表すことです．

例えば下図の赤線のようなデータを考えます．

ぐにゃぐにゃ曲がった関数ですが，実はこれは点線に示したsin波とcos波の重ね合わせで作られた関数です．具体的には，

f(x) = 0.5\cdot\sin4x + \cos x

です．つまり，赤線の関数というのは周期π/2と周期2πの成分が1：2の比率で混ざった関数ともいえるわけです．

今回の例では，人工的に作ったデータなので「周期π/2と周期2πの成分が1：2の比率で混ざった関数」と分かるわけですが，実際のデータではどの周期成分がどれだけ入っているかなんてすぐにはわかりません．

しかし，フーリエ変換を使えばこれが可能です．

Pythonでは，例えばnumpy.fft.rfftを使えます¹．コード例を示します．

import numpy as np

# データ
x = np.linspace(0, 2*np.pi, 500)
y1 = 0.5*np.sin(4*x)
y2 = np.cos(x)
y = y1 + y2  

# サンプリングサイズ
N = len(x)
# 実フーリエ変換
f = np.fft.rfft(y) / (N/2)
# 振幅
f = np.abs(f)
# 周波数軸のリスト
freq = np.fft.rfftfreq(len(x), d=1/len(x))

np.fft.rfft(y)はフーリエ変換の結果を返しますが，numpy内部では

A_k = \sum_{m=0}^{n-1} a_m \exp\left\{-2\pi i{mk \over n}\right\} \qquad k = 0,\ldots,n-1.

で波数成分$A_k$を求めているので，cos，sinの係数そのものを得たい場合は，これに$\frac{2}{n}$をかけてやる必要があります（cos成分の係数は実部，sin成分の係数は虚部にマイナスをかけたものからわかる）．

結果を可視化すると次のようになります．

可視化コード

# グラフ表示
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.plot(freq*np.pi, f.real, marker='o')
plt.xticks(range(5), labels=[f'{i}π' for i in range(5)])
plt.xlim(0, 5)
plt.title('cos')
plt.xlabel('frequency')
plt.grid()

plt.subplot(2, 1, 2)
plt.plot(freq*np.pi, -f.imag, marker='o')
plt.xticks(range(5), labels=[f'{i}π' for i in range(5)])
plt.xlim(0, 5)
plt.title('sin')
plt.xlabel('frequency')
plt.grid()
plt.show()

cosの周期$\frac{\pi}{2}$成分の係数が1，sinの周期$2\pi$成分の係数が0.5という結果になり，きちんと合致しています．

気圧データに対してフーリエ変換してみる

次は，気象庁の「過去の気象データダウンロード」から取得した現地気圧データに対してフーリエ変換を行い，どういう周期で変動しているのか調べてみましょう．

地点：京都
期間：2022年11月26日～2023年11月26日
データの種類：時別値

基本的には，先ほどと同じように実フーリエ変換をnp.fft.rfftで実行するだけです．

気を付けるポイントとしては，サンプリング周波数の指定です（グラフ表示に影響）．
今回使うデータは1時間に1個のデータなので，サンプリング周波数は 24[/日] となります．

また，どのような周期成分が含まれているかを調べたいのであって，sin・cosの区別には興味がないのでFの絶対値を見ればよいです．ただし，そのままだと値が極端で見づらいので常用対数を取っています．

# サンプリング数
N = len(df['pres'])

# 実フーリエ変換
F = np.fft.rfft(df['pres']) * (2/N)

# 周波数軸の値を計算
freq = np.fft.rfftfreq(N, d=1/24)  # 1day当たり24サンプル

# 周波数スペクトルの複素数を絶対値に変換
F_abs = np.abs(F)
# 見やすいように常用対数に変換
F_log = np.log10(F_abs)

# グラフ表示
plt.figure(figsize=(10, 5))
plt.plot(freq, F_log)
plt.xticks(np.arange(0, 12))
plt.xlabel('frequency [/day]')
plt.ylabel('log power spectrum')
plt.grid()
plt.show()

見ると1，2[/day]くらいの周期成分が特に目立ちます．このことから，気圧は半日周期で変動しているといえそうですね．

ソースコード

Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．

本記事のコードは01-JMA_data_analyticsの中に入っています．

1. これは実フーリエ変換であり，データが実数のときに使えます．データが複素数の場合はnumpy.fft.fftを使いましょう． ↩