高速フーリエ変換 その１
高速フーリエ変換に関することを、自分の勉強及び今後の記事のために書きます。本記事（その1）では基礎的なことやアルゴリズムのことを、その2では群の表現論との関係を書こうと思います。
参考文献は記事最後の文献欄にあるRef.[1][2]です。
表記に関して2つほど:
$${\{ \ \}}$$は集合を表す
$${\vec x}$$は$${x_i}$$を並べて作ったベクトルを表す
離散フーリエ変換
離散フーリエ変換(DFT)は$${\{x\}}$$から$${\{y\}}$$への以下の変換です:
$$
\begin{aligned}
y_j&=\sum_{k=0}^{N-1}\exp(-2\pi i j k/N)x_k\\
&(i\text{は虚数単位})
\end{aligned}
\tag{1}
$$
これは例えば時間の関数を周波数成分に分解するような変換です。「音をドレミで表す」ような操作に対応します。
フーリエ変換は実に様々な場面で使われます。音声、音楽、そして画像の圧縮やフィルターに使われるのはよく知られていると思います。物理学だと、偏微分方程式の解の構成、特に境界値問題を解くのに重要です。また光の回折現象にもフーリエ変換は現れます。情報・通信技術や、医学だとMRI・CTスキャン・エコーなどの画像構成にも使われています。
以下、DFTを速く計算する方法である高速フーリエ変換（Fast Fourier Transform、FFT）に関して書きます。
要素数が8の場合の離散フーリエ変換
最初に、$${N=8}$$の場合のDFTをFFTにより具体的に計算することで、素朴にEq.(1)を計算するより計算量が減らせることを確認します。
Eq.(1)を書き下すと以下のようになります：
$$
\begin{aligned}
y_0= &x_0& + &x_1 &+ &x_2 &+ &x_3 &+ &x_4 &+ &x_5 &+ &x_6 &+ &x_7&\\
y_1 =& x_0& + &x_1W_8 &+ &x_2W_8^2 &+ &x_3W_8^3 &+ &x_4W_8^4 &+ &x_5W_8^5& + &x_6W_8^6 &+ &x_7W_8^7&\\
y_2=& x_0& + &x_1W_8^2 &+ &x_2W_8^4 &+ &x_3W_8^6 &+ &x_4 &+ &x_5W_8^2& + &x_6W_8^4 &+ &x_7W_8^6&\\
y_3=& x_0& + &x_1W_8^3 &+ &x_2W_8^6 &+ &x_3W_8 &+ &x_4W_8^4 &+ &x_5W_8^7& + &x_6W_8^2 &+ &x_7W_8^5&\\
y_4=& x_0& + &x_1W_8^4 &+ &x_2 &+ &x_3W_8^4 &+ &x_4 &+ &x_5W_8^4 &+ &x_6 &+ &x_7W_8^4&\\
y_5=& x_0& + &x_1W_8^5 &+ &x_2W_8^2 &+ &x_3W_8^7 &+ &x_4W_8^4 &+ &x_5W_8 &+ &x_6W_8^6 &+ &x_7W_8^3&\\
y_6=& x_0& + &x_1W_8^6 &+ &x_2W_8^4 &+ &x_3W_8^2 &+ &x_4 &+ &x_5W_8^6 &+ &x_6W_8^4 &+ &x_7W_8^2&\\
y_7=& x_0& + &x_1W_8^7&+ &x_2W_8^6 &+ &x_3W_8^5 &+ &x_4W_8^4 &+ &x_5W_8^3 &+ &x_6W_8^2 &+ &x_7W_8&
\end{aligned}
\\
({\rm y1})
$$
ここで$${W_8:=\exp(-2\pi i/8)}$$であり、また$${W_8^8=1}$$であることを用いて変形しています。
この式は大変きれいな形をしています。そしてよく見ると至るところに同じような計算が存在しています。ここで$${W_8^4=-W_8}$$を用いると
$$
x_0\pm x_4, x_1\pm x_5,x_2\pm x_6,x_3\pm x_7
$$
の組み合わせのみが登場していることがわかると思います。そこで、$${x}$$から作られる新たな変数として
$$
\begin{aligned}
x^{(1)}_0=&x_0&+&x_4,& \\
x^{(1)}_1=&x_1&+&x_5,& \\
x^{(1)}_2=&x_2&+&x_6,& \\
x^{(1)}_3=&x_3&+&x_7,& \\
x^{(1)}_4=&x_0&-&x_4,& \\
x^{(1)}_5=&x_1&-&x_5,& \\
x^{(1)}_6=&x_2&-&x_6,& \\
x^{(1)}_7=&x_3&-&x_7,&
\tag{x1}
\end{aligned}
$$
を定義します。するとEq.(y1)は以下のようになります：
$$
\begin{aligned}
y_0=& x_0^{(1)} &+& x_1^{(1)} &+& x_2^{(1)} &+& x_3^{(1)} ,&\\
y_1=& x_4^{(1)} &+& x_5^{(1)}W_8 &+& x_6^{(1)}W_8^2 &+& x_7^{(1)}W_8^3,&\\
y_2=& x_0^{(1)} &+& x_1^{(1)}W_8^2 &+& x_2^{(1)}W_8^4 &+& x_3^{(1)}W_8^6,&\\
y_3=& x_4^{(1)} &+& x_5^{(1)}W_8^3 &+& x_6^{(1)}W_8^6 &+& x_7^{(1)}W_8,&\\
y_4=& x_0^{(1)} &+& x_1^{(1)}W_8^4 &+& x_2^{(1)} &+& x_3^{(1)}W_8^4,&\\
y_5=& x_4^{(1)} &+& x_5^{(1)}W_8^5 &+& x_6^{(1)}W_8^2 &+& x_7^{(1)}W_8^7,&\\
y_6=& x_0^{(1)} &+& x_1^{(1)}W_8^6 &+& x_2^{(1)}W_8^4 &+& x_3^{(1)}W_8^2,&\\
y_7=& x_4^{(1)} &+& x_5^{(1)}W_8^7 &+& x_6^{(1)}W_8^6 &+& x_7^{(1)}W_8^5,&
\tag{y2}
\end{aligned}
$$
これをみると、Eq.(y1)と比較して右辺の和の数が半分になっています。
更に同様の変換を施すことができます。$${W_8^0=-W_8^4, \ W_8^2=-W_8^6}$$より、次の変数の組
$$
\begin{aligned}
x^{(2)}_0=&x_0^{(1)}&+&x_2^{(1)},&\\
x^{(2)}_1=&x_1^{(1)}&+&x_3^{(1)},&\\
x^{(2)}_2=&x_0^{(1)}&-&x_2^{(1)},&\\
x^{(2)}_3=&x_1^{(1)}&-&x_3^{(1)},&\\
x^{(2)}_4=&x_4^{(1)}&+&x_6^{(1)}W_8^2,&\\
x^{(2)}_5=&x_5^{(1)}&+&x_7^{(1)}W_8^2,&\\
x^{(2)}_6=&x_4^{(1)}&-&x_6^{(1)}W_8^2,&\\
x^{(2)}_7=&x_5^{(1)}&-&x_7^{(1)}W_8^2&
\tag{x2}
\end{aligned}
$$
を定義すると
$$
\begin{aligned}
y_0=& x_0^{(1)} &+& x_1^{(1)},&\\
y_1=& x_4^{(1)} &+& x_5^{(1)}W_8,&\\
y_2=& x_3^{(1)} &+& x_4^{(1)}W_8^2,&\\
y_3=& x_6^{(1)} &+& x_7^{(1)}W_8^3,&\\
y_4=& x_0^{(1)} &+& x_1^{(1)}W_8^4,&\\
y_5=& x_4^{(1)} &+& x_5^{(1)}W_8^5,&\\
y_6=& x_3^{(1)} &+& x_4^{(1)}W_8^6,&\\
y_7=& x_6^{(1)} &+& x_7^{(1)}W_8^7&
\tag{y3}
\end{aligned}
$$
となります。最後に
$$
\begin{aligned}
x^{(3)}_0=&x_0^{(2)}&+&x_1^{(2)},& \\
x^{(3)}_1=&x_0^{(2)}&-&x_1^{(2)},& \\
x^{(3)}_2=&x_2^{(2)}&+&x_3^{(2)}W_8^2,& \\
x^{(3)}_3=&x_2^{(2)}&-&x_3^{(2)}W_8^2,&\\
x^{(3)}_4=&x_4^{(2)}&+&x_5^{(2)}W_8^1,& \\
x^{(3)}_5=&x_4^{(2)}&-&x_5^{(2)}W_8^1,& \\
x^{(3)}_6=&x_6^{(2)}&+&x_7^{(2)}W_8^3,& \\
x^{(3)}_7=&x_6^{(2)}&-&x_7^{(2)}W_8^3&
\tag{x3}
\end{aligned}
$$
と変換すると
$$
\begin{aligned}
y_0=&x^{(3)}_0,&\\
y_1=&x^{(3)}_4,&\\
y_2=&x^{(3)}_2,&\\
y_3=&x^{(3)}_6,&\\
y_4=&x^{(3)}_1,&\\
y_5=&x^{(3)}_5,&\\
y_6=&x^{(3)}_3,&\\
y_7=&x^{(3)}_7&
\tag{y4}
\end{aligned}
$$
を得ます。これでDFTが計算できました。つまりEq.(x1-x3)の変数変換を施すことはDFTを実行することに対応しています。これがFFTです。
最後の表式Eq.(y4)を得るのにどれだけの計算量が必要か見積もります。Eq.(y4)を得るにはEq.(x1-x3)を計算すればよいです。これらは位相をかける操作を除けば、それぞれ8個の和が存在します。よって$${3\times 8=24}$$回の計算が必要です。
この議論を$${N=2^n}$$のときに一般化します。1回のEq.(x1-x3)のような変数変換で、$${\vec y}$$の表式の和はそれぞれ半分になります。$${n=\log_2 N}$$回の変数変換により$${\vec y}$$の表式の和はなくなり、Eq.(x3)のように最終的な変数そのものになります。どの変数変換でも、1回の和を含む$${N}$$個の変数の取り直しが必要だから、計算量は$${N\log_2 N}$$回程度になります。
元の表式だと、要素数$${N}$$のときには、縦$${N}$$個の足し算が$${N}$$個あるので、およそ$${N\times N = N^2}$$回の計算が必要です。よって計算量が本質的に違います。
Ref.[1]に、FFTにおける掛け算と足し算の回数を、素朴なDFTのそれと比較したテーブルが記載されているので引用しておきます（素朴なDFTでは$${{\mathcal O}(N^2)}$$、FFTでは$${{\mathcal O}(N\log_2N)}$$のアルゴリズムを使用）：
$$
\begin{array}{ccccc}
\hline
\hline
\text{data size } N & \text{DFTの積} & \text{DFTの和}& \text{FFTの積} &\text{FFTの和}\\
\hline
8 & 64& 56 & 12 & 24\\
16 & 256 & 240 & 32 & 64\\
32 & 1,024 & 992 & 80 & 160\\
64 & 4,096 & 4,032 & 192 & 384\\
\hline
\end{array}
\\
{}\\
(\text{from Ref.[1] P55 Table 3.2})
$$
$${N}$$が大きい時、FFTは素朴なDFTと比較して圧倒的に計算回数が少ないことがわかります。
疎行列による高速フーリエ変換の表現
前章で行った計算は、行列の積で以下のように表すことができます（Ref.[1]）。
$$
\vec y_{\rm BRO} = A_1 A_2 A_3\vec x
$$
ここで$${A_1, A_2, A_3}$$はそれぞれ以下の行列です：
$$
\begin{aligned}
A_1&=
\begin{aligned}
{\rm diag}
\left(
\begin{pmatrix}
1 & W_8^0\\
1 & -W_8^0
\end{pmatrix}
,
\begin{pmatrix}
1 & W_8^2\\
1 & -W_8^2
\end{pmatrix}
,
\begin{pmatrix}
1 & W_8^1\\
1 & -W_8^1
\end{pmatrix}
,\begin{pmatrix}
1 & W_8^3\\
1 & -W_8^3
\end{pmatrix}
\right)
\end{aligned}
,\\
A_2&=
\begin{aligned}
{\rm diag}
\left(
\begin{pmatrix}
I_2 & W_8^0 I_2\\
I_2 & -W_8^0 I_2
\end{pmatrix}
,
\begin{pmatrix}
I_2 & W_8^2 I_2\\
I_2 & -W_8^2 I_2
\end{pmatrix}
\right)
\end{aligned}
,\\
A_3&=
\begin{pmatrix}
I_4 & I_4\\
I_4 & -I_4
\end{pmatrix}
\end{aligned}
$$
$${I_n}$$は$${n\times n}$$の単位行列です。$${{\rm diag}}$$はカッコの中の要素を行列の対角部分に並べたものです。BROの添字は「ビットを逆に並べた順序」（Bit Reversed Order）の略です。整数$${b}$$がバイナリ（=2進数）表示で$${b=b_{n-1}2^{n-1}+b_{n-2}2^{n-2}+\ldots +b_12^1 +b_0 2^0}$$であるとき、これを$${b_{n-1}b_{n-2}\ldots b_1b_0}$$と表すことにします。これを逆順序に並べると$${b_0b_1\ldots b_{n-1}}$$になります。$${{b_{n-1}b_{n-2}\ldots b_1b_0}\rightarrow b_0b_1\ldots b_{n-1}}$$の対応をBROといいます。上記$${N=8}$$ならば以下のようになります：
$$
\begin{array}{c|c|c|c}
元の数の並び & 2進表示 & 2進表示の逆並び & {\rm BRO} \\
0 & 000 & 000 & 0\\
1 & 001 & 100 & 4\\
2 & 010 & 010 & 2\\
3 & 011 & 110 & 6\\
4 & 100 & 001 & 1\\
5 & 101 & 101 & 5\\
6 & 110 & 011 & 3\\
7 & 111 & 111 & 7\\
\end{array}
\\
\text{ (from Ref.[1] P46)}
$$
BROは前章のEq.(y4)の$${y}$$と$${x^{(3)}}$$の対応に現れています。
$${A_1, A_2, A_3}$$を順に$${\vec x}$$に作用させたものと、前章の変数変換Eq.(x1)(x2)(x3)がそれぞれ対応していることはすぐに確認できます。
これは一般の$${N}$$に拡張できて
$$
\vec y_{\rm BRO}={W^{nk}_N}_{\rm BRO}\vec x_{\rm BRO}
$$
のように$${\vec y}$$と$${\vec x}$$を結ぶ$${{W^{nk}_N}_{\rm BRO}}$$は
$$
{W^{nk}_N}_{\rm BRO}=A_1A_2\ldots A_{\log_2 N}
$$
と書けます(Ref.[1] (3.4)の上部参照)。ただし$${A_i \ (i=1,2,\ldots \log_2 N)}$$は単位行列、対角行列、またはブロック対角行列のような疎行列（多くの部分が0である行列）です。疎行列ではない行列を、何らかの変換を施すなどして疎行列の積に書き換えると、一般に計算量は減ります。
このような分解は、群の表現論的な視点から言うと、有限アーベル群上のフーリエ変換を剰余類の直和に分解することに対応します（次の記事で言及します）。表現行列を直和に分解すると、行列はブロック対角の形になり、そのため計算が減ります。
バタフライ演算
前章の$${N=8}$$のFFTを図にすると以下のように表現できます。
この図は一見難しそうですが、そうでもないです。図は左から右に見ていきます。最初$${x_i}$$が一番左の白丸にそれぞれ存在していると思ってください。これら$${x_i}$$がひとつ右の白丸に移るとき、以下の2つの操作が施されます：
線上に色付き太線があれば、その太線に対応したファクターをその線上の$${x}$$にかける。色付き太線とファクターの関係は図の下部分参照(オレンジ一本が$${-1}$$、青1本が$${W_8}$$に対応する。複数存在するときにはその本数分かける)。
他の場所からある線上の白丸に線がつながっている場合、他の場所から来た線の元にある$${x}$$を線上の$${x}$$に足す
例えば上から5番目の線$${x_4}$$の位置にある、左から2番目の白丸には
$$
-1\times x_4 + x_0= x_0-x_4
$$
が存在します。$${A_i}$$の文字の下の操作が、前章の$${A_i}$$の行列に対応します。最後の出力はBROであることに注意してください。
この図に表される演算はバタフライ演算と呼ばれます。
図から以下のことがわかります：
操作の塊が3つ$${=\log_2 N}$$個ある
$${A_2}$$から右側は上4つと下4つが分離している。$${ A_1}$$から右側は$${(x_0,x_1),(x_2,x_3),(x_4,x_5),(x_6,x_7)}$$ がそれぞれ分離している。
1.の3つの部分はEq.(x1-x3)の変数変換に対応しています。2.は$${N=8}$$のDFTが、4要素のDFT→2要素のDFTのように分解されることを示しています（次章で説明します）。
Radix-2における一般の高速フーリエ変換
ここまで主に$${N=8}$$の場合のFFTについて考えてきました。以下ではもっと一般の、要素数が$${2}$$のべき乗:$${N=2^n}$$の場合を考えます（本章のタイトルにあるRadix-2とは要素数が2のべき乗のことを表します）。
$${W_N:=\exp(-2\pi i/N)}$$とすると、$${W_N}$$は以下の性質を持ちます:
$$
\begin{aligned}
W_N^{N/2}&=-W_N, \ W_N^N=1,\\
W_N^m&=W_{N/m} \ \ \ (m\in {\mathbb N})
\end{aligned}
$$
これがFFTに重要です。
DFT:$${y_j=\sum_{k=0}^{N-1}x_k W_N^{jk}}$$において、$${k}$$の和をe偶数と奇数に分解します。$${k=2r}$$と$${k=2r+1}$$にわけ、それぞれの和を書き直すと以下のようになります：
$$
\begin{aligned}
y_j &= \sum_{k=0}^{N-1}W_N^{jk}x_k\\
&=\sum_{r=0}^{N/2-1}W_N^{j(2r)} x_{2r}
+\sum_{r=0}^{N/2-1}W_N^{j(2r+1)} x_{2r+1}\\
&=\sum_{r=0}^{N/2-1}W_{N/2}^{jr} x_{2r}
+W_N^j\sum_{r=0}^{N/2-1}W_{N/2}^{jr} x_{2r+1}
\end{aligned}
$$
これは、$${x_k}$$を偶奇にわけると、それぞれが独立なDFTに帰着することを意味します。そしてDFTの要素数は元のDFTの半分になります。
この作業は再帰的に行うことができます。すなわち、偶奇にわけたDFTそれぞれに存在する和をさらに偶奇にわければ、更に要素数が半分のDFTに帰着します。そして長さ1のDFTは恒等変換、すなわち$${x_0,y_0}$$に対して
$$
y_j=\sum_{k=0}^{1-1}\exp(2\pi j k/1)x_k=x_0 \ \ \ (j=0)
$$
なので、上記操作を繰り返し最終的に要素数1のDFTに辿りつけば作業は終わります。
FFTをフローチャートにすると以下のようになります：
図の$${\{x\}}$$の下に記したe,oの並びは、上で説明した偶(e)奇(o)の分解に対応するインデックスです。インデックスのe,oの並びを逆順にし、$${e\rightarrow 0, o\rightarrow 1}$$に対応させます。例えば$${{x}_{ooe}}$$なら$${eoo\rightarrow 011}$$とします。$${\{x\}_{oee}}$$は、$${x_i}$$の$${i}$$をバイナリ表示したとき、下位3ビットが$${011}$$となる2進数$${n}$$ケタのすべての$${i}$$の集合に対応する$${x_i}$$の集合です。一般には
$$
\begin{aligned}
&\hspace{1cm} \{x\}_{a_0a_1\ldots a_m}:=\{x_i\}_{i\in X}\\
&ただし\\
&\hspace{1cm} X=\{[j_{n-1}j_{n-2}\ldots j_{m+1}f(a_m)f(a_{m-1})\ldots f(a_{1})f(a_{0})] \ | \ j_{m+1},j_{m+2},\ldots,j_{n-2},j_{n-1}\in \{0,1\}\}\\
&ここで\\
& \hspace{1cm} \bullet [\hspace{0.2cm}]はバイナリ表示を表す。すなわち\\
& \hspace{1cm} \ \ \ [b_{n-1}b_{n-2}\ldots b_1b_0]:=b_{n-1}2^{n-1}+b_{n-2}2^{n-2} + \ldots + b_02^0 \ \ (b_i\in \{0,1\})\\
&\hspace{1cm} \bullet a_k \in \{e,o\} \ \ \ (k=0,1,\ldots,m-1)\\
&\hspace{1cm} \bullet f(a)=
\begin{cases}
0 & a=e\\
1 & a=o
\end{cases}
\end{aligned}
$$
で定義されます。$${ j_{m+1},j_{m+2},\ldots,j_{n-2},j_{n-1}\in \{0,1\}}$$に関しては、可能な0,1の組すべてに関して集合を作ります。図にある下線は、eoの並びと、$${x_{i}}$$の$${i}$$のビットの並びとの対応を示す補助的なものです。
Radix-2における一般の変数変換の構成
FFTの変数変換を、一般の$${N=2^n}$$において具体的に構成します。
改めて、サイズ$${N=2^n}$$の以下のDFTは以下で定義されます:
$$
y_j=\sum_{k=0}^{N-1}\exp(-2\pi i jk/N)x_{k} \ \ \ \ (j=0,1,\ldots,N-1)
$$
以下
$$
W_N:=\exp(-2\pi i/N)
$$
とします。また$${[\hspace{0.2cm}]}$$は前章と同様、数$${b}$$のバイナリ表示であり、
$$
\begin{aligned}
b&=b_{n-1}2^{n-1}+b_{n-2}2^{n-2}+\ldots+b_02^0
\rightarrow [b_{n-1}b_{n-2}\ldots b_0]\\
&(b_0,b_1,\ldots, b_{n-1}\in \{0,1\})
\end{aligned}
$$
です。このときFFTの変数変換は以下のように表されます(脚注)：
$$
\begin{aligned}
&x^{(0)}_{[k_0k_1\ldots k_{n-1}]}
=x_{[k_{n-1}k_{n-2}\ldots k_1k_0]}
\\
&\begin{cases}
x^{(1)}_{[k_0\ldots k_{n-2}|0]}
=x^{(0)}_{[k_0\ldots k_{n-2}0]}
+x^{(0)}_{[k_0\ldots k_{n-2}1]}
\\
x^{(1)}_{[k_0\ldots k_{n-2}|1]}
=x^{(0)}_{[k_0\ldots k_{n-2}0]}
-x^{(0)}_{[k_0\ldots k_{n-2}1]}
\end{cases}
\\
&\begin{cases}
x^{(2)}_{[k_0\ldots k_{n-3}|0j_0]}
=x^{(1)}_{[k_0\ldots k_{n-3}0|j_0]}
+W_N^{[j_0]N/2^2}
x^{(1)}_{[k_0\ldots k_{n-3}1|j_0]}
\\
x^{(2)}_{[k_0\ldots k_{n-3}|1j_0]}
= x^{(1)}_{[k_0\ldots k_{n-3}0|j_0]}
-W_N^{[j_0]N/2^2}
x^{(1)}_{[k_0\ldots k_{n-3}1|j_0]}
\end{cases}
\\
&\begin{cases}
x^{(3)}_{[k_0\ldots k_{n-4}|0j_1j_0]}
=x^{(2)}_{[k_0\ldots k_{n-4}0|j_1j_0]}
+W_N^{[j_1j_0]N/2^3}
x^{(2)}_{[k_0\ldots k_{n-4}1|j_1j_0]}
\\
x^{(3)}_{[k_0\ldots k_{n-4}|1j_1j_0]}
=x^{(2)}_{[k_0\ldots k_{n-4}0|j_1j_0]}
-W_N^{[j_1j_0]N/2^3}
x^{(2)}_{[k_0\ldots k_{n-4}1|j_1j_0]}
\end{cases}
\\
&\hspace{4cm}\vdots
\\
&\begin{cases}
x^{(n)}_{[0j_{n-2}\ldots j_1j_0]}
=x^{(n-1)}_{[0|j_{n-2}\ldots j_1j_0]}
+W_N^{[j_{n-2}\ldots j_1j_0]N/2^n}
x^{(n-1)}_{[1|j_{n-2}\ldots j_1j_0]}
\\
x^{(n)}_{[1j_{n-2}\ldots j_1j_0]}
=x^{(n-1)}_{[0|j_{n-2}\ldots j_1j_0]}
-W_N^{[j_{n-2}\ldots j_1j_0]N/2^n}
x^{(n-1)}_{[1|j_{n-2}\ldots j_1j_0]}
\end{cases}
\\
& y_{[j_{n-1}j_{n-2}\ldots j_1j_0]}=
x^{(n)}_{[j_{n-1}j_{n-2}\ldots j_1j_0]}
\end{aligned}
\tag{2}
$$
バイナリ表記に含まれる縦線$${|}$$は操作の区切りを表すものであり、２進数として読む際には無視してください。
上式Eq.(2)が$${N=2^1,2^2,2^3}$$で正しいことを確認するのは難しくないですが、表記の意味がわかりにくいかもしれないので、$${N=2^3}$$の場合の計算をAppedixに記しておきます。$${x_k^{(0)}}$$と$${x_k}$$の対応はBROであることに注意してください。
まとめ
本記事では高速フーリエ変換(FFT)を、$${N=8}$$の場合の具体例と共に説明しました。位相$${W_N:=\exp(-2\pi i/N)}$$がもつ対称性により、FTの長さを半分にする操作を再帰的に行うことができて、これがFFTに対応します。またFFTは疎行列による表現やバタフライ図による図示が可能であり、これらはFFTの理論的・直感的理解にそれぞれ役立ちます。最後の章「Radix-2における一般の変数変換の構成」では、一般の$${N=2^n}$$における具体的な変数変換の構成を示しました。
実際のコンピュータコードを作成するときには、和を偶奇にわける作業を再帰的に行うのが、一番楽にFFTを実装する方法かもしれません。
FFTはフーリエ変換の群論・表現論的側面と関係します。次の記事：「高速フーリエ変換のこと その２」では、これに関して説明する予定です。
おしまい。$${{}_\blacksquare}$$
(脚注) この表式はRef.[2]に記載されているものと基本的に同じですが、誤植と思われる部分を修正してあります。$${N=2^1,2^2,2^3}$$では式が正しいことを確かめていますが、一般に正しいかは保証しかねますので（たぶん大丈夫だとは思うのですが…）ご留意ください。
続きの記事：
References
[1] Rao, K. R., Kim, D. N. & Hwang, J.-J. "Fast Fourier Transform - Algorithms and Applications" (Springer Netherlands, 2010).
[2] 宮下精二「数値計算」応用数学基礎講座7, 朝倉書店, 2002.
Appendix: Eq.(2)を確認する
以下本文のEq.(2)が正しいことを$${N=2^3}$$の場合に確認します。$${\log_2 N=3}$$なので、$${x^{(1)},x^{(2)},x^{(3)}}$$の3つの変数変換（およびBROに対応する$${x^{(0)}}$$）があります。これを書き下すと以下のようになります：
$$
\begin{aligned}
&x^{(0)}_{[k_0k_1k_2]}=x_{[k_2k_1k_0]}\\
&\begin{cases}
x^{(1)}_{[k_0k_10]}&=x^{(0)}_{[k_0k_10]}+x^{(0)}_{[k_0k_11]},\\
x^{(1)}_{[k_0k_11]}&=x^{(0)}_{[k_0k_10]}-x^{(0)}_{[k_0k_11]}
\end{cases}
\\
&\begin{cases}
x^{(2)}_{[k_00j_0]}&=x^{(1)}_{[k_00j_0]}+W_8^{[j_0]}x^{(1)}_{[k_01j_0]},\\
x^{(2)}_{[k_01j_0]}&=x^{(1)}_{[k_00j_0]}-W_8^{[j_0]}x^{(1)}_{[k_01j_0]}
\end{cases}
\\
&\begin{cases}
x^{(3)}_{[0j_1j_0]}&=x^{(2)}_{[0j_1j_0]}+W_8^{[j_1j_0]}x^{(2)}_{[1j_1j_0]},\\
x^{(3)}_{[1j_1j_0]}&=x^{(2)}_{[0j_1j_0]}-W_8^{[j_1j_0]}x^{(2)}_{[1j_1j_0]}
\end{cases}
\end{aligned}
$$
ここでは操作の区切りを示す縦線$${|}$$は省略しました。$${j_i,k_i}$$に可能な$${0,1}$$の組をすべて代入し、それをバイナリ表示から10進表示に直します。これを$${x^{(3)}}$$と$${x^{(2)}}$$に関して行ってみます：
$$
\text{左辺の最上位ビットが0}\\
\begin{aligned}
j_0=0, j_1=0&: \ x^{(3)}_{[000]}=x^{(2)}_{[000]}+W_8^{[00]}x^{(2)}_{[100]}\leftrightarrow
x^{(3)}_{0}=x^{(2)}_{0}+W_8^{0}x^{(2)}_{4}\\
j_0=1, j_1=0&: \ x^{(3)}_{[001]}=x^{(2)}_{[001]}+W_8^{[01]}x^{(2)}_{[101]}\leftrightarrow
x^{(3)}_{1}=x^{(2)}_{1}+W_8^{1}x^{(2)}_{5}\\
j_0=0, j_1=1&: \ x^{(3)}_{[010]}=x^{(2)}_{[010]}+W_8^{[10]}x^{(2)}_{[110]}\leftrightarrow
x^{(3)}_{2}=x^{(2)}_{2}+W_8^{2}x^{(2)}_{6}\\
j_0=1, j_1=1&: \ x^{(3)}_{[011]}=x^{(2)}_{[011]}+W_8^{[11]}x^{(2)}_{[111]}\leftrightarrow
x^{(3)}_{3}=x^{(2)}_{3}+W_8^{3}x^{(2)}_{7}
\end{aligned}
\\{}
\\
\text{左辺の最上位ビットが1}\\
\begin{aligned}
j_0=0, j_1=0&: x^{(3)}_{[100]}=x^{(2)}_{[000]}-W_8^{[00]}x^{(2)}_{[100]}\leftrightarrow
x^{(3)}_{4}=x^{(2)}_{0}-W_8^{0}x^{(2)}_{4}\\
j_0=1, j_1=0&: x^{(3)}_{[101]}=x^{(2)}_{[001]}-W_8^{[01]}x^{(2)}_{[101]}\leftrightarrow
x^{(3)}_{5}=x^{(2)}_{1}-W_8^{1}x^{(2)}_{5}\\
j_0=0, j_1=1&: x^{(3)}_{[110]}=x^{(2)}_{[010]}-W_8^{[10]}x^{(2)}_{[110]}\leftrightarrow
x^{(3)}_{6}=x^{(2)}_{2}-W_8^{2}x^{(2)}_{6}\\
j_0=1, j_1=1&: x^{(3)}_{[111]}=x^{(2)}_{[011]}-W_8^{[11]}x^{(2)}_{[111]}\leftrightarrow
x^{(3)}_{7}=x^{(2)}_{3}-W_8^{3}x^{(2)}_{7}
\end{aligned}
$$
同様に$${x^{(2)}}$$と$${x^{(1)}}$$の関係は以下のようになります：
$$
\text{左辺の中位ビットが0}\\
\begin{aligned}
j_0=0, k_0=0&: \ x^{(2)}_{[000]}=
x^{(1)}_{[000]}
+W_8^{2[0]}x^{(1)}_{[010]}\leftrightarrow
x^{(2)}_{0}=
x^{(1)}_{0}
+W_8^{0}x^{(1)}_{2}\\
j_0=1, k_0=0&: \ x^{(2)}_{[001]}=
x^{(1)}_{[001]}
+W_8^{2[1]}x^{(1)}_{[011]}\leftrightarrow
x^{(2)}_{1}=x^{(1)}_{1}
+W_8^{2}x^{(1)}_{3}\\
j_0=0, k_0=1&: \ x^{(2)}_{[100]}=
x^{(1)}_{[100]}
+W_8^{2[0]}x^{(1)}_{[110]}\leftrightarrow
x^{(2)}_{4}=
x^{(1)}_{4}
+W_8^{0}x^{(1)}_{6}\\
j_0=1, k_0=1&: \ x^{(2)}_{[101]}=
x^{(1)}_{[101]}
+W_8^{2[1]}x^{(1)}_{[111]}\leftrightarrow
x^{(2)}_{5}=
x^{(1)}_{5}
+W_8^{2}x^{(1)}_{7}
\end{aligned}
\\{}
\\
\text{左辺の中位ビットが1}\\
\begin{aligned}
j_0=0, k_0=0&: x^{(2)}_{[010]}=
x^{(1)}_{[000]}
-W_8^{2[0]}x^{(1)}_{[010]}\leftrightarrow
x^{(2)}_{2}=
x^{(1)}_{0}
-W_8^{0}x^{(1)}_{2}\\
j_0=1, k_0=0&: x^{(2)}_{[011]}=
x^{(1)}_{[001]}
-W_8^{2[1]}x^{(1)}_{[011]}\leftrightarrow
x^{(2)}_{3}=
x^{(1)}_{1}
-W_8^{2}x^{(1)}_{3}\\
j_0=0, k_0=1&: x^{(2)}_{[110]}=
x^{(1)}_{[100]}-W_8^{2[0]}
x^{(1)}_{[110]}\leftrightarrow
x^{(2)}_{6}=
x^{(1)}_{4}
-W_8^{0}x^{(1)}_{6}\\
j_0=1, k_0=1&: x^{(2)}_{[111]}=
x^{(1)}_{[101]}
-W_8^{2[1]}
x^{(1)}_{[111]}\leftrightarrow
x^{(2)}_{7}=
x^{(1)}_{5}
-W_8^{2}x^{(1)}_{7}
\end{aligned}
$$
更に$${x^{(1)}}$$と$${x^{(0)}}$$の関係は以下のようになります：
$$
\text{左辺の最下位ビットが0}\\
\begin{aligned}
k_1=0, k_0=0&: \ x^{(1)}_{[000]}=
x^{(0)}_{[000]}
+x^{(0)}_{[001]}
\leftrightarrow
x^{(1)}_{0}=
x^{(0)}_{0}
+x^{(0)}_{1}\\
k_1=1, k_0=0&: \ x^{(1)}_{[010]}=
x^{(0)}_{[010]}
+x^{(0)}_{[011]}
\leftrightarrow
x^{(1)}_{2}=x^{(0)}_{2}
+x^{(0)}_{3}\\
k_1=0, k_0=1&: \ x^{(1)}_{[100]}=
x^{(0)}_{[100]}
+x^{(0)}_{[101]}
\leftrightarrow
x^{(1)}_{4}=
x^{(0)}_{4}
+x^{(0)}_{5}\\
k_1=1, k_0=1&: \ x^{(1)}_{[110]}=
x^{(0)}_{[110]}
+x^{(0)}_{[111]}
\leftrightarrow
x^{(1)}_{6}=
x^{(0)}_{6}
+x^{(0)}_{7}
\end{aligned}
\\{}
\\
\text{左辺の最下位ビットが1}\\
\begin{aligned}
k_1=0, k_0=0&: x^{(1)}_{[001]}=
x^{(0)}_{[000]}
-x^{(0)}_{[001]}
\leftrightarrow
x^{(1)}_{1}=
x^{(0)}_{0}
-x^{(0)}_{1}\\
k_1=1, k_0=0&: x^{(1)}_{[011]}=
x^{(0)}_{[010]}
-x^{(0)}_{[011]}
\leftrightarrow
x^{(1)}_{3}=
x^{(0)}_{2}
-x^{(0)}_{3}\\
k_1=0, k_0=1&: x^{(1)}_{[101]}=
x^{(0)}_{[100]}-
x^{(0)}_{[101]}
\leftrightarrow
x^{(1)}_{5}=
x^{(0)}_{4}
-x^{(1)}_{5}\\
k_1=1, k_0=1&: x^{(1)}_{[111]}=
x^{(0)}_{[110]}
-x^{(0)}_{[111]}
\leftrightarrow
x^{(1)}_{7}=
x^{(0)}_{6}
-x^{(0)}_{7}
\end{aligned}
$$
$${x^{(0)}_i}$$と$${x_i}$$の関係はBROであり（本文のBROの表参照）、$${y_i}$$と$${x^{(3)}_i}$$は足がそのまま対応します：$${y_i=x^{(3)}_i \ \ \ (i=0,1,\ldots 7)}$$。
これらよりFFTは以下のようになります：
$$
\begin{aligned}
&\begin{aligned}
y_0&=x_0^{(3)}\\
&=x_0^{(2)}+x_4^{(2)}\\
&=(x_0^{(1)}+x_2^{(1)})+(x_4^{(1)}+x_6^{(1)})\\
&=(x_0^{(0)}+x_1^{(0)})+(x_2^{(0)}+x_3^{(0)})
+(x_4^{(0)}+x_5^{(0)})+(x_6^{(0)}+x_7^{(0)})\\
&=(x_0+x_4)+(x_1+x_5)+(x_2+x_6)+(x_3+x_7)
\end{aligned}
\\
&\begin{aligned}
y_1&=x_1^{(3)}\\
&=x_1^{(2)}+W_8 x_5^{(2)}\\
&=(x_1^{(1)}+W_8^2x_3^{(1)})+W_8(x_5^{(1)}+W_8^2x_7^{(1)})\\
&=(x_0^{(0)}-x_1^{(0)})+W_8^2(x_2^{(0)}-x_3^{(0)})
+W_8(x_4^{(0)}-x_5^{(0)})+W_8^3(x_6^{(0)}-x_7^{(0)})\\
&=(x_0-x_4)+W_8^2(x_2-x_6)+W_8(x_1-x_5)+W_8^3(x_3-x_7)
\end{aligned}
\\
&\begin{aligned}
y_2&=x_2^{(3)}\\
&=x_2^{(2)}+W_8^2 x_6^{(2)}\\
&=(x_0^{(1)}-x_2^{(1)})+W_8^2(x_4^{(1)}-x_6^{(1)})\\
&=(x_0^{(0)}+x_1^{(0)})-(x_2^{(0)}+x_3^{(0)})
+W_8^2(x_4^{(0)}+x_5^{(0)})-W_8^2(x_6^{(0)}+x_7^{(0)})\\
&=(x_0+x_4)-(x_2+x_6)
+W_8^2(x_1+x_5)-W_8^2(x_3+x_7)
\end{aligned}
\\
&\begin{aligned}
y_3&=x_3^{(3)}\\
&=x_3^{(2)}+W_8^3 x_7^{(2)}\\
&=(x_1^{(1)}-W_8^2x_3^{(1)})+W_8^2(x_5^{(1)}-W_8^2x_7^{(1)})\\
&=(x_0^{(0)}-x_1^{(0)})-W_8^2(x_2^{(0)}-x_3^{(0)})
+W_8^3(x_4^{(0)}-x_5^{(0)})-W_8^5(x_6^{(0)}-x_7^{(0)})\\
&=(x_0-x_4)-W_8^2(x_2-x_6)
+W_8^3(x_1-x_5)-W_8^5(x_3-x_7)
\end{aligned}
\\
&\begin{aligned}
y_4&=x_4^{(3)}\\
&=x_0^{(2)}-x_4^{(2)}\\
&=(x_1^{(1)}+x_2^{(1)})-(x_4^{(1)}+x_6^{(1)})\\
&=(x_0^{(0)}+x_1^{(0)})+(x_2^{(0)}+x_3^{(0)})
-(x_4^{(0)}+x_5^{(0)})-(x_6^{(0)}+x_7^{(0)})\\
&=(x_0+x_4)+(x_2+x_6)
-(x_1+x_5)-(x_3+x_7)
\end{aligned}
\\
&\begin{aligned}
y_5&=x_5^{(3)}\\
&=x_1^{(2)}-W_8x_5^{(2)}\\
&=(x_1^{(1)}+W_8^2x_3^{(1)})-W_8(x_5^{(1)}+W_8^2x_7^{(1)})\\
&=(x_0^{(0)}-x_1^{(0)})+W_8^2(x_2^{(0)}-x_3^{(0)})
-W_8(x_4^{(0)}-x_5^{(0)})-W_8^3(x_6^{(0)}-x_7^{(0)})\\
&=(x_0-x_4)+W_8^2(x_2-x_6)
-W_8(x_1-x_5)-W_8^3(x_3-x_7)
\end{aligned}
\\
&\begin{aligned}
y_6&=x_6^{(3)}\\
&=x_2^{(2)}-W_8^2x_6^{(2)}\\
&=(x_0^{(1)}-x_2^{(1)})-W_8^2(x_4^{(1)}-x_6^{(1)})\\
&=(x_0^{(0)}+x_1^{(0)})-(x_2^{(0)}+x_3^{(0)})
-W_8^2(x_4^{(0)}+x_5^{(0)})+W_8^2(x_6^{(0)}+x_7^{(0)})\\
&=(x_0+x_4)-(x_2+x_6)
-W_8^2(x_1+x_5)+W_8^2(x_3+x_7)
\end{aligned}
\\
&\begin{aligned}
y_7&=x_7^{(3)}\\
&=x_3^{(2)}-W_8^3x_7^{(2)}\\
&=(x_1^{(1)}-W_8^2x_3^{(1)})-W_8^3(x_5^{(1)}-W_8^2x_7^{(1)})\\
&=(x_0^{(0)}-x_1^{(0)})-W_8^2(x_2^{(0)}-x_3^{(0)})
-W_8^3(x_4^{(0)}-x_5^{(0)})+W_8^5(x_6^{(0)}-x_7^{(0)})\\
&=(x_0-x_4)-W_8^2(x_2-x_6)
-W_8^3(x_1-x_5)+W_8^5(x_3-x_7)
\end{aligned}
\end{aligned}
$$
これは本文のEq.(y1)を正しく再現しています。$${{}_\blacksquare}$$