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

レーダーのアルゴリズムとは、MITのリンカーン研究所のチャールズ・M・レーダーにより考案された高速フーリエ変換のアルゴリズムである 。このアルゴリズムではサイズが素数の離散フーリエ変換（DFT）を巡回畳み込みに置き換えることで計算コストを減らす（ のアルゴリズムもまたDFTを巡回畳み込みと置き換えるアルゴリズムである)。
レーダーのアルゴリズムはDFTカーネルの周期性
のみを利用しているので、このアルゴリズムは同様の特徴をもつ（素数サイズの） やに対しても適用することができる。
このアルゴリズムはデータの並び変えを少し変更し、実数データに対して更なる高速化が可能である 。また、実数データについてに対する別の最適化として、 によって示された離散ハートレイ変換を用いたアルゴリズムがある。
ウィノグラードはレーダーのアルゴリズムを拡張し、素数の冪乗のサイズのDFTに適用できるアルゴリズムを考案した (; ）。このため、レーダーのアルゴリズムは、ウィノグラードのアルゴリズムの特別なケースとみなされる。ウィノグラードのアルゴリズムは乗算的FFTアルゴリズム とも言われる。このアルゴリズムは素数の冪乗以外のサイズにも適用できるが、合成数のサイズのDFTについては、クーリー・ターキーのFFTアルゴリズムを用いるほうがより単純で実装も容易であるため、このアルゴリズムは通常、クーリー・ターキーのアルゴリズムの再帰分割により得られたDFTのうち、大きな素数サイズのDFTにのみ適用される 。
アルゴリズム.
"N" が素数の場合、"n" =1...,"N"–1の添字はNを法とする乗算について群をなす。数論によるとこのような群については生成元（原始根）が存在し、これを "gとすると、すべての n は" q = 0...,N–2 を用いて "n" = "g""q" (mod "N") と表される。(qとnの間に全単射が存在する。) n と k の添字をこの置き換えを用いて新しい添字 q と p により表すと、定式は以下のように置き換えられる。
上式は以下に示される長さ N–1(q =0...,N–2) の数列 "a""q" と "b""q" の巡回畳み込み積分となっている。
畳み込み積分の計算.
 畳み込み定理 により、FFTを用いて計算することができる。また、"N"−1 は合成数のため、より高速なクーリーターキー型のFFTアルゴリズムが適用できる。しかしながら、"N"−1 が大きな素数の素因数をもつ場合、レーダーのアルゴリズムを再帰的に使う必要がある。その代わりに、長さ "N"−1 のDFTをゼロ埋めによって2の冪乗にすることで、高速なクーリーターキーのアルゴリズムを用いることができる。ゼロ埋め後のデータサイズが最低 2("N"−1)−1 以上あれば元のデータを完全に復元することができ、ゼロ埋めによる結果の変化は生じない。
畳み込み積分の計算にゼロ埋めを用いたFFTを用いることでのこのアルゴリズムは加算について O("N") の計算時間とそれに加えと畳み込み積分についてと O("N" log "N") の計算時間で実行できる。しかしながら、このアルゴリズムは付近の合成数のサイズの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"2) の計算時間がかかる可能性がある。このような性質をもつ "N" は ソフィー・ジェルマン素数と呼ばれ、上記の数列は一次の Cunningham(ビル-カニンガム)チェーンと呼ばれる。しかしながら、これまでの研究ではカニンガムチェーンの成長は log2("N") よりも遅いことが分かっているため、レーダーのアルゴリズムの再帰によりかかる計算時間は O("N"2) よりかは速いと思われる。幸いにも、畳み込み計算にゼロ埋めを用いたFFTを使えば計算時間は O("N" log "N") のオーダーになることが保証されている。
プログラムの例.
以下はC言語による実装例である。
int powerMod(int a, int b, int m);
int primePrimitiveRoot(int p);
void DFT(double **x, int n);
void IDFT(double **x, int n);
void FFT_prime(double **x, int n){
 int i,len,flag,pad;
 int g,ig;
 double X[2] = {0,0};
 double **b,**a,**c;
 /*ゼロ埋めのサイズの計算*/
 flag = 0;
 len=1;
 for(i=n-1;i>1;i»=1){
 len «= 1;
 
 if(flag){
 len «= 2; /*2(n-1) - 1 以上の2の冪乗*/
 }else{
 len = n-1; /*n-1が2の冪乗の場合ゼロ埋めはしない*/
 pad = len - (n-1); /*ゼロ埋めの数*/
 g = primePrimitiveRoot(n); /*nの原始根*/
 ig = powerMod(g,n-2,n); /*nの原始根の逆数*/
 /*数列 a_q の作成*/
 a = malloc(len*sizeof(double*));
 for(i=0;i<len;i++){
 a[i] = malloc(2*sizeof(double));
 /*1番目と2番目の要素の間をゼロ埋めする。*/
 if(i == 0){
 a[0][0] = x[1][0];
 a[0][1] = x[1][1];
 }else if(i <= pad){
 a[i][0] = 0;
 a[i][1] = 0;
 }else{
 a[i][0] = x[powerMod(g,i-pad,n)][0];
 a[i][1] = x[powerMod(g,i-pad,n)][1];
 /*数列 b_q の作成*/
 b = malloc(len*sizeof(double*));
 for(i=0;i<len;i++){
 b[i] = malloc(2*sizeof(double));
 b[i][0] = cos(2*M_PI*powerMod(ig,i,n)/n);
 b[i][1] = -sin(2*M_PI*powerMod(ig,i,n)/n);
 /*離散フーリエ変換*/
 DFT(a,len);
 DFT(b,len);
 /*離散フーリエ変換されたaとbの要素ごとの積を計算し、配列cに格納する。*/
 c = malloc(len*sizeof(double*));
 for(i=0;i<len;i++){
 c[i] = malloc(2*sizeof(double));
 /*複素数の積の計算*/
 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(c,len);
 /*x0をXg^-qに足す*/
 for(i=0;i<len;i++){
 c[i][0] += x[0][0];
 c[i][1] += x[0][1];
 /*X0の計算*/
 for(i=0;i<n;i++){
 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(i=0;i<n-1;i++){
 x[powerMod(ig,i,n)][0] = c[i][0];
 x[powerMod(ig,i,n)][1] = c[i][1];
 
 free(a);
 free(b);
 free(c);
/*冪剰余を求める関数*/
int powerMod(int a, int b, int m){
 int i;
 for(i=1;b;b»=1){
 if(b & 1){
 i = (i*a)%m;
 a = (a*a)%m;
 return i;
/*素数の最小の原始根を見つける関数*/
int primePrimitiveRoot(int p){
 int i,j,size=0;
 int *list;
 list = malloc((p-1)*sizeof(int));
 /*p-1の素因数を見つけ、listに格納する。*/
 for(i=2,j=p-1;i*i<=j;i++){
 if(j % i == 0){
 list[size] = i;
 size++;
 do{
 j /= i;
 }while(j%i==0);
 for(i=0;i<size;i++){
 list[i] = (p-1)/list[i];
 for(i=2;i<=p-1;i++){
 for(j=0;j<size;j++){
 if(powerMod(i,list[j],p)==1){
 goto next;
 break;
 next: ;
 free(list);
 return i;
/*離散フーリエ変換の関数*/
void DFT(double **x, int n){
 int i,j;
 double **c;
 c = malloc(n*sizeof(double*));
 for(i=0;i<n;i++){
 c[i] = malloc(sizeof(double));
 c[i][0] = 0;
 c[i][1] = 0;
 
 for(i=0;i<n;i++){
 for(j=0;j<n;j++){
 c[i][0] += x[j][0]*cos(2*M_PI*i*j/n) + x[j][1]*sin(2*M_PI*i*j/n);
 c[i][1] += -x[j][0]*sin(2*M_PI*i*j/n) + x[j][1]*cos(2*M_PI*i*j/n);
 for(i=0;i<n;i++){
 x[i][0] = c[i][0];
 x[i][1] = c[i][1];
 
 free(c);
/*逆離散フーリエ変換の関数*/
void IDFT(double **x, int n){
 int i,j;
 double **c;
 c = malloc(n*sizeof(double*));
 for(i=0;i<n;i++){
 c[i] = malloc(sizeof(double));
 c[i][0] = 0;
 c[i][1] = 0;
 
 for(i=0;i<n;i++){
 for(j=0;j<n;j++){
 c[i][0] += x[j][0]*cos(2*M_PI*i*j/n) - x[j][1]*sin(2*M_PI*i*j/n);
 c[i][1] += x[j][0]*sin(2*M_PI*i*j/n) + x[j][1]*cos(2*M_PI*i*j/n);
 for(i=0;i<n;i++){
 x[i][0] = c[i][0]/n;
 x[i][1] = c[i][1]/n;
 
 free(c);