バタフライ図も行列も見たくない人のための高速フーリエ変換

高速フーリエ変換（FFT）を実装しようと思って、まず理論を知らねばと調べたところ、まったくチンプンカンプンだったのですが、とりあえず、なんとか証明を追うことができました。というわけで、バタフライ演算の図を見てもなんもわからん人が FFT を実装するための資料として、この記事を書き残したいと思います。

参考にした資料は、富永氏による GNU Scientific Library に関するドキュメントの翻訳の「FFT アルゴリズム」の「基数2の時間空間法」です。周波数空間については追っていません（どっちでやっても計算量は同じなので）。

なぜ離散フーリエ変換（DFT）を高速化できるのか？

この記事内の式で使用する文字を定義しつつ、 FFT の概要を紹介していきます。

とりあえず、 DFT の式を出しておきましょう（フーリエ変換自体を解説できるほどの理解はないので、ごめんなさい）。

$\displaystyle DFT_a(\vec{x}) = \sum_{b=0}^{N-1} x_b e^{-2 \pi i \frac{ab}{N}}$

記号	意味
$\vec{x}$	入力データ（実数列）
$N$	入力データ数
$a$	出力のインデックス（周波数ビン） $0, 1, \cdots, N-1$
$b$	入力のインデックス（時刻） $0, 1, \cdots, N-1$
$x_b$	入力データの $b$ 番目（0オリジン）
$i$	虚数単位

周波数ビン $a$ に対するフーリエ係数 $DFT_a(\vec{x})$ を求めるには、データ数だけループを回す必要があります。これをすべての周波数ビンについて求めるならば、二重ループが必要になります。つまり、計算量は $O(N^2)$ であり、とても重い処理だということです。 FFT では、これを頑張って $O(N \log N)$ まで下げます。これでなんとか大量のデータを実用的に解析することができるというわけです。

次の説明に行く前に、指数関数の部分を文字で置き換えましょう。いわゆる回転因子というやつです。

$\displaystyle W_N = e^{-2 \pi i \frac{1}{N}}$

指数法則が効くので、例えば $W_N^2 = e^{-2 \pi i \frac{2}{N}}$ となります。

回転因子の意味ですが、複素平面上の単位円を何個の点に分割するか、ということです。例えば、 $N=4$ ならば、次の図のように分割され、 0 乗のとき 1 で、それから 1 乗されるごとに時計回りに移動していきます。 1 周まわり終わると、また同じところをぐるぐる回るだけなので、 $W_N^N=W_N^0=1, W_N^{N+1} = W_N^1$ だったりします（もっと一般化すると $W_N^x = W_N^{x \bmod N}$ ）。

f:id:azyobuzin:20190602024727p:plain — 回転因子 N=4

今導入した $W_N$ を使って DFT の式を書き直すと、次のようになります。

$\displaystyle DFT_a(\vec{x}) = \sum_{b=0}^{N-1} x_b W_N^{ab}$

さて、新たな文字の導入も終わったので、計算量を減らすトリックの紹介です。 DFT の計算は、 $b$ を偶数と奇数に分けて考えると、ループ回数を半分にできるという性質があります*1。というわけで、式変形で、偶数と奇数に分けてみましょう。ただし、 $N$ が 2 の倍数でないと添え字が合わなくなってしまいます（これを繰り返すので、データ数は 2 のべき乗である必要があります）。第1項が偶数のときの計算、第2項が奇数のときの計算になっています。

$\displaystyle DFT_a(\vec{x}) = \sum_{b=0}^{\frac{N}{2} - 1} x_{2b} W_N^{2ab} + \sum_{b=0}^{\frac{N}{2} - 1} x_{2b+1} W_N^{a(2b+1)}$

ここからさらに式変形をしていきます。 $W_N^{2ab}$ は指数法則で $\left( W_N^2 \right) ^{ab}$ と分解できます。そして、 $W_N^2$ のほうを変形していきます。

$\displaystyle \begin{eqnarray} W_N^2 &=& e^{-2 \pi i \frac{2}{N}} = e^{-2 \pi i \frac{1}{N/2}} \\ &=& W_{N/2} \end{eqnarray}$

といった感じで $W_N^2 = W_{N/2}$ ということが判明したので、 $W_N^{2ab} = W_{N/2}^{ab}$ です。

さらに $W_N^{a(2b+1)}$ のほうも、指数法則で変形していきます。

$\displaystyle W_N^{a(2b+1)} = W_N^{2ab} W_N^a = W_{N/2}^{ab} W_N^a$

これで、もう一度 DFT の式を書き直します。

$\displaystyle \begin{split} DFT_a(\vec{x}) &= \sum_{b=0}^{\frac{N}{2} - 1} x_{2b} W_N^{2ab} + \sum_{b=0}^{\frac{N}{2} - 1} x_{2b+1} W_N^{a(2b+1)} \\ &= \sum_{b=0}^{\frac{N}{2} - 1} x_{2b} W_{N/2}^{ab} + W_N^a \sum_{b=0}^{\frac{N}{2} - 1} x_{2b+1} W_{N/2}^{ab} \end{split}$

これを、最初の DFT の式 $\sum_{b=0}^{N-1} x_b W_N^{ab}$ と見比べてみてください。長さが半分になった DFT を 2 個足し合わせた形 $DFT_a(偶数) + W_N^a DFT_a(奇数)$ になっています！各項の DFT の式もさらに半分に分けることができるので、どんどん分割していけます。例として、入力データを 8 個与えた場合を考えてみましょう。

$\displaystyle \begin{split} &DFT_a(x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7) \\ &= DFT_a(x_0, x_2, x_4, x_6) + W_8^a DFT_a(x_1, x_3, x_5, x_7) \\ &= DFT_a(x_0, x_4) + W_4^a DFT_a(x_2, x_6) \\&\quad + W_8^a \left\{ DFT_a(x_1, x_5) + W_4^a DFT_a(x_3, x_7) \right\} \\ &= x_0 + W_2^a x_4 + W_4^a (x_2 + W_2^a x_6) \\&\quad + W_8^a \left\{ x_1 + W_2^a x_5 + W_4^a (x_3 + W_2^a x_7) \right\} \end{split}$

いい感じに分割できました。この分割する作業を $\log_2 N$ 回行い、すべての周波数ビン $a$ に対して計算するので、計算量が $O(N \log N)$ になります。

分割とインデックスのビットとの関係

今の分割を図に表すと、次のような感じになっています。

f:id:azyobuzin:20190602052931p:plain — 分割の様子

説明の都合上、一番下を0段目と呼ぶことにしました。では、0段目の並び順には、どのような規則があるのでしょうか？それは、インデックスを2進法で表すとわかります。

10進	0	4	2	6	1	5	3	7
2進	000	100	010	110	001	101	011	111

これはちょうど、0～7を2進法で表したときのビットを、左右反転させた並びになっています。

0段目から1段目を求める

ここからは、分割した結果このようにな並びになったのではなく、最初からこの並び順になっているものとして、話を進めます（つまり事前にこのように並べ替えます）。この並び順から、 $DFT_a$ を求めるには、分割結果の式から、次のようにすればよいことがわかります。

f:id:azyobuzin:20190602050730p:plain — 分割結果から計算する様子

求める結果は、すべての周波数ビン $a$ についてですので、そこに注目しながら考えていきます。

まずは、0段目のデータを使って1段目を求めることを考えます。いったん元の順番のことは忘れて、 $g_m(b, a)$ とは、周波数ビン $a$ について計算するときに、 $m$ 段目の左から $b$ 番目（0オリジン）のデータを表すものとします。図通りに式を書くならこんな感じになりそうですね。

$\displaystyle g_1(b, a) = g_0(2b) + W_2^a g_0(2b+1)$

（ $g_0(b, a)$ は、 $a$ にまったく関係なく同じ値を返すので $g_0(b)$ と表しています。）

ここで、 $W_2^a$ について観察すると、 $W_2^a = W_2^{a \bmod 2}$ なので、 $a$ を2進法で表して $a = 2^2 a_2 + 2^1 a_1 + 2^0 a_0 = [ a_2 a_1 a_0 ]$ としたとき（以降、2進数を [ ] で表します）、 $a \bmod 2 = a_0$ より、 $W_2^a = W_2^{a_0}$ ということがわかります。つまり、すべての $a$ について求めるとしても、 $a_0$ しか関係ありません。

$\displaystyle g_1(b, a) = g_1(b, [a_0])$

このことから、すべての $a$ について $g_1(b, a)$ を求めた結果は、次のように格納することができます。ここでは、 $W_2^0 = 1, W_2^1 = -1$ を利用しています。

$b=[b_1 b_0]$ → ↓ $a=[a_0]$	0 = [00]	1 = [01]	2 = [10]	3 = [11]
0	$g_0(0)+g_0(1)$	$g_0(2)+g_0(3)$	$g_0(4)+g_0(5)$	$g_0(6)+g_0(7)$
1	$g_0(0)-g_0(1)$	$g_0(2)-g_0(3)$	$g_0(4)-g_0(5)$	$g_0(6)-g_0(7)$

必要な容量は、表の大きさから 8 個、すなわち入力データ数と同じです。そして、2段目を計算するとき、0段目のデータははもう必要ありません。というわけで、0段目のデータが入っていた領域を、1段目の計算結果を入れておく場所として再利用しようということを考えます。

ある $b = [b_1 b_0]$ について、 $a_0 = 0, 1$ それぞれについて求めるとき、表内の式に表れるインデックスを2進法で表すと、次の法則があります。

$\displaystyle \begin{eqnarray} g_1([b_1 b_0], [0]) &=& g_0([b_1 b_0 0]) + g_0([b_1 b_0 1]) \\ g_1([b_1 b_0], [1]) &=& g_0([b_1 b_0 0]) - g_0([b_1 b_0 1]) \end{eqnarray}$

この式から次のことが言えます:

$g_0([b_1 b_0 0]), g_0([b_1 b_0 1])$ にアクセスするのは $g_1([b_1 b_0], [0]), g_1([b_1 b_0], [1])$ を求めるときだけ。

ここで、作業用のメモリを $g(b)$ としましょう。1段目を計算する前は $g(b) = g_0(b)$ です。 $g_1([b_1 b_0], [0]), g_1([b_1 b_0], [1])$ を求めた後、もう $g([b_1 b_0 0]), g([b_1 b_0 1])$ にアクセスすることはないということは、これと同じ領域 $g([b_1 b_0 0]), g([b_1 b_0 1])$ に結果を保存しておけばよいということになります。といっても、式で書いていると、メモリの状況をうまく表現できないので、手続き的に書き表してみます。一時変数 $t_0, t_1$ を用いて

$t_0 \leftarrow g_1([b_1 b_0], [0]) = g([b_1 b_0 0]) + g([b_1 b_0 1])$
$t_1 \leftarrow g_1([b_1 b_0], [1]) = g([b_1 b_0 0]) - g([b_1 b_0 1])$
$g([b_1 b_0 0]) \leftarrow t_0$
$g([b_1 b_0 1]) \leftarrow t_1$

といった感じの手順で、うまくメモリを上書きできます。図で表すとこんな感じでしょうか。

f:id:azyobuzin:20190602165828p:plain — $g_1(0, [a_0])$ を求めるときのメモリ

1段目から2段目を求める

次の段も同じようにできそうなので、やっていきましょう。図から式を抽出すると、こんな感じです。

$\displaystyle g_2(b, a) = g_1(2b, a) + W_4^a g_1(2b+1, a)$

2段目は2個しかないので、 $b=0,1$ です。データ8個だと簡単になりすぎてしまって、一般化まで持っていくのが大変ですが、データ数を増やすと図がつらくなるので許してください。

1段目と同様に、 $W_4^a$ は、回転因子の周期によって、次の変形ができます。

$\displaystyle W_4^a = W_4^{a \bmod 4} = W_4^{[a_1 a_0]}$

さらにもう一歩、指数法則で分解しておきます。

$\displaystyle W_4^{[a_1 a_0]} = W_4^{2a_1 + [a_0]} = W_4^{2a_1} W_4^{[a_0]}$

というわけで、これと、 $g_1(b, a)=g_1(b, a_0)$ を利用すると、 $g_2(b,a)$ は、次のようになります。

$\displaystyle g_2(b,a) = g_2(b, [a_1 a_0]) = g_1(2b, a_0) + W_4^{2a_1} W_4^{[a_0]} g_1(2b+1, a_0)$

さらに、 $a_1$ で場合分けすると、 $W_4^0 = 1, W_4^2 = -1$ より、次のようになります。

$\displaystyle g_2(b, [a_1 a_0]) = \begin{cases} g_1(2b, a_0) + W_4^{a_0} g_1(2b+1, a_0) & (a_1 = 0) \\ g_1(2b, a_0) - W_4^{a_0} g_1(2b+1, a_0) & (a_1 = 1) \end{cases}$

では、すべての $a,b$ について計算してみましょう。

$b$ → ↓ $a=[a_1 a_0]$	0	1
0 = [00]	$g_1(0, 0) + W_4^0 g_1(1, 0)$	$g_1(2, 0) + W_4^0 g_1(3, 0)$
1 = [01]	$g_1(0, 1) + W_4^1 g_1(1, 1)$	$g_1(2, 1) + W_4^1 g_1(3, 1)$
2 = [10]	$g_1(0, 0) - W_4^0 g_1(1, 0)$	$g_1(2, 0) - W_4^0 g_1(3, 0)$
3 = [11]	$g_1(0, 1) - W_4^1 g_1(1, 1)$	$g_1(2, 1) - W_4^1 g_1(3, 1)$

1段目の結果はひとつの配列（作業用メモリ）に保存されているのでした。したがって、この表から、同じメモリへのアクセスする計算が2つずつあることがわかります。まとめると、こういうことです。

$g_1([b_0 0], [a_0]), g_1([b_0 1], [a_0])$ にアクセスするのは、 $g_2([b_0], [0 a_0]), g_2([b_0], [1 a_0])$ を求めるときだけ。

メモリのインデックスの関係から

$\displaystyle \begin{eqnarray} g_1(2b, [a_0]) = g_1([b_0 0], [a_0]) &=& g([b_0 0 a_0]) \\ g_1(2b+1, [a_0]) = g_1([b_0 1], [a_0]) &=& g([b_0 1 a_0]) \end{eqnarray}$

なので、 $g([b_0 0 a_0]), g([b_0 1 a_0])$ のデータを使って計算し、その結果をまた $g([b_0 0 a_0]), g([b_0 1 a_0])$ に保存すればやりくりできそうです。ちょうど、 $a_1$ について、場合分けしてあり、 $a_1=0,1$ のそれぞれについて一緒に計算できるよう準備をしておいたので、あとは結果を $g([b_0 a_1 a_0]) \leftarrow g_2([b_0], [a_1 a_0])$ のように保存すればよいわけです。手続きで考えると、次のようになります。

$t_0 \leftarrow g_2([b_0], [0 a_0]) = g([b_0 0 a_0]) + W_4^{[a_0]} g([b_0 1 a_0])$
$t_1 \leftarrow g_2([b_0], [1 a_0]) = g([b_1 1 a_0]) - W_4^{[a_0]} g([b_0 1 a_0])$
$g([b_0 0 a_0]) \leftarrow t_0$
$g([b_0 1 a_1]) \leftarrow t_1$

f:id:azyobuzin:20190602170522p:plain — $g_2(0, [a_1 0])$ を求めるときのメモリ

m-1段目からm段目を求める

さて、大体一般化する準備はできたと思います。まず計算の図から明らかっぽいのは、次の式です。

$\displaystyle g_m(b, a) = g_{m-1}(2b, a) + W_{2^m}^a g_{m-1}(2b+1, a)$

ここで、 $b$ は、1段進むごとに半分になるので、 $0,1,\cdots,\frac{N}{2^m} - 1$ です。

そして、回転因子の周期を使った変形ができるのでした。

$\displaystyle W_{2^m}^a = W_{2^m}^{a \bmod 2^m} = W_{2^m}^{[a_{m-1} a_{m-2} \cdots a_0]}$

さらに、 $a_{m-1}$ を外に出します。

$\displaystyle \begin{split} W_{2^m}^{[a_{m-1} a_{m-2} \cdots a_0]} &= W_{2^m}^{2^{m-1} a_{m-1} + [a_{m-2} \cdots a_0]} \\ &= W_{2^m}^{2^{m-1} a_{m-1}} W_{2^m}^{[a_{m-2} \cdots a_0]} \end{split}$

このようにすることで、次のような場合分けができるようになります。

$\displaystyle W_{2^m}^{2^{m-1} a_{m-1}} = \begin{cases} W_{2^m}^0=1 & (a_{m-1}=0) \\ W_{2^m}^{2^{m-1}}=W_2^1=-1 & (a_{m-1}=1) \end{cases}$

ここまでの変形をまとめると、次のようになります。

$\displaystyle g_m(b, a) = g_{m-1}(2b, a) + W_{2^m}^{2^{m-1} a_{m-1}} W_{2^m}^{[a_{m-2} \cdots a_0]} g_{m-1}(2b+1, a)$

このとき、 $g_m$ は $a$ について、下位 $m$ ビットしか利用しないということがわかります。ということは、 $g_{m-1}$ では、下位 $m-1$ ビットしか利用しないですね。というわけで、ビットごとに考えるとこうなります。

$\displaystyle \begin{split} g_m(b, [a_{m-1} \cdots a_0]) &= g_{m-1}(2b, [a_{m-2} \cdots a_0]) \\ & + W_{2^m}^{2^{m-1} a_{m-1}} W_{2^m}^{[a_{m-2} \cdots a_0]} g_{m-1}(2b+1, [a_{m-2} \cdots a_0]) \end{split}$

さらに、 $a_{m-1}$ の場合分けによって、関数を分けます。

$\displaystyle \begin{eqnarray} g_m(b, [0 a_{m-2} \cdots a_0]) &=& g_{m-1}(2b, [a_{m-2} \cdots a_0]) \\&+& W_{2^m}^{[a_{m-2} \cdots a_0]} g_{m-1}(2b+1, [a_{m-2} \cdots a_0]) \\ g_m(b, [1 a_{m-2} \cdots a_0]) &=& g_{m-1}(2b, [a_{m-2} \cdots a_0]) \\&-& W_{2^m}^{[a_{m-2} \cdots a_0]} g_{m-1}(2b+1, [a_{m-2} \cdots a_0]) \end{eqnarray}$

また、 $b$ のほうもビットで考えると、 $N=2^n$ のとき、 $b=[b_{n-m-1} b_{n-m-2} \cdots b_0]$ とします。

それでは、またメモリアクセスについて考えます。 $g_m$ の式から、 $2b=[b_{n-m-1} \cdots b_0 0], 2b+1=[b_{n-m-1} \cdots b_0 1]$ を用いて、次のことが言えます。

$g_{m-1}([b_{n-m-1} \cdots b_0 0], [a_{m-2} \cdots a_0]), g_{m-1}([b_{n-m-1} \cdots b_0 1], [a_{m-2} \cdots a_0])$ にアクセスするのは、 $g_m([b_{n-m-1} \cdots b_0], [0 a_{m-2} \cdots a_0]), g_m([b_{n-m-1} \cdots b_0], [1 a_{m-2} \cdots a_0])$ を求めるときだけ。

2段目までの議論から、天下り式に $g_m(b,a)$ の結果を保存する場所 $g(b)$ を設定します*2。上位ビットが $b$ 、下位ビットが $[a_{m-1} \cdots a_0]$ になります。

$\displaystyle g([b_{n-m-1} \cdots b_0 a_{m-1} \cdots a_0]) \leftarrow g_m([b_{n-m-1} \cdots b_0], [a_{m-1} \cdots a_0])$

この結果から、どの $g_m$ においても、 $g$ の大きさは $n$ ビット、すなわち $N$ であればよいということがわかります。

では、 $g_{m-1}(2b, [a_{m-2} \cdots a_0]), g_{m-1}(2b+1, [a_{m-2} \cdots a_0])$ をメモリアクセスとして書き直してみます。

$\displaystyle \begin{eqnarray} g_{m-1}(2b, [a_{m-2} \cdots a_0]) &=& g([b_{n-m-1} \cdots b_0 0 a_{m-2} \cdots a_0]) \\ g_{m-1}(2b+1, [a_{m-2} \cdots a_0]) &=& g([b_{n-m-1} \cdots b_0 1 a_{m-2} \cdots a_0]) \end{eqnarray}$

この作業用メモリのインデックスは、ちょうど $g_m(b, [0 a_{m-2} \cdots a_0]), g_m(b, [1 a_{m-2} \cdots a_0])$ に対応するインデックスと同じところです。したがって、2段目までと同じように、2つ取り出し、計算して、格納し直す方法が使えそうです。そうすると、もう見慣れた手続きですが、このようになりますね（あまりに横に長いので、一時変数を分けました）。

$t_0 \leftarrow g([b_{n-m-1} \cdots b_0 0 a_{m-2} \cdots a_0])$
$t_1 \leftarrow g([b_{n-m-1} \cdots b_0 1 a_{m-2} \cdots a_0])$
$t_2 \leftarrow g_m(b, [0 a_{m-2} \cdots a_0]) = t_0 + W_{2^m}^{[a_{m-2} \cdots a_0]} t_1$
$t_3 \leftarrow g_m(b, [1 a_{m-2} \cdots a_0]) = t_0 - W_{2^m}^{[a_{m-2} \cdots a_0]} t_1$
$g([b_{n-m-1} \cdots b_0 0 a_{m-2} \cdots a_0]) \leftarrow t_2$
$g([b_{n-m-1} \cdots b_0 1 a_{m-2} \cdots a_0]) \leftarrow t_3$

というわけで、1段前の結果を使って、次の段を計算することを一般化することができました。 $m=n$ まで達すると、 $g$ のインデックスは、 $g(a_{m-1} \cdots a_0)$ すなわち $g(a)$ となり、すべての周波数ビンについての DFT 結果が $g$ に入っている状態になるので、これで処理は完了です。

どうやってプログラムとして実現するか

1からn段目まで順番に計算していくわけですが、今までの話だと、まだプログラムとしてどう書くかが見えていない気がするので、もう少し話をします。

重要なポイントは $a_{m-1}=0,1$ それぞれの計算を一緒に行ってしまうことで、メモリをうまくやりくりできることでした。逆に言うと、それ以外の数値 $[b_{n-m-1} \cdots b_0], [a_{m-2} \cdots a_0]$ は、すべてのパターンについて計算を行う必要があります。したがって、 a と b の二重ループを行います。

$a_{m-1}$ は同時に計算を行うので、1段の計算において、 a は $a_{m-2}$ までのビットを使って $0,1,\cdots,2^{m-1}-1$ の範囲で動きます。 b は $0,1,\cdots,2^{n-m} - 1$ ですね。あとは、 b をビットシフトしたりすれば、インデックスが合います。

そもそも最初に、ビット反転順に並べ替えなきゃいけないじゃん、というのもあります。入力と出力の配列を別々にするならば、簡単なビット演算でビットを反転して、そのインデックスに詰め替えることで実現できます。入力と出力に同じ配列を使う場合は、計算自体と同じようなループを回して値を交換していけばできるのですが、もうちょっと計算量的にマシな方法があるらしいのです。しかし、解説できるほど理解していません……。

解説は以上です。ちゃんとした疑似コードは、最初に挙げた参考資料をご覧ください。「ビットシフトしたり」が掛け算で実現されているくらいで、あまり驚かずに読めるようになっているのではないでしょうか。ただし、文字 b の扱い方は結構違うので、注意してください。この記事では最初からビット反転してあるものとして扱ってしまいましたが、参考資料では、最初の並び順で導出したのち、反転するとメモリの格納の都合上便利といった風に解説されています。

*1:周波数のほうを分けてもいいです。あと、偶数と奇数ではなく、割ると割り切れるもの、 1 余るもの……のように一般化することもできますが、割る数は 2 か 4 が効率が良いようです。そして僕の頭では 2 しかわかりません！！

*2:1段目では、ビット反転順のデータを相手に計算したときに、うまくデータを置き換えられるということで説明しました。2段目では、1段目と同じようにインデックスを定めると、うまくいきそうだという話をしました。