離散コサイン変換(DCT)の理論とPython実装:JPEG・音声圧縮の基礎

scipy.fft.dct・scipy.fft.idct・scipy.fftpack.dct・scipy.fft.dctn で離散コサイン変換(DCT)の理論・Python 実装・JPEG/MP3 圧縮応用を体系解説。DCT-I/II/III/IV と MDCT(変形 DCT)の使い分け、直交変換、エネルギー集中(energy compaction)、FFT を介した高速計算、JPEG 8×8 ブロック量子化、MP3/AAC のフレーム重畳までカバー。

はじめに

デジタルカメラで撮影した写真をJPEGで保存するとき、スマートフォンでMP3を再生するとき、あるいはH.264でエンコードされた動画を視聴するとき——これらすべての背後に**離散コサイン変換(DCT: Discrete Cosine Transform)**が動いています。

DCTは1974年にNasir AhmedらによってIEEE Transactions on Computersに発表された変換であり、信号をコサイン基底関数の重みつき和として表現します。 DFT(離散フーリエ変換) と密接に関係しながらも、実数演算のみで完結し、エネルギー集中性に優れるという特性から、データ圧縮の世界で事実上の標準となっています。

本記事では、DCTの数学的定義からエネルギー集中性の原理、FFTを利用した高速計算法、そしてJPEG圧縮・MP3音声符号化への応用まで、Pythonコードを交えて体系的に解説します。

DCT-IIの定義

DCTにはDCT-I〜DCT-IVまでの変種がありますが、データ圧縮で「DCT」と呼ぶ場合は通常DCT-IIを指します。

長さ \(N\) の実数列 \(x[n]\) (\(n = 0, 1, \ldots, N-1\) )に対して、DCT-IIは次のように定義されます。

\[X[k] = \sum_{n=0}^{N-1} x[n] \cos\left(\frac{\pi k (2n+1)}{2N}\right), \quad k = 0, 1, \ldots, N-1 \tag{1}\]

コサイン引数 \(\frac{\pi k (2n+1)}{2N}\) の形状が重要です。半整数シフト \((2n+1)/2\) により、各基底関数は \(n = -0.5\) と \(n = N - 0.5\) で偶対称となり、これがエネルギー集中性の源泉になります(後述)。

DFTとの比較

DFTは複素指数関数 \(e^{-j2\pi kn/N}\) を基底とします。

\[X_{\text{DFT}}[k] = \sum_{n=0}^{N-1} x[n]\, e^{-j2\pi kn/N} \tag{2}\]

DCTとDFTの主な違いを以下に整理します。

特性DCT-IIDFT
基底関数コサイン(実数)複素指数
出力実数複素数
対称性偶対称拡張周期拡張
エネルギー集中性高い中程度
ギブス現象抑制発生しやすい
主な応用JPEG, MP3, H.264スペクトル解析, 通信

正規化形式

実用上は直交正規化形式が用いられます。

\[X[k] = w(k) \sum_{n=0}^{N-1} x[n] \cos\left(\frac{\pi k (2n+1)}{2N}\right) \tag{3}\] \[w(k) = \begin{cases} \sqrt{1/N} & (k = 0) \\ \sqrt{2/N} & (k \geq 1) \end{cases} \tag{4}\]

この正規化のもとでパーセバルの等式が成立します。

\[\sum_{n=0}^{N-1} x[n]^2 = \sum_{k=0}^{N-1} X[k]^2 \tag{5}\]

時間領域の総エネルギーと変換領域の総エネルギーが等しくなり、可逆変換であることが保証されます。

逆DCT(DCT-III)

DCT-IIの逆変換はDCT-IIIです。正規化形式では、DCT-IIIはDCT-IIの転置行列と一致します。

\[x[n] = \sum_{k=0}^{N-1} w(k)\, X[k] \cos\left(\frac{\pi k (2n+1)}{2N}\right) \tag{6}\]

変換行列を \(\mathbf{C}\) とすると \(\mathbf{C}\mathbf{C}^T = \mathbf{I}\) が成り立ち、これが逆変換の一意性を保証します。

エネルギー集中性

DCTが圧縮に強い最大の理由は**エネルギー集中性(Energy Compaction)**です。自然信号(画像・音声など)のDCT係数は、低次数(低周波)の少数の係数に大半のエネルギーが集中します。

ギブス現象との比較

DFT は信号の両端を暗黙的に周期接続します。信号の先端と末端に不連続がある場合(矩形波など)、スペクトルに高周波成分が現れるギブス現象が発生し、エネルギーが広い周波数帯域に分散します。

DCT-IIは、信号 \(x[n]\) を境界で偶対称に拡張(\(\tilde{x}[-n-1] = x[n]\) )してからDFTを適用するのと等価です(次節参照)。この偶拡張により、境界での不連続がなくなり、ギブス現象が抑制されます。その結果、エネルギーが低周波に集中します。

エネルギー集中度の定量的比較

自然画像・音声の代表的な統計モデルである1次マルコフ過程(隣接サンプルの相関係数 \(\rho\) が高い自己回帰過程)を使って、DCTとDFTのエネルギー集中度を比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import dct

# 1次マルコフ過程(自然画像・音声の統計モデル)で自然信号を模擬
np.random.seed(0)
N = 64
rho = 0.95
noise = np.random.randn(N)
x = np.zeros(N)
x[0] = noise[0]
for i in range(1, N):
    x[i] = rho * x[i - 1] + noise[i]

print(f"境界値: x[0]={x[0]:.4f}, x[N-1]={x[-1]:.4f}, 周期接続時の段差={x[-1] - x[0]:.4f}")

X_dct = dct(x, type=2, norm='ortho')
X_dft = np.fft.rfft(x)  # 実信号の非冗長な片側スペクトル(N//2+1 = 33ビン)

energy_dct = X_dct ** 2       # 64個、各1自由度(実数)
energy_dft = np.abs(X_dft) ** 2  # 33個、DC・ナイキスト以外は複素数=2自由度

# (a) 素朴な「係数の本数」で比較
order_dct = np.argsort(energy_dct)[::-1]
order_dft = np.argsort(energy_dft)[::-1]
cum_dct = np.cumsum(energy_dct[order_dct]) / energy_dct.sum()
cum_dft = np.cumsum(energy_dft[order_dft]) / energy_dft.sum()
n95_dct = np.searchsorted(cum_dct, 0.95) + 1
n95_dft = np.searchsorted(cum_dft, 0.95) + 1
print(f"[素朴比較] 95%エネルギーに必要な本数 — DCT: {n95_dct}, DFT(rfft): {n95_dft}")

# (b) 実自由度(DOF)で正規化した公平な比較
dof_weight = np.full(len(X_dft), 2.0)
dof_weight[0] = 1.0            # 直流成分は実数=1自由度
if N % 2 == 0:
    dof_weight[-1] = 1.0        # ナイキスト成分も実数=1自由度

cum_dof_dct = np.cumsum(np.ones_like(order_dct, dtype=float))
cum_dof_dft = np.cumsum(dof_weight[order_dft])
dof95_dct = cum_dof_dct[np.searchsorted(cum_dct, 0.95)]
dof95_dft = cum_dof_dft[np.searchsorted(cum_dft, 0.95)]
print(f"[DOF正規化] 95%エネルギーに必要な実自由度数 — DCT: {dof95_dct:.0f}, DFT: {dof95_dft:.0f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.6))
axes[0].plot(np.arange(1, len(cum_dct) + 1), cum_dct, color='#2a78d6', lw=2, label='DCT-II (coefficient count)')
axes[0].plot(np.arange(1, len(cum_dft) + 1), cum_dft, color='#e34948', lw=2, label='DFT (rfft bin count)')
axes[0].axhline(0.95, color='#898781', linestyle='--', lw=1, label='95% threshold')
axes[0].set_xlabel('Number of coefficients / bins')
axes[0].set_ylabel('Cumulative energy fraction')
axes[0].set_title('(a) Naive count comparison\n(looks nearly equal)')
axes[0].legend(fontsize=8, loc='lower right')
axes[0].grid(True, alpha=0.3)

axes[1].plot(cum_dof_dct, cum_dct, color='#2a78d6', lw=2, label='DCT-II (real DOF)')
axes[1].plot(cum_dof_dft, cum_dft, color='#e34948', lw=2, label='DFT (real DOF)')
axes[1].axhline(0.95, color='#898781', linestyle='--', lw=1, label='95% threshold')
axes[1].axvline(dof95_dct, color='#2a78d6', linestyle=':', lw=1)
axes[1].axvline(dof95_dft, color='#e34948', linestyle=':', lw=1)
axes[1].set_xlabel('Real degrees of freedom used')
axes[1].set_ylabel('Cumulative energy fraction')
axes[1].set_title(f'(b) DOF-normalized comparison\nDCT: {dof95_dct:.0f} DOF vs DFT: {dof95_dft:.0f} DOF')
axes[1].legend(fontsize=8, loc='lower right')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('energy_compaction.png', dpi=150, facecolor='white')
境界値: x[0]=1.7641, x[N-1]=-6.0537, 周期接続時の段差=-7.8177
[素朴比較] 95%エネルギーに必要な本数 — DCT: 6, DFT(rfft): 6
[DOF正規化] 95%エネルギーに必要な実自由度数 — DCT: 6, DFT: 11

DCTとDFTのエネルギー集中性比較。左図は係数・ビンの本数だけで比較すると両者はほぼ同等に見えるが、右図のように実自由度(DOF)で正規化するとDCTは6自由度、DFTは11自由度を要し、DCTがおよそ2倍効率的であることが分かる

実行結果は一見意外です。素朴に「本数」で比較すると、DCTもDFT(rfftの片側スペクトル)も95%エネルギーに6本の係数・ビンしか要らず、差がないように見えます。しかしこれは比較のフェアネスを欠いています。rfft が返す33個のビンのうち、直流成分とナイキスト成分(\(N\) が偶数の場合)を除く残り31個は複素数であり、実部・虚部の2自由度を持ちます。一方DCTの64個の係数はすべて実数で1自由度です。したがって「本数」ではなく信号を再構成するのに必要な**実自由度(degrees of freedom)**で数え直すと、DCTは6自由度で足りるのに対し、DFTは11自由度を要し、DCTの方がおよそ2倍効率的にエネルギーを集約できていることが分かります。この自由度の数え間違いは、DCTとDFTの圧縮性能を比較する際によくある落とし穴なので注意してください(詳しくは後述の「エッジケースと注意点」を参照)。

DCTとDFTの関係:偶拡張トリック

DCT-IIとDFTの数学的関係を導きます。

長さ \(N\) の信号 \(x[n]\) に対して、長さ \(2N\) の偶拡張信号 \(\tilde{x}[n]\) を構成します。

\[\tilde{x}[n] = \begin{cases} x[n] & 0 \leq n \leq N-1 \\ x[2N-1-n] & N \leq n \leq 2N-1 \end{cases} \tag{7}\]

この \(\tilde{x}[n]\) の長さ \(2N\) のDFTを \(\tilde{X}[k]\) とすると、次の関係が成立します。

\[X_{\text{DCT-II}}[k] = \frac{1}{2}\text{Re}\left(e^{-j\pi k/(2N)}\, \tilde{X}[k]\right), \quad k = 0, 1, \ldots, N-1 \tag{8}\]

完全な証明: 定義どおりに \(\tilde{X}[k]\) を展開します。

\[\tilde{X}[k] = \sum_{n=0}^{2N-1} \tilde{x}[n]\, e^{-j\pi kn/N}\]

和を \(n = 0, \ldots, N-1\) の前半部分と \(n = N, \ldots, 2N-1\) の後半部分に分割し、偶拡張の定義 \(\tilde{x}[n] = x[n]\) (前半)、\(\tilde{x}[n] = x[2N-1-n]\) (後半)を代入します。

\[\tilde{X}[k] = \sum_{n=0}^{N-1} x[n]\, e^{-j\pi kn/N} + \sum_{n=N}^{2N-1} x[2N-1-n]\, e^{-j\pi kn/N}\]

後半の和で \(m = 2N - 1 - n\) (\(n = N \Rightarrow m = N-1\) 、\(n = 2N-1 \Rightarrow m = 0\) )と変数変換すると、\(n = 2N - 1 - m\) なので、

\[\sum_{n=N}^{2N-1} x[2N-1-n]\, e^{-j\pi kn/N} = \sum_{m=0}^{N-1} x[m]\, e^{-j\pi k(2N-1-m)/N} = e^{j\pi k/N} \sum_{m=0}^{N-1} x[m]\, e^{j\pi km/N}\]

ここで \(e^{-j\pi k (2N-1)/N} = e^{-j2\pi k}\, e^{j\pi k/N} = e^{j\pi k/N}\) (\(e^{-j2\pi k}=1\) )を使いました。両辺に位相補正因子 \(e^{-j\pi k/(2N)}\) を掛けます。

\[e^{-j\pi k/(2N)}\, \tilde{X}[k] = \sum_{n=0}^{N-1} x[n]\, e^{-j\pi k(2n+1)/(2N)} + \sum_{n=0}^{N-1} x[n]\, e^{j\pi k(2n+1)/(2N)}\]

(第2項では \(e^{-j\pi k/(2N)} \cdot e^{j\pi k/N} \cdot e^{j\pi km/N} = e^{j\pi k(2m+1)/(2N)}\) を用いています。)括弧内の指数の和はオイラーの公式 \(e^{j\theta} + e^{-j\theta} = 2\cos\theta\) によりちょうど打ち消しあって実数になります。

\[e^{-j\pi k/(2N)}\, \tilde{X}[k] = \sum_{n=0}^{N-1} x[n] \left(e^{-j\theta_n} + e^{j\theta_n}\right) = 2\sum_{n=0}^{N-1} x[n] \cos\left(\frac{\pi k(2n+1)}{2N}\right), \quad \theta_n = \frac{\pi k(2n+1)}{2N}\]

右辺は式 \((1)\) の \(X_{\text{DCT-II}}[k]\) の2倍そのものであり、かつ既に実数なので \(\text{Re}(\cdot)\) は恒等的に作用します。ゆえに式 \((8)\) が成り立ちます(\(1/2\) の係数を落とすと式を誤るので注意)。

この式は、DCT-IIの計算はFFTを用いて \(O(N \log N)\) で実行できることを意味します。位相補正因子 \(e^{-j\pi k/(2N)}\) の乗算は \(O(N)\) で済みます。

高速DCT計算アルゴリズム

FFTを利用した高速DCT

式 \((8)\) を直接実装した高速DCTアルゴリズムは次のとおりです。

import numpy as np
from scipy.fft import dct

def fast_dct2(x):
    """FFTを利用したDCT-II実装(O(N log N))"""
    N = len(x)

    # ステップ1: 偶拡張(長さ2Nの信号を構成)
    x_ext = np.concatenate([x, x[::-1]])

    # ステップ2: FFT(長さ2N)
    X_fft = np.fft.fft(x_ext)

    # ステップ3: 位相補正と実部の取得
    k = np.arange(N)
    phase = np.exp(-1j * np.pi * k / (2 * N))
    X_dct = np.real(phase * X_fft[:N])

    return X_dct

# 検証:SciPyのdctと比較
np.random.seed(1)
x = np.random.randn(256)
X_fast = fast_dct2(x)
X_scipy = dct(x, type=2, norm=None)

print(f"最大誤差: {np.max(np.abs(X_fast - X_scipy)):.2e}")
print(f"平均誤差: {np.mean(np.abs(X_fast - X_scipy)):.2e}")
最大誤差: 2.13e-14
平均誤差: 4.31e-15

誤差は倍精度の丸め誤差レベル(\(10^{-14}\) オーダー)に収まっており、式 \((8)\) の導出が数値的にも正しいことが確認できます。ここで fast_dct2 は式 \((8)\) の右辺 \(\text{Re}(e^{-j\pi k/(2N)}\tilde{X}[k])\) をそのまま(\(1/2\) を掛けずに)返している点に注意してください。これは誤りではありません。SciPyの dct(x, type=2, norm=None) 自体が式 \((1)\) の \(2\) 倍、すなわち \(2\sum_{n} x[n]\cos(\cdot)\) を返す規約になっているため、\(1/2\) を掛けない生の値どうしがちょうど一致します(この正規化規約の違いは次節「エッジケースと注意点」で詳しく扱います)。

計算量の比較

さらなる高速化として、偶拡張後のFFTの対称性を利用した最適化があります。元の信号を偶数インデックス \(y[n] = x[2n]\) と奇数インデックス \(z[n] = x[2n+1]\) に分割し、長さ \(N/2\) の2本のFFTとして実行できます(Chen et al., 1977)。

手法演算量
素朴なDCT(定義通り)\(O(N^2)\)
偶拡張 + FFT\(O(N \log N)\)
Chen (1977) 最適化版\(\approx \frac{3}{2} N \log_2 N\) 乗算

実用的には scipy.fft.dct を使用すれば、内部でこれらの最適化が適用されます。

2次元DCT:画像処理への応用

2D DCT-IIの定義と分離可能性

\(M \times N\) の画像ブロック \(x[m, n]\) に対する2次元DCT-IIは次の式で定義されます。

\[X[u, v] = w(u)\, w(v) \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} x[m,n] \cos\left(\frac{\pi u (2m+1)}{2M}\right) \cos\left(\frac{\pi v (2n+1)}{2N}\right) \tag{9}\]

2D DCTは分離可能であり、まず各行に1D DCTを適用し、次に各列に1D DCTを適用することで計算できます。計算量は \(O(MN \log(MN))\) です。

8×8 DCT基底関数の視覚化

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import dct, idct

def dct2d(block):
    """2次元DCT-II(行→列の順に適用)"""
    return dct(dct(block, axis=0, norm='ortho'), axis=1, norm='ortho')

def idct2d(block):
    """2次元逆DCT"""
    return idct(idct(block, axis=0, norm='ortho'), axis=1, norm='ortho')

# 8x8 DCT基底関数の可視化
fig, axes = plt.subplots(8, 8, figsize=(10, 10))
for u in range(8):
    for v in range(8):
        basis = np.zeros((8, 8))
        basis[u, v] = 1.0
        basis_img = idct2d(basis)
        axes[u, v].imshow(basis_img, cmap='RdBu', vmin=-0.5, vmax=0.5)
        axes[u, v].axis('off')

plt.suptitle('8x8 DCT Basis Functions (row=u, col=v)', fontsize=14)
plt.tight_layout()
plt.savefig('dct_basis_functions.png', dpi=150, facecolor='white')

8×8 DCT基底関数(u, v = 0..7)。左上(u=0, v=0)は直流成分(一様な平均輝度)で、右方向にvが増えると水平方向の空間周波数が、下方向にuが増えると垂直方向の空間周波数が上がる。右下の高周波基底ほど自然画像のエネルギーが小さく、量子化でゼロになりやすい

左上(\(u=0, v=0\) )の基底は直流成分(全体の平均輝度)を表し、右下ほど高周波の格子模様になります。自然画像ではほとんどのエネルギーが左上の係数に集中するため、右下の係数は小さく、量子化でゼロに丸めることができます。

JPEG圧縮でのDCT

JPEGの圧縮手順

JPEG(Joint Photographic Experts Group)は1992年に標準化された静止画像圧縮フォーマットです。

圧縮手順の概略:

  1. 色空間変換: RGBからYCbCrへ変換(輝度・色差分離)
  2. クロマサブサンプリング: 人間の目が色差の細部に鈍感であることを利用
  3. 8×8ブロック分割: 画像を8×8ピクセルのブロックに分割
  4. レベルシフト: 各ピクセル値から128を引いて \([-128, 127]\) の範囲にシフト
  5. 2D DCT: 各ブロックに8×8 DCT-IIを適用
  6. 量子化: DCT係数を量子化テーブルで割り、整数に丸める(非可逆圧縮の本質
  7. エントロピー符号化: ジグザグスキャン後にハフマン符号化

量子化ステップは次式で表されます。

\[\hat{X}[u, v] = \text{round}\left(\frac{X[u, v]}{Q[u, v]}\right) \tag{10}\]

量子化テーブル \(Q[u, v]\) は低周波係数を細かく、高周波係数を粗く量子化します。これにより、視覚的に重要な低周波成分を保持しながら、高周波の細部を積極的に切り捨てます。

JPEGライク圧縮のPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import dct, idct

# JPEG標準輝度量子化テーブル(品質50相当)
QUANT_TABLE = np.array([
    [16, 11, 10, 16, 24, 40, 51, 61],
    [12, 12, 14, 19, 26, 58, 60, 55],
    [14, 13, 16, 24, 40, 57, 69, 56],
    [14, 17, 22, 29, 51, 87, 80, 62],
    [18, 22, 37, 56, 68, 109, 103, 77],
    [24, 35, 55, 64, 81, 104, 113, 92],
    [49, 64, 78, 87, 103, 121, 120, 101],
    [72, 92, 95, 98, 112, 100, 103, 99],
], dtype=float)

def dct2d(block):
    return dct(dct(block, axis=0, norm='ortho'), axis=1, norm='ortho')

def idct2d(block):
    return idct(idct(block, axis=0, norm='ortho'), axis=1, norm='ortho')

def jpeg_compress_block(block, quality=50):
    """8x8ブロックをJPEGライク圧縮"""
    if quality < 50:
        scale = 5000 / quality
    else:
        scale = 200 - 2 * quality
    q_table = np.clip(np.floor((QUANT_TABLE * scale + 50) / 100), 1, 255)

    shifted = block.astype(float) - 128
    coeffs = dct2d(shifted)
    quantized = np.round(coeffs / q_table)
    return quantized, q_table

def jpeg_decompress_block(quantized, q_table):
    """8x8ブロックを復元"""
    dequantized = quantized * q_table
    restored = idct2d(dequantized) + 128
    return np.clip(restored, 0, 255).astype(np.uint8)

# テスト画像(輝度グラデーション)の生成
np.random.seed(0)
img_size = 64
img = np.zeros((img_size, img_size), dtype=np.uint8)
for i in range(img_size):
    img[i, :] = np.clip(128 + 80 * np.sin(2 * np.pi * i / img_size), 0, 255)
img = (img.astype(int) + np.random.randint(0, 20, img.shape)).clip(0, 255).astype(np.uint8)

# 品質ごとの圧縮・復元
fig, axes = plt.subplots(1, 4, figsize=(14, 4))
qualities = [90, 50, 20, 5]

for ax, quality in zip(axes, qualities):
    restored_img = np.zeros_like(img, dtype=float)
    total_nonzero = 0
    total_coeffs = 0

    for row in range(0, img_size, 8):
        for col in range(0, img_size, 8):
            block = img[row:row+8, col:col+8].astype(float)
            quantized, q_table = jpeg_compress_block(block, quality)
            restored_block = jpeg_decompress_block(quantized, q_table)
            restored_img[row:row+8, col:col+8] = restored_block
            total_nonzero += np.count_nonzero(quantized)
            total_coeffs += quantized.size

    psnr_val = 10 * np.log10(255**2 / np.mean((img.astype(float) - restored_img)**2))
    ratio = total_nonzero / total_coeffs * 100

    ax.imshow(restored_img, cmap='gray', vmin=0, vmax=255)
    ax.set_title(f'Quality={quality}\nPSNR={psnr_val:.1f}dB\nnon-zero={ratio:.1f}%')
    ax.axis('off')
    print(f"Quality={quality}: PSNR={psnr_val:.2f} dB, 非ゼロ係数率={ratio:.2f}%")

plt.suptitle('JPEG-like Compression at Different Quality Levels', fontsize=13, y=1.05)
plt.tight_layout()
plt.savefig('jpeg_quality_comparison.png', dpi=150, facecolor='white', bbox_inches='tight')
Quality=90: PSNR=37.08 dB, 非ゼロ係数率=40.36%
Quality=50: PSNR=33.43 dB, 非ゼロ係数率=9.20%
Quality=20: PSNR=32.37 dB, 非ゼロ係数率=3.54%
Quality=5: PSNR=29.65 dB, 非ゼロ係数率=2.51%

品質90/50/20/5でのJPEGライク圧縮結果。品質を下げるほどPSNRが低下し(37.1→29.6dB)、非ゼロ係数の割合も40.4%→2.5%まで減少する。品質5では8×8ブロック境界に明瞭なブロッキングアーティファクトが視認できる

品質を下げるほどPSNRが低下し、非ゼロ係数の割合が減少(=圧縮率向上)します。品質90ではPSNR 37.1dB・非ゼロ係数40.4%と高品質を維持しますが、品質5ではPSNR 29.6dBまで低下し、非ゼロ係数はわずか2.5%まで減少します。図の品質5の画像をよく見ると、8×8ブロックの境界に沿って明るさが不連続に変化するブロッキングアーティファクトが視認できます。これは各ブロックが独立にDCT・量子化されるため、隣接ブロック間でDC成分(平均輝度)の量子化誤差が揃わないことに起因します。H.264/HEVCなどの動画コーデックがブロック境界をなだらかにするデブロッキングフィルタを後処理として持つのはこのためです。

なぜDCTが画像圧縮に適しているか

自然画像の統計的性質がDCTに適している理由は次の2点です。

  1. 1/f特性: 自然画像のパワースペクトルは低周波ほど大きく(\(P(f) \propto 1/f^2\) )、この性質がDCT係数の低次集中と対応します。
  2. マルコフ性とKLT近似: 隣接ピクセル間の相関が強く、1次マルコフモデルで近似できる場合、DCT変換行列はKLT(Karhunen-Loève変換)行列に漸近します。KLTは任意の相関行列を完全に非相関化する理想的な変換であり、DCTはその優れた近似です。

SciPyによるPython実装:DCTの圧縮効果デモ

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import dct, idct

def dct_compress(x, keep_ratio):
    """DCT係数の上位 keep_ratio 割合を保持して再構成"""
    X = dct(x, type=2, norm='ortho')
    N = len(X)
    k = max(1, int(N * keep_ratio))

    # エネルギーの大きい係数のみ保持
    idx = np.argsort(np.abs(X))[::-1]
    X_sparse = np.zeros_like(X)
    X_sparse[idx[:k]] = X[idx[:k]]

    x_rec = idct(X_sparse, type=2, norm='ortho')
    mse = np.mean((x - x_rec) ** 2)
    snr = 10 * np.log10(np.mean(x**2) / (mse + 1e-12))
    return x_rec, snr

# テスト信号
np.random.seed(7)
N = 256
n = np.arange(N)
x = (np.cos(2 * np.pi * 10 * n / N)
     + 0.5 * np.cos(2 * np.pi * 30 * n / N)
     + 0.3 * np.random.randn(N))

ratios = [0.5, 0.2, 0.1, 0.05]
fig, axes = plt.subplots(len(ratios) + 1, 1, figsize=(12, 12), sharex=True)

axes[0].plot(n, x, color='steelblue', label='Original')
axes[0].set_title('Original Signal')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

for ax, ratio in zip(axes[1:], ratios):
    x_rec, snr = dct_compress(x, ratio)
    print(f"keep_ratio={ratio}: SNR={snr:.2f} dB")
    ax.plot(n, x, color='steelblue', alpha=0.4, label='Original')
    ax.plot(n, x_rec, color='tomato',
            label=f'Reconstructed ({ratio*100:.0f}% coeffs, SNR={snr:.1f}dB)')
    ax.set_title(f'{ratio*100:.0f}% coefficients - SNR: {snr:.1f} dB')
    ax.legend()
    ax.grid(True, alpha=0.3)

axes[-1].set_xlabel('Sample index')
plt.tight_layout()
plt.savefig('dct_threshold_compression.png', dpi=150, facecolor='white')
keep_ratio=0.5: SNR=21.05 dB
keep_ratio=0.2: SNR=13.54 dB
keep_ratio=0.1: SNR=11.34 dB
keep_ratio=0.05: SNR=10.15 dB

DCT係数の上位k%だけを保持した再構成波形(原信号は薄い青、再構成は赤)。50%保持でSNR 21.0dB、10%保持でも11.3dBを維持しており、係数を大幅に間引いても波形の概形は保たれる。5%まで減らすとSNRは10.2dBまで落ち、高周波のノイズ成分の再現性が失われ始める

保持する係数の割合を50%→5%まで減らすとSNRは21.05dB→10.15dBまで単調に低下しますが、10%(26係数程度)を保持するだけでもSNRは11.34dBあり、視覚的には波形の概形をほぼ再現できています。これがDCTベース圧縮の威力です。ただし本実験の信号にはノイズ成分(0.3 * np.random.randn(N))が含まれており、ノイズはDCT領域でも広い係数に分散するため、保持比率を上げてもSNRの伸びは鈍化します(ノイズフロアへの漸近)。この鈍化自体もエネルギー集中性の限界を示す良い例です。

音声圧縮におけるMDCT

MP3やAACなどの音声符号化では、DCT-IIの変形である**MDCT(Modified Discrete Cosine Transform、変形離散コサイン変換)**が使われます。

MDCTの定義

MDCTは長さ \(2N\) の入力から \(N\) 個の係数を生成します。

\[X[k] = \sum_{n=0}^{2N-1} x[n] \cos\left(\frac{\pi}{N}\left(n + \frac{1}{2} + \frac{N}{2}\right)\left(k + \frac{1}{2}\right)\right), \quad k = 0, \ldots, N-1 \tag{11}\]

MDCTの重要な性質:

  1. 50%オーバーラップ: 連続するフレームが半分ずつ重なる
  2. TDAC(Time-Domain Aliasing Cancellation): フレーム間のエイリアシングが相殺され、完全再構成が可能
  3. 臨界サンプリング: \(2N\) サンプルから \(N\) 係数を生成(50%削減)

MP3では、サブバンドフィルタリングと組み合わせてMDCTを使用し、さらに心理音響モデル(人間の聴覚の感度が周波数により異なることを利用)で不要な係数を積極的に量子化します。

ケプストラム分析 で扱う対数スペクトルのDCT(メル周波数ケプストラム係数: MFCC)も音声認識の特徴量として広く使われています。より深いオーバーラップ処理・TDACの数学的証明・窓関数設計は MDCT(修正離散コサイン変換)とフィルタバンク で専門的に扱っているので、本記事ではDCT-IIとの接続点にとどめます。

エッジケースと注意点

DCTを実務で使う際につまずきやすい落とし穴を、実際に検証した数値とともにまとめます。

1. エネルギー集中性の比較は「本数」ではなく「自由度」で行う

前述のとおり、DCT係数とDFTビンを単純な本数で比較すると差がないように見えることがあります(実測: いずれも6本)。しかしDFTの非DC・非ナイキストビンは複素数(2自由度)を持つため、実自由度で数え直すとDCTは6、DFTは11となり、DCTが約1.8倍効率的です。信号処理ライブラリのベンチマークやブログ記事で「DCTとDFTの圧縮率を比較した」と謳うものの中には、この自由度の数え間違いに気づかず本数だけで比較しているケースがあるため、一次情報(自由度・ビット数)に立ち返って検証する習慣が重要です。

2. 正規化規約の違い(norm=None vs norm='ortho'

DCT-IIの「非正規化」形式には歴史的に複数の流儀があります。SciPyの dct(x, type=2, norm=None) は式 \((1)\) の2倍、すなわち \(2\sum_{n} x[n]\cos(\cdot)\) を返す規約を採用しています(実測で比を取ると厳密に2.0倍でした)。一方、正規直交形式 norm='ortho' は式 \((3)(4)\) のとおりで、norm=None の結果に \(w(k)/2\) を掛けたものと一致します。この違いに気づかず異なる正規化のDCT出力を混在させて演算すると、パーセバルの等式(式5)によるエネルギー保存チェックが合わなくなったり、量子化テーブルのスケールが2倍ずれたりする不具合を招きます。実装をまたいでDCT係数を扱う際は、必ず正規化オプションを明示し、両者の変換係数を検算してから使うことを推奨します。

3. DCT-I/II/III/IVは境界条件が異なり、用途を取り違えると誤差が出る

「DCT」と一括りにされがちですが、境界での対称性の取り方が変種ごとに異なります。

変種境界の対称性主な用途
DCT-I両端のサンプル点で対称(\(\tilde x[-n] = x[n]\) )まれ。境界サンプルの重みが半分になる点に注意
DCT-II両端のサンプル間の中点で対称(\(n=-0.5\) , \(N-0.5\) )JPEG・一般的な圧縮(本記事の主題)
DCT-IIIDCT-IIの逆変換(転置)逆量子化後の復元(IDCT)
DCT-IV両端ともサンプル間の中点で対称、かつ自己逆変換に近い性質MDCT(オーバーラップ変換)の基礎

JPEGの実装でうっかりDCT-Iを使うと係数の本数が合わない(DCT-Iは長さ \(N\) の入力に対し実質 \(N-1\) 個の独立自由度しか持たない)、MDCTの実装でDCT-IIを使うとTDACの相殺条件が満たされない、といった不具合につながります。ライブラリの type= 引数を必ず確認してください。

4. 「品質100」は非可逆圧縮では真のロスレスではない

JPEGライク圧縮のPythonコードで quality=100 を指定すると、量子化テーブルは全エントリが最小値の 1 になります(実測で確認済み: np.unique(q_table)[1.] のみ)。しかしこれは量子化ステップが最小になるだけであり、round() による丸め・レベルシフト・uint8 へのクリップは残るため、真のロスレスにはなりません。実測では品質100でもPSNRは51.1dB、最大画素誤差は1(256階調中)でした。真のロスレス圧縮が必要な場合はDCTベースのJPEGではなく、PNG(可逆圧縮)やJPEG-LS、あるいはJPEG 2000のロスレスモードを検討すべきです。

5. ブロッキングアーティファクトとその対策

8×8ブロックごとに独立してDCT・量子化するため、低品質ではブロック境界に不連続が生じるブロッキングアーティファクトが発生します(前掲の品質5の画像で視認可能)。これはブロックDCTの本質的な限界であり、H.264/HEVC/AV1などの動画コーデックはブロック境界を平滑化するデブロッキングフィルタ(ループ内フィルタ)を標準搭載しています。またMDCT(後述)はオーバーラップ処理によりこの問題を根本的に回避します。

6. 非2のべき乗長でのFFT高速化

式 \((8)\) のFFTベース高速DCTは、内部の np.fft.fft が \(N\) の素因数分解に応じて速度が大きく変わります。\(N\) が2の累乗であれば \(O(N \log N)\) が保証されますが、\(N\) が大きな素数を含む場合はSciPyの pocketfft バックエンドでも著しく遅くなることがあります。実務では、ブロックサイズを2の累乗(8, 16, 32など)に固定するか、必要に応じてゼロパディングで長さを調整することが推奨されます。

発展:DCTを置き換える学習型画像符号化——JPEG AI

DCTは50年以上にわたり画像・音声圧縮の主役であり続けていますが、2025年2月5日、ISO/IEC/ITU-TはJPEG AI(ISO/IEC 6048-1:2025)を国際標準として正式に採択しました。これは固定のブロックDCT変換を、学習済みニューラルネットワーク(オートエンコーダ)による変換に置き換えた、初のエンドツーエンド学習ベース画像符号化の国際標準です。JPEG AIのテストモデルは、VVC(Versatile Video Coding)のイントラ符号化アンカーに対して最大28〜29%のビットレート削減(BD-rate改善)を達成したと報告されています。これは、DCTという固定の解析的変換に対して、データから学習した非線形変換が符号化効率で優位に立ちうることを示す具体例です。

ただし、JPEG AIは既存のJPEG(DCTベース)との後方互換性を持たず、エンコード・デコードの両方に高いオンライン計算コスト(ニューラルネット推論)を要するため、Webブラウザでの画像表示や組み込み機器のようにDCTベースJPEG・MPEG系コーデックが低コストで動作する用途では、当面DCTが主流であり続けると見られます。DCTは「枯れた技術」でありながら、今なお比較対象として最新の学習型符号化研究のベースラインであり続けている点は興味深い事実です。

まとめ:DCT vs DFT vs ウェーブレット

特性DCT(DCT-II)DFTウェーブレット変換
出力実数複素数実数(マルチスケール)
計算量\(O(N \log N)\)\(O(N \log N)\)\(O(N)\)
エネルギー集中性高(自然信号に最適)高(時間局在も考慮)
境界処理偶対称拡張(ギブスなし)周期拡張(ギブスあり)反射・ゼロパディング等
時間周波数局在なし(大域的)なし(大域的)あり(マルチスケール)
主な応用JPEG, MP3, H.264スペクトル解析, 通信JPEG2000, 医用画像
逆変換DCT-III(実数)IDFT(複素数)逆ウェーブレット変換

本記事では、離散コサイン変換(DCT-II)の数学的定義から実用的な応用まで体系的に解説しました。

  • DCT-IIの定義: コサイン基底の重みつき和。正規化形式でパーセバルの等式が成立する
  • エネルギー集中性: 偶対称拡張によりギブス現象を回避し、低周波への高い集中性を実現
  • DFTとの関係: 偶拡張 + FFT + 位相補正で \(O(N \log N)\) の高速計算が可能
  • 2D DCT: 分離可能性により行・列に順次1D DCTを適用。8×8ブロックDCTがJPEGの中核
  • JPEG量子化: 視覚的な重要性に応じて高周波係数を粗く量子化し、高圧縮率を達成
  • MDCT: オーバーラップ処理とTDACにより、MP3/AACで完全再構成を実現
  • エッジケース: エネルギー集中性比較は自由度ベースで行う、正規化規約の違い(norm=Noneは式(1)の2倍)、DCT変種ごとの境界条件の違い、品質100でも非可逆であることなど、実務上の落とし穴に注意
  • 学習型符号化の台頭: 2025年に国際標準化されたJPEG AI(ISO/IEC 6048-1:2025)はDCTを学習済み変換で置き換え、VVCイントラ比で約28〜29%のビットレート削減を報告している

関連記事

参考文献

  • Ahmed, N., Natarajan, T., & Rao, K. R. (1974). “Discrete cosine transform”. IEEE Transactions on Computers, C-23(1), 90-93.
  • Wallace, G. K. (1991). “The JPEG still picture compression standard”. Communications of the ACM, 34(4), 30-44.
  • Chen, W. H., Smith, C. H., & Fralick, S. C. (1977). “A fast computational algorithm for the discrete cosine transform”. IEEE Transactions on Communications, 25(9), 1004-1009.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • ISO/IEC 6048-1:2025. “Information technology — JPEG AI learning-based image coding system — Part 1: Core coding system”. International Organization for Standardization(2025年2月採択、DCTを学習型変換に置き換える初の国際標準)。
  • SciPy DCT documentation