N = p (素数)の場合の FFT

自分で使うため，そして勉強のため，FFT を計算するコードを書いていたが，長さ $N=p$ (素数)の場合に高速化する手法(Rader の FFT)を知らなかった．いろいろと調べた結果，何とか理解することができたのでまとめておこうと思う．

準備

$x(n)\;\;n=0,\cdots,N-1$ のスペクトルを $X(k)\;\;k=0,\cdots,N-1$ とする．つまり
$\;\;\;\;X(k)=\sum_{n=0}^{N-1} x(n) W_N^{kn}\;\;{\rm where}\;W_N=\exp \left( -\frac{2 \pi j}{N} \right)$

行列表現では

$\left[ \begin{array} X(0) \\ X(1) \\ X(2) \\ \vdots \\ X(N-1) \end{array} \right] = \left[ \begin{array} 1 & 1 & 1 & \cdots & 1 \\ 1 & W_N & W_N^2 & \cdots & W_N^{N-1} \\ 1 & W_N^2 & W_N^4 & \cdots & W_N^{N-2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & W_N^{N-1} & W_N^{N-2} & \cdots & W_N \end{array} \right] \left[ \begin{array} x(0) \\ x(1) \\ x(2) \\ \vdots \\ x(N-1) \end{array} \right]$

変換行列

簡単のため，変換行列の $W_N$ の0以外の指数のみを取りだした行列 $F$ を考えてみる．
$F=\left[ \begin{array} 1 & 2 & 3 & \cdots & N-1 \\ 2 & 4 & 6 & \cdots & N-2 \\ 3 & 6 & 9 & \cdots & N-3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ N-1 & N-2 & N-3 & \cdots & 1 \end{array} \right]$
この行列の $i$ 行 $j$ 列の要素 $f_{ij}$ は
$\;\;\;\;f_{ij} = (i \times j)\;{\rm mod}\;N$
となる．modulo-N であるのは，回転子 $W_N$ の性質による．

ガロア体

ガロア体についての詳しい説明は避けるが，ガロア体 $GF(p)$ は $p$ 個( $N$ =p (素数))の要素を持ち，四則演算について閉じているものであると考えてよい．今， $a$ を $GF(p)$ の原始根の1つであるとする．つまり， $GF(p)$ 上で $a^1,\cdots,a^{p-1}$ を計算すると， $a^{p-1}$ のみが1となるような要素である．

例えば， $p=5$ の場合， $GF(5)$ の原始根は2と3の2つが考えられる．
$\;\;\;\;2^0=1,\;2^1=2,\;2^2=4,\;2^3=3,\;2^4=1$
$\;\;\;\;3^0=1,\;3^1=3,\;3^2=4,\;3^3=2,\;3^4=1$

$a^m$ の逆元

$a$ が原始根ならば， $a^{p-1}=1$ が成り立っているから， $a^m$ の逆元 $a^{-m}$ は
$\;\;\;\;a^m a^{-m} = 1 = a^{p-1}$
より
$\;\;\;\;a^{-m} = a^{p-1-m}$

行列の変形

行列 $F$ を行については原始根 $a$ の数列 $a^0,a^1,\cdots,a^{N-2}$ ，列についてはそれらの逆元の数列 $a^0,a^{N-2},\cdots,a$ の順番に入れ替え，その行列を $F'$ とする．
$F'=\left[ \begin{array} 1 & a^{N-2} & a^{N-3} & \cdots & a \\ a & 1 & a^{N-2} & \cdots & a^2 \\ a^2 & a & 1 & \cdots & a^3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a^{N-2} & a^{N-3} & a^{N-4} & \cdots & 1 \end{array} \right]$
行列 $F'$ の要素は対角方向に同じものが並んでいる．

FFTの計算式の変形

以上のような変形を最初のFFTの計算式(行列表現)に持ち込むと，

$\left[ \begin{array} X(0) \\ X(1) \\ X(a) \\ X(a^2) \\ \vdots \\ X(a^{N-2}) \end{array} \right] = \left[ \begin{array} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & W_N & W_N^{a^{N-2}} & W_N^{a^{N-3}} & \cdots & W_N^a \\ 1 & W_N^a & W_N & W_N^{a^{N-2}} & \cdots & W_N^{a^2} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & W_N^{a^{N-2}} & W_N^{a^{N-3}} & W_N^{a^{N-4}} & \cdots & W_N \end{array} \right] \left[ \begin{array} x(0) \\ x(1) \\ x(a^{N-2}) \\ x(a^{N-3}) \\ \vdots \\ x(a) \end{array} \right]$

のようになる． $X(0)$ は単純に和を取ればよいだけなので，それ以外について考える．以下のような置き換えを考え，
$\;\;\;X(a^i)=x(0)+X'(i)\;\;\;{\rm where}\;i=0,\cdots,N-2$
$\;\;\;x(a^i)=x'(i)\;\;\;{\rm where}\;i=0,\cdots,N-2$
$\;\;\;W_N^{a^{-i}}=W(i)\;\;\;{\rm where}\;i=0,\cdots,N-2$
これ用いてを行列を数式に分解すれば，
$\begin{array} X'(0) &=& W(0) & x'(0)+ & W(1) & x'(N-2)+ & W(2) & x'(N-3)+ & \cdots & + W(N-2) & x'(1) \\ X'(1) &=& W(0) & x'(1)+ & W(1) & x'(0)+ & W(2) & x'(N-2)+ & \cdots & + W(N-2) & x'(2) \\ X'(2) &=& W(0) & x'(2)+ & W(1) & x'(1)+ & W(2) & x'(0)+ & \cdots & + W(N-2) & x'(3) \\ & \vdots & & & & & & & & & & \\ X'(N-2) &=& W(0) & x'(N-2)+ & W(1) & x'(N-3)+ & W(2) & x'(N-4) & \cdots & + W(N-2) & x'(0) \end{array}$
となる．これは $\{X'(i)\}$ が $\{W(i)\}$ と $\{x'(i)\}$ の巡回畳み込みで計算出来ることを意味している．この巡回畳み込みの計算には長さ $N-1$ の FFT および IFFT が使用できる．ここがこの手法のミソである．

まとめ

長さ $N=p$ (素数)の場合のFFTの計算方法

$X(0)=x(0)+x(1)+\cdots+x(N-1)$
ガロア体 $GF(p)$ の原始根 $a$ を1つ見つける．
を用いて，データの並べ替えを行う．
- $x'(i)=x(a^i)\;\;{\rm where}\;i=0,\cdots,N-2$
- $W(i)=W_N^{a^{-i}}\;\;{\rm where}\;i=0,\cdots,N-2$
長さ $N-1$ の FFT および IFFT で $\{W(i)\}$ と $\{x'(i)\}$ の巡回畳み込みを計算し，その結果を $\{X'(i)\}$ とする．
$X(a^i)=x(0)+X'(i)\;\;{\rm where}\;i=0,\cdots,N-2$
$X(0),\cdots,X(N-1)$ を得る．