DFT と FFT

離散フーリエ変換（DFT）と高速フーリエ変換（FFT）については、 参考資料にて挙げているWebサイトに大変詳しくかつ丁寧に 説明がありますので、私が改めて説明を加えるよりは、 そちらを参照いただくのが良いと思います。 ここでは、単に、DFTにおける回転子の符号についての考察と、 十進BASIC にて FFT を実現するためのサブルーチンの作成を 行います。 プログラムとしてはあまり美しくないかもしれませんがあしからず。 可能ならば、将来十進BASICに FFT が組み込まれるのが一番楽なのですが．．．

離散フーリエ変換（DFT）はN個のデータの要素番号 k を、 k=1,2,3,...,N とすると、以下の式で定義されます。

変換

逆変換

ここで、W は以下で定義され、回転子と呼ばれます。

また、この回転子には以下のような重要な特性があります。


回転子は一回転をＮ等分しているということがわかれば、 どちらの特性もおのずとわかります。 ２番目の特性は時計の針のようなもので、一回転加えても（引いても） またもとの場所にもどります。 この２番目の特性を利用すると、DFTの周期性が導けます。

さて、ここで DFT には実はいくつかの定義の仕方があります。 基本的に、変換・逆変換をあわせた演算は同じ結果になるのですが、 次のパラメータに違いがあります。

• 変換・逆変換をあわせると全体として 1/N の係数がかかりますが、 この係数を変換・逆変換に分配するやり方
  （例、変換 1、逆変換 1/N や、変換 1/sqrt(N)、逆変換 1/sqrt(N) など）
• 回転子の EXP の肩にのっている項の符号の取り方（＋/−）

ここで、大きさに関する項 1/N をどのようにかけるかは、 スペクトルの絶対値を問題にしない限りどのようにとっても、 ほぼ問題ありません。逆にスペクトルの絶対値にどう影響が でるかというと、もっとも極端な例として、直流 1 をサンプルした 場合を考えます。DFT変換後の値は、最初の X1 を除いて すべて 0 となりますが、X1 の値は上記の式ですとサンプル数 N と 等しくなります。これは言い換えるとスペクトル強度（振幅?）が サンプル数に比例することになります。 同じ信号をサンプルしたときに、サンプリング間隔の違いにより 強度が違うというのはあまり気持ちよくありません。 これを避けるには単純に変換時に 1/N をかければよいということに なりますが、反対にスペクトルのプロファイルが与えられている ときに、周波数サンプル間隔をかえてサンプルしたデータを逆変換 すると強度がサンプル数に依存するということにもなります。 いずれにせよ、このパラメータはシミュレーションの対象に依存 しますが、変換・逆変換を常にペアで行う場合はほとんど気にしなく ても問題ありません。
さて、回転子の符号の取り方についてですが、これもスペクトル 強度を見るだけであれば特にどちらでも問題はありません。 しかし、通常多くの技術計算では、 周波数領域においてフィルタリング、 位相特性の確認や群遅延の評価などを行います。 この場合、符号のとり方が問題となります。
結論としては、S = j ω なる置き換えを行う場合 （もしくは、群遅延＝-dθ/dωにて定義する場合） はマイナスにとるということになります。
以下、10進BASICのプログラムを用いて、上記を確認します。 サンプルファイルのプログラムファイル名は dft.bas です。 キーとなる DFT 計算については、外部サブルーチン（外部副関数）にて 上記定義式をほぼそのまま記述しています。 （したがって遅いです．．．）

サブルーチンの内容表示 （テキストエディタによるファイル表示で、行番号はエディタによるものです）

この例では、既知の2次LPFフィルタインパルスレスポンス

を時間軸上でサンプリングし、 これに DFT をかけることで得られるスペクトル強度および 位相特性を、ラプラス逆変換前の式から得られる値と比較しています。
入力となるインパルスレスポンスは、


これをサンプリングして DFT をかけたスペクトル強度は、 DFT の対象性から後半半分のデータを無視すると、


（青い線がラプラス逆変換前の式から得られた値です。）
位相特性は、


（水色がラプラス逆変換前の式から得られた値です。）
これを逆変換すると、元の波形が得られます。


これを、DFT サブルーチン内の回転子の符号を逆にして、
exp( complex(0, 2 * pi / nx * dir * (k-1) * (n-1)) )
とプログラムを書き直して、試してみると、 位相特性のみがきれいに正負逆になります。 当然ですが、フィルタリングなどで S = jω と置いて計算した、 周波数領域での係数をこのデータにかけるとおかしな結果が得られます。 ここではあえて結果を示しませんが、 ご興味のある方は、いろいろ変更して試してみてください。
結局、回転子は位相の進行方向の定義に依存していて、あまり意識せずに 行っている、S = jω なる置き換えもまた、実は位相の進行方向を 定義しているため、両者が一致（位相遅れが負の位相になる）して いないと計算がおかしなことになる、ということになります。

基数２の周波数間引き Cooley-Tukey 型 FFT を用いたサブルーチンの例です。 直接サンプルファイルをご確認いただければ良いかと思います。 FFT のアルゴリズムについては、参考文献に挙げたWebサイトを 参照ください。
サブルーチンのみのサンプルファイル fft_sub.bas と、 検算のため DFT のときに用いたものと同じインパルスレスポンス を計算するプログラム、fft.bas の２つがあります。 サブルーチンは、もっとも単純に再帰呼び出しを利用して周波数間引きを 実現しています。残念ながら 2 の倍数ではないデータ数に対しては、 単純に DFT を実行しますので、非常に遅いです。 しかし、2 のべき乗のデータ数で FFT をかけるとなんとか使えると思います。 大きなデータをいきなり計算する前に、少し小さめのデータ数で試される ことをお勧めいたします。 ちなみに、2.4GHz Pentium4、DDR333 のマシンで、 2^17(131072)個のデータを FFT変換・逆変換すると 15秒程度ですが、 2^11-1(2047)個のデータを変換・逆変換すると、これまた15秒程度かかります．．．