|
|
●概要 FFTすることは決まったが、演算方式は種々提案されており、選択する必要がある。 ●データ数 ○理論値 当然、2N とする。既に、基数は10,000 、最大配列数は10,002,000 以内と決まっているので、実際の最大の2の冪数 L は、以下のように求められる。 2L ≦ 10,002,000 なので、 L ≦ 23 Nmax = 223 = 8,388,608 となる。つまり、800万データとなる。 ○設計値 高速化のために三角関数やビット逆順は、最初に全域のテーブルを生成し、実行時は参照のみで、値を得る方式とするので、最大データ配列長に影響される。三角関数(Sin)は、Double(8 バイト)、逆順はInteger(4 バイト)で持つので、Nmax * 12 が必要となる。理論値では、 Nmax * 12 = 8,388,608 * 12 = 100,663,296 = 100MB 100MBもメモリが必要になる。マシンによっては負担が大きい。今回は、L = 21 とした。これで、メモリは 25MB となる。これでも、FFTは、200万データ、つまり乗算では、400万桁は可能となる。但し、この値は自由に変えられる。ライブラリの仕様上、最大は100万桁としている。 余裕として、200万桁を線形畳込みできるようにする。 つまり、200万桁 * 2 / 4 = 100万データまでのFFTが必要となる。 ●演算方式 DFTは、複素数演算で、実は、乗算は、実数の場合の4倍となる。加減算は2倍。FFT乗算の場合は、結果として実数部分しか使わないので、虚数部は無駄のように思える（実際には無駄ではないが）。そこで、今回、2組の実数データを一度のDFTで同時に求める方式を採用する。また、コンボリューションも、データの 対称性を利用して、半分の乗算で済ませる方法とする。ここは、FFTそのもを紹介する場ではないので、なぜそうなるかは他書に依られたい。 ○同時FFT 今、N個の２組の実数データを、a()、b() とすれば、これらより、新しい複素数c()を、 c() = a() + j・b() (j は虚数単位) と、a()を実数部、b()を虚数部にそれぞれあてがう。このフーリエ変換を、C() とすれば、a()、b() のフーリエ変換、A()、B() のn番目の要素は、 A(n) =1/2((Cr(n) + Cr(N - n)) + j・(Cj(n) - Cj(N -n))) B(n) =1/2((Cr(n) + Cr(N - n)) - j・(Cj(n) - Cj(N -n))) 但し、 となる。つまり、一回のDFTで、2組のフーリエ変換が求まる。 ○コンボリューション 更に、フーリエ空間のコンボリューション、F(n) = A(n) ・ B(n) では、データの対称性から、 ・n ≦ N/2 のとき、 F(n) = A(n) ・ B(n) ・n ＞ N/2 のとき、 F(n) = Fr(N -n) - j・Fj(N - n) となるので、コンボリューションの複素数乗算が半分になる。が、実際には採用していない。 ●桁上げ処理 F() をＩFFT したものを、f() とすれば、 下の桁から、以下のようにする。最終結果配列、m() のサイズは、N+1 する(最上位に1つ付加する)。演算は、全てLong にて行う。 Redim m(N) 結果は、m() の、0 〜 N-2 に格納される。 ●実際のプログラムは？ Webを見ると、結構Cライブラリなどがソースで提供されており、その中から、適当にDLLとして利用させてもらうつもりでいたが、結構マシンやオーナ依存性が高く、.NET環境特に、VB.NET環境に馴染まないものであった。また。徹頭徹尾、VB.NET で統一するポリシー(あんまり意味はないが)を貫く意味でも、Webのライブラリを参考にさせて貰い自分で開発することにした。 |