C言語によるMulti-RadixのFFTプログラム

FFTのCooley-Turkeyアルゴリズムでは、データを分割してそれらのフーリエ変換の組み合わせとして元のフーリエ変換を表すことで、計算量を削減する。データの分割の単位をRadixというが、広く知られているFFTはRadix-2のFFTである。この場合、データ数は2の冪乗に限られる。
Cooley-Turkeyアルゴリズムは実際にはどんなRadixでも動く。この場合データ数は2の冪乗に限られない。ただし、サイズが素数の場合は通常のDFTになる。もちろんRadix-2は素数データ数は扱えないので、これはMulti-Radixだけの問題というわけではない。

ここではc言語によるMulti-Radixのプログラムを紹介する。

Cooley-Turkey型アルゴリズムの原理

この図はサイズ18のフーリエ変換の図である。行は入力データに対応し、列は出力データに対応する。この図の右にいくにつれ時計は回転しており、下に行くほど時計の回転が早くなっている。この時計に対応する位相をもつ三角関数をデータの各点にかけ、それを横に進みながら足していけば離散フーリエ変換が得られる。
別の言い方をすれば、元データをベクトルとみなし、この図に対応するような1のN乗根からなる行列を掛けると離散フーリエ変換をしたことになる。

Cooley-Turkeyのアルゴリズムは、この行列の列をうまく並び替えると元の行列が小さなサイズのフーリエ変換の行列に分解できることを利用している。

この図の最初の並び替えで現れた4つのパターンのうち、左側はサイズ9のフーリエ変換のパターンとまったく同じである。しかし、右側はそうなってはいない。しかし、よく注意してみると、右側のパターンは左側のパターンの時計を下に行くにつれ徐々にずらしていくことで作れることが分かる。これが回転因子をかけることに対応する。

Multi-Radix Cooley-Turkeyアルゴリズムのステップ

• データサイズ$N$を$N_1\times N_2$に分解する。
• データ列を並べなおす。（$N_1$の剰余が同じグループに分ける）
• 各グループ毎にそれぞれ離散フーリエ変換を施す。
• 上で得られた結果を回転因子をかけながら足し合わせる。

Multi-Radix Cooley-Turkeyアルゴリズムのc言語による実装
実装には色々工夫が必要である。
再帰を単純にしたいなら、小さなサイズの離散フーリエ変換に分けるたびに、新しくそのサイズの配列を生成すればいいが、それは効率が悪いため、in-placeで元の配列を更新していくのが望ましい。しかし、そのためには元の配列の分割された領域を操作する関数を再帰させないといけないがこの関数の設計は骨が折れた。
Radixが2や4のときはhard-codeにしたり、
SinやCosが0や1になるときはピッタリそうなるようにもした。

Data.h

typedef struct st_data{
 double **elem;
 int row;
 int column;
}Data_Sub;
typedef struct st_data *Data;
Data Data_create(int row, int column);
void Data_delete(Data data);

Data.c

Data Data_create(int row, int column){
 Data data = malloc(sizeof(Data_Sub));
 data->elem = malloc(row*sizeof(double*));
 data->elem[0] = malloc(row*column*sizeof(double));
 for(i=1;i<row;i++){data->elem[i] = data->elem[0] + i*column;}
 data->row = row;
 data->column = column;
 return data;
}
void Data_delete(Data data){
 if(data){
  free(data->elem[0]);
  free(data->elem);
  free(data);
 }
}

Data.h Data.c は動的な二次元配列のクラスである。

Data_FFT.h

# include <Data.h>
void Data_FFT(Data data);

Data_FFT.c

# include <Data_FFT.h>
# define Data_FFT_prime_table_size 100
int Data_FFT_prime_table[Data_FFT_prime_table_size]=
{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
211,223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277,
281, 283,293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359,
367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439,
443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521,
523, 541
};
double Data_FFT_cos(int n, int m){
	if(n == 0){return 1;}
	if(n == m){/*2 Pi*/return 1;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return -1;}
		if(m/n == 4){/*Pi/2*/return 0;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return 0;}
	}
	return cos(2*M_PI*n/m);
}
double Data_FFT_sin(int n, int m){
	if(n == 0){return 0;}
	if(n == m){return 0;}
	if(m%n == 0){
		if(m/n == 2){/*Pi*/return 0;}
		if(m/n == 4){/*Pi/2*/return 1;}
	}
	if(m%(m-n)==0){
		if(m/(m-n)==4){/*3 Pi/2*/return -1;}
	}
	return sin(2*M_PI*n/m);
}
void Data_FFT_mod(Data data, Data temp, int n, int sign){
	int i,j,k,l,m,ip,s,t;
	m = data->row/n;
	/*hard-code for size 1,2,4 DFT*/
	switch(m){
	  case 1:
		return;
	  case 2:
		for(i=0;i<n;i++){
			temp->elem[i*2][0] = data->elem[i*2][0] + data->elem[i*2+1][0];
			temp->elem[i*2+1][0] = data->elem[i*2][0] - data->elem[i*2+1][0];
			temp->elem[i*2][1] = data->elem[i*2][1] + data->elem[i*2+1][1];
			temp->elem[i*2+1][1] = data->elem[i*2][1] - data->elem[i*2+1][1];
		}
		for(i=0;i<n;i++){
			data->elem[i*2][0] = temp->elem[i*2][0];
			data->elem[i*2][1] = temp->elem[i*2][1];
			data->elem[i*2+1][0] = temp->elem[i*2+1][0];
			data->elem[i*2+1][1] = temp->elem[i*2+1][1];
		}
		return;
	  case 4:
		for(i=0;i<n;i++){
			temp->elem[i*4][0] = data->elem[i*4][0] + data->elem[i*4+1][0] + data->elem[i*4+2][0] + data->elem[i*4+3][0];
			temp->elem[i*4+1][0] = data->elem[i*4][0] - sign*data->elem[i*4+1][1] - data->elem[i*4+2][0] + sign*data->elem[i*4+3][1];
			temp->elem[i*4+2][0] = data->elem[i*4][0] - data->elem[i*4+1][0] + data->elem[i*4+2][0] - data->elem[i*4+3][0];
			temp->elem[i*4+3][0] = data->elem[i*4][0] + sign*data->elem[i*4+1][1] - data->elem[i*4+2][0] - sign*data->elem[i*4+3][1];
			temp->elem[i*4][1] = data->elem[i*4][1] + data->elem[i*4+1][1] + data->elem[i*4+2][1] + data->elem[i*4+3][1];
			temp->elem[i*4+1][1] = data->elem[i*4][1] + sign*data->elem[i*4+1][0] - data->elem[i*4+2][1] - sign*data->elem[i*4+3][0];
			temp->elem[i*4+2][1] = data->elem[i*4][1] - data->elem[i*4+1][1] + data->elem[i*4+2][1] - data->elem[i*4+3][1];
			temp->elem[i*4+3][1] = data->elem[i*4][1] - sign*data->elem[i*4+1][0] - data->elem[i*4+2][1] + sign*data->elem[i*4+3][0];
		}
		for(i=0;i<n;i++){
			data->elem[i*4][0]   = temp->elem[i*4][0];
			data->elem[i*4+1][0] = temp->elem[i*4+1][0];
			data->elem[i*4+2][0] = temp->elem[i*4+2][0];
			data->elem[i*4+3][0] = temp->elem[i*4+3][0];
			data->elem[i*4][1]   = temp->elem[i*4][1];
			data->elem[i*4+1][1] = temp->elem[i*4+1][1];
			data->elem[i*4+2][1] = temp->elem[i*4+2][1];
			data->elem[i*4+3][1] = temp->elem[i*4+3][1];
		}
		return;
	}
	/*Try factorization using prime table*/
	for(ip=0;ip<Data_FFT_prime_table_size;ip++){
		if(m % Data_FFT_prime_table[ip] == 0){
			s = Data_FFT_prime_table[ip];
			t = m/s;
			/*permute*/
			for(i=0;i<n;i++){
				for(j=0;j<s;j++){
					for(k=0;k<t;k++){
						temp->elem[i*m + j*t + k][0] = data->elem[i*m + s*k + j][0];
						temp->elem[i*m + j*t + k][1] = data->elem[i*m + s*k + j][1];
					}
				}
			}
			for(i=0;i<data->row;i++){
				data->elem[i][0] = temp->elem[i][0];
				data->elem[i][1] = temp->elem[i][1];
			}
            /*recursion*/
			Data_FFT_mod(data,temp,n*s,sign);
			for(i=0;i<data->row;i++){
				temp->elem[i][0] = temp->elem[i][1] = 0;
			}
			/*weave in twiddle factors*/
			for(i=0;i<n;i++){
				for(j=0;j<s;j++){
					for(l=0;l<t;l++){
						for(k=0;k<s;k++){
							temp->elem[i*m + j*t + l][0] 
							+= data->elem[i*m + k*t + l][0]*temp->elem[(k*(j*t+l)%m)*n][2] 
							- sign*data->elem[i*m + k*t + l][1]*temp->elem[(k*(j*t+l)%m)*n][3];
							temp->elem[i*m + j*t + l][1] 
							+= sign*data->elem[i*m + k*t + l][0]*temp->elem[(k*(j*t+l)%m)*n][3] 
							+ data->elem[i*m + k*t + l][1]*temp->elem[(k*(j*t+l)%m)*n][2];
						}
					} 
				}
			}
			for(i=0;i<data->row;i++){
				data->elem[i][0] = temp->elem[i][0];
				data->elem[i][1] = temp->elem[i][1];
			}
			return;
		}
	}
	/*Normal DFT*/
	for(i=0;i<data->row;i++){
		temp->elem[i][0] = temp->elem[i][1] = 0;
	}
	for(i=0;i<n;i++){
		for(j=0;j<m;j++){
			for(k=0;k<m;k++){
				temp->elem[i*m+j][0] += data->elem[i*m+k][0]*temp->elem[j*k*n%m][2] 
				- sign*data->elem[i*m+k][1]*temp->elem[j*k*n%m][3];
				temp->elem[i*m+j][1] += sign*data->elem[i*m+k][0]*temp->elem[j*k*n%m][3] 
				+ data->elem[i*m+k][1]*temp->elem[j*k*n%m][2];
			}
		}
	}
	for(i=0;i<data->row;i++){
		data->elem[i][0] = temp->elem[i][0];
		data->elem[i][1] = temp->elem[i][1];
	}
	return;
}
void Data_FFT(Data data){
	Data temp;
	int i;
	if(data->column!=2){
		fprintf(stdout,"Data_FFT : data column length must be 2\n");
		return;
	}
	temp = Data_create(data->row,4);
	for(i=0;i<data->row;i++){
		temp->elem[i][2] = Data_FFT_cos(i,data->row);
		temp->elem[i][3] = Data_FFT_sin(i,data->row);
	}
	Data_FFT_mod(data,temp,1,-1);
	Data_delete(temp);
}

Data_FFT_mod(Data data,Data temp,int n, int sign)が再帰の関数である。
dataは入出力データで実数部と虚数部の2列のデータからなり、関数の実行後、中身がFFTの結果に置き換わる。
tempという二次元配列は4列からなり、行数は入力データ数と同じである。
最初の2列は作業用で行列演算の結果を一時的に記憶するために使い、
3列目と4列目には予め使う三角関数の値を入れておく。
n は現在の分割数を表している。再帰はn=1から始まって、再帰のたびに素数がかかっていく。
signは指数関数の符号でフーリエ変換のとき-1、逆フーリエ変換のときに1を入れる。