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

レーダーのアルゴリズムとは、MITのリンカーン研究所のチャールズ・M・レーダーにより考案された高速フーリエ変換のアルゴリズムである (Rader 1968)。このアルゴリズムではサイズが素数の離散フーリエ変換（DFT）を巡回畳み込みに置き換えることで計算コストを減らす（Bluestein のアルゴリズムもまたDFTを巡回畳み込みと置き換えるアルゴリズムである)。

レーダーのアルゴリズムはDFTカーネルの周期性のみを利用しているので、このアルゴリズムは同様の特徴をもつ（素数サイズの）数論的変換（英語版）や離散ハートレー変換（英語版）に対しても適用することができる。

このアルゴリズムはデータの並び変えを少し変更し、実数データに対して更なる高速化が可能である (Chu & Burrus 1982)。また、実数データについてに対する別の最適化として、Johnson & Frigo (2007)^{[要文献特定詳細情報]} によって示された離散ハートレイ変換を用いたアルゴリズムがある。

ウィノグラードはレーダーのアルゴリズムを拡張し、素数の冪乗のサイズのDFTに適用できるアルゴリズムを考案した (Winograd 1976; Winograd 1978）。このため、レーダーのアルゴリズムは、ウィノグラードのアルゴリズムの特別なケースとみなされる。ウィノグラードのアルゴリズムは乗算的FFTアルゴリズム (Tolimieri, An & Lu 1997) とも言われる。このアルゴリズムは素数の冪乗以外のサイズにも適用できるが、合成数のサイズのDFTについては、クーリー・ターキーのFFTアルゴリズムを用いるほうがより単純で実装も容易であるため、このアルゴリズムは通常、クーリー・ターキーのアルゴリズムの再帰分割により得られたDFTのうち、大きな素数サイズのDFTにのみ適用される (Frigo & Johnson 2005)。

アルゴリズム

X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\qquad k=0,\dots ,N-1.

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

X_{0}=\sum _{n=0}^{N-1}x_{n},

X_{g^{-p}}=x_{0}+\sum _{q=0}^{N-2}x_{g^{q}}e^{-{\frac {2\pi i}{N}}g^{-(p-q)}}\qquad p=0,\dots ,N-2.

上式は以下に示される長さ N–1(q =0,...,N–2) の数列 a_q と b_q の巡回畳み込み積分となっている。

a_{q}=x_{g^{q}}

b_{q}=e^{-{\frac {2\pi i}{N}}g^{-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₂ により N−1= 2N₂ と表され、また N₂−1 が素数 N₃ により N₂−1 = 2N₃ と表され、以下同様に続いていく場合である。このような場合、レーダーのアルゴリズムの再帰が連続することになり、O(N²) の計算時間がかかる可能性がある。このような性質をもつ N はソフィー・ジェルマン素数と呼ばれ、上記の数列は一次の Cunningham(ビル-カニンガム)チェーンと呼ばれる。しかしながら、これまでの研究ではカニンガムチェーンの成長は log₂(N) よりも遅いことが分かっているため、レーダーのアルゴリズムの再帰によりかかる計算時間は O(N²) よりかは速いと思われる。幸いにも、畳み込み計算にゼロ埋めを用いたFFTを使えば計算時間は O(N log N) のオーダーになることが保証されている。

プログラムの例

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

C言語による実装例

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
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){
		if(i & 1){flag=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];
	}
	
	for(i=0;i<len;i++){free(a[i]);}
	free(a);
	for(i=0;i<len;i++){free(b[i]);}
	free(b);
	for(i=0;i<len;i++){free(c[i]);}
	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];
	}
	
	for(i=0;i<n;i++){free(c[i]);}
	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;
	}
	
	for(i=0;i<n;i++){free(c[i]);}
	free(c);
}

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

参考文献

Rader, C. M. (1968年6月). “Discrete Fourier transforms when the number of data samples is prime”. Proceedings of the IEEE 56 (6): 1107–1108. doi:10.1109/PROC.1968.6477. ISSN 0018-9219.
Chu, Shuni; Burrus, C. (1982年4月). “A prime factor FTT algorithm using distributed arithmetic”. IEEE Transactions on Acoustics, Speech, and Signal Processing 30 (2): 217–227. doi:10.1109/TASSP.1982.1163875. ISSN 0096-3518.
Frigo, M.; Johnson, S. G. (2005年2月). “The Design and Implementation of FFTW3” (PDF). Proceedings of the IEEE 93 (2): 216–231. doi:10.1109/JPROC.2004.840301. ISSN 0018-9219. NAID 80017181219. http://fftw.org/fftw-paper-ieee.pdf.
Winograd, Shmuel (Apr 1976). “On computing the Discrete Fourier Transform”. Proc Natl Acad Sci USA 73 (4): 1005–1006. ISSN 0027-8424.
Winograd, S. (1978). “On Computing the Discrete Fourier Transform”. Mathematics of Computation 32 (141): 175–199. doi:10.2307/2006266. ISSN 00255718.
Tolimieri, Richard; An, Myoung; Lu, Chao (1997). Algorithms for Discrete Fourier Transform and Convolution (2nd ed.). Springer-Verlag. doi:10.1007/978-1-4757-2767-8. ISBN 978-0-387-98261-8. ISSN 1431-7893