Pythonプログラミング(ステップ7・リスト・データの周期性)
このページでは、系列データの中の周期性の検出方法について考えてみる。
周期関数のテンプレート
以下では、簡単のため、データの平均は0と仮定する。そうでないデータの場合、全体から平均を差し引いておいてから考えればよい。
等しい間隔
データ点は
ここで、
三角関数の総和が増大しないことの考察
ここで考えているのは
ここで、
ところが、テスト波形の角振動数
以上の考察から、時系列データから角振動数と振幅を推定するアルゴリズムとして、以下が考えられる:
を色々と変えながら以下の量を計算する:- 振幅
を計算する: - 振幅
を計算する: が最大となる がデータの角振動数、 を振幅とする。
ここまでの考えは、元々の時系列が単一の角振動数のみで構成されていることを何ら前提にしていないので、
周期と振幅の解析
上記のアルゴリズムを使って、東北電力が公開している東北・新潟地区の電力消費量の時系列データ(データはからダウンロードできる)の周期性を検出してみよう。 このデータはCSV形式で、データ点は1時間おきに与えられている。
時系列データをCSVファイルから読み込み、周期
データ: juyo_2018_tohoku.txt
# coding: utf-8 import pandas as pd import math import matplotlib.pyplot as plt data = pd.read_csv('juyo_2018_tohoku.txt', encoding='Shift_JIS', skiprows=1) date = data['DATE'].str.cat(data['TIME'],sep=" ") datetime = pd.to_datetime(date,format='%Y/%m/%d %H:%M') unixtime = datetime.map(pd.Timestamp.timestamp) x = data['実績(万kW)'].tolist() N = len(x) AMP2 = [ ] PERIOD = [] for period in range(2,200,1): omega = 2*math.pi/period s=0 c=0 for i in range(N): s = s + math.sin(omega*i) * x[i] c = c + math.cos(omega*i) * x[i] AMP2.append( (s*2/N)**2 + (c*2/N)**2) PERIOD.append(period) plt.plot(PERIOD,AMP2, color=(1.0,0.0,0.0), linewidth=1.0) plt.xlabel('PERIOD [HOURS]') plt.ylabel('AMPLITUDE^2') plt.grid(True) plt.show()
以下が解析結果のプロットで、横軸は周期(時間)、縦軸は振幅の2乗である。
明らかに24時間と1週間(
実行結果
フーリエ変換とFFT
ここで注意したいのは、(11)式のように、離散フーリエ変換では着目する角振動数を
データ点が多くなると、全ての
うまいことに、中間的な計算結果を巧妙に使い回すことによって、これを
以下のコードは、NumPyのFFT関数を使って、時系列配列 x
のフーリエ変換を求め、
それから、データの絶対値の二乗をプロットする例である。
# coding: utf-8 import pandas as pd import math import numpy as np import matplotlib.pyplot as plt data = pd.read_csv('juyo_2018_tohoku.txt', encoding='Shift_JIS', skiprows=1) date = data['DATE'].str.cat(data['TIME'],sep=" ") datetime = pd.to_datetime(date,format='%Y/%m/%d %H:%M') unixtime = datetime.map(pd.Timestamp.timestamp) x = data['実績(万kW)'].tolist() N = len(x) freq = np.linspace(0,1.0,N) y = np.fft.fft(x, norm=None) amp2 = np.absolute(y*2/N)**2 plt.plot(freq[10:N//2], amp2[10:N//2], color=(1.0,0.0,0.0), linewidth=1.0) plt.xlabel('FREQ [1/HOUR]') plt.ylabel('AMPLITUDE^2') plt.grid(True) plt.show()
上記のコードを実行すると、横軸に周波数(1/時間)、縦軸のその周波数成分の振幅がプロットされる。
もっとも大きなピークは24時間、すなわち
実行結果
解説: 離散フーリエ変換についての注意点
データの前処理
分析の上で興味の無い、周波数が0およびそれに近い成分の振幅が相対的に大きいため、上記のコードでは
サンプリング間隔
連続的な時間現象を、一定間隔
高い周波数成分を持つ現象を、それには不十分なサンプリング速度で計測すると、折返し雑音と呼ばれる、本来は存在しない見かけ上の周期性が検知されてしまうので、注意が必要である。
周期
窓関数
離散フーリエ変換では、
そこで、窓関数と呼ばれる、両端でなめらかに零になるような関数とデータをかけて、両端が連続的に繋がるような前処理を行なっておく。 窓関数としては、細かな形状の違いにより、Hann窓、Hamming窓、Blackman窓などが用いられる。
以上を考慮して、若干手直ししたコードの例を示す:
# coding: utf-8 import pandas as pd import math import numpy as np import matplotlib.pyplot as plt data = pd.read_csv('juyo_2018_tohoku.txt', encoding='Shift_JIS', skiprows=1) date = data['DATE'].str.cat(data['TIME'],sep=" ") datetime = pd.to_datetime(date,format='%Y/%m/%d %H:%M') unixtime = datetime.map(pd.Timestamp.timestamp) x = data['実績(万kW)'].tolist() N = len(x) x0 = x - np.mean(x) # 直流成分除去 window = np.hanning(N) # Hann窓関数定義 ACF=1.63 # 振幅調整 x1 = (x0 * window) * ACF # 窓関数を適用 time = np.linspace(0,N,N) freq = np.linspace(0,1.0,N) y = np.fft.fft(x1, norm=None) amp2 = np.absolute(y*2/N)**2 fig = plt.figure() ax1 = fig.add_subplot(3, 1, 1) ax2 = fig.add_subplot(3, 1, 2) ax3 = fig.add_subplot(3, 1, 3) ax1.plot(time, x0, color=(0.0,0.0,1.0), linewidth=0.2) ax1.set_ylabel('10^4 KWH') ax2.plot(time, x1, color=(0.0,0.0,1.0), linewidth=0.2) ax2.set_xlabel('TIME [HOUR]') ax2.set_ylabel('10^4 KWH') ax3.plot(freq[0:N//2], amp2[0:N//2], color=(1.0,0.0,0.0), linewidth=1.0) ax3.set_xlabel('FREQ [1/HOUR]') ax3.set_ylabel('AMPLITUDE^2') ax3.grid(True) fig.tight_layout() plt.show()
実行結果
上から、解析対象の時系列データ(電力消費量の推移)、窓関数を適用後の時系列、周波数と振幅の二乗の関係。