高速フーリエ変換（FFT）ライブラリの紹介

２次元高速フーリエ変換（FFT）を実行するライブラリと利用方法を解説します。 本項では、ガウス関数と呼ばれる次の関数のフーリエ変換による展開係数の計算と、展開係数から元の関数へ逆変換を実行する方法を解説します。

手順１：関数の準備

関数fと項数Nと実空間の区間Lを用意します。

//項数と空間距離
var N = Math.pow(2, 7); //128
var L = 10;
//元の関数
function f( x, y ){
	var sigma = L/10;
	x = x - L/2;
	y = y - L/2;
	return Math.exp( -(x*x+y*y)/(sigma*sigma)  );
}

手順２：高速フーリエ変換関数を実行

FFT2D_f関数の第一引数と第二引数に、フーリエ変換対象関数の実数部と虚数部、第三引数に区間の長さ、第四引数に項数（空間分割数）を与えます。 戻り値のresults.anとresults.bnに展開係数の実数部と虚数部が格納されます。結果を取得後グラフ描画するだけです。

//展開係数の計算
var results = FFT2D_f( f, null, L, N, false );

手順３：高速逆フーリエ変換関数を実行

IFFT2D関数の第一引数と第二引数に、展開係数の実部と虚部をそれぞれ２次元配列で渡します。 第三引数と第四引数に逆変換後の結果を受け取る２次元配列を渡します。

//展開係数の取得
var an = results.an;
var bn = results.bn;
//結果受取用の配列の準備
var fr = [];
var fi = [];
for( var i=0; i<N; i++ ){
	fr[i] = []; 
	fi[i] = [];
}
//展開係数の計算
IFFT2D( an, bn, fr, fi );

２次元フーリエ変換の結果

実際に計算した結果をグラフ描画（グラフ描画方法）します。

高速フーリエ変換（FFT）ライブラリ

１次元高速フーリエ変換

function FFT( an, bn, N, Inverse ){
	/////////////////////////////////////////
	//参考：Javaで学ぶシミュレーションの基礎
	/////////////////////////////////////////
	// 入力
	// N  ： 項数（2のべき乗）
	// an : 実数配列（順変換：実数空間データを項数Nで指定、逆変換：展開係数a(n)）
	// bn : 実数配列（順変換：虚数空間データを項数Nで指定、逆変換：展開係数b(n)）
	// Inverse : 逆変換の場合に true
	/////////////////////////////////////////
	// 出力
	// an : 実数配列（順変換：展開係数a(n)、逆変換：実数空間データ）
	// bn : 実数配列（順変換：展開係数b(n)、逆変換：虚数空間データ）
	/////////////////////////////////////////
	var ff = Inverse ? 1 : -1;
	var rot = new Array(N);
	for( var i = 0; i < rot.length; i++ ) rot[ i ] = 0;
	var nhalf = N/2, num = N/2;
	var sc = 2 * Math.PI / N;
	while( num >= 1 ){
		for(var j = 0; j < N; j += 2 * num ){
			var phi = rot[j] / 2;
			var phi0 = phi + nhalf;
			var c = Math.cos( sc * phi );
			var s = Math.sin( sc * phi * ff );
			for( var k = j; k < j + num; k++ ){
				var k1 = k + num;
				var a0 = an[ k1 ] * c - bn[ k1 ] *s;
				var b0 = an[ k1 ] * s + bn[ k1 ] *c;
				an[ k1 ] = an[ k ] - a0;
				bn[ k1 ] = bn[ k ] - b0;
				an[ k ] = an[ k ] + a0;
				bn[ k ] = bn[ k ] + b0;
				rot[ k ] = phi;
				rot[ k1 ] = phi0;
			}
		}
		num = num / 2;
	}
	for( var i = 0; i < N ; i++ ){
		var j = rot[ i ]; 
		if( j > i ){
			var tmp = an[ i ];
			an[ i ] = an[ j ];
			an[ j ] = tmp;
			var tmp = bn[ i ];
			bn[ i ] = bn[ j ];
			bn[ j ] = tmp;	
		}
	}
	for( var i = 0; i < N ; i++ ){
		an[ i ] = an[ i ] / Math.sqrt(N);
		bn[ i ] = bn[ i ] / Math.sqrt(N);
	}
}

１次元高速フーリエ変換（関数形指定型）

function FFT_f( Fr, Fi, L, N ){
	/////////////////////////////////////////
	// 入力
	// Fr : フーリエ変換を実行する関数の実部（function(x)）
	// Fi : フーリエ変換を実行する関数の虚部（function(x)）
	// L  : 実空間の区間 0からL             （float）
	// N  : 項数（2のべき乗）               （int）
	/////////////////////////////////////////
	// 出力 return results;
	// results.an : フーリエ変換後の展開係数の実部
	// results.bn : フーリエ変換後の展開係数の虚部
	// results.replace.an : フーリエ変換後の展開係数の実部（n=0を配列要素の中心に）
	// results.replace.bn : フーリエ変換後の展開係数の虚部（n=0を配列要素の中心に）
	/////////////////////////////////////////
	if( Fr == undefined ) Fr = function(){ return 0; };
	if( Fi == undefined ) Fi = function(){ return 0; };

	//展開係数
	var fr = [];
	var fi = [];
	var an = [];
	var bn = [];

	for( var i=0; i<N; i++ ){
		var x = L/N * i;
		fr[ i ] = Fr( x );
		fi[ i ] = Fi( x );
	}

	//フーリへ変換の実行
	FFT( fr, fi, N, false);

	//結果を格納するオブジェクト
	var results = {};
	results.an = [];
	results.bn = [];
	results.replace = {};
	results.replace.an = [];
	results.replace.bn = [];
	for (var i = 0; i < N; i++) {
		results.an[ i ] = fr[ i ];
		results.bn[ i ] = fi[ i ];
		var ii;
		if( i < N/2 ) ii = N/2 + i;
		else ii = i - N/2; 
		results.replace.an[ i ] = fr[ ii ];
		results.replace.bn[ i ] = fi[ ii ];
	}

	return results;
}

２次元高速フーリエ変換（関数形指定型）

function FFT2D_f( Fr, Fi, L, N ){
	/////////////////////////////////////////
	// 入力
	// Fr : フーリエ変換を実行する関数の実部（function(x,y)）
	// Fi : フーリエ変換を実行する関数の虚部（function(x,y)）
	// L  : 実空間の区間 0からL             （float）
	// N  : 項数（2のべき乗）               （int）
	/////////////////////////////////////////
	// 出力 return results;
	// results.an : フーリエ変換後の展開係数の実部
	// results.bn : フーリエ変換後の展開係数の虚部
	// results.replace.an : フーリエ変換後の展開係数の実部（n=0を配列要素の中心に）
	// results.replace.bn : フーリエ変換後の展開係数の虚部（n=0を配列要素の中心に）
	/////////////////////////////////////////
	
	if( Fr == undefined ) Fr = function(){ return 0; };
	if( Fi == undefined ) Fi = function(){ return 0; };
	if( Fr == undefined ) Fr = function(){ return 0; };
	if( Fi == undefined ) Fi = function(){ return 0; };

	//展開係数
	var fr = [];
	var fi = [];
	var an = [];
	var bn = [];
	for( var i=0; i<N; i++ ){
		fr[i] = []; 
		fi[i] = [];
		an[i] = []; 
		bn[i] = [];
	}
	for( var i=0; i<N; i++ ){
		for( var j=0; j<N; j++ ){
			var x = L/N * i;
			var y = L/N * j;
			fr[ i ][ j ] = Fr( x, y );
			fi[ i ][ j ] = Fi( x, y );
		}
	}

	//実空間格子上値から展開係数を計算
	FFT2D( fr, fi, an, bn );

	//結果
	var results = {};
	results.an = [];
	results.bn = [];
	results.replace = {};
	results.replace.an = [];
	results.replace.bn = [];
	for (var i = 0; i < N; i++) {
		results.an[ i ] = [];
		results.bn[ i ] = [];
		results.replace.an[ i ] = [];
		results.replace.bn[ i ] = [];

		for (var j = 0; j < N; j++) {
			results.an[ i ][ j ] = an[ i ][ j ];
			results.bn[ i ][ j ] = bn[ i ][ j ];

			var ii = (i < N/2)? N/2 + i : i - N/2;
			var jj = (j < N/2)? N/2 + j : j - N/2;
			results.replace.an[ i ][ j ] = an[ ii ][ jj ];
			results.replace.bn[ i ][ j ] = bn[ ii ][ jj ];
		}
	}
	return results;
}

２次元高速フーリエ変換（実空間格子上値から展開係数を計算）

function FFT2D( fr, fi, an, bn ){
	/////////////////////////////////////////
	// 入力
	// fr : フーリエ変換を実行する関数の実部（[][]）
	// fi : フーリエ変換を実行する関数の虚部（[][]）
	/////////////////////////////////////////
	// 出力 
	// an : フーリエ変換後の展開係数の実部（[][]）
	// bn : フーリエ変換後の展開係数の虚部（[][]）
	/////////////////////////////////////////

	var N = fr.length;
	for( var i=0; i<N; i++ ){
		//x軸に平行なFFTの実行
		FFT( fr[ i ], fi[ i ], N );
	}
	for( var i=0; i<N; i++ ){
		for( var j=0; j<N; j++ ){
			an[ i ][ j ] = fr[ j ][ i ];
			bn[ i ][ j ] = fi[ j ][ i ];
		}
	}	
	for( var i=0; i<N; i++ ){
		//y軸に平行なFFTの実行
		FFT( an[ i ], bn[ i ], N);
	}

}

２次元高速逆フーリエ変換（展開係数から実空間格子上値を計算）

function IFFT2D( an, bn, fr, fi ){
	/////////////////////////////////////////
	// 入力
	// an : 逆フーリエ変換対象の展開係数の実部（[][]）
	// bn : 逆フーリエ変換対象の展開係数の虚部（[][]）
	/////////////////////////////////////////
	// 出力 
	// fr : 逆フーリエ変換後の関数の実部（[][]）
	// fi : 逆フーリエ変換後の関数の虚部（[][]）
	/////////////////////////////////////////

	var N = an.length;
	for( var i=0; i<N; i++ ){
		//kxに平行なFFTの実行
		FFT( an[ i ], bn[ i ], N, true );
	}
	for( var i=0; i<N; i++ ){
		for( var j=0; j<N; j++ ){
			fr[ i ][ j ] = an[ j ][ i ];		
			fi[ i ][ j ] = bn[ j ][ i ];		
		}
	}
	for( var i=0; i<N; i++ ){
		//ky軸に平行なFFTの実行
		FFT( fr[ i ], fi[ i ], N, true);
	}

}