Complex FFT Example — Complex Exponential Signal Analysis

This example builds a two-tone complex exponential signal \(z(t) = e^{2\pi i \cdot 5t} + 0.5 \cdot e^{2\pi i \cdot 12t}\) with N = 64 samples and demonstrates:

The round-trip error \(\max|z - \text{IFFT}(\text{FFT}(z))|\) is verified to be at floating-point machine precision, and the point-wise difference between fft_complex_forward and fast_dft confirms numerical equivalence.

Example Code

"""IMSL complex FFT example: complex exponential signal analysis.

Signal: z(t) = exp(2πi·5·t) + 0.5·exp(2πi·12·t), N=64 points

Covers:
  - fft_complex_forward / fft_complex_backward  (IMSL FFTCF / FFTCB)
  - fft_complex_init / fft_real_init            (IMSL FFTCI / FFTRI)
  - fast_dft                                    (IMSL FAST_DFT)

Outputs:
  - Summary printed to stdout showing FFT vs direct-DFT max error
  - SVG plot saved to test_output/example_imsl_complex_fft.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 (
    fast_dft,
    fft_complex_backward,
    fft_complex_forward,
    fft_complex_init,
    fft_real_init,
)


def run_demo_imsl_complex_fft() -> Dict[str, object]:
    """Run IMSL complex FFT example: complex exponential signal analysis.

    Builds a two-tone complex exponential z(t) = exp(2πi·5·t) +
    0.5·exp(2πi·12·t) sampled at N=64 points, computes its spectrum via
    fft_complex_forward, recovers the signal via fft_complex_backward, and
    cross-checks the result against fast_dft (direct O(N²) DFT).  Init
    metadata from fft_complex_init and fft_real_init is also demonstrated.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``max_diff_fft_dft`` (float),
            ``roundtrip_error`` (float), and ``plot_path`` (str).
    """
    N = 64
    dt = 1.0 / N
    t = np.arange(N) * dt

    # Complex signal: two complex exponentials
    z = np.exp(2j * np.pi * 5 * t) + 0.5 * np.exp(2j * np.pi * 12 * t)

    # --- init metadata (compatibility stubs) ---
    c_params = fft_complex_init(N)
    r_params = fft_real_init(N)

    # --- forward complex FFT ---
    Z_fft = fft_complex_forward(z)

    # --- inverse complex FFT: round-trip ---
    z_recovered = fft_complex_backward(Z_fft)
    roundtrip_error = float(np.max(np.abs(z_recovered - z)))

    # --- direct DFT for comparison ---
    Z_dft = fast_dft(z)
    max_diff = float(np.max(np.abs(Z_fft - Z_dft)))

    freqs = np.fft.fftfreq(N, d=dt)

    print("\nIMSL Complex FFT Example: z(t) = exp(2πi·5t) + 0.5·exp(2πi·12t), N=64")
    print("-" * 62)
    print(f"{'Init metadata (complex) n':<42} {c_params['n']:>10}")
    print(f"{'Init metadata (real) type':<42} {r_params['type']:>10}")
    print(f"{'Round-trip max |z - IFFT(FFT(z))|':<42} {roundtrip_error:>10.2e}")
    print(f"{'Max |fft_complex_forward - fast_dft|':<42} {max_diff:>10.2e}")
    print("-" * 62)

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

    fig, axes = plt.subplots(3, 1, figsize=(10, 10))

    # Signal real and imag parts
    axes[0].plot(t, z.real, color="#6d28d9", linewidth=1.5, label="Re[z(t)]")
    axes[0].plot(t, z.imag, color="#f59e0b", linewidth=1.5, label="Im[z(t)]", alpha=0.8)
    axes[0].set_xlabel("Time (s)")
    axes[0].set_ylabel("Amplitude")
    axes[0].set_title("Complex Signal: z(t) = e^{2πi·5t} + 0.5·e^{2πi·12t}")
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)

    # Magnitude spectrum (fft_complex_forward)
    mag_fft = np.abs(Z_fft)
    axes[1].vlines(freqs, 0, mag_fft, color="#7c3aed", linewidth=1.2)
    axes[1].plot(freqs, mag_fft, "o", color="#7c3aed", markersize=3)
    axes[1].axvline(5, color="#dc2626", linestyle="--", linewidth=1.2, label="f=+5 Hz")
    axes[1].axvline(12, color="#d97706", linestyle="--", linewidth=1.2, label="f=+12 Hz")
    axes[1].set_xlabel("Frequency (Hz)")
    axes[1].set_ylabel("|Z(f)|")
    axes[1].set_title("Magnitude Spectrum (fft_complex_forward)")
    axes[1].set_xlim(-N // 2, N // 2)
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)

    # Point-wise difference: FFT vs direct DFT
    diff = np.abs(Z_fft - Z_dft) + 1e-18
    axes[2].semilogy(np.arange(N), diff, color="#059669", linewidth=1.5)
    axes[2].set_xlabel("Bin index k")
    axes[2].set_ylabel("|FFT[k] − DFT[k]| (log scale)")
    axes[2].set_title("Numerical Difference: fft_complex_forward vs fast_dft")
    axes[2].grid(True, alpha=0.3)

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

    return {
        "max_diff_fft_dft": max_diff,
        "roundtrip_error": roundtrip_error,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_complex_fft()

Plot Output

Three-panel plot of complex signal, magnitude spectrum, and FFT vs DFT difference

Top: real and imaginary parts of the complex exponential signal. Middle: two-sided magnitude spectrum with peaks at ±5 Hz and ±12 Hz. Bottom: point-wise numerical difference between fft_complex_forward and fast_dft (log scale), confirming machine-precision agreement.

Console Output

IMSL Complex FFT Example: z(t) = exp(2πi·5t) + 0.5·exp(2πi·12t), N=64
--------------------------------------------------------------
Init metadata (complex) n                          64
Init metadata (real) type                        real
Round-trip max |z - IFFT(FFT(z))|            4.97e-16
Max |fft_complex_forward - fast_dft|         0.00e+00
--------------------------------------------------------------