2.5（発展）スペクトル解析¶
最後に時系列データの解析手法の代表であるスペクトル解析を実践しよう。 用語解説を含めた詳細の説明は、別途こちら(B.フーリエ解析とスペクトル解析）にまとめておいたので、興味があれば参照されたい。 ここでは実例を示すにとどめる。
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import copy
from scipy.fftpack import fft,ifft,fftfreq
from scipy import signal
file = "../AMeDAS/data/Kyoto_1hr/Kyoto_2017_T.csv"
df = pd.read_csv(file,encoding = "shift-jis",header=2,index_col=0,parse_dates=[0])#,skiprows=[1,2,3])
df = df.iloc[1:,:] # ヘッダーの余計な部分を除く
df.head()
Out[1]:
In [2]:
tmp = df['気温(℃)']
#tmp = df['風速(m/s)']
#tmp = df['現地気圧(hPa)']
tmp = tmp.interpolate()
dT = 1 # データ間隔-> 1 (hour)
# 元データからトレンドを差し引いて、偏差時系列を作っておく([関数] signal.detrend)
tmp_a = signal.detrend(tmp)
N = tmp_a.shape[0] # 時間方向の配列の大きさ
pi = np.pi
domg = 2.0*pi/(N*dT) #最小周波数
# FFT
yfft = fft(tmp_a)/N
# パワースペクトル密度に変換(PSD)
pw = abs(yfft)**2
psd = pw/domg
# 片側スペクトルに変換
psd_oneside = psd[0:int(N/2)+1]*2 # 最初の半分をとってきて２倍する
psd_oneside[int(N/2)] /= 2.0 ## 0, N/2の成分(N:偶数のとき）は二倍しなくて良いので元に戻す（N/2の成分はほとんど寄与しないが）
psd_oneside[0] /= 2.0
#### 周波数軸の準備---------
# 周波数
xfft = np.linspace(0.0,domg*(N/2),int(N/2)+1)
# "[関数] fftfreq"を用いて以下のようにもできる
#xfft2 = 2.0*pi*fftfreq(N,d=dT) # 両側スペクトルの周波数（負の値もある）
#xfft = xfft2[0:int(N/2)+1] # 片側だけ（正の値だけ）
# 周期（hr）に変換
xfft_T = 2*pi/xfft
##-------------------------
# 描画
fig = plt.figure(figsize=(10,3))
ax = fig.add_subplot(111)
# PSD
ax.plot(xfft_T,psd_oneside,label='PSD')
##（注意）a[start:end]としたとき、endは含まないことに注意
# omega*PSD (x軸をlog座標で表示する際にはこうすると各周波数帯のエネルギーの寄与が見やすい)
# ax.plot(xfft_T[1:int(N/2)+1],psd_oneside[1:int(N/2)+1]*xfft[1:int(N/2)+1],label='ω*PSD')
#設定
ax.set_yscale('log') # 軸を対数表示に設定
ax.set_xscale('log')
ax.invert_xaxis() # X軸を反転
ax.set_title('Power Spectrum',size='x-large') # タイトルとか
ax.legend(loc=2) # 凡例
plt.show()
時間スケールが大きいほどエネルギーが大きい右肩下がりの分布をしており、赤色ノイズ (red noise)とも言われる。
ただしその中にも、日周期(24 hr, 12 hr,....)に明瞭なピークが見られる。
各周波数帯からの寄与¶
まず、周波数についてのPSDの積分が全分散に等しい(パーセバルの式）ことを確かめておく。
In [3]:
#元の時系列の分散がPSDの全積分に等しいことを確認する（パーセバルの式）
# (1)分散
var = tmp_a.var()
# (2) PSD積分
var2 = (psd_oneside[0:int(N/2)+1]*domg).sum()
print('variance=', var)
print('Integral of PSD=',var2)
従って、特定の周期帯にわたって積分したパワースペクトルを分散と比べることで、その周期帯の寄与率を見積もることができる。
例として、総観規模擾乱(3-10日程度；高・低気圧スケールの変動）の寄与を見積もってみよう。
In [4]:
xfft_Ts = xfft_T[0:int(N/2)+1]
psd_oneside_sub = psd_oneside[0:int(N/2)+1]
mask = (xfft_Ts < 240) & (xfft_Ts > 72)
psd_synoptic = psd_oneside_sub[mask] #synoptic（総観規模）
print(psd_synoptic.shape)
var_synoptic = (psd_synoptic*domg).sum()
print('variance for synoptic-scale dist. = ', var_synoptic)
全分散に対する寄与は、かなり小さい（<3%)。 これは、少なくとも気温については季節変化の影響が非常に大きいためだと考えられる（他の変数ではどうかも試してみるとよい）。