1.一次元高速フーリエ変換

　一次元関数 f(x) を離散化したうえで数値的にフーリエ変換する方法とそのプログラムを以下で 紹介します。
　なおここでの計算方法は「光とフーリエ変換」の「5.高速フーリエ変換」を参考にしています。
　(ただしこのページに載せたプログラムは独自に作成したものです)
　2007年1月24日 説明を一部追加するとともに「6.プログラム使用上の注意」を追加しました

1.離散化フーリエ変換の式

　まず f(x) を間隔Tに離散化して x= n*T (nは整数)の値のみを考えるようにし、さらに n を0からN-1 までの N個の点に制限します。

　この離散化関数を、フーリエ変換したものをF(ν) とします。(引数 νは周波数になります)
　この段階で　F(ν)は周波数に対して連続関数となっていますが、次にこれを離散化します。
元の関数 f(x) を間隔TのN個の点に離散化したことにより、フーリエ変換した式は間隔 1/(N*T) で離散化することに なります(詳しい事は上の参考文献や、高速フーリエ変換に関する本を参照してください。)。

　離散化した結果の式を書くと以下の様になります。なおこの式で k は0からN-1の値をとる整数で、引数の k/( N*T) は 離散化した周波数のｋ-1番目(k =0 が1番目だから)の値を表します。

...(1)

　注意 :以下の式では(1)式のΣの前に掛けてある間隔T が省略されています。また例として 載せているプログラムでもこのTが掛けてないので、実際に使用するときは間隔Tを掛けるように修正してください。

2.高速フーリエ変換とは

　(1)式で　F( k/NT) をF(k)に、f( nT )をf(n)に、exp(-i2π/N) を W と書き換えて下の(2)式のように 簡単に書くことにします。

...(2)

3.高速フーリエ変換のアルゴリズム

　まず離散点の数Nを、N = r₁*r₂ と二つの数の積に分解します。
　つぎに以下に示すようにk を二つの数k₀、k₁で、n も同じように二つの数n₀、n₁を 使って書き直します。

...(3)

...(4-2)

...(4-1)

4.具体的な計算の手順

　ここではN=2^m の場合のより具体的な手順について記します。

�@r₁=2、r₂= N /2 の場合

...(5)

	...(5-1)
	...(5-2)

	...(6-1)
	...(6-2)

	...(7-1)
	...(7-2)

5.プログラムの例

　プログラムは C++ で書いてあります。ソースファイルの拡張子は cpp にしておいて下さい。
　以下に各部分について簡単な説明を入れておきます。

　関数 DecimFrq_FFTCall が計算の本体で第一引数に計算したい関数を離散化した値を収めた配列を与えます。
この配列に計算したフーリエ変換値をいれて送り返すので、元の関数の値を保持しておきたいときには注意してください。
なおこのフーリエ変換値には式(1)に掛けてある間隔 T が掛けてないので注意してください。
この引数は complex型になっています。
第二引数は第一引数の配列要素数を与えてください。

　DecimFrq_FFT は繰り返し計算をおこなう部分です。

　CalcRev_BinaryExp は計算結果の順番を入れ替えている部分です。

宣言ファイルfftdef.h

/*     高速フーリエ変換　周波数間引き法プログラム     関数宣言ファイル       関数の分割数が2^N のもの */ #ifndef ___fastfrt  #define ___fastfrt #include <complex> using namespace std;  //introduces namespace std /*      関数宣言 */ //高速フーリエ変換 int    DecimFrq_FFTCall(complex<double>* ,unsigned long ); //高速フーリエ変換計算部分 int    DecimFrq_FFT( complex<double>* ,unsigned long ,unsigned long ,complex<double> ) ; //二進表記で逆の順にした数値を作成 unsigned long  CalcRev_BinaryExp( unsigned long ,unsigned long ); //与えた数値が2^N の形をしているか調べる int    Check_PowerOf2( unsigned long , unsigned long* ) ;   #endif

定義ファイルfftdef.cpp

#include "fftdef.h" #include <cmath> /*      func_array ==>計算終了時はフーリエ変換したものが入っている                    元の関数は、別に保管しておく事 */ int    DecimFrq_FFTCall(complex<double> *func_array ,//対象離散化関数                         unsigned long div_numb )//関数分割数 を与える {     //div_numb のチェック     unsigned long pow_num ; //関数分割数2^N の Nを与える          if( !Check_PowerOf2( div_numb ,&pow_num ) )          return 0 ; //2^N の形になっていない場合 //    unsigned long div_numb = pow(2,pow_num ) ;       //高速フーリエ変換の計算     double pi = acos(-1.0) ;     complex<double> exp_arg = complex<double>(0.0 ,-2.0*pi/div_numb ) ;     complex<double> the_factor = exp(exp_arg) ;          DecimFrq_FFT(func_array ,0 ,div_numb ,the_factor) ;       //フーリエ変換 func_array は順番が変わっているので元の順に並び替える     int the_roop = 0 ;     unsigned long *flag_array = new unsigned long[div_numb] ;        for(the_roop = 0 ;the_roop < div_numb ;the_roop++)          flag_array[the_roop] = 0 ;       unsigned long rebny_num = 0 ;     complex<double> temp_value ;          for(the_roop = 0 ;the_roop < div_numb ;the_roop++ ){            if( flag_array[the_roop] == 0 ){//まだ交換されていない            rebny_num = CalcRev_BinaryExp((unsigned long)the_roop ,pow_num ) ;            //入れ替え            temp_value = func_array[rebny_num] ;            func_array[rebny_num] =func_array[the_roop] ;            func_array[the_roop] = temp_value ;              flag_array[the_roop] = flag_array[rebny_num] = 1 ;          }          }          delete[] flag_array ;//動的配列解放       return 1 ; } /*     高速フーリエ変換　周波数間引き法 計算本体       関数分割数 2^N */ int    DecimFrq_FFT( complex<double> *func_array ,//対象離散化関数                      unsigned long step_num , //繰り返し回数 root で０                      unsigned long div_numb , //関数分割数(rootで2^N )                      complex<double> the_factor ) //exp(-i*2*PI/(2^N) ) {     unsigned long nwdiv_numb = div_numb/2 ;     complex<double> *fun_pre = &(func_array[0]) ;     complex<double> *fun_aft = &(func_array[nwdiv_numb]) ;     unsigned long stp_factor= pow(2,step_num) ;  //2^step_num     complex<double>  tmp_sav1 ,tmp_sav2 ; //値の一時保管用          for( int the_roop=0 ;the_roop < nwdiv_numb ;the_roop++ ){          tmp_sav1 = func_array[the_roop] ;          tmp_sav2 = func_array[ the_roop + nwdiv_numb] ;            fun_pre[the_roop] = tmp_sav1 + tmp_sav2 ;          fun_aft[the_roop] =(tmp_sav1 - tmp_sav2 )*pow(the_factor ,the_roop*stp_factor) ;          }          if( nwdiv_numb != 1 ){//再帰呼び出し            step_num++ ; //繰り返し回数増やす          DecimFrq_FFT(fun_pre ,step_num ,nwdiv_numb ,the_factor ) ;          DecimFrq_FFT(fun_aft ,step_num ,nwdiv_numb ,the_factor ) ;          }       return 1 ; }   /*    二進表示で逆に読む数値を作成    (unsigned long とする) */ unsigned long  CalcRev_BinaryExp( unsigned long the_source , //対象となる数字                                   unsigned long the_num ) //pow(2,N ) の N {     unsigned long num_revby = 0UL ;          for( int the_roop = 0 ;the_roop <= the_num-2 ;the_roop++ ){            if( the_source & 1UL ){            num_revby += 1UL ;          }            the_source >>= 1 ;          num_revby <<= 1 ;          }          //the_roop = the_num -1        if( the_source & 1UL ){          num_revby += 1UL ;        }       return num_revby ; } /*     与えた数値が2^N の形をしているかチェックするとともに、Nを獲得する       返し値         0 :2^N の形をしていない         1 :2^N の形をしている */ int    Check_PowerOf2( unsigned long the_number , //対象となる数値                        unsigned long *the_power )//2^N のN {     bool the_flag = false ;          *the_power = 0 ;//初期化          while( the_number && !the_flag ){//the_number 二進表示で0以外の数値が含まれる          if( the_number & 1UL )//二進表示の末尾が1の場合            the_flag = true ;            the_number>>=1 ;//ビットシフトして二進表示を左へずらす          (*the_power)++ ;        }          if( the_number==0 && the_flag ){//2^N の形をしている場合のみここに来る          (*the_power)-- ;//最後にひとつ余計にカウントするのでひとつ戻す            return 1 ;        }       return 0 ; }

---