FFT — 高速フーリエ変換

概要

任意長の離散フーリエ変換 (DFT) を高速に計算する。

• 混合基数 — radix-2, 3, 5 のバタフライで $N = 2^a \cdot 3^b \cdot 5^c$ を $O(N \log N)$ で処理
• Bluestein — 上記で分解できない任意長を Chirp-Z 変換で処理
• SIMD 最適化 — 2 のべき乗サイズで AVX2/SSE 自動ディスパッチ
• 実数 FFT — 共役対称性を利用して $N/2+1$ ビンに圧縮
• 2D / 3D FFT — 行・列・奥行き方向の逐次 1D FFT
• NTT — ModularInt 向け数論変換
• ユーティリティ — DCT, DST, Hilbert 変換, Walsh-Hadamard, Hartley, Hankel
• バッチ FFT — 複数信号の同時処理 (AVX2 インターリーブ), ステレオ FFT

全 API は calx 名前空間 (バッチのみ calx::simd_fft)。

ヘッダ

#include <math/fft/fft.hpp>         // FFT, RealFFT, FFTEngine, FFTDivMode
#include <math/fft/ntt.hpp>         // Ntt (数論変換)
#include <math/fft/fft2d.hpp>       // fft2d, ifft2d, fft2d_real
#include <math/fft/fft3d.hpp>       // Tensor3D, fft3d, ifft3d, fft3d_real, fftshift3d
#include <math/fft/fft_utils.hpp>   // DCT, DST, Hilbert, Walsh-Hadamard, Hartley, Hankel
#include <math/fft/fft_batch.hpp>   // batchFFT, batchIFFT, stereoFFT

リンク

// CMake
target_link_libraries(your_target PRIVATE calx_fft calx_ntt)

FFTDivMode

FFTEngine の正規化モードを指定する列挙型。

enum class FFTDivMode {
    None,       // 正規化なし
    Forward,    // 順変換で 1/N 正規化
    Inverse,    // 逆変換で 1/N 正規化 (既定)
    Both        // 両方で 1/sqrt(N) 正規化
};

注: FFT<T> と RealFFT<Real> は常に Inverse 相当 (逆変換で $1/N$ 正規化)。 FFTEngine のみ任意のモードを選択可能。

FFT<T, Real> クラス

複素数および ModularInt 向けの静的 FFT API。

template<typename T, typename Real = double>
class FFT;

T は std::complex<float>, std::complex<double>, または ModularInt<P>。 Real は回転子の実数型。

is_valid_size / next_valid_size

static bool is_valid_size(int n);
static int  next_valid_size(int n);

注: 有効でないサイズも fft() に渡せる (Bluestein にフォールバック)。

fft

static void fft(std::span<T> data);

順変換 (DFT)。正規化なし。サイズに応じて自動ディスパッチ:

• $2^k$ — SIMD 最適化 radix-2 (float/double), スカラー radix-2 (ModularInt)
• $2^a \cdot 3^b \cdot 5^c$ — 混合基数
• その他 — Bluestein (Chirp-Z)

ifft

static void ifft(std::span<T> data);

逆変換 (IDFT)。$1/N$ 正規化済み。複素数型では共役 + FFT + 共役 + $1/N$ で実装。

convolve

static std::vector<T> convolve(std::span<const T> a, std::span<const T> b);

返り値: サイズ a.size() + b.size() - 1 の畳み込み結果。内部で FFT サイズを 2 のべき乗に切り上げ。

用途: 型 T は複素数型 (std::complex<double> 等) 専用。 多項式乗算や複素係数の畳み込みに用いる。 実数波形 (float / double) の畳み込みには、 実数専用の RealFFT::convolve を使うこと (FIR フィルタ, インパルス応答適用等)。

ntt / intt / convolve_ntt (deprecated)

非推奨: これらのメソッドは非推奨です。代わりに Ntt クラスを使用してください。 Ntt は 5-smooth サイズ対応、AVX2 Montgomery 最適化、Bluestein 任意長など、上位互換の機能を提供します。

template<int P>
[[deprecated]] static void ntt(std::span<ModularInt<P>> data, ModularInt<P> primitive_root);

template<int P>
[[deprecated]] static void intt(std::span<ModularInt<P>> data, ModularInt<P> primitive_root);

template<int P>
[[deprecated]] static std::vector<ModularInt<P>> convolve_ntt(
    std::span<const ModularInt<P>> a, std::span<const ModularInt<P>> b,
    ModularInt<P> primitive_root);

Ntt クラス — 数論変換 (NTT)

64-bit 整数に対する Number Theoretic Transform (NTT)。 内部で AVX2 Montgomery NTT パイプラインを使用し、 5-smooth サイズ ($N = 2^a \cdot 3^b \cdot 5^c$, $b,c \in \{0,1\}$) および任意サイズ (Bluestein) に対応する。

#include <math/fft/ntt.hpp>
// リンク: calx_ntt

class Ntt;

Prime 構造体

struct Prime {
    uint64_t p;          // NTT-friendly 素数
    uint64_t g;          // 原始根
    size_t maxNttSize;   // 最大 NTT 長 = 2^47
};

推奨素数

3 素数とも $2^{47} \mid (p - 1)$ を満たし、最大 $2^{47}$ 点の NTT に対応。 係数の範囲は $[0, p)$ ($p \approx 2^{62}$)。

forward / inverse

static void forward(uint64_t* data, size_t N, const Prime& prime);
static void inverse(uint64_t* data, size_t N, const Prime& prime);

5-smooth サイズは混合基数 NTT で直接処理。それ以外は Bluestein (Chirp-Z) アルゴリズムで 5-smooth 畳み込みに帰着する。 inverse は $1/N$ スケーリング込み。

例外: std::invalid_argument — $N \nmid (p - 1)$ の場合。

pointwiseMul

static void pointwiseMul(uint64_t* r, const uint64_t* a, const uint64_t* b,
                          size_t N, const Prime& prime);

要素ごとの剰余乗算: $r[i] = (a[i] \times b[i]) \bmod p$。NTT 領域での畳み込みに使用。

polyMul

static void polyMul(uint64_t* r, const uint64_t* a, size_t an,
                     const uint64_t* b, size_t bn, const Prime& prime);

$\bmod p$ での多項式乗算。小サイズ ($an + bn \leq 16$) は schoolbook、それ以上は NTT 経由。 内部で Montgomery 変換を自動的に適用する。

サイズユーティリティ

使用例

#include <math/fft/ntt.hpp>
#include <vector>
#include <iostream>

int main() {
    using calx::Ntt;
    auto prime = Ntt::prime1();

    // 多項式乗算: (1 + 2x) × (3 + 4x) = 3 + 10x + 8x²
    uint64_t a[] = {1, 2};
    uint64_t b[] = {3, 4};
    uint64_t r[4] = {};
    Ntt::polyMul(r, a, 2, b, 2, prime);

    for (int i = 0; i < 3; ++i)
        std::cout << r[i] << " ";   // 3 10 8
    std::cout << std::endl;

    // 手動 NTT: forward → pointwise multiply → inverse
    size_t N = 4;  // zero-padded power of 2
    std::vector<uint64_t> fa(N, 0), fb(N, 0);
    fa[0] = 1; fa[1] = 2;
    fb[0] = 3; fb[1] = 4;

    Ntt::forward(fa.data(), N, prime);
    Ntt::forward(fb.data(), N, prime);
    Ntt::pointwiseMul(fa.data(), fa.data(), fb.data(), N, prime);
    Ntt::inverse(fa.data(), N, prime);

    for (size_t i = 0; i < 3; ++i)
        std::cout << fa[i] << " ";   // 3 10 8
    std::cout << std::endl;
}

RealFFT<Real> クラス

実数データに特化した静的 FFT API。共役対称性を利用して出力を $N/2+1$ ビンに圧縮。

template<typename Real>
class RealFFT;

fft

static std::vector<std::complex<Real>> fft(std::span<const Real> data);

返り値: サイズ $N/2+1$ の複素数ベクトル。

ifft

static std::vector<Real> ifft(std::span<const std::complex<Real>> spectrum,
                               int output_size);

返り値: サイズ $N$ の実数ベクトル。$1/N$ 正規化済み。

例外: MathError — spectrum.size() != output_size / 2 + 1 の場合。

convolve

static std::vector<Real> convolve(std::span<const Real> a, std::span<const Real> b);

実数データの畳み込み。内部で RealFFT::fft を使用。

FFTEngine<Real> クラス

インスタンスベースの FFT ラッパー。正規化モード (FFTDivMode) を保持。

template<typename Real = double>
class FFTEngine;

コンストラクタ

explicit FFTEngine(int nfft, FFTDivMode mode = FFTDivMode::Forward);

例外: std::invalid_argument — nfft が 2 のべき乗でない場合。

メンバ関数

入力サイズが nfft と一致しない場合は std::invalid_argument。

2D FFT

行方向・列方向に 1D FFT を逐次適用する 2 次元フーリエ変換。Matrix<complex> ベース。

#include <math/fft/fft2d.hpp>

fft2d

template<typename Real = double>
Matrix<std::complex<Real>> fft2d(const Matrix<std::complex<Real>>& input);

2D FFT (行方向 → 列方向)。入出力は同サイズ。

ifft2d

template<typename Real = double>
Matrix<std::complex<Real>> ifft2d(const Matrix<std::complex<Real>>& input);

2D IFFT。$1/N$ 正規化は各軸の FFT::ifft で自動適用。

fft2d_real

template<typename Real = double>
Matrix<std::complex<Real>> fft2d_real(const Matrix<Real>& input);

実数行列を複素行列に変換してから fft2d を適用。

3D FFT

$x \to y \to z$ 方向に 1D FFT を逐次適用する 3 次元フーリエ変換。

#include <math/fft/fft3d.hpp>

Tensor3D<T>

template<typename T>
class Tensor3D;

3 次元テンソルコンテナ。row-major 配置: data[z * ny * nx + y * nx + x]。

fft3d / ifft3d

template<typename Real = double>
Tensor3D<std::complex<Real>> fft3d(const Tensor3D<std::complex<Real>>& input);

template<typename Real = double>
Tensor3D<std::complex<Real>> ifft3d(const Tensor3D<std::complex<Real>>& input);

fft3d_real

template<typename Real = double>
Tensor3D<std::complex<Real>> fft3d_real(const Tensor3D<Real>& input);

fftshift3d

template<typename T>
Tensor3D<T> fftshift3d(const Tensor3D<T>& input);

DC 成分をテンソルの中心に移動。各軸方向に半分ずらす。

ユーティリティ

#include <math/fft/fft_utils.hpp>

スペクトル解析

template<typename Real>
std::vector<Real> amplitude_spectrum(std::span<const std::complex<Real>> spectrum);

template<typename Real>
std::vector<Real> phase_spectrum(std::span<const std::complex<Real>> spectrum);

template<typename Real>
std::vector<Real> unwrap_phase(std::span<const Real> phase,
                                Real tol = std::numbers::pi_v<Real>);

template<typename Real>
std::vector<std::complex<Real>> reconstruct_spectrum(
    std::span<const Real> amplitude, std::span<const Real> phase);

位相アンラップ: phase_spectrum は std::arg() の主値を返すため、 位相応答が $\pm\pi$ でラップする。群遅延 $-d\phi/d\omega$ の計算、minimum-phase 推定、 complex cepstrum などで連続位相が必要な場合は unwrap_phase を適用する。 tol は隣接サンプル間の位相差がこの値を超えたら $2\pi$ の倍数を加減算するしきい値 (既定: $\pi$)。

auto spec  = RealFFT<double>::fft(signal);
auto phase = phase_spectrum<double>(spec);       // 主値 (ラップあり)
auto cont  = unwrap_phase<double>(phase);        // 連続位相
// 群遅延: 差分で近似 -dφ/dω

多項式

template<typename T>
std::vector<T> polynomial_multiply(std::span<const T> a, std::span<const T> b);

template<typename T, typename U>
auto evaluate_polynomial(std::span<const T> poly, const U& x);

polynomial_multiply は FFT::convolve への委譲。evaluate_polynomial はホーナー法。

DCT / IDCT

template<typename Real>
std::vector<Real> dct(std::span<const Real> data);    // DCT-II

template<typename Real>
std::vector<Real> idct(std::span<const Real> coeffs);  // DCT-III (IDCT)

template<typename Real>
std::vector<Real> fast_dct(std::span<const Real> data); // dct の別名

正規化 DCT-II: $X[k] = c(k) \sum_{i=0}^{N-1} x[i] \cos\bigl(\pi(i+0.5)k/N\bigr)$, $c(0) = \sqrt{1/N}$, $c(k) = \sqrt{2/N}$ ($k > 0$)。

$N$ が FFT 対応サイズ ($2^a \cdot 3^b \cdot 5^c$) なら Makhoul 法 $O(N \log N)$。それ以外は直接計算 $O(N^2)$。

DST / IDST

template<typename Real>
std::vector<Real> dst(std::span<const Real> data);     // DST-II

template<typename Real>
std::vector<Real> idst(std::span<const Real> coeffs);   // DST-III (IDST)

正規化 DST-II。現在は直接計算 $O(N^2)$。

Hilbert 変換

template<typename Real>
std::vector<std::complex<Real>> analyticSignal(std::span<const Real> data);

template<typename Real>
std::vector<Real> hilbertTransform(std::span<const Real> data);

template<typename Real>
std::vector<Real> instantaneousAmplitude(std::span<const Real> data);

template<typename Real>
std::vector<Real> instantaneousFrequency(std::span<const Real> data,
                                          Real sampleRate = Real(1));

Walsh-Hadamard 変換

template<typename Real>
std::vector<Real> walshHadamard(std::span<const Real> data);

template<typename Real>
std::vector<Real> iwalshHadamard(std::span<const Real> coeffs);

template<typename Real>
void walshHadamardInplace(std::vector<Real>& data);

高速 Walsh-Hadamard 変換。サイズは 2 のべき乗。$O(N \log N)$ バタフライ。自己逆変換 ($1/N$ 正規化)。

Hartley 変換

template<typename Real>
std::vector<Real> hartley(std::span<const Real> data);

template<typename Real>
std::vector<Real> ihartley(std::span<const Real> coeffs);

$H[k] = \sum_{n=0}^{N-1} x[n] \cdot \mathrm{cas}(2\pi nk/N)$, $\mathrm{cas}(\theta) = \cos\theta + \sin\theta$。 実数入力 → 実数出力。FFT 経由で計算。自己逆変換 ($1/N$ 正規化)。

Hankel 変換

template<typename Real>
struct HankelTransform {
    HankelTransform(int n, Real radius, Real order = Real(0));

    std::vector<Real> r;    // サンプル点
    std::vector<Real> k;    // 周波数点

    std::vector<Real> transform(std::span<const Real> f) const;
    std::vector<Real> inverse(std::span<const Real> F) const;
};

離散 Hankel 変換 (GSL gsl_dht 相当)。任意次数 $\nu$ 対応。$\nu = 0$ は McMahon + $J_0$/$J_1$ 近似で高速パス。

NTT ユーティリティ

template<int P>
std::optional<ModularInt<P>> find_ntt_primitive_root(int n);

サイズ $n$ (2 のべき乗) に適した原始根を探索。見つからなければ std::nullopt。

バッチ FFT

#include <math/fft/fft_batch.hpp>
namespace calx::simd_fft { ... }

同一サイズの複数信号を一括処理する。4 信号インターリーブで AVX2 バタフライを適用。

batchFFT / batchIFFT

void batchFFT(std::complex<float>** signals, int count, int fftSize);
void batchIFFT(std::complex<float>** signals, int count, int fftSize);

4 信号ずつ AVX2 バッチ処理し、端数は通常の simd_fft::fft() で処理。回転子テーブルは thread_local でキャッシュ。

stereoFFT

void stereoFFT(
    const float* interleaved, int fftSize,
    const float* window,
    std::complex<float>* outL,
    std::complex<float>* outR);

ステレオ WAV のインターリーブデータを 1 回の複素 FFT で L/R スペクトルに分離。Hermitian 対称性を利用。

使用例

基本 FFT

#include <math/fft/fft.hpp>
#include <iostream>
#include <complex>
#include <vector>

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

    // 8-point FFT
    std::vector<C> signal = {1, 2, 3, 4, 4, 3, 2, 1};
    calx::FFT<C>::fft(signal);

    std::cout << "FFT result:" << std::endl;
    for (auto& x : signal)
        std::cout << "  " << x << std::endl;

    // Inverse FFT to recover original
    calx::FFT<C>::ifft(signal);

    std::cout << "IFFT result:" << std::endl;
    for (auto& x : signal)
        std::cout << "  " << x.real() << std::endl;
}

スペクトル解析

#include <math/fft/fft.hpp>
#include <math/fft/fft_utils.hpp>
#include <vector>
#include <cmath>
#include <iostream>

int main() {
    constexpr int N = 1024;
    constexpr double fs = 44100.0;
    constexpr double f0 = 440.0;

    // Generate 440 Hz sine wave
    std::vector<double> signal(N);
    for (int i = 0; i < N; ++i)
        signal[i] = std::sin(2.0 * M_PI * f0 * i / fs);

    // Real FFT
    auto spectrum = calx::RealFFT<double>::fft(signal);

    // Amplitude spectrum
    auto amp = calx::amplitude_spectrum<double>(spectrum);

    // Find peak frequency
    int peak = 0;
    for (int k = 1; k < (int)amp.size(); ++k)
        if (amp[k] > amp[peak]) peak = k;

    double peakFreq = peak * fs / N;
    std::cout << "Peak frequency: " << peakFreq << " Hz" << std::endl;
}

畳み込み (多項式乗算)

#include <math/fft/fft.hpp>
#include <iostream>

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

    // p(x) = 1 + 2x + 3x^2
    // q(x) = 4 + 5x
    std::vector<C> p = {1, 2, 3};
    std::vector<C> q = {4, 5};

    auto result = calx::FFT<C>::convolve(p, q);
    // result = {4, 13, 22, 15}  (4 + 13x + 22x^2 + 15x^3)

    for (auto& c : result)
        std::cout << c.real() << " ";
    std::cout << std::endl;
}

2D FFT

#include <math/fft/fft2d.hpp>
#include <math/core/matrix.hpp>
#include <iostream>

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

    // 4x4 matrix
    calx::Matrix<C> img(4, 4);
    for (int r = 0; r < 4; ++r)
        for (int c = 0; c < 4; ++c)
            img(r, c) = C(r + c, 0);

    // Forward 2D FFT
    auto freq = calx::fft2d(img);

    // Inverse 2D FFT
    auto recovered = calx::ifft2d(freq);

    std::cout << "Recovered (0,0): " << recovered(0, 0).real() << std::endl;
}

DCT (離散コサイン変換)

#include <math/fft/fft_utils.hpp>
#include <vector>
#include <iostream>

int main() {
    std::vector<double> data = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};

    auto coeffs = calx::dct<double>(data);
    auto recovered = calx::idct<double>(coeffs);

    std::cout << "DCT coefficients:";
    for (auto v : coeffs) std::cout << " " << v;
    std::cout << std::endl;

    std::cout << "Recovered:";
    for (auto v : recovered) std::cout << " " << v;
    std::cout << std::endl;
}

Hilbert 変換 (包絡線検出)

#include <math/fft/fft_utils.hpp>
#include <vector>
#include <cmath>
#include <iostream>

int main() {
    constexpr int N = 256;

    // AM signal: (1 + 0.5*cos(2pi*5*t)) * cos(2pi*50*t)
    std::vector<double> signal(N);
    for (int i = 0; i < N; ++i) {
        double t = double(i) / N;
        signal[i] = (1.0 + 0.5 * std::cos(2 * M_PI * 5 * t))
                   * std::cos(2 * M_PI * 50 * t);
    }

    // Instantaneous amplitude (envelope)
    auto envelope = calx::instantaneousAmplitude<double>(signal);

    std::cout << "Envelope[0]: " << envelope[0] << std::endl;
}

関連する数学的背景

以下の記事では、FFT モジュールの基盤となる数学的概念を解説している。

• 高速フーリエ変換 — Cooley-Tukey アルゴリズム・計算量解析
• 高速乗算アルゴリズム — NTT・Schönhage-Strassen 法