離散フーリエ変換の基礎

Note

本節のポイント

• 離散フーリエ変換(DFT)の数学的定義を理解する。
• ndarrayを用いて、DFTを定義通りに実装する方法を学ぶ。
• DFTを行列演算（線形変換）として解釈し、線形代数で学んだ知識と結びつける。
• 直接計算の計算量 という課題を理解する。

コンピュータで扱うデータは、一定の時間間隔でサンプリングされた離散的な数値の列です。このような離散データに対してフーリエ変換を行う手法を離散フーリエ変換 (Discrete Fourier Transform, DFT) と呼びます。

DFTの定義

長さ の複素数データ系列 に対して、その離散フーリエ変換 は次のように定義されます。

ここで は虚数単位です。 は、元の信号に含まれる周波数成分（振幅と位相）を表します。

逆に、 から元の信号 を復元する逆離散フーリエ変換 (IDFT) は以下の通りです。

DFTの実装 (ndarray)

ndarray入門で学んだndarrayを使用して、定義式に基づいてDFTを計算するコードを実装してみましょう。複素数の扱いはnum-complexクレート（Complex64）を使用します。

[dependencies]
ndarray = "0.17"
num-complex = "0.4"

use ndarray::{Array1, ArrayView1};
use num_complex::Complex64;
use std::f64::consts::PI;

/// 定義式に基づいたDFTの直接計算
fn dft(x: ArrayView1<Complex64>) -> Array1<Complex64> {
    let n = x.len();
    let mut x_k = Array1::zeros(n);

    for k in 0..n {
        let mut sum = Complex64::new(0.0, 0.0);
        for n_idx in 0..n {
            let angle = -2.0 * PI * (k as f64) * (n_idx as f64) / (n as f64);
            // オイラーの公式 exp(iθ) = cos θ + i sin θ を用いた計算
            let exponent = Complex64::from_polar(1.0, angle);
            sum += x[n_idx] * exponent;
        }
        x_k[k] = sum;
    }
    x_k
}

fn main() {
    // 例：4点のデータ
    let data = Array1::from(vec![
        Complex64::new(1.0, 0.0),
        Complex64::new(2.0, 0.0),
        Complex64::new(3.0, 0.0),
        Complex64::new(4.0, 0.0),
    ]);

    let result = dft(data.view());

    for (k, val) in result.iter().enumerate() {
        println!("X[{}] = {:.3} + {:.3}i", k, val.re, val.im);
    }
}

X[0] = 10.000 + 0.000i
X[1] = -2.000 + 2.000i
X[2] = -2.000 + -0.000i
X[3] = -2.000 + -2.000i

線形代数としての解釈：DFT行列

線形代数で学んだように、DFTはベクトルに対する線形変換とみなすことができます。 とおくと、DFTは以下のような行列とベクトルの積で表現できます。

この行列を DFT行列 と呼びます。ndarrayの行列積dotを利用すれば、DFTはさらに簡潔に記述可能です。

fn dft_matrix_method(x: &Array1<Complex64>) -> Array1<Complex64> {
    let n = x.len();
    // DFT行列の作成
    let mut w = Array2::<Complex64>::zeros((n, n));
    for k in 0..n {
        for n_idx in 0..n {
            let angle = -2.0 * PI * (k as f64) * (n_idx as f64) / (n as f64);
            w[[k, n_idx]] = Complex64::from_polar(1.0, angle);
        }
    }
    // 行列とベクトルの積
    w.dot(x)
}

計算量と課題

いずれの実装でも、全ての周波数 に対して 回の和をとるため、計算量は となります。

データ数 が大きくなると計算時間が爆発的に増加します。

• のとき、 回の複素数演算。
• のとき、 回の複素数演算（一般的なPCでは実用困難）。

この の壁を打破し、実用的な速度でフーリエ変換を行う手法が、次節で学ぶ高速フーリエ変換 (FFT) です。

次節では、計算効率を劇的に改善する FFTのアルゴリズムとライブラリの使い方について学びます。