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

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

# 典型的な自然信号（徐々に変化する成分 + 小さなノイズ）
np.random.seed(42)
N = 64
n = np.arange(N)
x = (np.cos(2 * np.pi * 2 * n / N)
     + 0.5 * np.cos(2 * np.pi * 5 * n / N)
     + 0.1 * np.random.randn(N))

# DCT-II（正規化あり）
X_dct = dct(x, type=2, norm='ortho')

# DFT（片側スペクトル）
X_dft = np.fft.rfft(x)
energy_dft = np.abs(X_dft) ** 2
energy_dct = X_dct ** 2

# 上位k個の係数で占めるエネルギー累積分布
energy_dct_sorted = np.sort(energy_dct)[::-1]
energy_dft_sorted = np.sort(energy_dft)[::-1]

cumsum_dct = np.cumsum(energy_dct_sorted) / np.sum(energy_dct_sorted)
cumsum_dft = np.cumsum(energy_dft_sorted) / np.sum(energy_dft_sorted)

# プロット
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].bar(range(N), energy_dct / energy_dct.sum(),
            alpha=0.7, label='DCT', color='steelblue')
axes[0].bar(range(len(energy_dft)), energy_dft / energy_dft.sum(),
            alpha=0.5, label='DFT', color='tomato')
axes[0].set_xlabel('Coefficient index')
axes[0].set_ylabel('Normalized energy')
axes[0].set_title('Energy Distribution: DCT vs DFT')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(cumsum_dct, label='DCT', color='steelblue')
axes[1].plot(cumsum_dft, label='DFT', color='tomato')
axes[1].axhline(0.95, color='gray', linestyle='--', label='95% threshold')
axes[1].set_xlabel('Number of coefficients (sorted by energy)')
axes[1].set_ylabel('Cumulative energy fraction')
axes[1].set_title('Cumulative Energy Compaction')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 95%エネルギーに必要な係数数
n95_dct = np.searchsorted(cumsum_dct, 0.95) + 1
n95_dft = np.searchsorted(cumsum_dft, 0.95) + 1
print(f"95%エネルギーに必要な係数数 — DCT: {n95_dct}, DFT: {n95_dft}")

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と比較
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}")
# 出力例: 最大誤差: 2.13e-12

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('8×8 DCT Basis Functions (row=u, col=v)', fontsize=14)
plt.tight_layout()
plt.show()

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\n非ゼロ={ratio:.1f}%')
    ax.axis('off')

plt.suptitle('JPEG-like Compression at Different Quality Levels', fontsize=12)
plt.tight_layout()
plt.show()

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

# テスト信号
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)
    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.show()