Quarter-Wave Transforms Example — DCT/DST Signal Compression

This example demonstrates JPEG-style 1-D lossy compression using the quarter-wave cosine (type-2 DCT) and quarter-wave sine (type-1 DST) transforms. A piecewise sinusoidal signal of N = 128 samples is transformed, the smallest 75% of coefficients are zeroed out, and the signal is reconstructed. The signal-to-noise ratio (SNR) quantifies the compression quality.

Functions demonstrated:

Example Code

"""IMSL quarter-wave transforms example: DCT/DST signal compression.

Demonstrates:
  - fft_quarter_cosine_forward / fft_quarter_cosine_backward  (IMSL QCOSF/QCOSB)
  - fft_quarter_sine_forward   / fft_quarter_sine_backward    (IMSL QSINF/QSINB)
  - fft_sine_init / fft_cosine_init                           (IMSL FSINI/FCOSI)
  - fft_quarter_sine_init / fft_quarter_cosine_init           (IMSL QSINI/QCOSI)

Application: JPEG-style 1-D DCT compression — keep top 25% of coefficients,
reconstruct, and report signal-to-noise ratio.

Outputs:
  - Compression statistics printed to stdout
  - SVG plot saved to test_output/example_imsl_quarter_transforms.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from transforms import (
    fft_cosine_init,
    fft_quarter_cosine_backward,
    fft_quarter_cosine_forward,
    fft_quarter_cosine_init,
    fft_quarter_sine_backward,
    fft_quarter_sine_forward,
    fft_quarter_sine_init,
    fft_sine_init,
)


def run_demo_imsl_quarter_transforms() -> Dict[str, object]:
    """Run IMSL quarter-wave DCT/DST compression example.

    Builds a piecewise sinusoidal signal, applies the type-2 DCT (quarter-wave
    cosine forward), zeros out the smallest 75% of coefficients to simulate
    lossy compression, then reconstructs via the inverse DCT (quarter-wave
    cosine backward).  The same signal is processed with the type-1 DST
    (quarter-wave sine) for comparison.  Init stubs for all four transform
    variants are also exercised.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``kept_pct`` (float),
            ``dct_snr_db`` (float), ``dst_snr_db`` (float), and
            ``plot_path`` (str).
    """
    N = 128
    t = np.linspace(0, 1, N, endpoint=False)
    signal = (
        np.sin(2 * np.pi * 3 * t)
        + 0.5 * np.sin(2 * np.pi * 9 * t)
        + 0.25 * np.sin(2 * np.pi * 18 * t)
    )

    # --- init stubs ---
    s_init  = fft_sine_init(N)
    c_init  = fft_cosine_init(N)
    qs_init = fft_quarter_sine_init(N)
    qc_init = fft_quarter_cosine_init(N)

    # --- DCT compression (quarter cosine, type-2) ---
    dct_coeffs = fft_quarter_cosine_forward(signal)
    keep = int(N * 0.25)
    idx_sorted = np.argsort(np.abs(dct_coeffs))[::-1]
    dct_compressed = np.zeros(N)
    dct_compressed[idx_sorted[:keep]] = dct_coeffs[idx_sorted[:keep]]
    dct_reconstructed = fft_quarter_cosine_backward(dct_compressed)
    dct_error = float(np.sqrt(np.mean((dct_reconstructed - signal) ** 2)))
    signal_rms = float(np.sqrt(np.mean(signal ** 2)))
    dct_snr = float(20 * np.log10(signal_rms / (dct_error + 1e-15)))

    # --- DST round-trip (quarter sine, type-1) ---
    dst_coeffs = fft_quarter_sine_forward(signal)
    dst_compressed = np.zeros(N)
    dst_compressed[idx_sorted[:keep]] = dst_coeffs[idx_sorted[:keep]]
    dst_reconstructed = fft_quarter_sine_backward(dst_compressed)
    dst_error = float(np.sqrt(np.mean((dst_reconstructed - signal) ** 2)))
    dst_snr = float(20 * np.log10(signal_rms / (dst_error + 1e-15)))

    kept_pct = 100.0 * keep / N

    print("\nIMSL Quarter-Wave Transforms Example (N=128)")
    print("-" * 55)
    print(f"{'Init type (sine)':<38} {s_init['type']:>12}")
    print(f"{'Init type (cosine)':<38} {c_init['type']:>12}")
    print(f"{'Init type (quarter_sine)':<38} {qs_init['type']:>12}")
    print(f"{'Init type (quarter_cosine)':<38} {qc_init['type']:>12}")
    print(f"{'Coefficients kept (%)':<38} {kept_pct:>12.1f}")
    print(f"{'DCT SNR (dB) — 25% kept':<38} {dct_snr:>12.2f}")
    print(f"{'DST SNR (dB) — 25% kept':<38} {dst_snr:>12.2f}")
    print("-" * 55)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_quarter_transforms.svg"

    fig, axes = plt.subplots(2, 2, figsize=(12, 8))

    # Original vs DCT reconstruction
    axes[0, 0].plot(t, signal, color="#6d28d9", linewidth=1.5, label="Original")
    axes[0, 0].plot(t, dct_reconstructed, "--", color="#dc2626",
                    linewidth=1.5, label=f"DCT 25% ({dct_snr:.1f} dB SNR)")
    axes[0, 0].set_title("Original vs DCT Reconstruction (25% coefficients)")
    axes[0, 0].set_xlabel("t")
    axes[0, 0].set_ylabel("Amplitude")
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)

    # Original vs DST reconstruction
    axes[0, 1].plot(t, signal, color="#6d28d9", linewidth=1.5, label="Original")
    axes[0, 1].plot(t, dst_reconstructed, "--", color="#059669",
                    linewidth=1.5, label=f"DST 25% ({dst_snr:.1f} dB SNR)")
    axes[0, 1].set_title("Original vs DST Reconstruction (25% coefficients)")
    axes[0, 1].set_xlabel("t")
    axes[0, 1].set_ylabel("Amplitude")
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)

    # DCT coefficient magnitudes
    axes[1, 0].bar(np.arange(N), np.abs(dct_coeffs), color="#7c3aed", width=1.0, alpha=0.85)
    axes[1, 0].axvline(keep - 0.5, color="#dc2626", linestyle="--", linewidth=1.2,
                       label=f"Keep threshold (top {keep})")
    axes[1, 0].set_title("DCT Coefficient Magnitudes (type-2)")
    axes[1, 0].set_xlabel("Coefficient index (sorted by magnitude)")
    axes[1, 0].set_ylabel("|C_k|")
    axes[1, 0].set_xlim(0, N)
    axes[1, 0].grid(True, alpha=0.3)

    # DST coefficient magnitudes
    axes[1, 1].bar(np.arange(N), np.abs(dst_coeffs), color="#f59e0b", width=1.0, alpha=0.85)
    axes[1, 1].axvline(keep - 0.5, color="#dc2626", linestyle="--", linewidth=1.2,
                       label=f"Keep threshold (top {keep})")
    axes[1, 1].set_title("DST Coefficient Magnitudes (type-1)")
    axes[1, 1].set_xlabel("Coefficient index (sorted by magnitude)")
    axes[1, 1].set_ylabel("|S_k|")
    axes[1, 1].set_xlim(0, N)
    axes[1, 1].grid(True, alpha=0.3)

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {
        "kept_pct": kept_pct,
        "dct_snr_db": dct_snr,
        "dst_snr_db": dst_snr,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_quarter_transforms()

Plot Output

Four-panel plot showing DCT/DST compression results

Top-left: original signal vs DCT reconstruction (25% of coefficients kept). Top-right: original signal vs DST reconstruction (25% of coefficients kept). Bottom-left: DCT coefficient magnitudes (type-2). Bottom-right: DST coefficient magnitudes (type-1).

Console Output

IMSL Quarter-Wave Transforms Example (N=128)
-------------------------------------------------------
Init type (sine)                               sine
Init type (cosine)                           cosine
Init type (quarter_sine)               quarter_sine
Init type (quarter_cosine)             quarter_cosine
Coefficients kept (%)                          25.0
DCT SNR (dB) ΓÇö 25% kept                       34.11
DST SNR (dB) ΓÇö 25% kept                       20.79
-------------------------------------------------------