実際のディジタル信号処理では、信号は離散的にサンプリングされます。離散フーリエ変換（DFT）は、有限長の離散データを周波数領域で解析する手法であり、高速フーリエ変換（FFT）アルゴリズムにより高速に計算できます。
本記事の内容
- DFTの定義と逆DFT
- 連続フーリエ変換との関係
- 周波数分解能とサンプリング定理
- FFTの計算量
- Pythonでの実装（numpy.fft）
前提知識
DFTの定義
$N$ 点の離散データ $x[0], x[1], \dots, x[N-1]$ のDFTは:
$$ \boxed{X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-i2\pi kn/N}, \quad k = 0, 1, \dots, N-1} $$
逆DFT（IDFT）:
$$ \boxed{x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k] \, e^{i2\pi kn/N}, \quad n = 0, 1, \dots, N-1} $$
$W_N = e^{-i2\pi/N}$（回転因子）とおくと、$X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn}$ と簡潔に書けます。
DFTの行列表現
DFTは行列-ベクトル積として表せます。
$$ \bm{X} = \bm{W} \bm{x} $$
ここで $W_{kn} = e^{-i2\pi kn/N}$ はDFT行列です。直接計算では $O(N^2)$ の演算量が必要です。
周波数分解能
サンプリング周波数 $f_s$、データ長 $N$ のとき:
- 周波数分解能: $\Delta f = f_s / N$
- 最大周波数（ナイキスト周波数）: $f_{\max} = f_s / 2$
- $X[k]$ が対応する周波数: $f_k = k \cdot f_s / N$
FFT（高速フーリエ変換）
FFTは分割統治法によりDFTの計算量を $O(N^2)$ から $O(N\log N)$ に削減するアルゴリズムです。
$N$ が2の冪のとき、$x[n]$ を偶数番目と奇数番目に分割:
$$ X[k] = \underbrace{\sum_{m=0}^{N/2-1}x[2m]W_{N/2}^{km}}_{\text{偶数項のDFT}} + W_N^k \underbrace{\sum_{m=0}^{N/2-1}x[2m+1]W_{N/2}^{km}}_{\text{奇数項のDFT}} $$
これを再帰的に適用することで $O(N\log_2 N)$ を達成します。
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
# --- サンプル信号の生成 ---
fs = 1000 # サンプリング周波数 [Hz]
T = 1.0 # 信号長 [s]
N = int(fs * T)
t = np.arange(N) / fs
# 3つの正弦波の合成 + ノイズ
f1, f2, f3 = 50, 120, 300 # Hz
signal = 1.0*np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t) + 0.3*np.sin(2*np.pi*f3*t)
signal += 0.2 * np.random.randn(N) # ノイズ
# --- DFT（numpy.fft） ---
X = np.fft.fft(signal)
freqs = np.fft.fftfreq(N, d=1/fs)
# 正の周波数のみ
positive_mask = freqs >= 0
freqs_pos = freqs[positive_mask]
amplitude = 2 * np.abs(X[positive_mask]) / N # 振幅スペクトル
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
# 時間領域信号
axes[0, 0].plot(t[:200], signal[:200], 'b-', linewidth=0.8)
axes[0, 0].set_title('時間領域信号')
axes[0, 0].set_xlabel('時間 [s]')
axes[0, 0].set_ylabel('振幅')
axes[0, 0].grid(True, alpha=0.3)
# 振幅スペクトル
axes[0, 1].plot(freqs_pos, amplitude, 'r-', linewidth=0.8)
axes[0, 1].set_title('振幅スペクトル（DFT）')
axes[0, 1].set_xlabel('周波数 [Hz]')
axes[0, 1].set_ylabel('振幅')
axes[0, 1].set_xlim(0, 400)
axes[0, 1].grid(True, alpha=0.3)
# パワースペクトル（dB）
power_db = 10 * np.log10(amplitude**2 + 1e-20)
axes[1, 0].plot(freqs_pos, power_db, 'g-', linewidth=0.8)
axes[1, 0].set_title('パワースペクトル [dB]')
axes[1, 0].set_xlabel('周波数 [Hz]')
axes[1, 0].set_ylabel('パワー [dB]')
axes[1, 0].set_xlim(0, 400)
axes[1, 0].grid(True, alpha=0.3)
# スクラッチDFT vs FFT（計算時間比較）
from time import perf_counter
sizes = [2**k for k in range(5, 13)]
times_fft = []
for n in sizes:
x_test = np.random.randn(n)
start = perf_counter()
for _ in range(100):
np.fft.fft(x_test)
times_fft.append((perf_counter() - start) / 100)
axes[1, 1].loglog(sizes, times_fft, 'bo-', label='FFT (numpy)')
axes[1, 1].loglog(sizes, [s**2 * 1e-9 for s in sizes], 'r--', label='$O(N^2)$')
axes[1, 1].loglog(sizes, [s*np.log2(s) * 1e-8 for s in sizes], 'g--', label='$O(N\\log N)$')
axes[1, 1].set_title('FFTの計算量スケーリング')
axes[1, 1].set_xlabel('N')
axes[1, 1].set_ylabel('計算時間 [s]')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('dft_fft.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"検出されたピーク周波数: {freqs_pos[np.argsort(amplitude)[-3:][::-1]]} Hz")
スペクトルから50Hz、120Hz、300Hzの3成分が明確に検出されます。
まとめ
- DFT: $X[k] = \sum x[n]e^{-i2\pi kn/N}$（離散→離散の変換）
- 周波数分解能 $\Delta f = f_s/N$、ナイキスト周波数 $f_s/2$
- FFTにより $O(N^2) \to O(N\log N)$ に高速化
numpy.fft.fftで簡単に実装可能
次のステップとして、以下の記事も参考にしてください。