離散フーリエ変換(DFT)の理論と実装

コンピュータは連続関数を直接扱えません。実際のデータは有限個のサンプル点（離散データ）として記録されます。音声データなら1秒間に44,100個（CD品質）、センサーデータなら1秒間に数百〜数千個のサンプルです。

これらの離散データに対してフーリエ変換を適用するにはどうすればよいでしょうか。連続フーリエ変換の積分を離散データの和に置き換えたものが離散フーリエ変換（Discrete Fourier Transform, DFT）です。

DFTの基本的なアイデアは直感的です。楽器のチューナーを想像してください。ギターの弦を弾いたとき、チューナーは「この音の中にどの周波数がどれだけ含まれているか」を判定して、たとえば「440Hz（A音）にほぼ一致」と表示します。DFTはまさにこれと同じことを数学的に行います。入力データの中に含まれる各周波数成分の強さと位相を計算するのです。

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

この式の意味を噛み砕くと、「データ $x[n]$ と周波数 $k$ の波 $e^{-i2\pi kn/N}$ の相関を計算する」ということです。データの中に周波数 $k$ の成分が多く含まれていれば $|X[k]|$ が大きくなり、含まれていなければ小さくなります。

DFTを理解すると、以下のような応用が開けます。

音声分析: 音声信号のスペクトログラム（周波数-時間表示）
振動解析: 機械の振動データから周波数成分を特定
画像処理: 2次元DFTによる空間周波数フィルタリング
通信工学: OFDM（直交周波数分割多重）方式の基盤
データ圧縮: MP3やJPEGの周波数変換

本記事の内容

連続フーリエ変換から離散化への移行
DFTの数学的定義と逆変換
DFT行列の性質（ユニタリ性、直交性）
周波数分解能とサンプリング定理
Pythonによるスクラッチ実装と応用

前提知識

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

フーリエ変換の定義と性質 — 連続フーリエ変換
フーリエ級数の導出と収束条件 — 周期関数の展開
行列の基本演算 — 行列表現に必要

連続から離散へ — DFTの動機

サンプリングと離散化

連続フーリエ変換から離散フーリエ変換への移行を、一歩ずつ追っていきましょう。この過程を理解することで、DFTが「どこからやってきたのか」が明確になり、DFTの適用範囲と限界もわかるようになります。

連続信号 $f(t)$ をサンプリング周期 $T_s$ で $N$ 点サンプリングしたデータ $x[n] = f(nT_s)$（$n = 0, 1, \ldots, N-1$）を考えます。ここで、サンプリングとは連続信号を一定間隔で「つまみ食い」する操作です。たとえばCDでは $T_s = 1/44100$ 秒（約22.7マイクロ秒）ごとに音声波形の値を記録しています。

連続フーリエ変換 $\hat{f}(\omega) = \int f(t) e^{-i\omega t}\, dt$ の積分を台形則で近似すると

$$ \hat{f}(\omega) \approx T_s \sum_{n=0}^{N-1} x[n]\, e^{-i\omega n T_s} $$

周波数も離散化して $\omega_k = 2\pi k / (NT_s)$（$k = 0, 1, \ldots, N-1$）とおくと

$$ \hat{f}(\omega_k) \approx T_s \sum_{n=0}^{N-1} x[n]\, e^{-i2\pi kn/N} $$

右辺の和（$T_s$ を除く）がまさにDFTです。つまりDFTは、連続フーリエ変換を離散データで近似したものと解釈できます。

この離散化の過程で2つの重要な制約が生まれます。

サンプリング定理（ナイキストの定理）: 元の信号を正しく再現するためには、信号に含まれる最高周波数 $f_{\max}$ の2倍以上のサンプリング周波数 $f_s \geq 2f_{\max}$ が必要です。これを満たさないと「エイリアシング」と呼ばれる周波数の折り返し現象が起きます。たとえば、映画で車のホイールがゆっくり逆回転して見えるのはエイリアシングの身近な例です。フレームレート（サンプリング周波数）がホイールの回転速度（信号の周波数）に対して不十分なために、本来の高い周波数が低い周波数として誤って知覚されるのです。
有限データ長の制約: $N$ 個のデータからは $N$ 個の周波数成分しか得られません。これは周波数分解能 $\Delta f = f_s / N$ に制約を与えます。周波数分解能を高めたければ、より多くのデータ点（より長い時間の信号）が必要です。

なぜ $e^{-i2\pi kn/N}$ なのか

指数項 $e^{-i2\pi kn/N}$ は回転因子（twiddle factor）と呼ばれ、$W_N = e^{-i2\pi/N}$ とおくと $W_N^{kn}$ と書けます。$W_N$ は複素平面上の単位円を $N$ 等分する原始 $N$ 乗根です。

回転因子の幾何学的な意味を考えてみましょう。$e^{-i\theta}$ は複素平面上で単位円上の角度 $\theta$ の点を表します。$W_N^{kn} = e^{-i2\pi kn/N}$ では、$n$ が1増えるごとに角度が $2\pi k/N$ ずつ進みます。つまり、$k = 1$ のときは $N$ ステップで単位円を1周し、$k = 2$ のときは2周します。これは周波数 $k$ の「デジタル版の波」を表しているのです。

$W_N$ の基本的な性質として

周期性: $W_N^{k+N} = W_N^k$（$N$ ステップで元に戻る）
対称性: $W_N^{k+N/2} = -W_N^k$（$N$ が偶数のとき、半周で符号反転）
共役性: $(W_N^k)^* = W_N^{-k} = W_N^{N-k}$（複素共役は逆回転に対応）

があり、これらの性質がFFTアルゴリズムの高速化の鍵となります。特に対称性は、$N$ 点のDFTを $N/2$ 点の2つのDFTに分割する「バタフライ演算」の根拠となっています。

DFTの動機と回転因子の意味を理解したところで、正式な定義に進みましょう。

DFTの数学的定義

定義

長さ $N$ の系列 $\{x[0], x[1], \ldots, x[N-1]\}$ の離散フーリエ変換（DFT）は

$$ \begin{equation} X[k] = \sum_{n=0}^{N-1} x[n]\, W_N^{kn} = \sum_{n=0}^{N-1} x[n]\, e^{-i2\pi kn/N}, \quad k = 0, 1, \ldots, N-1 \end{equation} $$

逆離散フーリエ変換（IDFT）は

$$ \begin{equation} x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k]\, e^{i2\pi kn/N}, \quad n = 0, 1, \ldots, N-1 \end{equation} $$

$X[k]$ は第 $k$ 周波数ビン（bin）の複素振幅であり、$|X[k]|$ が振幅、$\angle X[k]$ が位相を表します。

DFTの式を「内積」として解釈すると理解が深まります。$X[k]$ は、入力信号ベクトル $\bm{x} = (x[0], x[1], \ldots, x[N-1])$ と、周波数 $k$ の基底ベクトル $\bm{w}_k = (1, e^{-i2\pi k/N}, e^{-i2\pi k \cdot 2/N}, \ldots)$ の内積です。内積が大きければ信号と基底が「似ている」（その周波数成分を多く含んでいる）ことを意味し、小さければ「似ていない」ことを意味します。これはまさにパターンマッチングの一種です。

逆DFTは「分析した周波数成分を元に戻す合成」に対応します。各周波数 $k$ の成分 $X[k]$ に基底波 $e^{i2\pi kn/N}$ を掛けて足し合わせることで、元の信号を完全に再構成できます。DFTとIDFTの対（ペア）は、信号の分解と合成の完全な枠組みを提供しているのです。

各ビンの物理的意味

ビン番号 $k$	対応する周波数	物理的意味
$k = 0$	$0$ Hz（DC成分）	信号の平均値 × $N$
$k = 1$	$f_s / N$	最低周波数成分
$k = N/2$	$f_s / 2$（ナイキスト周波数）	最高周波数成分
$k > N/2$	負の周波数に対応	$X[k] = X^*[N-k]$（実数信号の場合）

ここで $f_s = 1/T_s$ はサンプリング周波数です。

$k = 0$ のDC成分は特に重要です。定義式に $k = 0$ を代入すると

$$ X[0] = \sum_{n=0}^{N-1} x[n] \cdot e^0 = \sum_{n=0}^{N-1} x[n] $$

つまり、DC成分は信号の全サンプルの単純な和であり、$X[0]/N$ が信号の平均値に対応します。音声信号であればDC成分は音の「基準レベルのずれ」を表し、通常はゼロに近い値になります。

また、$k > N/2$ のビンが「負の周波数」に対応するという点は注意が必要です。実数信号の場合、$X[k] = X^*[N-k]$ という共役対称性が成り立つため、実質的に独立な周波数情報は $k = 0$ から $k = N/2$ までの $N/2 + 1$ 個しかありません。これが「片側スペクトル」を表示する理由です。

DFTの定義と各ビンの意味を理解したところで、次にDFTを線形代数の枠組みで表現してみましょう。行列表現を導入することで、DFTの構造がより明確になります。

行列表現

DFTは行列とベクトルの積として書けます。

$$ \begin{equation} \bm{X} = \bm{F}_N \bm{x} \end{equation} $$

ここで $\bm{F}_N$ は $N \times N$ のDFT行列です。

$$ (\bm{F}_N)_{kn} = W_N^{kn} = e^{-i2\pi kn/N} $$

$N = 4$ の場合、$W_4 = e^{-i\pi/2} = -i$ なので

$$ \bm{F}_4 = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix} $$

DFT行列の性質

DFT行列は以下の重要な性質を持ちます。

直交性: $\bm{F}_N \bm{F}_N^* = N\bm{I}_N$

ここで $\bm{F}_N^*$ は共役転置です。よって

$$ \bm{F}_N^{-1} = \frac{1}{N}\bm{F}_N^* $$

これは逆DFTの公式 $x[n] = \frac{1}{N}\sum_k X[k] e^{i2\pi kn/N}$ に対応します。

証明: $(\bm{F}_N \bm{F}_N^*)_{kl}$ を計算すると

$$ \sum_{n=0}^{N-1} W_N^{kn} (W_N^{ln})^* = \sum_{n=0}^{N-1} W_N^{(k-l)n} = \sum_{n=0}^{N-1} e^{-i2\pi(k-l)n/N} $$

$k = l$ のとき、全ての項が1なので和は $N$ です。

$k \neq l$ のとき、等比級数の和の公式を使うと

$$ \sum_{n=0}^{N-1} (W_N^{k-l})^n = \frac{1 – W_N^{(k-l)N}}{1 – W_N^{k-l}} = \frac{1 – 1}{1 – W_N^{k-l}} = 0 $$

$W_N^{(k-l)N} = e^{-i2\pi(k-l)} = 1$ を使いました。$\square$

この直交性の証明の核心は等比級数の和の公式です。$k \neq l$ のとき各項が単位円上をぐるぐる回って打ち消し合い、合計がゼロになります。一方 $k = l$ のとき全ての項が $1$ になって合計が $N$ になる — これが直交性の本質です。

直交性の意味を考えると、DFT行列の各行は互いに直交する「基底」であり、DFTはこの直交基底への射影（正射影分解）として理解できます。これは正規直交基底への座標変換と同じ構造であり、逆変換が共役転置で与えられるのも自然なことです。

DFT行列がユニタリ行列（$1/\sqrt{N}$ で正規化した場合）であることは、DFTが「情報を失わない変換」であることを保証しています。時間領域の情報は周波数領域に完全に保存され、いつでも元に戻せるのです。

パーシバルの定理（DFT版）

$$ \begin{equation} \sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N}\sum_{k=0}^{N-1} |X[k]|^2 \end{equation} $$

時間領域のエネルギーと周波数領域のエネルギーが（$1/N$ の定数を除いて）等しいことを表しています。

この定理の物理的な意味は深いです。信号が持つ「パワー」は表現方法（時間で見るか周波数で見るか）によらず一定であるということです。音楽に例えると、CDに録音された曲の総エネルギーは、時間波形で計算しても周波数スペクトルで計算しても同じ値になります。DFTがエネルギーを「作り出したり消したりしない」という意味で、物理法則（エネルギー保存則）と整合しています。

パーシバルの定理はDFTの実装の検証にも使えます。自作のDFT関数が正しく動作しているかどうかを、時間領域と周波数領域のエネルギーを比較することで確認できるのです。後のPython実装で実際にこの検証を行います。

周波数分解能

DFTの周波数分解能 $\Delta f$ は

$$ \Delta f = \frac{f_s}{N} = \frac{1}{NT_s} = \frac{1}{T} $$

ここで $T = NT_s$ はデータの総収録時間です。周波数分解能を高めるには、データの収録時間を長くする（$N$ を増やす、または $T_s$ を大きくする）必要があります。

この周波数分解能の制約は、不確定性原理のフーリエ変換版とも言えます。量子力学の不確定性原理 $\Delta x \cdot \Delta p \geq \hbar/2$ と同様に、信号処理にも $\Delta t \cdot \Delta f \geq 1$ という「時間-周波数の不確定性」が存在します。短い時間窓で信号を切り出すと時間分解能は上がりますが周波数分解能は下がり、逆もまた然りです。

たとえば、楽器の演奏をスペクトル分析する場合を考えましょう。サンプリング周波数 $f_s = 44100$ Hz で $N = 4410$ 点（0.1秒分）を使うと、周波数分解能は $\Delta f = 10$ Hz です。ピアノの中央C（約262Hz）とその半音上のC#（約277Hz）の差は約15Hzなので、これらを分離するには十分です。しかし、$N = 441$ 点（0.01秒分）にすると $\Delta f = 100$ Hz となり、隣接する音を分離できなくなります。

DFTの計算量

DFTの定義式をそのまま実装すると、各 $X[k]$ の計算に $N$ 回の乗算と加算が必要で、$k$ が $N$ 通りあるため、全体の計算量は $O(N^2)$ です。$N = 10^6$（100万点）のデータでは約 $10^{12}$ 回の演算が必要になり、現代のコンピュータでも数秒〜数十秒を要します。

この計算量の壁を打ち破ったのがCooley-Tukey（クーリー・テューキー）の高速フーリエ変換（FFT）アルゴリズム（1965年）です。FFTは回転因子の周期性と対称性を巧みに利用して計算量を $O(N\log N)$ に削減します。$N = 10^6$ の場合、$N^2 = 10^{12}$ が $N\log_2 N \approx 2 \times 10^7$ に削減されるため、約5万倍の高速化が実現します。FFTの登場はデジタル信号処理の実用化に決定的な役割を果たしました。

理論的な準備が整ったところで、DFTをPythonでスクラッチ実装し、動作を確認してみましょう。

Pythonでの実装

DFTのスクラッチ実装

import numpy as np
import matplotlib.pyplot as plt

def dft(x):
    """離散フーリエ変換のスクラッチ実装 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

def idft(X):
    """逆離散フーリエ変換のスクラッチ実装"""
    N = len(X)
    x = np.zeros(N, dtype=complex)
    for n in range(N):
        for k in range(N):
            x[n] += X[k] * np.exp(2j * np.pi * k * n / N)
    return x / N

# テスト信号: 2つの正弦波の重ね合わせ + ノイズ
np.random.seed(42)
N = 128
fs = 128  # サンプリング周波数
t = np.arange(N) / fs
f1, f2 = 10, 25  # 周波数成分
x = 1.0 * np.sin(2*np.pi*f1*t) + 0.5 * np.sin(2*np.pi*f2*t) + 0.3 * np.random.randn(N)

# DFT計算
X_my = dft(x)
X_np = np.fft.fft(x)

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

# 逆変換の検証
x_recovered = idft(X_my)
print(f"IDFT復元誤差: ||x - IDFT(DFT(x))|| = {np.linalg.norm(x - x_recovered.real):.2e}")

# 可視化
freq = np.arange(N) * fs / N

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

# (a) 時間波形
ax = axes[0, 0]
ax.plot(t, x, "b-", linewidth=1)
ax.set_xlabel("Time (s)", fontsize=12)
ax.set_ylabel("Amplitude", fontsize=12)
ax.set_title("Time-domain Signal", fontsize=13)
ax.grid(True, alpha=0.3)

# (b) 振幅スペクトル
ax = axes[0, 1]
amplitude = np.abs(X_my) / N * 2  # 片側スペクトルの正規化
ax.plot(freq[:N//2], amplitude[:N//2], "r-", linewidth=1.5)
ax.set_xlabel("Frequency (Hz)", fontsize=12)
ax.set_ylabel("Amplitude", fontsize=12)
ax.set_title("Amplitude Spectrum |X[k]|", fontsize=13)
ax.grid(True, alpha=0.3)

# ピーク位置をマーク
peaks = np.argsort(amplitude[:N//2])[-2:]
for p in peaks:
    ax.annotate(f"{freq[p]:.0f} Hz", (freq[p], amplitude[p]),
                textcoords="offset points", xytext=(10, 10),
                fontsize=10, fontweight="bold",
                arrowprops=dict(arrowstyle="->", color="black"))

# (c) 位相スペクトル
ax = axes[1, 0]
phase = np.angle(X_my)
ax.plot(freq[:N//2], phase[:N//2], "g.", markersize=3)
ax.set_xlabel("Frequency (Hz)", fontsize=12)
ax.set_ylabel("Phase (rad)", fontsize=12)
ax.set_title("Phase Spectrum", fontsize=13)
ax.grid(True, alpha=0.3)

# (d) パーシバルの定理
ax = axes[1, 1]
energy_time = np.sum(np.abs(x)**2)
energy_freq = np.sum(np.abs(X_my)**2) / N

ax.bar(["Time domain\n$\\sum|x[n]|^2$", "Freq domain\n$\\frac{1}{N}\\sum|X[k]|^2$"],
       [energy_time, energy_freq], color=["blue", "red"], alpha=0.7)
ax.set_ylabel("Energy", fontsize=12)
ax.set_title(f"Parseval: {energy_time:.4f} ≈ {energy_freq:.4f}", fontsize=13)
ax.grid(True, alpha=0.3, axis="y")

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

このグラフから、DFTの基本的な動作が確認できます。

左上（時間波形）: 10Hzと25Hzの正弦波にノイズを加えた合成信号です。時間領域では2つの周波数成分を視覚的に分離するのは困難です。
右上（振幅スペクトル）: DFTにより10Hzと25Hzに明確なピークが検出されています。10Hzのピークが25Hzより高いのは、元の振幅が1.0 vs 0.5であることを正しく反映しています。ノイズ成分はスペクトル全体に均等に分布しています。
左下（位相スペクトル）: 各周波数ビンの位相情報です。信号のピーク周波数以外ではノイズの影響で位相がランダムに分布しています。
右下（パーシバルの定理）: 時間領域と周波数領域のエネルギーが完全に一致しており、DFTがエネルギーを保存することが確認されています。

窓関数とスペクトルリーケージ

import numpy as np
import matplotlib.pyplot as plt

N = 256
fs = 256
t = np.arange(N) / fs
f_signal = 10.5  # 周波数ビンの中間（リーケージが発生）
x = np.sin(2*np.pi*f_signal*t)

# 各種窓関数
windows = {
    "Rectangle": np.ones(N),
    "Hann": np.hanning(N),
    "Hamming": np.hamming(N),
    "Blackman": np.blackman(N)
}

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

for ax, (name, w) in zip(axes.flat, windows.items()):
    xw = x * w
    X = np.fft.fft(xw)
    freq = np.arange(N) * fs / N
    amplitude_db = 20 * np.log10(np.abs(X[:N//2]) / np.max(np.abs(X[:N//2])) + 1e-10)

    ax.plot(freq[:N//2], amplitude_db, linewidth=1.5)
    ax.axvline(f_signal, color="red", linestyle="--", alpha=0.5, label=f"True freq = {f_signal} Hz")
    ax.set_xlabel("Frequency (Hz)", fontsize=11)
    ax.set_ylabel("Amplitude (dB)", fontsize=11)
    ax.set_title(f"{name} Window", fontsize=12)
    ax.set_ylim(-80, 5)
    ax.set_xlim(0, 30)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

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

このグラフから、窓関数がスペクトルリーケージに与える影響が明確に読み取れます。

矩形窓（Rectangle、左上）では、信号の周波数（10.5Hz）がDFTの周波数ビン（整数Hz）に一致しないため、ピークの周辺に大きなサイドローブ（副ローブ）が発生しています。サイドローブのレベルは最大ピークから約-13dBと高く、他の周波数成分があるかのような誤った印象を与えます。これが「スペクトルリーケージ」です。

Hann窓（右上）を適用すると、サイドローブが-31dB程度まで抑制されます。ただし、メインローブ（主ローブ）の幅が矩形窓の約2倍に広がっており、これが周波数分解能の低下に対応します。Hamming窓（左下）はサイドローブを-42dB程度に抑制しますが、最も近いサイドローブが完全にはゼロにならないという特徴があります。Blackman窓（右下）はサイドローブを-58dB以下に抑制しますが、メインローブはさらに広くなります。

窓関数の選択は用途に応じて行います。近接した2つの周波数を分離したい場合は矩形窓（高い周波数分解能）が有利で、微弱な信号成分を強い信号の近くで検出したい場合はBlackman窓（低いサイドローブ）が有利です。実用的には、Hann窓が分解能とサイドローブ抑制のバランスが良く、汎用的に広く使われています。

具体的な計算例

DFTの理論を手計算で確認してみましょう。小さなデータで具体的に計算することで、DFTの仕組みが実感を伴って理解できます。

例1: $N = 4$ の手計算

データ $x = [1, 2, 3, 4]$ の4点DFTを計算します。

$W_4 = e^{-i2\pi/4} = e^{-i\pi/2} = -i$ ですから、$W_4^0 = 1$, $W_4^1 = -i$, $W_4^2 = -1$, $W_4^3 = i$ です。

$k = 0$ の場合：全てのデータの和

$$ X[0] = 1 \cdot 1 + 2 \cdot 1 + 3 \cdot 1 + 4 \cdot 1 = 10 $$

$k = 1$ の場合：

$$ X[1] = 1 \cdot W_4^0 + 2 \cdot W_4^1 + 3 \cdot W_4^2 + 4 \cdot W_4^3 $$

各項を計算すると $1 \cdot 1 + 2 \cdot (-i) + 3 \cdot (-1) + 4 \cdot i = 1 – 2i – 3 + 4i = -2 + 2i$ です。

$k = 2$ の場合：

$$ X[2] = 1 \cdot W_4^0 + 2 \cdot W_4^2 + 3 \cdot W_4^4 + 4 \cdot W_4^6 $$

周期性 $W_4^4 = W_4^0 = 1$, $W_4^6 = W_4^2 = -1$ を使うと、$1 \cdot 1 + 2 \cdot (-1) + 3 \cdot 1 + 4 \cdot (-1) = 1 – 2 + 3 – 4 = -2$ です。

$k = 3$ の場合：同様の計算で $X[3] = -2 – 2i$ です。

結果をまとめると $X = [10, -2+2i, -2, -2-2i]$ です。$X[1]$ と $X[3]$ が互いに共役（$X[3] = X[1]^*$）になっていることに注目してください。これは入力 $x$ が実数であることの帰結です。

パーシバルの定理も確認しましょう。時間領域のエネルギーは $|1|^2 + |2|^2 + |3|^2 + |4|^2 = 1 + 4 + 9 + 16 = 30$ です。周波数領域のエネルギーは $(|10|^2 + |-2+2i|^2 + |-2|^2 + |-2-2i|^2)/4 = (100 + 8 + 4 + 8)/4 = 120/4 = 30$ です。確かに等しくなっています。

例2: 正弦波のDFT

$N = 8$ 点で、周波数ビン $k = 2$ にぴったり一致する正弦波 $x[n] = \sin(2\pi \cdot 2 \cdot n / 8)$ のDFTを考えます。

オイラーの公式より $\sin\theta = (e^{i\theta} – e^{-i\theta})/(2i)$ ですから

$$ x[n] = \frac{1}{2i}\left(e^{i2\pi \cdot 2n/8} – e^{-i2\pi \cdot 2n/8}\right) $$

DFTの直交性を使うと、$X[k]$ の計算は

$$ X[k] = \sum_{n=0}^{7} \frac{1}{2i}\left(e^{i2\pi \cdot 2n/8} – e^{-i2\pi \cdot 2n/8}\right) e^{-i2\pi kn/8} $$

第1項は $k = 2$ のときのみ和が $N = 8$ となり、それ以外はゼロです。第2項は $e^{-i2\pi(k+2)n/8}$ の形をしているため、$k + 2 = 8$ すなわち $k = 6$ のときのみ和が $N = 8$ となります。

よって $X[2] = 8/(2i) = -4i$ で、$X[6] = -8/(2i) = 4i$ です。それ以外の $X[k]$ はすべてゼロです。

振幅スペクトルは $|X[2]|/N \cdot 2 = 4/8 \cdot 2 = 1.0$ となり、元の正弦波の振幅1.0を正しく反映しています。位相は $\angle X[2] = \angle(-4i) = -\pi/2$ であり、$\sin$ 関数が $\cos$ 関数に対して $-\pi/2$ の位相差を持つことと整合します。

例3: 周波数分解能の限界

サンプリング周波数 $f_s = 100$ Hz で $N = 100$ 点（収録時間 $T = 1$ 秒）のデータがあるとします。このとき周波数分解能は $\Delta f = f_s/N = 1$ Hz です。

もし信号が $f_1 = 10$ Hz と $f_2 = 10.5$ Hz の2つの正弦波を含んでいたら、周波数差 $0.5$ Hz は分解能 $1$ Hz よりも小さいため、これら2つのピークをDFTで分離することはできません。スペクトル上では1つの太いピークとして現れてしまいます。

2つのピークを分離するためには $\Delta f \leq 0.5$ Hz が必要で、これには $N \geq f_s / 0.5 = 200$ 点（収録時間 $T \geq 2$ 秒）が必要です。このように、周波数分解能の制約は実際のデータ分析で頻繁に問題となります。

まとめ

本記事では、離散フーリエ変換（DFT）の理論と実装について解説しました。

DFTは連続フーリエ変換を離散化したものであり、$N$ 個のサンプルから $N$ 個の周波数成分を計算する
DFT行列 $\bm{F}_N$ はユニタリ行列（$1/\sqrt{N}$ で正規化した場合）であり、逆変換は共役転置で与えられる
周波数分解能 $\Delta f = f_s/N = 1/T$ はデータの総収録時間の逆数で決まる
パーシバルの定理により、時間領域と周波数領域のエネルギーが保存される
窓関数はスペクトルリーケージを軽減するが、周波数分解能とのトレードオフがある
スクラッチ実装は $O(N^2)$ の計算量であり、次の記事で紹介するFFTにより $O(N\log N)$ に高速化される

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

FFT(高速フーリエ変換)のアルゴリズムと実装 — $O(N\log N)$ の高速化
畳み込み定理とフーリエ変換の関係 — 離散畳み込み
ウェーブレット変換の基礎理論と実装 — 時間-周波数解析