DFT解析プログラム dft.py 技術ドキュメント

プログラムの動作

dft.py は、離散フーリエ変換（DFT）および高速フーリエ変換（FFT）の異なる実装を比較するためのPythonスクリプトです。このプログラムは、特定の時間領域信号を内部で生成し、以下の3つの方法で周波数領域に変換します。

dft 関数: 回転因子を事前に計算し、乗算で再帰的に利用することで計算を最適化した手動実装のDFT。
dft2 関数: 各ステップで三角関数（sin, cos）を直接計算する、最適化されていない手動実装のDFT。
fft 関数: numpy ライブラリの np.fft.fft を使用したFFT。

プログラムは、これらの各変換の実行時間を計測し、その差を標準出力（コンソール）に出力します。最終的に、生成された時間領域データと各変換方法による周波数領域の結果（実数部と虚数部）をCSVファイルに保存します。これにより、DFTの最適化手法（回転因子の利用）が計算速度に与える影響や、numpyによるFFTの高速性を比較・理解できます。コマンドライン引数でデータの点数を指定することが可能です。

原理

このプログラムは、時間領域の離散信号を周波数領域の離散信号に変換する離散フーリエ変換 (DFT) の原理に基づいています。

DFTは、以下の数式で定義されます。 $$X_k = \sum_{n=0}^{N-1} x_n e^{-j 2\pi k n / N}$$ ここで、$X_k$ は$k$番目の周波数成分、$x_n$ は$n$番目の時間領域サンプル、$N$ はサンプル数、$j$ は虚数単位です。

プログラム内の dft および dft2 関数では、入力 y が複素数であることを考慮し、y[i].real と y[i].imag を使用して計算しています。プログラムが内部で生成する信号 y は実数ですが、計算の内部表現は複素数として扱われます。指数部分の符号を制御する定数 $C$ を用いていますが、main 関数では常に $C=1.0$ と設定されています。このため、実質的には以下の形式で変換が行われています。

\[X_k = \sum_{n=0}^{N-1} x_n e^{j 2\pi k n / N}\]

この形式では、$e^{j\theta} = \cos(\theta) + j\sin(\theta)$ であるため、各周波数成分の実数部 $X_{k,real}$ と虚数部 $X_{k,imag}$ は以下のように計算されます。 $$X_{k,real} = \sum_{n=0}^{N-1} (x_{n,real} \cos(\phi_{nk}) - x_{n,imag} \sin(\phi_{nk}))$$X_{k,imag} = \sum_{n=0}^{N-1} (x_{n,real} \sin(\phi_{nk}) + x_{n,imag} \cos(\phi_{nk}))$$ ここで、$\phi_{nk} = 2\pi k n / N$ です。

dft2 関数 は、この計算を直接的に実装しており、内側のループで毎回 $\cos(\phi_{nk})$ と $\sin(\phi_{nk})$ を計算します。この方法は計算コストが高いです。

dft 関数 は、回転因子の概念を用いて最適化されています。回転因子 $W_N = e^{j 2\pi / N}$ を利用し、そのべき乗 $W_N^{ik} = e^{j 2\pi ik / N}$ の実数部と虚数部を事前に計算して配列 wr と wi に格納します。これにより、内側のループでの高価な三角関数の計算を減らし、乗算と加算で処理できるため、dft2 よりも高速に動作します。

fft 関数 は、numpy.fft.fft を利用しています。これは、Cooley-Tukeyアルゴリズムなどの高速フーリエ変換 (FFT) アルゴリズムの実装であり、DFTを大幅に高速化するものです。FFTは、DFTの計算量を $O(N^2)$ から $O(N \log N)$ へと削減します。なお、numpy.fft.fftは標準的なDFTの定義である $e^{-j 2\pi k n / N}$ を使用しているため、dftおよびdft2の結果とは虚数部分の符号が反転します。

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

このプログラムは、以下の非標準ライブラリに依存しています。

numpy: 数値計算全般に使用されます。特に高速フーリエ変換（FFT）の実装（np.fft.fft）、三角関数（np.sin, np.cos, np.tan）、および配列操作に利用されます。

インストールは、Pythonのパッケージマネージャー pip を使用して行います。

pip install numpy

必要な入力ファイル

dft.py は、外部の入力ファイルを必要としません。プログラムが内部で時間領域の信号データを生成します。

生成される信号は func(x) 関数によって定義され、以下の3つのサイン波の重ね合わせにわずかなランダムノイズが加えられたものです。 $$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) + \text{random\_noise}$$ ここで、振幅 $A_i$、周波数 $f_i$、位相 $p_i$ はそれぞれ以下の通りです。

$A_1=1.0, f_1=1.5, p_1=\pi/4.0$
$A_2=0.3, f_2=3.0, p_2=\pi/3.0$
$A_3=0.5, f_3=10.0, p_3=\pi/6.0$

ランダムノイズは $0.03 \times \text{np.random.rand()}$ によって導入されます。時間軸 x の範囲はデフォルトで $0.0$ から $1.0$ までです。

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

プログラムは、実行結果を dft.csv というCSVファイルに保存します。

ファイル名: dft.csv
内容:
- ヘッダー行: x, y, FT(DFT).r, FT(DFT).i, FT(DFT2).r, FT(DFT2).i, FT(FFT).r, FT(FFT).i
- データ行: 各行が時間領域のデータ点および対応する周波数領域の変換結果を示します。
  - x: 時間軸上のサンプリング点 ($0.0$ から $1.0$ までの範囲)。
  - y: func(x) によって生成された元の時間領域信号の振幅。
  - FT(DFT).r: dft 関数によるDFT結果の実数部。
  - FT(DFT).i: dft 関数によるDFT結果の虚数部。
  - FT(DFT2).r: dft2 関数によるDFT結果の実数部。
  - FT(DFT2).i: dft2 関数によるDFT結果の虚数部。
  - FT(FFT).r: np.fft.fft 関数によるFFT結果の実数部。
  - FT(FFT).i: np.fft.fft 関数によるFFT結果の虚数部。

出力される周波数成分は、データ点数 ndata の半分 ($N/2$) までの正の周波数成分のみです。なお、dftとdft2の周波数インデックスは1から始まるため、0Hz（DC成分）は含まれず、最低周波数は$\Delta f$となります。

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

基本的なコマンドライン引数は以下の通りです。

python dft.py [ndata]

[ndata]: データ点数を指定するオプションの整数引数です。指定しない場合、プログラムはデフォルト値の 1024 を使用します。

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

例1: デフォルトのデータ点数で実行

ndata を指定せずに実行する場合、デフォルト値である 1024 のデータ点数が使用されます。

python dft.py

実行結果の例:

Make data by func(x) with random scattering
x = (0.0, 1.0, 0.0009765625)
ndata=1024
Perform DFT (using rotation factor, faster): 0.0875248908996582 sec

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

Perform FFT(fast): 0.00010972023010253906 sec

Save FFTed data to [dft.csv]

この実行結果では、デフォルトの ndata=1024 が使用され、dft、dft2、fft の順に処理が行われ、それぞれの実行時間が表示されます。dft は dft2 よりも高速であり、fft が圧倒的に高速であることがわかります。最後に、結果が dft.csv に保存されたことが示されます。

例2: データ点数を `2048` に指定して実行

データ点数を 2048 に増やして実行する場合。

python dft.py 2048