FFT（高速フーリエ変換）は化学のみならず音声処理（雑音除去など）、画像処理（MRIなど）等々身の回りで活躍しています。
しかし、「高速って、どれだけ早いんだ？」「何をしているかは全く知らないよ。」という方も多いかと思います。
このページではFFTがどれくらい早くて、なぜ早いのかを解説し、さらにpythonによるFFTのプログラミング方法も紹介します。
FFT（高速フーリエ変換）とは
FFTとは高速フーリエ変換（Fast Fourier Transform）のことで、文字通りフーリエ変換を高速に行うアルゴリズムのことです。
通常のフーリエ変換の場合、$\{x_k\}$を$\{X_k\}$へフーリエ変換するには、$X_l=\sum_k^Nx_k*\{\cos(2\pi{kl/N})-j\sin(2\pi{kl/N})\}$を計算すればよく（$N$個はデータ点数）、pythonだとわずか30行ほどでプログラミングできます（いわゆるDFT（離散フーリエ変換）アルゴリズムです）。
DFTのプログラミングコードはこちら
$2^m$のデータ（$[1,1,\cdots,1,0,0,\cdots,0]$）をDFTでフーリエ変換します。($m$はプログラミングコード内で変更してください)。
import math
import numpy as np
import time
m = 8
n = 2**m
xr = np.zeros((n,m+1),dtype = complex)
a = np.zeros(2**m,dtype = complex)
for i in range(n):
xr[i][0] = complex(1,0)
for i in range(n//2,n):
xr[i][0] = complex(0,0)
pi = math.acos(-1.0)
theta = -2*pi/(2**m)
w = complex(math.cos(theta),math.sin(theta))
start = time.perf_counter()
for j in range(2**m):
for k in range(2**m):
wp = w ** (j * k)
a[j] = a[j] + xr[k][0] * wp
end = time.perf_counter()
for j in range(2**m):
print('{0.real:.2f}+{0.imag:.2f}i'.format(a[j]))
print('{:.2f}'.format(end-start))しかし、このDFTアルゴリズムは遅いという致命的な問題があります。例えば、近年のFT-IRではシャープなスペクトルを得るために、データ点数$2^{16}$以上の測定を$2^8$回程度行っています。（日本分光学会 測定法シリーズ 「フーリエ変換赤外分光法」より）。しかし、先ほどのDFTアルゴリズムではデータ点数$2^{16}$個のフーリエ変換$1$回にノートパソコンでは$40$分以上要してしまいます。
もし、FT-IRが測定毎に$40$分かかるフーリエ変換を$2^8=216$回もしてたら使えたものじゃないね
FFTはこのDFTの欠点を改善した素晴らしいアルゴリズムで、$2^{18}$個のデータ点でも$15$秒ほどでフーリエ変換ができます。もちろんDFTとFFTで得られる結果自体は変わりません。
- FFTはフーリエ変換を高速に行うアルゴリズム
- FFTとDFTで得られる結果自体は変わらない
FFTのプログラミングコードはこちら
$2^m$のデータ（$[1,1,\cdots,1,0,0,\cdots,0]$）をDFTでフーリエ変換します。($\text{m}$はプログラミングコード内で変更してください)。$\text{ir}=\pm1$でフーリエ変換と逆フーリエ変換を切り替え可能です。$\text{ir}=1$での出力結果’input.csv’を読み込めば、$\text{ir}=-1$の逆フーリエ変換で元の関数に戻せます。（その際はコメントアウト部分を有効化してください）
def calc_p(n,m,r):
n_b = [0]*m
n_b_s = [0]*m
n_b_i = [0]*m
for i in range(m):
n_b[m-1-i] = n % 2
n = n //2
for i in range(0,m-r):
n_b_s[i+r] = n_b[i]
n_b_i = n_b_s[::-1]
p = 0
for i in range(m):
p = p + n_b_i[i] * (2**(m-i-1))
return p
import math
import numpy as np
import time
import csv
m=3
n=2**m
ir = 1
pi = math.acos(-1.0)
xr = np.zeros((n,m+1),dtype = complex)
#with open('input.csv') as f:
# reader = csv.reader(f, quoting=csv.QUOTE_NONNUMERIC)
# l_f = [ row for row in reader]
#for i in range(n):
# xr[i][0] = complex(l_f[i][0],l_f[i][1])
for i in range(n):
xr[i][0] = complex(1,0)
for i in range(n//2,n):
xr[i][0] = complex(0,0)
theta = -2*pi/(2**m)
w = complex(math.cos(theta),ir * math.sin(theta))
#w = complex(math.cos(theta),0)
start = time.perf_counter()
for r in range(m):
for j in range(2**r):
for i in range(2**(m-(r+1))):
p = calc_p(i+(2**(m-r))*j,m,m-(r+1))
wp = w ** p
k = i + (2**(m-r))*j
l = k + 2**(m-(r+1))
xr[k][r+1] = xr[k][r] + wp * xr[l][r]
xr[l][r+1] = xr[k][r] - wp * xr[l][r]
end = time.perf_counter()
if ir == -1:
for i in range(2**m):
xr[i][m] = xr[i][m] / (2**m)
for j in range(2**m):
k = calc_p(j,m,0)
rep = xr[k][m].real
imp = xr[k][m].imag
print("{:.8f}".format(rep),", ","{:.8f}".format(imp))
print("time ",'{:.2f}'.format(end-start))なぜFFTは高速なアルゴリズムなのか
FFTでは、なぜフーリエ変換が高速になるのでしょうか？
その理由は、「途中計算をうまくまとめて行うから」です。
簡単な例として、$f=1*4+2*4+3*3+5*3$を一つずつ計算すると掛け算と足し算が4回ずつ必要ですが、$f=(1+2+3+5)*4=40$とまとめると計算回数を減らせます。FFTでも計算をまとめることで計算回数を減らすことで高速化を行っています。
もう少しだけ理論的に説明すると
$N=4$のフーリエ変換を使って少し理論的に説明します。
DFTでは$X_l=\sum_{k=0}^3x_kW^{kl}$なので、$\{X_l\}$を全て得るのに$4\times4$回の足し算が必要です。これは$\begin{pmatrix}X_0\\X_1\\X_2\\X_3\end{pmatrix}$$=\begin{pmatrix}1&1&1&1\\1&-1&1&-1\\1&-j&-1&j\\1&j&-1&-j\end{pmatrix}\begin{pmatrix}x_0\\x_1\\x_2\\x_3\end{pmatrix}$と書くことができます。
FFTではこの演算を$\begin{pmatrix}X_0\\X_1\\X_2\\X_3\end{pmatrix}$$=\begin{pmatrix}1&1&0&0\\1&-1&0&0\\0&0&1&-j\\0&0&1&j\end{pmatrix}\begin{pmatrix}1&0&1&0\\0&1&0&1\\1&0&-1&0\\0&1&0&-1\end{pmatrix}\begin{pmatrix}x_0\\x_1\\x_2\\x_3\end{pmatrix}$と整理し、$0$の部分を増やすことでトータルの演算を減らしています。
【N=8で具体化】FFTアルゴリズムは3ステップでできる
実際に$N=8$のデータの場合を使ってFFTのアルゴリズムを解説します。このページでは理論的な背景は説明せず、pythonコードからFFTのアルゴリズムの説明を行います。
ここから少し話がややこしくなりますので注意です
バタフライ演算で処理を進める
$8=2^3$なので、FFTでは$3$回に処理を分けてフーリエ変換を行います。FFTのアルゴリズムを図に置き換えたのが「バラフライ演算」と呼ばれる表です。
実際のプログラミングコードを使ってバタフライ図の見方を把握しましょう
例えば、「1回目処理」で$\{x_i^1\}$から$\{x_i^2\}$を作る部分はプログラミングでは次のように書くことができます。wp2はバタフライ図の$W^2$を表します。
xr[0][2] = xr[0][1] + xr[2][1]
xr[1][2] = xr[1][1] + xr[3][1]
xr[2][2] = xr[0][1] - xr[2][1]
xr[3][2] = xr[1][1] - xr[3][1]
xr[4][2] = xr[4][1] + wp2* xr[6][1]
xr[5][2] = xr[5][1] + wp2* xr[7][1]
xr[6][2] = xr[4][1] - wp2* xr[6][1]
xr[7][2] = xr[5][1] - wp2* xr[7][1]このようにバタフライ図では、必要な足し算の組み合わせと重みづけの係数（コードではwp2）を視覚的に表現しています。
データ数$2^M$、$R=1$での処理に一般化した場合、バタフライ図は次のように書くことができ、プログラミングも数行で書くことができます。
theta = -2*pi/(2**m) # データ数が細かくなると角度も細かくなる
w = complex(math.cos(theta),ir * math.sin(theta))
# ir = 1はフーリエ変換、 ir = -1は逆フーリエ変換
for r in range(m): # 処理はm回行う
for j in range(2**r): # r回目の処理ではブロックは2**r個
for i in range(2**(m-(r+1))): # バタフライは2**(m-(r+1))組
p = calc_p(i+(2**(m-r))*j,m,m-(r+1))
# pは別関数で作成する
wp = w ** p # w を p乗して wpを作る
k = i + (2**(m-r))*j # j番目ブロックの初数はi+2**(m-r)*j
l = k + 2**(m-(r+1)) # kとペアを作るのはl = k + 2**(m-(r+1))
xr[k][r+1] = xr[k][r] + wp * xr[l][r]
xr[l][r+1] = xr[k][r] - wp * xr[l][r]重み係数Wpを得る方法
次に、重み係数$wp=w^p$の値が分かれば$r$回目の処理を実際に実行することができます。
このパラメーター$p$は、次の3ステップの通り、$K$を2進数化させた後にビットシフトと反転処理を行うことで得ることができます。
例えば、$K=5$であれば$5=101_{(2)}$、$K=6$であれば$6=110_{(2)}$になります。
例えば、データ数$2^3=8$での$K=5$の$R=1$回目処理であれば、$3-1-1=1$なので$5=101_{(2)}$を右に$1$つシフトさせ、$010_{(2)}$とします。$K=6$の$R=2$回目処理であれば、$3-2-1=0$なのでシフトさせず、$110_{(2)}$のままです。
例えば、$K=5$の$R=1$回目処理であれば、ビット反転後$010_{(2)}$なので$p=2^1=2$です。$K=6$の$R=2$回目処理であれば、ビット反転後$011_{(2)}$なので、$p=2^1+2^0=3$です。
このデータ数$2^m$の場合、$r$回目処理における$K=n$での$p$を得る手続きをpythonでプログラミングすると次のような関数にできます。
def calc_p(n,m,r):
n_b = [0]*m
n_b_s = [0]*m
n_b_i = [0]*m
for i in range(m):
n_b[m-1-i] = n % 2
n = n //2
for i in range(0,m-r):
n_b_s[i+r] = n_b[i]
n_b_i = n_b_s[::-1]
p = 0
for i in range(m):
p = p + n_b_i[i] * (2**(m-i-1))
return p最終的な並び替え
無事に$m$回の処理を終えた後に得られる$\{X_i^m\}$は実は求めたい$\{x_i\}$のフーリエ変換$\{X_i\}$の順番に並んでいません。
そのため、$\{X_i^m\}$を正しい順番に並べる際にもビット反転処理を行う必要があります。
例えば、データ数$2^3$の場合における$X_6^3$は$6=110_{(2)}$なので、ビット反転させると$011_{(2)}=3$となり、正しい順番としては$X_3$に対応することになります。同様に、$X_1^3$は$1=001_{(2)}$なので、ビット反転させると$X_4$に対応することがわかります。
ビット反転処理が必要な理由について
この理由はFFTではフーリエ変換$X_l=\sum_k^Nx_k*W^{kp}$に対して$k=0,1,2,\cdots$と順番通りに足し算を行わずに、効率よく計算するために順番を入れ替えて計算を行っているためです。下記文献によると、順番通りに出力されるとFFTのアルゴリズムとならないことが書かれています。
ビットの逆転をさけるために式(4.32)の和の順番を逆に$n_0$から行うと、（中略）、FFTのアルゴリズムにならない。
東京電機大学出版局 中村尚五著 ビギナーズデジタルフーリエ変換 P132より
まとめ
このページではN=8のデータ数を例にしてFFTのアルゴリズムを解説しました。
FFTは化学のみならず音声処理（雑音除去など）、画像処理（MRIなど）等々身の回りで活躍しています。
足し算の順番を工夫することで40分以上かかるようなフーリエ変換を数秒で終わらせるFFTのすばらしさを感じていただければと思います。
