Prime Factor 型 FFT の技法は Cooley-Tukey の FFT より以前から 知られていた方法です.しかし,この方法が注目され始めたのは, Cooley-Tukey の FFT から十年以上も後になってからです.この方法は, Rader による短い素数の長さの DFT の効率的な計算方法と組み合わせたときに 実用的な効果を発揮します.この Prime Factor 型 FFT の基本的な考え方は, Cooley-Tukey 型 FFT の基本的な考え方とほぼ同じで,N = N1*N2 の長さの DFT をネストされた N1 の長さと N2 の長さの DFT に分解することで 行います.ただ異なるのは,その分解の方法です.Prime Factor 型 FFT では, 分解する N1 と N2 は任意には選ばず,互いに素(最大公約数が 1)となるように 選びます.このとき,Cooley-Tukey 型 FFT で現れる回転因子の乗算が 除去できて,乗算回数を大幅に削減することが可能となります.
Prime Factor 型 FFT での分解はまず,長さ N = N1*N2 の DFT
N-1 jk -2πi / N A = Σ a W , W = e k j=0 j N N
の添字 j, k を,N1 を法とする添字 j1, k1 と,N2 を法とする添字 j2, k2 に 変換することで行います[参考文献]. この変換は,N1 と N2 が互いに素な数であれば,N より小さい数 j, k は,
j = N1 j2 + N2 j1 mod N k = N1 k2 + N2 k1 mod N
と一意に分解できるという事実を用いて行います.例えば,N = 15 で N1 = 3, N2 = 5 のときは次のようになります.
リスト1.3.1-1. j = 3*j2+5*j1 mod 15 の分解 |
---|
j | 0 3 6 9 12 5 8 11 14 2 10 13 1 4 7 -----+------------------------------------------------------------- j1 | 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 j2 | 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 |
因みに,j, k から j1, j2, k1, k2 の求め方は,まず N1 と N2 に対する ユークリッドの互除法などにより,
1 = N1 t1 + N2 t2 mod N
を満たす整数 t1, t2 を求め,これを j 倍または K 倍することで求めます. この添字の変換でもとの DFT は,N1 と N2 の長さの DFT
N1-1 N2-1 N2 j1 k1 N1 j2 k2 A = Σ Σ a W W N1 k2 + N2 k1 j1=0 j2=0 N1 j2 + N2 j1 N1 N2
に分解できます.これは,指数部が N2, N1 倍だけ回転した長さ N1 × N2 の 二次元 DFT に相当します.通常の二次元 DFT への分解は,例えば添字 k を
k = N1 t1 k2 + N2 t2 k1 mod N
で変換し直すことで
N1-1 N2-1 j1 k1 j2 k2 A = Σ Σ a W W N1 t1 k2 + N2 t2 k1 j1=0 j2=0 N1 j2 + N2 j1 N1 N2
と得られます.これで単純な二次元 DFT になります.さらに,添字 j で同じ 変換を行うと
N1-1 N2-1 t2j1k1 t1j2k2 A = Σ Σ a W W N1t1k2 + N2t2k1 j1=0 j2=0 N1t1j2 + N2t2j1 N1 N2
となります.この最後の分解は,添字が変換前と変換後で連続に変化するので, 実際の応用でよく用いられます.Prime Factor 型 FFT は,これらの分解を 繰り返し行い,長さ N の DFT を互いに素な長さの多次元 DFT に変換する ことで構成します.この分解は,Cooley-Tukey 型 FFT の分解と比較すると 回転因子が含まれておらず,分解に要する乗算回数はゼロになります.その 代償として,N1 と N2 を互いに素になるように選ばなければなりません.
次に,Cooley-Tukey 型 FFT の例と同様に,最初の分解式を素直に C に 直した例を示します.
リスト1.3.1-2. 再帰呼び出しによる PFA |
---|
/* F[k]=Σ_j=0^n-1 a[j]*exp(±2*pi*i*j*k/n),0<=k<n を計算する. n はデータ数で任意, theta は ±2*PI/n, ar[0...n-1], ai[0...n-1], はデータの実部, 虚部, tmpr[0...n-1], tmpi[0...n-1] は作業領域. */ #include <math.h> void fft(int n, double theta, double ar[], double ai[], double tmpr[], double tmpi[]) { int n1, n2, j, j1, j2, k, k1, k2; double xr, xi, wr, wi; if (n <= 1) return; /* ---- factorization ---- */ for (n1 = 2; n1 * n1 <= n; n1++) { if (n % n1 == 0) break; } if (n % n1 != 0) n1 = n; for (n2 = n; n2 % n1 == 0; n2 /= n1); n1 = n / n2; /* ---- short DFTs ---- */ for (j2 = 0; j2 < n2; j2++) { for (k1 = 0; k1 < n1; k1++) { xr = ar[n1 * j2]; xi = ai[n1 * j2]; for (j1 = 1; j1 < n1; j1++) { j = (n1 * j2 + n2 * j1) % n; wr = cos(theta * n2 * j1 * k1); wi = sin(theta * n2 * j1 * k1); xr += wr * ar[j] - wi * ai[j]; xi += wr * ai[j] + wi * ar[j]; } tmpr[n2 * k1 + j2] = xr; tmpi[n2 * k1 + j2] = xi; } } for (j = 0; j < n; j += n2) { fft(n2, theta * n1, &tmpr[j], &tmpi[j], ar, ai); } /* ---- unscrambler ---- */ for (k1 = 0; k1 < n1; k1++) { j1 = n2 * k1 % n1; for (k2 = 0; k2 < n2; k2++) { j2 = n1 * k2 % n2; k = (n2 * k1 + n1 * k2) % n; ar[k] = tmpr[n2 * j1 + j2]; ai[k] = tmpi[n2 * j1 + j2]; } } } |
これは,Prime Factor Algorithm (PFA) ルーチンの一例で,多次元の DFT は,単純な一次元の DFT の直積で 計算されています.このルーチンは,考え方を示すもので,最適化は 全くされていません.実用的な PFA ルーチンでは,再帰呼び出しは展開し, 除算や剰余の計算はループから取り除きます.さらに,素因数に分解された DFT の計算は N^2 算法でなく個別に設計された最適な Short DFT アルゴリズムを使います.具体的な設計方法は後に示します.
Prime Factor 型 FFT の基本的な考え方は,添字の変換により,大きな 一次元 DFT を小さな互いに素な長さの多次元 DFT に分解することです. 通常,Prime Factor 型 FFT の長さ N は,このような短い長さの因数に 分解できるように選ばれます.なぜならば,短い長さの DFT は,演算回数が 最小になるように設計することが可能だからです.このような最適な Short DFT アルゴリズムは,データフローの構造や,乗算を重視するか 加算を重視するかなどでいくつか存在します.さらに,Prime Factor 型 FFT は,この Short DFT をどのように組み合わせて多次元 DFT を構成する かで,いくつかのアルゴリズムに分けられます.