自分用作業メモ
概要
以前からFFTを時系列データ(たとえば、地震の発生回数)の『予測』に使えないかと勉強してたのですけれども、無理だとわかったので、次にFFTを使いたくなった時用に経緯を一通りメモ。
でも、補間アルゴリズムとして使う分には良さげです。
FFTとは?
まずは(1次元)FFTのざっくりした説明。
音声に代表されるような時系列データを、周波数ごとに分解するのがフーリエ変換。
その一種、高速フーリエ変換=FFT。
逆FFTで元のデータに戻すことができます。
元データ → FFT → 周波数によるフィルタリング (ローパスフィルタ、ハイパスフィルタ、などなど) → 逆FFT → 時系列データへ戻す
みたいな流れが一般的です。
なぜFFTはデータ予測に使えないか
大抵のフーリエ変換の入門書に書かれていますが、、、フーリエ変換には
◎ 入力データは、ちょうど1周期分である
という前提をおいたアルゴリズムです。
なので、
音声データ (時系列データ) X[t]を入力に使うとして
・X[t] は時刻 tのときの音声の大きさ
・X[0] ~ X[n-1] を入力として与え、フーリエ変換する
・フーリエ変換=サインカーブ、コサインカーブの和(線形結合)として表現
・これをもとに、入力データの続き すなわち、
X[n] X[n+1] X[n+2] ……
を求めてみると、、、
X[n] == X[0]
X[n+1] == X[1]
X[n+2] == X[2]
……
となってしまう。
(入力データの最初に戻ってしまう)
すなわち、X[0]~X[n-1]の補間には使えるが、残念ながら、 X[n]以降を予測するには使えない。
上のグラフの赤い線が入力データ。
x[t]=sin(2*pi*t*5.47)+sin(2*pi*t*3.48+sin(2*pi*t*7.7)
入力の0~1.0の部分をフーリエ変換に用いている。
下はフーリエ変換後=周波数成分に分解されたもの。
入力と同じ0~1.0 の部分では青い丸が線の上に乗っているが、 1.0~2.0では、大きく位置がズレている。
(0~1.0 と同じ位置に青丸が位置している)
FFT (フーリエ変換)→逆フーリエ変換は、補間にはよいが、予測にはまったく使えないことがわかる。
cos, sinなど 正弦波関数や、expを基底関数において、線形回帰すれば予測にも使えるのかもしれないのだけれど、それは今後の勉強の課題、ということで。
今回のコード
# -*- coding: utf8 -*- import matplotlib.pyplot as plt import numpy as np import numpy.fft as fft #単一正弦波を位相をずらしつつFFTしてみる plt.ion() #interactive plt.clf() plt.cla() n=256 n2=n*2 y=np.zeros(n2) flist=np.random.random(3)*8 for f in flist: #for f in [2.3,1.8,1.04]: x=np.ogrid[0:2*float(n2)/float(n)*np.pi*f:(n2)*1j] y+=np.sin(x) z=fft.fft(y[:n]) plt.subplot(211) px=np.arange(n)/float(n) px2=np.arange(n2)/float(n) plt.plot(px2,y,"r-") y2=np.zeros(n2) for f in range(n//2): x=np.ogrid[0:2*float(n2)/float(n)*np.pi*f:(n2)*1j] #y2+= z[f]*(np.cos(x)+1j*np.sin(x))/n y2+= (z[f].real*np.cos(x)- z[f].imag*np.sin(x))/n #y2+= (z[f]*np.exp(x*1j)/n) plt.plot(px2[::8],2*y2.real[::8],"bo") plt.title("%s"%str(flist)) plt.subplot(212) plt.plot(px,z.real,"r-") plt.plot(px,z.imag,"b-")