ゆるふわブログ

プログラミング,大学の勉強,日常生活で感じたことをゆるふわに書いていけたらなと思います.技術的に拙いところがあっても温かい目で見守っていただければ幸いです.

高速フーリエ変換 (FFT)

お久しぶりです,Ysmr-Ry です.
試験範囲なので高速フーリエ変換 (FFT) について僕の理解を記しておこうと思います.
高速フーリエ変換とは離散フーリエ変換を高速に計算するためのアルゴリズムで,コンピューターで (離散) フーリエ変換をするために重要な手法です.
フーリエ変換・離散フーリエ変換については,下記のようなサイトで調べればわかると思います.たぶん.僕読んでないけど.

www.ic.is.tohoku.ac.jp

成果物

見て楽しくないとあれかなぁと思ったので一応 Visualizer を作ってみました.
矩形波 (もどき) を Fourier 変換して表示してます.Enter を押すとマイクから音を拾ってきます.
もう一度 Enter を押すとその波形で固定して Fourier 変換してくれます.

CodePen - Fourier Visualizer

また,高速フーリエ変換の理解に欠かせないバタフライダイアグラムを自動で作ってくれるのも作ったので記事と合わせてみると良いかもしれません.↑↓ で数を増減できます.

CodePen - Butterfly Diagram Maker

4 つのポイント

以下の 4 点を押さえていればおそらく FFT を理解できるであろうと思います.

  • 分配法則
  • 約分
  • 二進数
  • 再帰的 (漸化式的) 定義

上の 3 つはすぐわかると思います.4 つ目が曲者なのでそのつもりで.
この 4 点を順に確認した後,離散フーリエ変換の高速化を考えます.

分配法則

ご存知だと思いますが一応.分配法則とは,
 2 \times (1+3) = 2\times1+2\times3 = 2+6 = 8 という風に,かっこの中の各項に掛け算を分配することを言います.

約分

これも大丈夫でしょう.
 \frac{2}{8} = \frac{1}{4} のようにすることです.

二進数

コンピューターと相性のいい数の表し方です.
我々が普段使っているのは十進数で,桁が上がるごとに桁の担っている数が 10 倍になっていきます.一の位,十の位,百の位,千の位…という風に.
二進数はこれが 2 倍になっていきます.一の位,ニの位,四の位,八の位…という風です.
例えば,11 を表したいなら,1+2+8 より,1011 となります.

再帰的 (漸化式的) 定義

ここまでは「何を当たり前のことを」という感じだったかもしれませんが,これはちょっと難しいです.
これは高校で習ったであろう漸化式を思い出してもらえるとわかりやすいです.
ご存知フィボナッチ数列が良さそうですね.
フィボナッチ数列とは,1, 1 から始めて直前の 2 つの数を足し合わせて次を作っていく数列です.
 1, 1, 2, 3, 5, 8, 13, \cdots
左から数えて n 番目のフィボナッチ数列の数を  f_n と書くと,「直前の 2 つの数を足し合わせて次を作っていく」は,
 f_n = f_{n-1} + f_{n-2}\quad(n\geq3)
と表せます.これが漸化式です.
フィボナッチ数列の作り方のように,最初 ( f_1, f_2) から順に作っていけば  f_n を得られますが,今は何もない状態からぽっと  f_n を聞かれた状況を考えてみましょう.最初から数列を作っていないので何もわかりません.
その時この漸化式を見ると,「n 番目を求めるには n-1 番目と n-2 番目を求めて足せばいいのさ!」と語りかけてきます.
そう,左から n 番目を求める問題を,n-1 番目と n-2 番目を求めると言い換えられるのです.
n-1 番目は n-2 番目と n-3 番目がわかればよく,n-2 番目は…と繰り返していけば,いつか 1 番目と 2 番目がわかればわかるというところまできます.1, 2 番目は 1 と決まっていますから,問題が解けたことになります.
このように,より小さな (この場合は左側) な問題の組み合わせに言い換えていくことができるのです.
よって,問題を解くのに必要なのは,より 小さな問題に言い換えるルール (漸化式)もっとも単純な問題の答え (1, 2 番目が 1) だけです.
これを覚えておいてください.

離散フーリエ変換 (DFT)

 W_N = e^{-i\frac{2\pi}{N}} とし,これを回転子と言います.(いきなり何)
難しそうですがなんのことはなく,名前の通り 360° を N 等分した分時計回りに回ることを表します.
例えば  W_4 によって円周を一周すると図のようになります.

f:id:Ysmr_Ry:20171108232907p:plain

はい.ちなみに円の右端を 0° として時計回りに回ってます.番号の順に回ってます.
 W_N^k W_N の回転を  k 回連続で行うことを表します.
ここでは簡単のため, W_N^k を円の中心から回転で行き着いた点を結ぶ矢印で表します.
例えば, W_4^2 ならば,円の右端から半周するので左端にいます.図の 3 の位置です.
よって,矢印 ← で表します.
すると,大きさ 8 の離散フーリエ変換は与えられた  x_i\,(i=1,2,\cdots,8) に対し,次の行列計算を行うことを指します.

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_1 \\ X_2 \\ X_3 \\ X_4 \\ X_5 \\ X_6 \\ X_7
\end{pmatrix} = 
\begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{↘} & \mbox{↓} & \mbox{↙} & \mbox{←} & \mbox{↖} & \mbox{↑} & \mbox{↗} \\
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} & \mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{↙} & \mbox{↑} & \mbox{↘} & \mbox{←} & \mbox{↗} & \mbox{↓} & \mbox{↖} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\mbox{→} & \mbox{↖} & \mbox{↓} & \mbox{↗} & \mbox{←} & \mbox{↘} & \mbox{↑} & \mbox{↙} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} & \mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\mbox{→} & \mbox{↗} & \mbox{↑} & \mbox{↖} & \mbox{←} & \mbox{↙} & \mbox{↓} & \mbox{↘} \\
\end{pmatrix}\begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \\ x_7
\end{pmatrix}
\end{align*}

行列の 2 行目は 8 等分を 1 個飛ばしで,k 行目は 8 等分を k-1 個飛ばしで回っていることがわかります.
これは通常  O(n^2) の計算量になるのですが,隠れた対称性を見出すことでこの計算量を  O(n\log n) まで落とすのが高速フーリエ変換です.
バブルソートクイックソートぐらい違いますね.

高速フーリエ変換 (FFT)

高速フーリエ変換を一言でいうと,「 \frac{2}{8} って約分して  \frac{1}{4} にできんじゃん!」です.たぶん.
つまり, W_8^2 で 8 等分を 2 個飛ばしで回るのと, W_4^1 で 4 等分を順に回るのとでは区別がつきません.

f:id:Ysmr_Ry:20171109002725p:plain

f:id:Ysmr_Ry:20171108232907p:plain

これです.
これによって,行を半分に分けた時,回転していく様子が "対等" となり分配法則によってくくれるのです.
サイズの小さい場合から具体的に見てみましょう.

n = 2 の場合

離散フーリエ変換は,

 \begin{align*}
\begin{pmatrix} X_0 \\ X_1 \end{pmatrix} = \begin{pmatrix}
\mbox{→} & \mbox{→} \\ \mbox{→} & \mbox{←}
\end{pmatrix}\begin{pmatrix}
  x_0 \\ x_1
\end{pmatrix}
\end{align*}

となります.
ここで,行を半分 (1 個ずつ) に分けてみましょう.
1 行目はどちらも同じです.
2 行目は右側は左側を半周 (←) 回したものです.
よって,
 \begin{align*}
\begin{pmatrix} X_0 \\ X_1 \end{pmatrix} &= \begin{pmatrix}
\mbox{→} x_0 + \mbox{→} x_1 \\
\mbox{→} x_0 + \mbox{←} x_1 \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→}(x_0 +x_1) \\
\mbox{→}(x_0 + \mbox{←} x_1) \\
\end{pmatrix} \\
&= \mbox{→}
\begin{pmatrix}
x_0 +x_1 \\ x_0 + \mbox{←} x_1
\end{pmatrix} \\
&=
\begin{pmatrix}
x_0 +x_1 \\ x_0 + \mbox{←} x_1
\end{pmatrix}
\end{align*}

 \mbox{→} = W_2^0 = 1 ですからこうなります.
これを図にしてみるとこうなります.

f:id:Ysmr_Ry:20171109004714p:plain

左から線をたどっていきます.合流したら足します.
三角があったら流れてきた数に傍に書いてある回転子を書きます.
この図は, X_0 = x_0 + W_2^0 x_1 = x_0 + x_1 X_1 = x_0 + W_2^1 x_1 = x_0 + \mbox{←} x_1 を表しています.

この図を,蝶のようなのでバタフライダイアグラムといいます.

n = 4 の場合

離散フーリエ変換は,

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_1 \\ X_2 \\ X_3
\end{pmatrix} = \begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix}
\end{align*}

今度は偶数行目と奇数行目で考えてみましょう.偶数行目は行を半分に分けた時に同じものの繰り返しですが,奇数行目は半周分ずれています.
よって.

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_2
\end{pmatrix} &= \begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} x_0 + \mbox{→} x_1 + \mbox{→} x_2 + \mbox{→} x_3 \\
\mbox{→} x_0 + \mbox{←} x_1 + \mbox{→} x_2 + \mbox{←} x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} (x_0+x_2) + \mbox{→} (x_1+x_3) \\
\mbox{→} (x_0+x_2) + \mbox{←} (x_1+x_3) \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{←} \\
\end{pmatrix} \begin{pmatrix}
x_0+x_2 \\ x_1+x_3
\end{pmatrix} \\
\begin{pmatrix}
X_1 \\ X_3
\end{pmatrix} &= \begin{pmatrix}
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} x_0 + \mbox{↓} x_1 + \mbox{←} x_2 + \mbox{↑} x_3 \\
\mbox{→} x_0 + \mbox{↑} x_1 + \mbox{←} x_2 + \mbox{↓} x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} (x_0+\mbox{←}x_2) + \mbox{↓} (x_1+\mbox{←}x_3) \\
\mbox{→} (x_0+\mbox{←}x_2) + \mbox{↑} (x_1+\mbox{←}x_3) \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} & \mbox{↓} \\
\mbox{→} & \mbox{↑} \\
\end{pmatrix} \begin{pmatrix}
x_0+\mbox{←}x_2 \\ x_1+\mbox{←}x_3
\end{pmatrix} \\
\end{align*}

ここまでが第 1 段階です.これは n = 2 の場合とやっていることが全く同じです.
n = 2 の場合の  x_0 +x_1 が偶数行目に  x_0 + \mbox{←} x_1 が奇数行目に対応しています.(n = 2 の場合はそれぞれ 1 個ずつだったのが,今は 2 個ずつ)
しかし,添字の番号が違っていますね. x_0 + x_1 に対応するのが  x_0+x_2, x_1+x_3 です.ここにどのような対応があるのでしょう?
答えは二進数にあります.添字を二進数で表してみると,(0,2) が (00,10), (1,3) が (01,11) です.二進数の一番上の桁に注目すると,どちらの組も (0,1) となっていることがわかります.これがルールです.
これが第 k 段階ならば上から k 桁だけに注目します.あとで詳しく考えます.
次に,
 \begin{align*}
\begin{pmatrix}
X_0 \\ X_2
\end{pmatrix} &= \begin{pmatrix}
x_0+x_2 + \mbox{→}(x_1+x_3) \\ x_0+x_2 + \mbox{←}(x_1+x_3)
\end{pmatrix} \\
\begin{pmatrix}
X_1 \\ X_3
\end{pmatrix} &= \begin{pmatrix}
x_0+\mbox{←}x_2 + \mbox{↓}(x_1+\mbox{←}x_3) \\ x_0+\mbox{←}x_2 + \mbox{↑}(x_1+\mbox{←}x_3)
\end{pmatrix}
\end{align*}

で完成です.(第 2 段階)
これをバタフライダイアグラムで表すと,

f:id:Ysmr_Ry:20171109093939p:plain

となります.左側の列の蝶の集まりが第 1 段階,右側が第 2 段階に対応しています.
第 1 段階は n = 2 の場合に対応すると言ったように,n = 2 の時のバタフライが 2 つ入っています.
サイズ n の FFT を実行することを, fft_n(x_i) と表すとすると, fft_4(x_0,x_2,x_1,x_3) = fft_2(x_0,x_2) + fft_2(x_1,x_3) + (\mbox{一番右のバタフライ}) です.これが漸化式です.

n = 8 の場合

これも同様です.バタフライダイアグラムだけ示すと,

f:id:Ysmr_Ry:20171109094547p:plain

となります.1 列目 (第 1 段階) には n = 2 の場合のバタフライが 4 つ,1, 2 列目を合わせると n = 4 の場合のバタフライが 2 つ隠れています.
この場合も, fft_8(x_0,x_4,x_2,x_6,x_1,x_5,x_3,x_7) = fft_4(x_0,x_4,x_2,x_6) + fft_4(x_1,x_5,x_3,x_7) + (\mbox{一番右のバタフライ}) が漸化式になります.

先ほど確認したように,問題を解くには漸化式の他にもっとも単純な問題の答えが必要です.
もっとも単純なのはこの場合 n = 1 の場合を考えると簡単です.n = 1 の場合は,フーリエ変換する前もなくそのままになります.
 fft_1(x_i) = x_i です.これで問題が解けます.

解けるんですが,あと問題なのは,バタフライダイアグラムの左側の  x_i の並びです.これにはある規則があります.

なぜビットリバースの順なのか

バタフライダイアグラムの左側の  x_i の右に書いてあるかっこの中の数は添字を二進数で表したものです.
これをそれぞれ逆から読んでみると,000, 001, 010, 011, 100, … で,0, 1, 2, 3, 4, … と順に増えていることがわかります.これをビットリバースの順と言います.
なぜこうなるのか考えてみましょう.漸化式に沿って考えてみると,右から 2 番目のバタフライが 2 つある段階 (第  \log_2 n-1 段階) において,添字を二進数で表した時の上から  \log_2 n-1 桁だけ見ています.言い換えると,一番右の桁だけ無視しています.なので,( \log_2 n-1 桁)0 と ( \log_2 n-1 桁)1 の 2 つのセットがあります.これが 2 つのバタフライに対応します.
漸化式から一番右のバタフライと 2 つの半分のサイズの  fft に言い換えられるため,半分のサイズのバタフライダイアグラムに今と同じことを順に適用していくと,n = 2 の場合,一番上の桁が 0, 1 なので,ビットリバースの順に増えていくことに対応します.
例えば,n = 4 の場合,第 1 段階は 0, 1 ですが,第 2 段階まで考えると,その後ろに 0 をつけるか,1 をつけるかの 2 セットあります.これを順に書くと,00, 10, 01, 11 でビットリバースの順になっています.
n = 8 ならば,00, 10, 01, 11 の後ろに 0 をつけるか 1 をつけるかで,000, 100, 010, 110, 001, 101, 011, 111 です.
ビットリバースの順に増えていくのは,そこまでの列の後ろに 0 を加えたものと 1 を加えたものをくっつけて拡張するということを繰り返していくと作れます.
ここにも再帰的定義が隠れているのです.
(伝わる?)

まとめ

個人的にネットに転がってる FFT の説明がわかりにくかったので自分なりにまとめてみました.
徐々に難しくしていったつもりなんですが,ビットリバースのところがうまく伝わったか心配です…
不明な点があったらコメントください.間違っていた場合もお願いします.