レーダーのFFTアルゴリズム

レーダーのアルゴリズムとは、MITのリンカーン研究所のチャールズ・M・レーダーにより考案された高速フーリエ変換のアルゴリズムである Rader。このアルゴリズムではサイズが素数の離散フーリエ変換を巡回畳み込みに置き換えることで計算コストを減らす。このため、レーダーのアルゴリズムは、ウィノグラードのアルゴリズムの特別なケースとみなされる。ウィノグラードのアルゴリズムは乗算的FFTアルゴリズム Tolimieri とも言われる。このアルゴリズムは素数の冪乗以外のサイズにも適用できるが、合成数のサイズのDFTについては、クーリー・ターキーのFFTアルゴリズムを用いるほうがより単純で実装も容易であるため、このアルゴリズムは通常、クーリー・ターキーのアルゴリズムの再帰分割により得られたDFTのうち、大きな素数サイズのDFTにのみ適用される Frigo。

1. アルゴリズム
---------

//N//  が素数の場合、//n// =1,...,//N//–1の添字はNを法とする乗算について群をなす。数論によるとこのような群については生成元が存在し、これを  //gとすると、すべての n は//  q = 0,...,N–2 を用いて  //n// = //g////q//  と表される。 n と k の添字をこの置き換えを用いて新しい添字 q と p により表すと、定式は以下のように置き換えられる。

上式は以下に示される長さ N–1 の数列  //a////q//  と  //b////q//  の巡回畳み込み積分となっている。

1.1. 畳み込み積分の計算
--------------

  畳み込み定理  により、FFTを用いて計算することができる。また、//N//−1 は合成数のため、より高速なクーリーターキー型のFFTアルゴリズムが適用できる。しかしながら、//N//−1 が大きな素数の素因数をもつ場合、レーダーのアルゴリズムを再帰的に使う必要がある。その代わりに、長さ //N//−1 のDFTをゼロ埋めによって2の冪乗にすることで、高速なクーリーターキーのアルゴリズムを用いることができる。ゼロ埋め後のデータサイズが最低 2−1 以上あれば元のデータを完全に復元することができ、ゼロ埋めによる結果の変化は生じない。

畳み込み積分の計算にゼロ埋めを用いたFFTを用いることでのこのアルゴリズムは加算について O の計算時間とそれに加えと畳み込み積分についてと O の計算時間で実行できる。しかしながら、このアルゴリズムは付近の合成数のサイズのFFTに比べ、およそ3から10倍の時間がかかる。

畳み込み積分の計算をゼロ埋めなしのFFTで行う場合、計算時間は //N// の性質に強く依存する。最悪の場合、//N//−1 が素数  //N//2  により N−1= 2//N//2  と表され、また  //N//2−1 が素数 //N//3  により //N//2−1 = 2//N//3 と表され、以下同様に続いていく場合である。このような場合、レーダーのアルゴリズムの再帰が連続することになり、O の計算時間がかかる可能性がある。このような性質をもつ //N// は  ソフィー・ジェルマン素数と呼ばれ、上記の数列は一次の Cunninghamチェーンと呼ばれる。しかしながら、これまでの研究ではカニンガムチェーンの成長は log2 よりも遅いことが分かっているため、レーダーのアルゴリズムの再帰によりかかる計算時間は O よりかは速いと思われる。幸いにも、畳み込み計算にゼロ埋めを用いたFFTを使えば計算時間は O のオーダーになることが保証されている。

1.1.1. プログラムの例
--------------

以下はC言語による実装例である。

include
include
include int powerMod; int primePrimitiveRoot; void DFT; void IDFT;

void FFT_prime{ int i,len,flag,pad; int g,ig; double X[2] = {0,0}; double b,a,c;

/*ゼロ埋めのサイズの計算*/ flag = 0; len=1; for{ if{flag=1;} len <<= 1; }

if{ len <<= 2; /*2 - 1 以上の2の冪乗*/ }else{ len = n-1; /*n-1が2の冪乗の場合ゼロ埋めはしない*/ } pad = len - ; /*ゼロ埋めの数*/ g = primePrimitiveRoot; /*nの原始根*/ ig = powerMod; /*nの原始根の逆数*/

/*数列 a_q の作成*/ a = malloc); for{ a[i] = malloc); /*1番目と2番目の要素の間をゼロ埋めする。*/ if{ a[0][0] = x[1][0]; a[0][1] = x[1][1]; }else if{ a[i][0] = 0; a[i][1] = 0; }else{ a[i][0] = x[powerMod][0]; a[i][1] = x[powerMod][1]; } }

/*数列 b_q の作成*/ b = malloc); for{ b[i] = malloc); b[i][0] = cos/n); b[i][1] = -sin/n); }

/*離散フーリエ変換*/ DFT; DFT;

/*離散フーリエ変換されたaとbの要素ごとの積を計算し、配列cに格納する。*/ c = malloc); for{ c[i] = malloc); /*複素数の積の計算*/ c[i][0] = a[i][0]*b[i][0] - a[i][1]*b[i][1]; c[i][1] = a[i][0]*b[i][1] + a[i][1]*b[i][0]; }

/*逆離散フーリエ変換*/ IDFT;

/*x0をXg^-qに足す*/ for{ c[i][0] += x[0][0]; c[i][1] += x[0][1]; }

/*X0の計算*/ for{ X[0] += x[i][0]; X[1] += x[i][1]; }

/*X0,Xg^-q,を元の配列xに格納する。*/ x[0][0] = X[0]; x[0][1] = X[1];

/*添字の並び変え*/ for{ x[powerMod][0] = c[i][0]; x[powerMod][1] = c[i][1]; }

for{free;} free; for{free;} free; for{free;} free; }

/*冪剰余を求める関数*/ int powerMod{ int i; for{ if{ i = %m; } a = %m; } return i; }

/*素数の最小の原始根を見つける関数*/ int primePrimitiveRoot{ int i,j,size=0; int *list;

list = malloc*sizeof); /*p-1の素因数を見つけ、listに格納する。*/ for{ if{ list[size] = i; size++; do{ j /= i; }while; } } for{ list[i] = /list[i]; } for{ for{ if==1){ goto next; } } break; next: ; } free; return i; }

/*離散フーリエ変換の関数*/ void DFT{ int i,j; double c;

c = malloc); for{ c[i] = malloc); c[i][0] = 0; c[i][1] = 0; }

for{ for{ c[i][0] += x[j][0]*cos + x[j][1]*sin; c[i][1] += -x[j][0]*sin + x[j][1]*cos; } }

for{ x[i][0] = c[i][0]; x[i][1] = c[i][1]; }

for{free;} free; }

/*逆離散フーリエ変換の関数*/ void IDFT{ int i,j; double c;

c = malloc); for{ c[i] = malloc); c[i][0] = 0; c[i][1] = 0; }

for{ for{ c[i][0] += x[j][0]*cos - x[j][1]*sin; c[i][1] += x[j][0]*sin + x[j][1]*cos; } }

for{ x[i][0] = c[i][0]/n; x[i][1] = c[i][1]/n; }

for{free;} free; }<WtXmlEndTag />

複素数データを表現するため、一列目を実数部、二列目を虚数部とした二次元配列を用いている。
レーダーのアルゴリズムではgの逆数が生成する数列も必要であるため、フェルマーの小定理からgの逆数を求めている。
畳み込み積分を不変に保つゼロ埋めの手法として、aについて先頭要素とその次の要素の間にゼロを埋め、bについてはgの生成する数列がn-1の周期をもつことを自然に用いて元の数列を巡回させている。ゼロ埋めの手法はこの方法以外にも考えられる。
畳込み積分の計算に通常の離散フーリエ変換を用いているが、これを高速フーリエ変換のルーチンで置き換えることができる。
原始根を求めるアルゴリズムは原始根を参照。

2. 参考文献
-------

journal=Proc Natl Acad Sci USA