2021-01-19に投稿

Numpyを使用してFFT&トレンド除去

はじめに

Pythonで,Numpyを使って時系列データをFFT(Fast Fourier Transform: 高速フーリエ変換)する方法と,時系列データのトレンドを除去する方法について紹介しようと思います.FFTとは,DFT(Discrete Fourier Transform:離散フーリエ変換)を高速処理する計算方法です.この記事では理論には触れず,FFTを実行する最低限のコードを示します.参考にした文献は「スペクトル解析 著:日野幹雄 (朝倉書店)」.フーリエ解析の基礎からFFTの理論まで,この本1冊で十分です.

目次

1.時系列データ
2.FFT実行
3.トレンド除去
4.フーリエ成分を周波数平滑化(スムージング)
5.Appendix

本題

1. 時系列データ

サンプリング周波数10Hzの30分間データ.オレンジの線は移動平均です.トレンドがあるのが分かりますね.
元時系列.png
図1.元の時系列データ

このデータをFFTしていきます.

2. FFT実行

N =len(X)      #データ長
fs=10          #サンプリング周波数
dt =1/fs       # サンプリング間隔
t = np.arange(0.0, N*dt, dt) #時間軸
freq = np.linspace(0, fs,N) #周波数軸
fn=1/dt/2     #ナイキスト周波数

FFTはデータ長が2のべき乗である場合に最も計算速度が速くなるような計算手法ですが,そうでない場合でも計算は可能です(処理時間が多少長くなりますが).ただし,データ長が素数の場合には,2のべき乗の時に比べて相対的にかなり処理時間がかかるので,2のべき乗になるよう0パディングした方が良さそうです.Appendixにそれを確かめるコードを載せたので是非確かめてみてください.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0 #ナイキスト周波数以降をカット

plt.plot(freq,np.abs(F))#
plt.xlabel("[Hz]")
plt.ylabel("Amp")
plt.xlim(-0.01,0.5)
plt.grid()
plt.show()

np.fft.fftでは複素フーリエ成分で与えられるので,その絶対値をとったものをプロットしたのが下図です.トレンド成分が0Hz付近にありますね.
次節ではトレンド除去してみましょう.
元時系列をFFT.png
図2.元の時系列データに対してFFT実行

3. トレンド除去

図1の移動平均(オレンジ線)を見ると,時系列データの軸がx軸と乖離しており,データにトレンドがあるのが分かります.次の操作で,0.03Hz以下を0にすることでトレンドを除去してみましょう.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0
F[(freq<=0.03)]=0 #0.03HZ以下を除去
X_1=np.real(np.fft.ifft(F))*N

plt.xlabel("Time [s]")
plt.ylabel("Signal")
plt.xlim(-50,1850)
plt.grid()
plt.show()

周波数カット時系列.png
図3.トレンド除去後の時系列データ

x軸を起点に周期関数となっていることが分かります.

4. フーリエ成分を周波数平滑化(スムージング)

図3のデータをFFTすると図4になります.
周波数カット時系列をFFT.png
図4.トレンド除去後の時系列データに対してFFT

ガタガタしてるので,平滑化ウィンドーを作用させてスムージングをしてみます.

window=np.ones(5)/5 #平滑化ウィンドー
F3=np.convolve(F2,window,mode='same') #畳み込み
F3=np.convolve(F3,window,mode='same') #畳み込み
F3=np.convolve(F3,window,mode='same') #畳み込み

plt.plot(freq,np.abs(F3))
plt.xlabel("[Hz]")
plt.ylabel("Amp")
plt.xlim(-0.01,0.5)
plt.grid()
plt.show()

周波数平滑化.png
図5.平滑化した

5. Appendix

データ長さを3種類用意して計算時間を比較してみます.
①2^19(2のべき乗) ②2^(19)-1 (素数) ③2^(19)-2 (素数でも2のいべき乗でもない)

import time

if __name__ == '__main__':
    start = time.time()
    x2 = np.random.uniform(size=2**19-2)#2**19 , 2**19-1
    print(np.fft.fft(x2))
    elapsed_time = time.time() - start
    print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

計算結果
①0.04197[sec]
②0.1679[sec]
③0.05796[sec]

データ長が素数の場合には0パディングを行って2のべき乗にした方が良さそうですね.

ツイッターでシェア
みんなに共有、忘れないようにメモ

kawai_mizugorou

修士課程2年.自分用のメモとして使ってます.

Crieitは誰でも投稿できるサービスです。 是非記事の投稿をお願いします。どんな軽い内容でも投稿できます。

また、「こんな記事が読みたいけど見つからない!」という方は是非記事投稿リクエストボードへ!

有料記事を販売できるようになりました!

こじんまりと作業ログやメモ、進捗を書き残しておきたい方はボード機能をご利用ください。
ボードとは?

コメント