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
Ａ ＝ Σ ａ Ｗ ， Ｗ ＝ ｅ
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
Ａ ＝ Σ Σ ａ Ｗ Ｗ
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
Ａ ＝ Σ Σ ａ Ｗ Ｗ
N1 t1 k2 + N2 t2 k1 j1=0 j2=0 N1 j2 + N2 j1 N1 N2
と得られます．これで単純な二次元 DFT になります．さらに，添字 j で同じ 変換を行うと
N1-1 N2-1 t2j1k1 t1j2k2
Ａ ＝ Σ Σ ａ Ｗ Ｗ
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 を構成する かで，いくつかのアルゴリズムに分けられます．