高速フーリエ変換（FFT）の仕組みとPython実装

離散フーリエ変換（DFT）の定義から高速フーリエ変換（FFT）のアルゴリズムまで解説し、NumPy/SciPyを使った実践的な周波数解析のPython実装を紹介します。

4 min read

はじめに

信号処理やデータ解析において、信号がどのような周波数成分から構成されているかを知ることは極めて重要です。この周波数分解を行うのがフーリエ変換であり、離散データに対して適用するのが**離散フーリエ変換（DFT: Discrete Fourier Transform）**です。

しかしDFTを素朴に計算すると計算量が \(O(N^2)\) となり、大規模データへの適用は現実的ではありません。この問題を解決したのが**高速フーリエ変換（FFT: Fast Fourier Transform）**であり、計算量を \(O(N \log N)\) に削減するアルゴリズムです。

本記事では、DFTの数学的定義からFFTアルゴリズムの仕組みを解説し、Pythonによるスクラッチ実装およびNumPy/SciPyを使った実践的な周波数解析を紹介します。

離散フーリエ変換（DFT）

DFTの定義

長さ \(N\) の離散信号 \(x[n]\)（\(n = 0, 1, \ldots, N-1\)）に対して、DFTは次のように定義されます。

\[X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N}, \quad k = 0, 1, \ldots, N-1 \tag{1}\]

ここで、

  • \(X[k]\): 周波数インデックス \(k\) におけるスペクトル（複素数）
  • \(j\): 虚数単位（\(j^2 = -1\)）
  • \(e^{-j2\pi kn/N}\): 回転因子（twiddle factor）

\(X[k]\) の絶対値 \(|X[k]|\) が各周波数成分の振幅を、偏角 \(\angle X[k]\) が位相を表します。

逆DFT（IDFT）

スペクトル \(X[k]\) から元の時間領域信号 \(x[n]\) を復元する逆変換は次のとおりです。

\[x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k] \cdot e^{j2\pi kn/N}, \quad n = 0, 1, \ldots, N-1 \tag{2}\]

DFTと逆DFTは可逆であり、情報の損失なく時間領域と周波数領域を行き来できます。

DFTの計算量

式 \((1)\) を素朴に計算する場合、各 \(k\) について \(N\) 回の複素乗算と加算が必要であり、\(k\) は \(N\) 個あるため、計算量は \(O(N^2)\) です。例えば \(N = 10^6\)（100万サンプル）の場合、\(10^{12}\) 回の演算が必要となり、実用的な速度では計算できません。

FFTアルゴリズム（Cooley-Tukey法）

基本的なアイデア

1965年にCooleyとTukeyが発表したFFTアルゴリズムは、DFTの計算を分割統治法（divide and conquer）で効率化します。基本的なアイデアは、長さ \(N\) のDFTを偶数インデックスと奇数インデックスの2つの長さ \(N/2\) のDFTに分割することです。

入力信号 \(x[n]\) を偶数番目と奇数番目に分けます。

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

すると、DFTは次のように分解できます。

\[X[k] = \underbrace{\sum_{m=0}^{N/2-1} x_{\text{even}}[m] \cdot e^{-j2\pi km/(N/2)}}_{E[k]} + e^{-j2\pi k/N} \underbrace{\sum_{m=0}^{N/2-1} x_{\text{odd}}[m] \cdot e^{-j2\pi km/(N/2)}}_{O[k]} \tag{4}\]

ここで \(E[k]\) と \(O[k]\) はそれぞれ偶数部分列・奇数部分列の長さ \(N/2\) のDFTです。回転因子 \(W_N^k = e^{-j2\pi k/N}\) を用いると、次のように書けます。

\[X[k] = E[k] + W_N^k \cdot O[k], \quad k = 0, 1, \ldots, N/2 - 1 \tag{5}\]\[X[k + N/2] = E[k] - W_N^k \cdot O[k], \quad k = 0, 1, \ldots, N/2 - 1 \tag{6}\]

式 \((6)\) は \(W_N^{k+N/2} = -W_N^k\) という回転因子の対称性から導かれます。

バタフライ演算

式 \((5)\) と \((6)\) の組み合わせはバタフライ演算と呼ばれます。1組の入力 \((E[k], O[k])\) から2つの出力 \((X[k], X[k+N/2])\) を生成する演算で、その信号流れ図の形状が蝶に似ていることから名付けられました。

バタフライ演算の重要な性質は次のとおりです。

  • 各バタフライ演算は1回の複素乗算と2回の複素加減算で済む
  • \(N/2\) 個のバタフライ演算で1段が構成される
  • 再帰的に分割することで全体で \(\log_2 N\) 段になる

計算量の削減

長さ \(N\) のDFTを2つの長さ \(N/2\) のDFTに分割し、これを再帰的に繰り返すと、全体の計算量は次の漸化式に従います。

\[T(N) = 2T(N/2) + O(N) \tag{7}\]

これを解くと \(T(N) = O(N \log N)\) となります。\(N = 10^6\) の場合、\(O(N^2) = 10^{12}\) に対して \(O(N \log N) \approx 2 \times 10^7\) となり、約5万倍の高速化が実現されます。

Pythonによる再帰的FFTの実装

FFTアルゴリズムの理解を深めるため、再帰的な実装をPythonで示します。この実装は \(N\) が2の冪乗であることを前提としています。

import numpy as np

def fft_recursive(x):
    """再帰的なFFT実装（Cooley-Tukey radix-2）"""
    N = len(x)
    if N == 1:
        return x

    # 偶数・奇数インデックスに分割して再帰
    even = fft_recursive(x[0::2])
    odd = fft_recursive(x[1::2])

    # 回転因子を計算
    T = [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)]

    # バタフライ演算
    return [even[k] + T[k] for k in range(N // 2)] + \
           [even[k] - T[k] for k in range(N // 2)]

# 検証: NumPyのFFTと結果を比較
x = np.random.randn(256)
X_custom = np.array(fft_recursive(list(x)))
X_numpy = np.fft.fft(x)
print(f"最大誤差: {np.max(np.abs(X_custom - X_numpy)):.2e}")
# 出力例: 最大誤差: 1.42e-12

この実装はアルゴリズムの理解には有用ですが、Pythonのリスト操作によるオーバーヘッドが大きいため、実用にはNumPyの np.fft.fft を使用してください。

NumPy/SciPyによる実践的な周波数解析

基本的な使い方

NumPyには np.fft モジュールが用意されており、効率的なFFT計算が可能です。

import numpy as np

# 信号の生成
fs = 1000          # サンプリング周波数 [Hz]
t = np.arange(0, 1, 1/fs)  # 1秒間の時間軸
signal = np.sin(2 * np.pi * 50 * t)  # 50Hzの正弦波

# FFTの実行
N = len(signal)
X = np.fft.fft(signal)           # 複素スペクトル
freqs = np.fft.fftfreq(N, 1/fs)  # 周波数軸 [Hz]

np.fft.fftfreq(N, d) は各周波数ビンに対応する周波数値を返します。\(d\) はサンプリング間隔（\(1/f_s\)）です。

複合信号の周波数解析

実際の解析では、複数の周波数成分にノイズが重畳した信号を扱います。以下の例では、50Hzと120Hzの正弦波にガウスノイズを加えた信号のスペクトルを可視化します。

import numpy as np
import matplotlib.pyplot as plt

# --- 信号の生成 ---
fs = 1000  # サンプリング周波数 [Hz]
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)
signal += 0.3 * np.random.randn(len(t))

# --- FFT ---
N = len(signal)
X = np.fft.fft(signal)
freqs = np.fft.fftfreq(N, 1/fs)

# --- 振幅スペクトルの可視化（正の周波数のみ） ---
plt.figure(figsize=(10, 4))
plt.plot(freqs[:N//2], 2/N * np.abs(X[:N//2]))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.title('FFT Amplitude Spectrum')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

振幅スペクトルで 2/N を掛けている理由は、FFTの出力が両側スペクトルであるため、片側のみ表示する場合に2倍する必要があるためです（DC成分とナイキスト成分を除く）。

このスペクトルから、50Hzと120Hzにピークが明確に現れ、120Hzのピークが50Hzの約半分の振幅であることが確認できます。

窓関数の適用

スペクトル漏れ（Spectral Leakage）

FFTは有限長の信号を暗黙的に周期的と仮定します。信号が解析窓の両端で不連続になる場合、本来存在しない周波数成分がスペクトルに現れる現象が起こります。これがスペクトル漏れです。

例えば、サンプリング周波数1000Hzで正確に50.0Hzの信号を1秒間（1000サンプル）取得した場合は漏れが生じませんが、50.5Hzの信号では窓の両端で波形が不連続になり、50.5Hz以外の周波数にもエネルギーが広がります。

窓関数による軽減

窓関数は信号の両端を滑らかにゼロへ減衰させることで、不連続性を軽減します。代表的な窓関数にHann窓があります。

\[w[n] = 0.5 \left(1 - \cos\left(\frac{2\pi n}{N-1}\right)\right) \tag{8}\]
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hann

# スペクトル漏れが発生する信号（50.5Hz）
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 50.5 * t)
N = len(signal)

# 窓関数なし
X_no_win = np.fft.fft(signal)

# Hann窓を適用
window = hann(N)
X_hann = np.fft.fft(signal * window)
freqs = np.fft.fftfreq(N, 1/fs)

# 比較プロット
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(freqs[:N//2], 20 * np.log10(np.abs(X_no_win[:N//2]) + 1e-12))
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Without Window (Rectangular)')
axes[0].set_xlim(0, 200)
axes[0].grid(True, alpha=0.3)

axes[1].plot(freqs[:N//2], 20 * np.log10(np.abs(X_hann[:N//2]) + 1e-12))
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Magnitude [dB]')
axes[1].set_title('With Hann Window')
axes[1].set_xlim(0, 200)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

窓関数を適用すると、メインローブの幅はやや広がりますが、サイドローブが大幅に抑制され、スペクトル漏れが軽減されます。

STFTによる時間-周波数解析

STFTの必要性

通常のFFTは信号全体に対して1つのスペクトルを算出するため、「いつ」その周波数成分が存在したかという時間情報が失われます。音声や振動信号のように周波数が時間とともに変化する信号を解析するには、**短時間フーリエ変換（STFT: Short-Time Fourier Transform）**が必要です。

STFTは信号を短い窓で切り出し、各窓に対してFFTを適用することで、時間-周波数表現（スペクトログラム）を得ます。

\[\text{STFT}\{x[n]\}(m, k) = \sum_{n=0}^{L-1} x[n + mH] \cdot w[n] \cdot e^{-j2\pi kn/L} \tag{9}\]

ここで、\(L\) は窓長、\(H\) はホップサイズ（窓の移動量）、\(w[n]\) は窓関数です。

チャープ信号の解析例

周波数が時間とともに線形に変化するチャープ信号を生成し、スペクトログラムで可視化します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp, spectrogram

# チャープ信号の生成（0Hz → 500Hz、2秒間）
fs = 2000  # サンプリング周波数 [Hz]
t = np.arange(0, 2, 1/fs)
signal = chirp(t, f0=0, f1=500, t1=2, method='linear')

# スペクトログラムの計算
f, t_spec, Sxx = spectrogram(signal, fs=fs, nperseg=256, noverlap=128)

# 可視化
plt.figure(figsize=(10, 5))
plt.pcolormesh(t_spec, f, 10 * np.log10(Sxx + 1e-12), shading='gouraud')
plt.colorbar(label='Power [dB]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.title('Spectrogram of Linear Chirp Signal')
plt.tight_layout()
plt.show()

スペクトログラムでは、周波数が時間とともに線形に増加する様子が斜めの線として明確に可視化されます。nperseg（窓長）を大きくすると周波数分解能が向上しますが時間分解能が低下し、小さくするとその逆になります。これは不確定性原理に基づく本質的なトレードオフです。

まとめ

本記事では、離散フーリエ変換（DFT）の定義からCooley-Tukey FFTアルゴリズムの仕組み、NumPy/SciPyを使った実践的な周波数解析までを解説しました。

  • DFTは信号を周波数成分に分解する基本的な変換であり、計算量は \(O(N^2)\)
  • FFTは分割統治法により計算量を \(O(N \log N)\) に削減するアルゴリズム
  • 窓関数はスペクトル漏れを軽減するために用いる
  • STFTは時間変化する信号の周波数解析に用いる

FFTは信号処理の根幹をなすアルゴリズムであり、フィルタ設計、音声解析、画像処理、通信工学など幅広い分野で活用されています。

関連記事

参考文献

  • Cooley, J. W., & Tukey, J. W. (1965). “An algorithm for the machine calculation of complex Fourier series”. Mathematics of Computation, 19(90), 297-301.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • NumPy FFT documentation
  • SciPy Signal Processing documentation

関連記事