FFT(高速フーリエ変換)のアルゴリズムと実装

あなたのスマートフォンが電話を受けるとき、1秒間に何万もの電波サンプルを周波数成分に分解しています。音楽のストリーミングも、医療画像の撮影も、同じ処理を行っています。これらすべてを支えているのが、今回解説するFFT（高速フーリエ変換）です。

DFTの計算量は $O(N^2)$ です。$N = 10^6$ のデータ（1秒間の44.1kHz音声データに近い）では $10^{12}$ 回の演算が必要であり、現代のコンピュータでも数秒かかります。リアルタイム処理にはとても使えません。たとえば、あなたが電話で話しているとき、音声を20ミリ秒ごとに区切って処理する必要がありますが、1回のDFTに数秒もかかっていては会話が成り立ちません。

1965年、ジェームズ・クーリーとジョン・テューキーは、DFTの計算量を $O(N\log N)$ に削減するアルゴリズムを発表しました。$N = 10^6$ なら約 $2 \times 10^7$ 回の演算で済み、5万倍もの高速化です。このアルゴリズムが高速フーリエ変換（Fast Fourier Transform, FFT）です。

FFTのアイデアを日常生活のたとえで説明しましょう。100人のテスト成績の合計を求めるとき、「全員を1人ずつ足す」方法は $O(N)$ です。しかし、もし「クラス平均」と「クラス人数」が分かっていれば、各クラスの平均×人数を足すだけで済みます。FFTも同じ発想で、$N$ 個のデータを丸ごと処理する代わりに、半分に分割して個別に処理し、最後に組み合わせることで計算量を劇的に削減します。これが「分割統治法」と呼ばれるアルゴリズム設計技法です。

FFTは20世紀で最も重要なアルゴリズムの一つとされ、以下の分野を根本的に変革しました。

• 通信: 携帯電話のOFDM方式、デジタルTV
• 音声処理: リアルタイム音声認識、ノイズキャンセリング
• 科学計算: 高速な畳み込み計算、偏微分方程式の数値解法
• 画像処理: MRI画像再構成、CT画像のフィルタリング

なぜこれほど多くの分野でFFTが使われるのでしょうか。それは、フーリエ変換が「信号を周波数成分に分解する」という普遍的な操作であり、FFTがその計算を現実的な時間で実行可能にしたからです。FFTなくして現代のデジタル技術は成り立たないと言っても過言ではありません。

本記事の内容

• クーリー・テューキーアルゴリズムの導出（分割統治法）
• バタフライ演算の解説
• 再帰的実装と反復的実装
• 計算量の解析と実測ベンチマーク
• Pythonによるスクラッチ実装

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 離散フーリエ変換(DFT)の理論と実装 — DFTの定義と性質
• フーリエ変換の定義と性質 — 連続フーリエ変換
• アルゴリズムの計算量 — $O$記法と分割統治法

FFTのアイデア — 分割統治

DFTの問題点

$N$ 点のDFTを直接計算すると

$$ X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn}, \quad k = 0, \ldots, N-1 $$

各 $k$ について $N$ 回の複素乗算と加算が必要なので、全体で $O(N^2)$ です。

具体的に計算量を見てみましょう。$N = 8$ の場合、DFTの直接計算では各 $k$ について8回の複素乗算が必要で、$k$ は0から7の8通りなので $8 \times 8 = 64$ 回の複素乗算が必要です。$N = 1024$ なら $1024^2 = 1{,}048{,}576$ 回、$N = 10^6$ なら $10^{12}$ 回です。現代のCPUは1秒間に約 $10^9$ 回の浮動小数点演算が可能ですから、$N = 10^6$ のDFTには約1000秒（約17分）かかる計算です。音声のリアルタイム処理では20ミリ秒以内に計算を終えなければなりませんから、これでは全く実用になりません。

偶数・奇数に分割する

FFTの核心的なアイデアは、$N$ 点のDFTを2つの $N/2$ 点のDFTに分割することです。これはちょうど、分厚い本の全ページを検索する代わりに、前半と後半に分けてそれぞれ検索し、結果を統合するようなものです。

入力系列を偶数インデックスと奇数インデックスに分けます。

$$ x_{\text{even}}[m] = x[2m], \quad x_{\text{odd}}[m] = x[2m+1], \quad m = 0, \ldots, N/2 – 1 $$

たとえば $N = 8$ で $x = [x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7]$ なら

$$ x_{\text{even}} = [x_0, x_2, x_4, x_6], \quad x_{\text{odd}} = [x_1, x_3, x_5, x_7] $$

と分けるということです。DFTの和を偶数項（$n = 0, 2, 4, \ldots$）と奇数項（$n = 1, 3, 5, \ldots$）に分けると

$$ X[k] = \sum_{n:\text{even}} x[n] W_N^{nk} + \sum_{n:\text{odd}} x[n] W_N^{nk} $$

偶数項では $n = 2m$、奇数項では $n = 2m+1$ と置き換えると

$$ X[k] = \sum_{m=0}^{N/2-1} x[2m] W_N^{2mk} + \sum_{m=0}^{N/2-1} x[2m+1] W_N^{(2m+1)k} $$

ここで、回転因子 $W_N = e^{-i2\pi/N}$ の重要な性質を使います。$W_N^{2mk}$ の指数部分を整理すると

$$ W_N^{2mk} = \left(e^{-i2\pi/N}\right)^{2mk} = e^{-i2\pi \cdot 2mk/N} = e^{-i2\pi \cdot mk/(N/2)} = W_{N/2}^{mk} $$

つまり $W_N^2 = W_{N/2}$ という関係が成り立ちます。これは「$N$ 点の回転因子を2乗すると $N/2$ 点の回転因子になる」という性質です。

第2項では $W_N^{(2m+1)k} = W_N^{2mk} \cdot W_N^k = W_{N/2}^{mk} \cdot W_N^k$ なので、共通因子 $W_N^k$ をくくり出すことができます。

$$ X[k] = \sum_{m=0}^{N/2-1} x[2m] W_{N/2}^{mk} + W_N^k \sum_{m=0}^{N/2-1} x[2m+1] W_{N/2}^{mk} $$

ここで $E[k]$ と $O[k]$ を導入します。

$$ X[k] = \underbrace{\sum_{m=0}^{N/2-1} x_{\text{even}}[m] W_{N/2}^{mk}}_{E[k]} + W_N^k \underbrace{\sum_{m=0}^{N/2-1} x_{\text{odd}}[m] W_{N/2}^{mk}}_{O[k]} $$

$$ \begin{equation} X[k] = E[k] + W_N^k \cdot O[k] \end{equation} $$

$E[k]$ と $O[k]$ はそれぞれ $N/2$ 点のDFTです。ここでの $W_N^k$ は「回転因子」（twiddle factor）と呼ばれ、偶数部分と奇数部分の結果を組み合わせる際の位相調整の役割を果たします。直感的には、奇数インデックスの信号は偶数インデックスに対して1サンプルだけずれているため、その位相のずれを $W_N^k$ で補正しているのです。

周期性の活用

$E[k]$ と $O[k]$ は $N/2$ 点のDFTですから、周期 $N/2$ で周期的です。つまり

$$ E[k + N/2] = E[k], \quad O[k + N/2] = O[k] $$

が成り立ちます。これはDFTの周期性（$X[k + N] = X[k]$）そのものです。

さらに、回転因子には重要な半周期性があります。$W_N^{N/2} = e^{-i2\pi \cdot (N/2)/N} = e^{-i\pi} = -1$ なので

$$ W_N^{k+N/2} = W_N^k \cdot W_N^{N/2} = W_N^k \cdot (-1) = -W_N^k $$

この性質を使って $X[k + N/2]$ を書き下すと

$$ X[k + N/2] = E[k + N/2] + W_N^{k+N/2} \cdot O[k + N/2] $$

$E$ と $O$ の周期性、および $W_N^{k+N/2} = -W_N^k$ を代入すると

$$ \begin{equation} X[k + N/2] = E[k] – W_N^k \cdot O[k] \end{equation} $$

つまり、$E[k]$ と $O[k]$ を計算すれば、$X[k]$ と $X[k + N/2]$ が1回の加算と1回の乗算で同時に得られます。これがバタフライ演算です。

具体例として $N = 8$、$k = 1$ の場合を見てみましょう。$E[1]$ と $O[1]$ が計算済みだとすると

$$ X[1] = E[1] + W_8^1 \cdot O[1], \quad X[5] = E[1] – W_8^1 \cdot O[1] $$

ここで $W_8^1 = e^{-i2\pi/8} = e^{-i\pi/4} = \frac{\sqrt{2}}{2}(1 – i)$ です。仮に $E[1] = 3 + 2i$、$O[1] = 1 – i$ とすると

$$ W_8^1 \cdot O[1] = \frac{\sqrt{2}}{2}(1-i)(1-i) = \frac{\sqrt{2}}{2}(1 – 2i + i^2) = \frac{\sqrt{2}}{2}(-2i) = -\sqrt{2}\,i $$

したがって

$$ X[1] = (3 + 2i) + (-\sqrt{2}\,i) = 3 + (2-\sqrt{2})i \approx 3 + 0.586i $$

$$ X[5] = (3 + 2i) – (-\sqrt{2}\,i) = 3 + (2+\sqrt{2})i \approx 3 + 3.414i $$

このように、$W_N^k O[k]$ を1回だけ計算すれば、加算と減算だけで $X[k]$ と $X[k+N/2]$ の両方が求まります。これがFFTの効率性の源泉です。

計算量の解析

$T(N)$ を $N$ 点FFTの計算量（複素乗算の回数）とすると

$$ T(N) = 2T(N/2) + O(N) $$

この式は「$N$ 点FFT = $N/2$ 点FFTを2つ + $N/2$ 回のバタフライ演算」という構造を反映しています。この漸化式を展開してみましょう。

$$ T(N) = 2T(N/2) + cN $$

$T(N/2) = 2T(N/4) + c(N/2)$ を代入すると

$$ T(N) = 2[2T(N/4) + cN/2] + cN = 4T(N/4) + 2cN $$

さらに $T(N/4) = 2T(N/8) + c(N/4)$ を代入すると

$$ T(N) = 8T(N/8) + 3cN $$

$k$ 回展開すると

$$ T(N) = 2^k T(N/2^k) + kcN $$

$N = 2^p$ とすると、$k = p = \log_2 N$ で $T(1) = O(1)$ に到達するので

$$ T(N) = N \cdot O(1) + cN\log_2 N = O(N\log N) $$

マスター定理を使えば一発で同じ結論が出ますが、展開して確認することで「各段で $O(N)$ の仕事を $\log N$ 段繰り返す」という構造が見えます。

具体的な数値で比較してみましょう。

$N$ が大きくなるほど、FFTの優位性が指数的に増します。$N = 10^6$ での5万倍の高速化は、「17分かかる計算が0.02秒で終わる」ことを意味します。

FFTのアイデアを理解したところで、次にバタフライ演算の具体的な構造と、それを実装する方法を見ていきましょう。

バタフライ演算

バタフライの構造

$k = 0, \ldots, N/2 – 1$ の各 $k$ に対して

$$ \begin{align} X[k] &= E[k] + W_N^k \cdot O[k] \\ X[k + N/2] &= E[k] – W_N^k \cdot O[k] \end{align} $$

この2つの式を図で表すと、蝶の羽のような形になるため「バタフライ」と呼ばれます。

入力が $E[k]$ と $O[k]$、出力が $X[k]$ と $X[k+N/2]$ で、「回転因子」$W_N^k$ が掛かる加減算の対です。バタフライ演算のポイントは、1回の複素乗算（$W_N^k \cdot O[k]$）と2回の加減算だけで、2つの出力を同時に得られることです。

バタフライ演算を擬似コードで書くと、以下のようにシンプルです。

temp = W * O[k]
X[k]       = E[k] + temp
X[k + N/2] = E[k] - temp

temp を一時変数に保持するため、乗算は1回で済みます。もし $X[k]$ と $X[k+N/2]$ を個別に計算していたら2回の乗算が必要でしたが、バタフライ演算では1回で済むのです。

回転因子の幾何学的意味

回転因子 $W_N^k = e^{-i2\pi k/N}$ は複素平面上の単位円上の点です。$k = 0$ のとき $W_N^0 = 1$（実軸上の点）、$k$ が増えるにつれて時計回りに等間隔で回転していきます。

$N = 8$ の場合、$W_8^k$ は単位円を8等分する点に対応します。

$$ W_8^0 = 1, \quad W_8^1 = \frac{\sqrt{2}}{2}(1-i), \quad W_8^2 = -i, \quad W_8^3 = \frac{\sqrt{2}}{2}(-1-i) $$

$$ W_8^4 = -1, \quad W_8^5 = \frac{\sqrt{2}}{2}(-1+i), \quad W_8^6 = i, \quad W_8^7 = \frac{\sqrt{2}}{2}(1+i) $$

このうち $W_8^4 = -W_8^0$、$W_8^5 = -W_8^1$、$W_8^6 = -W_8^2$、$W_8^7 = -W_8^3$ です。つまり、後半の回転因子は前半の回転因子の符号反転であり、これがバタフライの「加算と減算」の構造に対応しています。

全体の構造

$N = 8$ の場合、FFTは3段階（$\log_2 8 = 3$）のバタフライ演算で構成されます。

• 第1段: 4つの2点FFT（$N/4$ 個のバタフライ）
• 第2段: 2つの4点FFT（$N/2$ 個のバタフライ）
• 第3段: 1つの8点FFT（$N/2$ 個のバタフライ）

各段で $N/2$ 個のバタフライ演算が行われ、$\log_2 N$ 段あるので、総計算量は $O(N\log N)$ です。

第1段では、2つの要素を1組としてバタフライ演算を行います。2点FFTは最もシンプルで、回転因子は $W_2^0 = 1$ のみです。

$$ X[0] = x[0] + x[1], \quad X[1] = x[0] – x[1] $$

これは単なる加算と減算であり、乗算すら不要です。第2段では、第1段の結果を2つずつ組み合わせて4点FFTを構成します。第3段では、さらに2つの4点FFTの結果を組み合わせて最終的な8点FFTを完成させます。

このように、バタフライ演算は「小さなFFTの結果を統合して大きなFFTを構成する」という分割統治法のボトムアップ過程そのものです。次に、この構造をPythonで実装していきましょう。

ビットリバーサル並べ替え

反復的FFTでは、入力データをあらかじめ「ビットリバーサル」順に並べ替えておく必要があります。ビットリバーサルとは、インデックスの2進表現を左右反転させる操作です。

$N = 8$ の場合、インデックスは3ビットで表現されます。

つまり、入力配列 $[x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7]$ は $[x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7]$ に並べ替えられます。

なぜこの並べ替えが必要なのでしょうか。再帰的FFTでは、偶数インデックスと奇数インデックスに再帰的に分割していきます。最初の分割で $[x_0, x_2, x_4, x_6]$ と $[x_1, x_3, x_5, x_7]$ に分かれ、さらにそれぞれを偶数・奇数に分割すると $[x_0, x_4]$、$[x_2, x_6]$、$[x_1, x_5]$、$[x_3, x_7]$ になります。この最終的な順序がビットリバーサル順に一致するのです。反復的実装では、この並べ替えを最初にまとめて行うことで、以降はインプレースでバタフライ演算を順番に適用できます。

Pythonでの実装

再帰的FFT

まず、分割統治法をそのまま再帰で実装します。

import numpy as np
import matplotlib.pyplot as plt
import time

def fft_recursive(x):
    """再帰的FFT（クーリー・テューキー、基数2）"""
    N = len(x)
    if N == 1:
        return x.copy()

    # 偶数・奇数インデックスに分割
    E = fft_recursive(x[0::2])  # 偶数インデックス
    O = fft_recursive(x[1::2])  # 奇数インデックス

    # バタフライ演算
    X = np.zeros(N, dtype=complex)
    for k in range(N // 2):
        W = np.exp(-2j * np.pi * k / N)
        X[k] = E[k] + W * O[k]
        X[k + N//2] = E[k] - W * O[k]

    return X

# テスト
np.random.seed(42)
N = 1024
x = np.random.randn(N)

X_my = fft_recursive(x)
X_np = np.fft.fft(x)

print(f"再帰FFT vs NumPy: ||差|| = {np.linalg.norm(X_my - X_np):.2e}")

この再帰的実装はアルゴリズムの構造を直接反映しており理解しやすいですが、再帰呼び出しのオーバーヘッドとメモリの配列コピーが発生します。出力の ||差|| は $10^{-10}$ 程度（機械イプシロン $\approx 2.2 \times 10^{-16}$ の数百倍程度）であり、丸め誤差の範囲内で NumPy の FFT と完全に一致していることが確認できます。

再帰的実装のコードを振り返ると、アルゴリズムの骨格は非常にシンプルです。

1. 基底ケース（$N = 1$）ではそのまま返す
2. 偶数インデックス x[0::2] と奇数インデックス x[1::2] にスライスして再帰呼び出し
3. 戻ってきた $E$ と $O$ に対してバタフライ演算を適用

再帰の深さは $\log_2 N$ です。$N = 1024$ なら10段の再帰で到達します。しかし、各段で新しい配列が作られるため、メモリ効率は良くありません。実用的には、次に紹介する反復的実装が使われます。

反復的FFT（インプレース）

実用的にはメモリ効率の良い反復的実装が使われます。

import numpy as np

def fft_iterative(x):
    """反復的FFT（ビットリバーサル + インプレースバタフライ）"""
    N = len(x)
    log2N = int(np.log2(N))
    assert N == 2**log2N, "Nは2のべき乗でなければなりません"

    # ビットリバーサル並べ替え
    X = np.array(x, dtype=complex)
    for i in range(N):
        j = int(bin(i)[2:].zfill(log2N)[::-1], 2)
        if i < j:
            X[i], X[j] = X[j], X[i]

    # バタフライ演算（ボトムアップ）
    for stage in range(log2N):
        half_size = 2**stage
        full_size = 2**(stage + 1)
        W_base = np.exp(-2j * np.pi / full_size)

        for group_start in range(0, N, full_size):
            W = 1.0
            for k in range(half_size):
                idx_top = group_start + k
                idx_bot = group_start + k + half_size
                temp = W * X[idx_bot]
                X[idx_bot] = X[idx_top] - temp
                X[idx_top] = X[idx_top] + temp
                W *= W_base

    return X

# テスト
X_iter = fft_iterative(x)
print(f"反復FFT vs NumPy: ||差|| = {np.linalg.norm(X_iter - X_np):.2e}")

反復的実装では、まず入力データを「ビットリバーサル」順に並べ替えてから、ボトムアップでバタフライ演算を行います。追加のメモリ割り当てが不要（インプレース）であり、再帰よりも効率的です。

コードの各部分を詳しく見ていきましょう。

ビットリバーサル部分（7〜10行目）: インデックス $i$ の2進表現を bin(i) で求め、zfill(log2N) でゼロパディングし、[::-1] で反転させます。反転後のインデックス $j$ と $i$ が異なり、かつ $i < j$ の場合にのみ交換します（$i < j$ の条件により、同じペアの2回交換を防ぎます）。

バタフライ演算部分（12〜21行目）: stage は現在の段（0から $\log_2 N – 1$ まで）を表します。各段では

• half_size = 2^stage: バタフライの半分のサイズ
• full_size = 2^(stage+1): バタフライグループの全体サイズ
• W_base = e^{-2\pi i / \text{full\_size}}: この段の基本回転因子

を設定し、各グループ内でバタフライ演算をインプレースで実行します。temp = W * X[idx_bot] を計算し、上端と下端を同時に更新するのが、先ほど説明したバタフライ演算の実装そのものです。

再帰版と反復版の出力を比較すると、どちらもNumPyのFFTと一致することが確認できます。反復版は再帰のスタックオーバーヘッドがないため、大きな $N$ でもメモリ効率よく動作します。

次に、ナイーブDFTとFFTの実行時間を実際に測定し、理論通りの高速化が得られるか確認しましょう。

計算量のベンチマーク

import numpy as np
import matplotlib.pyplot as plt
import time

def dft_naive(x):
    """ナイーブDFT O(N²)"""
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X

sizes = [2**k for k in range(4, 14)]
times_naive = []
times_fft = []
times_numpy = []

for N in sizes:
    x = np.random.randn(N)

    if N <= 4096:
        t0 = time.perf_counter()
        dft_naive(x)
        times_naive.append(time.perf_counter() - t0)
    else:
        times_naive.append(np.nan)

    t0 = time.perf_counter()
    fft_iterative(x)
    times_fft.append(time.perf_counter() - t0)

    t0 = time.perf_counter()
    np.fft.fft(x)
    times_numpy.append(time.perf_counter() - t0)

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

ax = axes[0]
ax.loglog(sizes, times_naive, "ro-", linewidth=2, markersize=6, label="Naive DFT $O(N^2)$")
ax.loglog(sizes, times_fft, "bs-", linewidth=2, markersize=6, label="My FFT $O(N\\log N)$")
ax.loglog(sizes, times_numpy, "g^-", linewidth=2, markersize=6, label="NumPy FFT")

n_arr = np.array(sizes, dtype=float)
ax.loglog(n_arr, 1e-8 * n_arr**2, "r:", alpha=0.3, label="$O(N^2)$ ref")
ax.loglog(n_arr, 1e-7 * n_arr * np.log2(n_arr), "b:", alpha=0.3, label="$O(N\\log N)$ ref")
ax.set_xlabel("N", fontsize=12)
ax.set_ylabel("Time (sec)", fontsize=12)
ax.set_title("DFT vs FFT: Computation Time", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")

ax = axes[1]
speedup = [tn/tf if not np.isnan(tn) else np.nan
           for tn, tf in zip(times_naive, times_fft)]
valid_sizes = [s for s, sp in zip(sizes, speedup) if not np.isnan(sp)]
valid_speedup = [sp for sp in speedup if not np.isnan(sp)]
ax.bar(range(len(valid_sizes)), valid_speedup, color="steelblue", alpha=0.8)
ax.set_xticks(range(len(valid_sizes)))
ax.set_xticklabels([str(s) for s in valid_sizes], rotation=45)
ax.set_xlabel("N", fontsize=12)
ax.set_ylabel("Speedup (DFT/FFT)", fontsize=12)
ax.set_title("FFT Speedup Factor", fontsize=13)
ax.grid(True, alpha=0.3, axis="y")

for i, sp in enumerate(valid_speedup):
    ax.text(i, sp + 0.5, f"{sp:.0f}x", ha="center", fontsize=9, fontweight="bold")

plt.tight_layout()
plt.savefig("fft_benchmark.png", dpi=150, bbox_inches="tight")
plt.show()

このベンチマーク結果から、FFTの劇的な高速化効果が確認できます。

1. 左図（計算時間）: ナイーブDFT（赤丸）は $O(N^2)$ に沿って急激に増大しますが、FFT（青四角）は $O(N\log N)$ に沿って緩やかに増加しています。$N = 4096$ ではナイーブDFTに約1秒かかるのに対し、FFTは0.1秒未満です。NumPyのFFT（緑三角）はCで最適化されているためさらに数百倍高速です。両対数グラフ上で傾きを確認すると、ナイーブDFTは傾き2（$N^2$ に対応）、FFTは傾き約1（$N\log N$ は対数グラフ上でほぼ線形に見える）であることがわかります。

2. 右図（速度比）: $N$ が増えるほどFFTの優位性が顕著になります。$N = 4096$ ではFFTがDFTの数百倍高速であり、$N = 10^6$ レベルでは数万倍の高速化が期待されます。この高速化率は理論値 $N / \log_2 N$ にほぼ一致します。$N = 4096$ では $4096/12 \approx 341$ 倍です。

3. NumPyとの比較: NumPyの np.fft.fft はCやFortranで書かれたFFTWライブラリを利用しており、キャッシュ効率やSIMD命令の活用により、Pythonで書いた同じアルゴリズムよりも数百倍高速です。これはアルゴリズム自体の違いではなく、実装レベルの最適化の差です。

具体的な数値計算例 — $N = 8$ のFFT

理論の理解を深めるために、$N = 8$ の小さな例でFFTの全過程を手計算で追ってみましょう。

入力データを $x = [1, 2, 3, 4, 4, 3, 2, 1]$ とします。この信号は左右対称であり、DFT結果も特徴的な性質を持つことが期待されます。

ステップ1: ビットリバーサル並べ替え

ビットリバーサル順は $[0, 4, 2, 6, 1, 5, 3, 7]$ でした。並べ替え後のデータは

$$ x_{\text{BR}} = [x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7] = [1, 4, 3, 2, 2, 3, 4, 1] $$

ステップ2: 第1段（2点FFT × 4組）

回転因子は $W_2^0 = 1$ のみです。各ペアに対して加算と減算を行います。

• ペア $(x_0, x_4) = (1, 4)$: $1 + 4 = 5$、$1 – 4 = -3$
• ペア $(x_2, x_6) = (3, 2)$: $3 + 2 = 5$、$3 – 2 = 1$
• ペア $(x_1, x_5) = (2, 3)$: $2 + 3 = 5$、$2 – 3 = -1$
• ペア $(x_3, x_7) = (4, 1)$: $4 + 1 = 5$、$4 – 1 = 3$

第1段後: $[5, -3, 5, 1, 5, -1, 5, 3]$

ステップ3: 第2段（4点FFT × 2組）

回転因子は $W_4^0 = 1$ と $W_4^1 = e^{-i\pi/2} = -i$ です。

前半グループ $[5, -3, 5, 1]$:

• $k=0$: $5 + 1 \cdot 5 = 10$、$5 – 1 \cdot 5 = 0$
• $k=1$: $-3 + (-i) \cdot 1 = -3 – i$、$-3 – (-i) \cdot 1 = -3 + i$

第2段前半後: $[10, -3-i, 0, -3+i]$

後半グループ $[5, -1, 5, 3]$:

• $k=0$: $5 + 1 \cdot 5 = 10$、$5 – 1 \cdot 5 = 0$
• $k=1$: $-1 + (-i) \cdot 3 = -1 – 3i$、$-1 – (-i) \cdot 3 = -1 + 3i$

第2段後半後: $[10, -1-3i, 0, -1+3i]$

ステップ4: 第3段（8点FFT × 1組）

回転因子は $W_8^k$（$k = 0, 1, 2, 3$）です。$W_8^0 = 1$、$W_8^1 = \frac{\sqrt{2}}{2}(1-i)$、$W_8^2 = -i$、$W_8^3 = \frac{\sqrt{2}}{2}(-1-i)$ を使います。

• $k=0$: $X[0] = 10 + 1 \cdot 10 = 20$、$X[4] = 10 – 1 \cdot 10 = 0$
• $k=1$: $W_8^1 \cdot (-1-3i) = \frac{\sqrt{2}}{2}(1-i)(-1-3i) = \frac{\sqrt{2}}{2}(-1-3i+i+3i^2) = \frac{\sqrt{2}}{2}(-4-2i) = -2\sqrt{2} – \sqrt{2}i$

$$ X[1] = (-3-i) + (-2\sqrt{2} – \sqrt{2}i) = -3-2\sqrt{2} + (-1-\sqrt{2})i $$

$$ X[5] = (-3-i) – (-2\sqrt{2} – \sqrt{2}i) = -3+2\sqrt{2} + (-1+\sqrt{2})i $$

数値的に $X[1] \approx -5.828 – 2.414i$、$X[5] \approx -0.172 + 0.414i$ です。

この手計算の結果は NumPy の np.fft.fft([1,2,3,4,4,3,2,1]) と一致します。入力が実数かつ対称なので、$X[0] = 20$（全要素の和）、$X[4] = 0$ などの特徴的な値が現れています。

FFTの応用例

FFTは単に周波数スペクトルを求めるだけでなく、さまざまな応用があります。ここでは代表的な2つの応用を紹介します。

応用1: FFTによる高速畳み込み

2つの信号 $x[n]$（長さ $M$）と $h[n]$（長さ $L$）の線形畳み込みは

$$ y[n] = \sum_{m=0}^{M-1} x[m] h[n-m], \quad n = 0, \ldots, M + L – 2 $$

直接計算すると $O(ML)$ の計算量ですが、畳み込み定理を使うと

$$ y = \text{IFFT}(\text{FFT}(x) \cdot \text{FFT}(h)) $$

と書け、FFTを3回（順変換2回 + 逆変換1回）適用するだけで済みます。$N = M + L – 1$ として $O(N\log N)$ の計算量です。

これは信号処理でFIRフィルタを適用する際に多用される手法です。フィルタ長が長い場合、直接畳み込みよりFFTベースの方法が圧倒的に高速です。たとえば、音声信号（$M = 44100$）にリバーブ効果のインパルス応答（$L = 88200$）を畳み込む場合、直接計算では $44100 \times 88200 \approx 3.9 \times 10^9$ 回の乗算が必要ですが、FFTベースでは $N \approx 131072 = 2^{17}$ として $3 \times 131072 \times 17 \approx 6.7 \times 10^6$ 回で済み、約580倍の高速化です。

応用2: スペクトル解析

FFTを使った最も基本的な応用がスペクトル解析です。時間領域の信号から周波数成分を抽出し、どの周波数がどれだけのエネルギーを持つかを調べます。

パワースペクトル密度（PSD）は、FFT結果の絶対値の2乗で定義されます。

$$ P[k] = \frac{1}{N} |X[k]|^2 $$

ここで $X[k]$ は入力信号のFFTです。$k$ 番目の周波数 $f_k$ は

$$ f_k = \frac{k \cdot f_s}{N} $$

で与えられます。$f_s$ はサンプリング周波数です。

たとえば、$f_s = 1000$ Hz でサンプリングした $N = 1024$ 点のデータに対してFFTを適用すると、周波数分解能は $\Delta f = f_s / N = 1000/1024 \approx 0.977$ Hz となります。つまり、約1Hz刻みで周波数成分を識別できます。ナイキスト周波数 $f_s/2 = 500$ Hz までの成分を分析可能です。

数値精度に関する注意

FFTの計算では浮動小数点演算を繰り返すため、丸め誤差が蓄積します。ここでは、その影響と対策について述べます。

$N$ 点FFTは $\log_2 N$ 段のバタフライ演算で構成されます。各段で加減算と複素乗算が行われるため、1段あたり約1ULP（Unit in the Last Place）の丸め誤差が発生する可能性があります。したがって、FFT全体の相対誤差は

$$ \epsilon_{\text{FFT}} \approx O(\log_2 N) \cdot \epsilon_{\text{machine}} $$

のオーダーとなります。倍精度（64ビット）では $\epsilon_{\text{machine}} \approx 2.2 \times 10^{-16}$ なので、$N = 2^{20}$ でも相対誤差は $20 \times 2.2 \times 10^{-16} \approx 4.4 \times 10^{-15}$ 程度に収まります。これは実用上問題にならないレベルです。

一方、ナイーブDFTでは $N$ 回の加算を行うため、相対誤差は $O(N) \cdot \epsilon_{\text{machine}}$ となります。$N = 10^6$ では $10^{-10}$ 程度の誤差が生じ得ます。つまり、FFTは高速なだけでなく、数値精度の面でもナイーブDFTより優れているのです。

ただし、回転因子の計算精度が重要です。$W_N^k = e^{-i2\pi k/N}$ を都度 np.exp で計算するか、$W = W \cdot W_{\text{base}}$ と逐次的に掛け算で更新するかで精度に差が出ます。逐次乗算方式は高速ですが、$k$ が大きくなると丸め誤差が蓄積します。高精度が求められる用途では、各 $k$ ごとに直接 np.exp(-2j * np.pi * k / N) で計算する方が安全です。

まとめ

本記事では、FFT（高速フーリエ変換）のアルゴリズムと実装について、分割統治法の基本アイデアから具体的な数値計算例、Pythonによるスクラッチ実装、そして応用例まで一貫して解説しました。

• FFTはDFTの $O(N^2)$ を $O(N\log N)$ に高速化するアルゴリズムであり、分割統治法に基づく
• クーリー・テューキーアルゴリズムは入力を偶数・奇数インデックスに分割し、$N/2$ 点のDFT 2つに帰着させる
• バタフライ演算 $X[k] = E[k] + W_N^k O[k]$、$X[k+N/2] = E[k] – W_N^k O[k]$ が基本演算単位であり、回転因子の半周期性 $W_N^{k+N/2} = -W_N^k$ を利用する
• ビットリバーサルは再帰的な偶数・奇数分割を一括して行う前処理であり、反復的実装の鍵となる
• 反復的実装はビットリバーサルとインプレースのバタフライ演算により、メモリ効率が良い
• $N = 10^6$ では約5万倍の高速化が実現され、リアルタイム信号処理を可能にする
• FFTは高速畳み込みやスペクトル解析など多くの応用を持ち、通信・音声・画像・科学計算の基盤技術である
• FFTはナイーブDFTに比べて数値精度も優れており、相対誤差は $O(\log N) \cdot \epsilon_{\text{machine}}$ に抑えられる

本記事で紹介したのは基数2のクーリー・テューキーアルゴリズムですが、FFTにはさまざまな変種があります。基数4や分割基数アルゴリズムはバタフライ当たりの乗算数を削減し、Bluesteinアルゴリズムは $N$ が2のべき乗でなくても適用できます。また、実数入力に特化したRFFT（Real FFT）は対称性を利用して計算量を半減させます。これらの発展的トピックについては、今後の記事で取り上げる予定です。

次のステップとして、以下の記事も参考にしてください。

• 畳み込み定理とフーリエ変換の関係 — FFTによる高速畳み込み
• ウェーブレット変換の基礎理論と実装 — 時間-周波数解析
• 離散フーリエ変換(DFT)の理論と実装 — DFTの基本理論