FFT デモ — 高速フーリエ変換

calx の FFT モジュールは混合基数 (radix-2, 3, 5) + Bluestein で任意長に対応し、 SIMD (AVX2/SSE) による自動高速化を備える。このページでは 4 つのデモを通じて基本操作を紹介する。出力はすべて calx で実際に計算したものである。ソースコードはページ末尾に掲載している。

Demo 1: 複素 FFT (順変換と逆変換)

4 点の複素データに対して FFT を実行し、逆変換で元に戻ることを確認する。 FFT<C>::fft は in-place で動作し、FFT<C>::ifft は $1/N$ 正規化を含む。

入力: (1,0) (2,0) (3,0) (4,0) FFT: (10,0) (-2,2) (-2,0) (-2,-2) IFFT: (1,0) (2,0) (3,0) (4,0) ✓ 元に戻った

$X[0] = 1+2+3+4 = 10$ (DC 成分 = 全サンプルの和)。逆変換後の値が入力と一致し、FFT→IFFT のラウンドトリップが正しいことがわかる。

Demo 2: 実数 FFT とスペクトル解析

8 点の実数信号 $x = [1, 0, -1, 0, 1, 0, -1, 0]$ (周期 4 のコサイン波) を RealFFT で変換し、振幅スペクトルを取得する。出力は共役対称性を利用して $N/2+1 = 5$ ビンに圧縮される。

入力: 1 0 -1 0 1 0 -1 0 (N=8) スペクトル (N/2+1 = 5 ビン): X[0] = (0, 0) |X[0]| = 0 X[1] = (0, 0) |X[1]| = 0 X[2] = (4, 0) |X[2]| = 4 ← 周波数 2/8 = 0.25 にピーク X[3] = (0, 0) |X[3]| = 0 X[4] = (0, 0) |X[4]| = 0

入力は周期 4 のコサイン波なので、$N = 8$ サンプルにちょうど 2 周期入る。 $X[2]$ だけが非ゼロとなり、期待どおりの結果である。実信号の共役対称性 $X[N{-}k] = \overline{X[k]}$ により、$N/2+1$ ビンで完全に情報を保持する。

Demo 3: 畳み込み (FIR フィルタ)

2 つの実数配列の畳み込みを RealFFT::convolve で高速に計算する。信号にインパルス応答を畳み込む FIR フィルタ適用と等価である。出力サイズは $|a| + |b| - 1$。

a = [1, 2, 3] b = [4, 5] convolve(a, b) = [4, 13, 22, 15]

検算: $(1 + 2x + 3x^2)(4 + 5x) = 4 + 13x + 22x^2 + 15x^3$。畳み込みと多項式乗算は同じ演算であるが、 FFT::convolve は複素係数用、RealFFT::convolve は float / double の波形用に最適化されている。

直接計算は $O(nm)$ だが、FFT 経由の畳み込みは $O((n+m) \log(n+m))$。長い FIR フィルタほど FFT の優位性が大きくなる。

Demo 4: 位相スペクトルとアンラップ

スペクトルの位相を phase_spectrum で取得し、 unwrap_phase で $2\pi$ のジャンプを除去して連続位相にする。群遅延 $-d\phi/d\omega$ の計算や位相応答の可視化に必要な前処理である。

ビン: phase (wrapped) phase (unwrapped) [0] 0.000 0.000 [1] -2.356 -2.356 [2] 3.142 3.142 [3] -2.356 -2.356 [4] 0.000 -6.283

phase_spectrum は std::arg() の主値 $(-\pi, \pi]$ を返すため、位相が $\pm\pi$ でジャンプする。unwrap_phase を適用すると連続した位相応答が得られ、差分で群遅延を求められる。

ソースコードと実行方法

example_fft_demo.cpp (完全なソースコード)

// example_fft_demo.cpp — FFT demo
#include <math/fft/fft.hpp>
#include <math/fft/fft_utils.hpp>
#include <iostream>
#include <iomanip>
using namespace calx;

int main() {
    using C = std::complex<double>;

    // --- Demo 1: Complex FFT ---
    std::cout << "--- Complex FFT ---\n";
    std::vector<C> data = {1, 2, 3, 4};
    auto orig = data;
    FFT<C>::fft(data);
    std::cout << "FFT: ";
    for (auto& v : data) std::cout << v << " ";
    std::cout << "\n";
    FFT<C>::ifft(data);
    std::cout << "IFFT: ";
    for (auto& v : data) std::cout << v << " ";
    std::cout << "\n";

    // --- Demo 2: Real FFT ---
    std::cout << "\n--- Real FFT ---\n";
    std::vector<double> signal = {1, 0, -1, 0, 1, 0, -1, 0};
    auto spec = RealFFT<double>::fft(signal);
    auto amp = amplitude_spectrum<double>(spec);
    for (size_t i = 0; i < spec.size(); ++i)
        std::cout << "X[" << i << "] = " << spec[i]
                  << "  |X| = " << amp[i] << "\n";

    // --- Demo 3: Convolution ---
    std::cout << "\n--- Convolution ---\n";
    std::vector<double> a = {1, 2, 3}, b = {4, 5};
    auto c = RealFFT<double>::convolve(a, b);
    std::cout << "convolve = ";
    for (auto v : c) std::cout << v << " ";
    std::cout << "\n";

    // --- Demo 4: Phase unwrap ---
    std::cout << "\n--- Phase unwrap ---\n";
    std::vector<double> sig2 = {1, 0.5, 0, -0.5, -1};
    auto spec2 = RealFFT<double>::fft(sig2);
    auto phase = phase_spectrum<double>(spec2);
    auto uphase = unwrap_phase<double>(phase);
    for (size_t i = 0; i < phase.size(); ++i)
        std::cout << std::fixed << std::setprecision(3)
                  << "  [" << i << "] " << phase[i]
                  << "  ->  " << uphase[i] << "\n";

    return 0;
}

API の詳細は FFT API リファレンスを参照のこと。

ビルドと実行

cd calx
mkdir build && cd build
cmake .. -G "Visual Studio 17 2022" -A x64
cmake --build . --config Release --target example-fft-demo
examples\Release\example-fft-demo.exe