DFT/FFTデモンストレーションスクリプト: dft.py

プログラムの動作

dft.py は、離散フーリエ変換 (DFT) および高速フーリエ変換 (FFT) のさまざまな実装を比較デモンストレーションすることを目的としたPythonスクリプトです。

主な機能は以下の通りです。

データ生成: 乱数によってわずかに変調された3つの異なる周波数の正弦波を合成し、シミュレートされた時間領域データを作成します。
フーリエ変換の実行:
- dft: 回転因子を事前に計算して利用する、最適化されたDFT実装。
- dft2: 各計算ステップで三角関数を直接評価する、基本的なDFT実装。
- fft: NumPyライブラリの np.fft.fft 関数を利用した高速フーリエ変換 (FFT)。
パフォーマンス計測: 各フーリエ変換の実行時間を計測し、標準出力に表示することで、実装間の効率の違いを示します。
結果出力: 元の時系列データと、各フーリエ変換によって得られた周波数スペクトルの実部および虚部を、dft.csv というCSVファイルに保存します。

このスクリプトは、異なるフーリエ変換アルゴリズムの計算効率と結果を比較することで、それぞれの特性を理解する手助けをします。

原理

本スクリプトは、離散フーリエ変換 (DFT) とその高速アルゴリズムである高速フーリエ変換 (FFT) の原理に基づいています。

離散フーリエ変換 (DFT)

DFTは、離散的な時間領域の信号を、離散的な周波数領域の信号に変換する数学的手法です。N個のサンプル \(x_n\) (ここで \(n = 0, 1, \ldots, N-1\)) を持つ信号の場合、そのDFT \(X_k\) (ここで \(k = 0, 1, \ldots, N-1\)) は以下の式で定義されます。

\[X_k = \sum_{n=0}^{N-1} x_n e^{-i2\pi \frac{kn}{N}}\]

ここで、\(i\) は虚数単位、\(e\) は自然対数の底です。

本スクリプトの dft および dft2 関数は、mode 引数に応じて指数部の符号を制御します。main 関数から 1 が渡されるため、C = 1.0 が設定され、以下の形式で計算が行われます。

\[X_k = \sum_{n=0}^{N-1} x_n e^{i2\pi \frac{kn}{N}}\]

この式は、一般的なDFTの定義とは指数部の符号が逆ですが、振幅スペクトルは同等になります。結果として得られる虚部成分は符号が反転します。複素数 \(x_n = R_n + iI_n\) と回転因子 \(e^{i\phi} = \cos\phi + i\sin\phi\) の積は、以下のように実部と虚部に展開されます。

\[(R_n + iI_n)(\cos\phi + i\sin\phi) = (R_n\cos\phi - I_n\sin\phi) + i(R_n\sin\phi + I_n\cos\phi)\]

dft および dft2 関数はこの展開を用いてDFTを計算しています。

dft関数 (最適化DFT)

dft 関数は、回転因子 \(e^{i2\pi \frac{k}{N}}\) の実部と虚部を事前に計算し、配列 wr と wi に格納します。これにより、内側のループで三角関数を繰り返し計算するオーバーヘッドを削減し、テーブルルックアップで済ませることで高速化を図っています。具体的には、\(e^{i2\pi \frac{j}{N}}\) の実部と虚部を事前に計算して格納しておき、必要な \(e^{i2\pi \frac{kn}{N}}\) は \(j = (kn) \pmod N\) のインデックスで参照します。

dft2関数 (基本的なDFT)

dft2 関数は、DFTの定義式に忠実に、各ステップ (内側のループ) で回転因子の実部と虚部 (\( \cos(2\pi \frac{kn}{N})\) と \( \sin(2\pi \frac{kn}{N})\)) を直接計算します。これは実装が直感的である一方で、三角関数の計算コストが高いため、dft 関数と比較して低速になります。

高速フーリエ変換 (FFT)

fft 関数は、NumPyライブラリに実装されている np.fft.fft を利用します。FFTは、Cooley-TukeyアルゴリズムなどによってDFTの計算量を \(O(N^2)\) から \(O(N \log N)\) に大幅に削減したアルゴリズムです。NumPyのFFTは標準的なDFT定義 (\(e^{-i2\pi \frac{kn}{N}}\)) に従って計算されます。

データ生成 (func関数)

func 関数は、以下の形式で合成された正弦波に微小なランダムノイズを加えたデータを生成します。

\[y(x) = A_1 \sin(2\pi f_1 x + p_1) + A_2 \sin(2\pi f_2 x + p_2) + A_3 \sin(2\pi f_3 x + p_3) + \epsilon\]

ここで、\(A\) は振幅、\(f\) は周波数、\(p\) は位相、\(\epsilon\) はランダムノイズです。これにより、複数周波数成分を含む、より現実的なテストデータが作成されます。

必要な非標準ライブラリとインストール方法

このプログラムの実行には、以下の非標準ライブラリが必要です。

NumPy: 数値計算、特に配列操作やFFTに利用されます。

これらのライブラリは、以下の pip コマンドでインストールできます。

pip install numpy

csv, sys, time はPythonの標準ライブラリに含まれているため、追加のインストールは不要です。

必要な入力ファイル

このプログラムは、外部からの入力ファイルを必要としません。プログラム内で合成正弦波データを生成します。

生成される出力ファイル

プログラムを実行すると、以下のCSVファイルが生成されます。

dft.csv
- 内容: 生成された元の時系列データと、3種類のフーリエ変換（dft、dft2、fft）の結果が格納されます。各変換結果は実部と虚部に分けて出力されます。
- ヘッダー行: 以下の形式でデータ列の内容を示します。 x, y, FT(DFT).r, FT(DFT).i, FT(DFT2).r, FT(DFT2).i, FT(FFT).r, FT(FFT).i
  - x: 時系列データの独立変数（時間または空間座標）。
  - y: func(x) によって生成された元の時系列データ。
  - FT(DFT).r: dft 関数によるフーリエ変換結果の実部。
  - FT(DFT).i: dft 関数によるフーリエ変換結果の虚部。
  - FT(DFT2).r: dft2 関数によるフーリエ変換結果の実部。
  - FT(DFT2).i: dft2 関数によるフーリエ変換結果の虚部。
  - FT(FFT).r: np.fft.fft 関数によるフーリエ変換結果の実部。
  - FT(FFT).i: np.fft.fft 関数によるフーリエ変換結果の虚部。

コマンドラインでの使用例 (Usage)

dft.py スクリプトは、以下の形式でコマンドラインから実行できます。

python dft.py [ndata]

[ndata]: データ点の総数を指定するオプションの引数です。
- この引数を省略した場合、デフォルト値の 1024 が使用されます。
- 指定された場合、その整数値がデータ点数として使用されます。

使用例:

python dft.py 
# データ点数1024で実行

python dft.py 2048
# データ点数2048で実行

コマンドラインでの具体的な使用例

データ点数を 2048 と指定して dft.py を実行する例を示します。

python dft.py 2048

実行結果の説明

上記のコマンドを実行すると、以下のような情報が標準出力に表示されます（実行環境やndataの値によって具体的な時間は変動します）。

Make data by func(x) with random scattering
x = (0.0, 1.0, 0.00048828125)
ndata=2048
Perform DFT (using rotation factor, faster): 0.6543210987654321 sec

Perform DFT2 (calculate sin/cos everytime, slow): 1.2345678901234567 sec

Perform FFT(fast): 0.000012345678901234567 sec

Save FFTed data to [dft.csv]

この出力から、以下の点が観察できます。

ndata=2048 が適用されていることが確認できます。
各フーリエ変換処理の実行時間が表示されます。
- dft (回転因子事前計算) は、dft2 (三角関数毎回計算) よりも高速であることが示されます。
- fft (NumPy実装) が、圧倒的に高速であることが示されます。これは、FFTアルゴリズムがDFTよりも計算効率が非常に高い (\(O(N \log N)\) vs \(O(N^2)\)) ためです。
最後に、dft.csv というCSVファイルに結果が保存されたことが示され、プログラムが終了します。このCSVファイルには、元の時系列データと、それぞれのフーリエ変換結果の実部と虚部が記録されています。