レーダーのFFTアルゴリズム
アルゴリズム > 畳み込み積分の計算 > プログラムの例

以下は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のFFT algoritmの視覚的表現。色つきの時計はDFT行列の位相を示す。11の原始根である2が生成する数列(2,4,8,5,10,9,7,3,6,1)を用いて行と列を先頭を除いて入れ替えることで、DFT行列は巡回行列となる。巡回行列をデータ列にかけることは巡回畳み込みと同値である。