自分用作業メモ
概要
以前からFFTを時系列データ（たとえば、地震の発生回数）の『予測』に使えないかと勉強してたのですけれども、無理だとわかったので、次にFFTを使いたくなった時用に経緯を一通りメモ。
でも、補間アルゴリズムとして使う分には良さげです。
FFTとは？
まずは（１次元）FFTのざっくりした説明。
音声に代表されるような時系列データを、周波数ごとに分解するのがフーリエ変換。
その一種、高速フーリエ変換＝FFT。
逆FFTで元のデータに戻すことができます。
元データ → FFT → 周波数によるフィルタリング （ローパスフィルタ、ハイパスフィルタ、などなど） → 逆FFT → 時系列データへ戻す
みたいな流れが一般的です。
なぜFFTはデータ予測に使えないか
大抵のフーリエ変換の入門書に書かれていますが、、、フーリエ変換には
◎ 入力データは、ちょうど１周期分である
という前提をおいたアルゴリズムです。
なので、
音声データ （時系列データ） 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-")