Power Spectral Density Example — Multi-Tone Signal

This example generates a noisy signal composed of sinusoids at 10 Hz, 50 Hz, and 120 Hz, computes the one-sided power spectral density (PSD) using transforms.fft_real_forward() (magnitude squared normalised by sample rate and length), and identifies the three dominant frequency peaks.

Example Code

"""IMSL Power Spectral Density example: PSD of a multi-tone signal with noise.

Signal: sum of sinusoids at 10 Hz, 50 Hz, 120 Hz embedded in white noise.

Outputs:
- Table printed to stdout showing dominant frequencies
- SVG plot saved to test_output/example_imsl_power_spectrum.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict, List

import matplotlib.pyplot as plt
import numpy as np

from transforms import fft_real_forward


def run_demo_imsl_power_spectrum() -> Dict[str, object]:
    """Run IMSL power spectral density example.

    Generates a noisy multi-tone signal with components at 10 Hz, 50 Hz, and
    120 Hz, computes the one-sided power spectral density (PSD) using
    :func:`transforms.fft_real_forward`, then identifies the dominant peaks.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``dominant_freqs``
            (list of float), ``psd_peak_values`` (list of float), and
            ``plot_path`` (str).
    """
    rng = np.random.default_rng(42)
    fs = 1024.0
    N = 2048
    dt = 1.0 / fs
    t = np.arange(N) * dt

    signal = (
        2.0 * np.sin(2 * np.pi * 10 * t)
        + 1.5 * np.sin(2 * np.pi * 50 * t)
        + 1.0 * np.sin(2 * np.pi * 120 * t)
        + 0.5 * rng.standard_normal(N)
    )

    X = fft_real_forward(signal)
    freqs = np.fft.rfftfreq(N, d=dt)
    # Power spectral density: |X|^2 / (fs * N), scaled to one-sided
    psd = (np.abs(X) ** 2) / (fs * N)
    psd[1:-1] *= 2.0  # double non-DC, non-Nyquist bins for one-sided

    # Identify top-3 peaks
    sorted_idx = np.argsort(psd)[::-1]
    top_idx = sorted_idx[:3]
    top_freqs: List[float] = freqs[top_idx].tolist()
    top_psd: List[float] = psd[top_idx].tolist()

    print("\nIMSL Power Spectral Density Example")
    print("Signal: 2·sin(2π·10t) + 1.5·sin(2π·50t) + sin(2π·120t) + noise")
    print(f"  Sample rate: {fs:.0f} Hz   N={N}")
    print("-" * 55)
    print(f"{'Rank':<6} {'Frequency (Hz)':>16} {'PSD (V²/Hz)':>14}")
    print("-" * 55)
    for i, (f, p) in enumerate(zip(top_freqs, top_psd), 1):
        print(f"{i:<6} {f:>16.2f} {p:>14.6f}")
    print("-" * 55)

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

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

    axes[0].plot(t[:512], signal[:512], color="#6d28d9", linewidth=0.9)
    axes[0].set_xlabel("Time (s)")
    axes[0].set_ylabel("Amplitude")
    axes[0].set_title("Time-Domain Signal (first 0.5 s)")
    axes[0].grid(True, alpha=0.3)

    axes[1].semilogy(freqs, psd, color="#7c3aed", linewidth=1.2)
    peak_colors = ["#dc2626", "#d97706", "#059669"]
    labels = ["10 Hz", "50 Hz", "120 Hz"]
    for freq, color, label in zip([10, 50, 120], peak_colors, labels):
        axes[1].axvline(freq, color=color, linestyle="--", linewidth=1.2, label=label)
    axes[1].set_xlabel("Frequency (Hz)")
    axes[1].set_ylabel("PSD (V²/Hz, log scale)")
    axes[1].set_title("One-Sided Power Spectral Density")
    axes[1].set_xlim(0, 200)
    axes[1].legend()
    axes[1].grid(True, alpha=0.3, which="both")

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

    return {
        "dominant_freqs": top_freqs,
        "psd_peak_values": top_psd,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_power_spectrum()

Plot Output

Top: time-domain signal showing the noisy composite waveform. Bottom: one-sided PSD (log scale) with vertical markers at 10 Hz, 50 Hz, and 120 Hz — the three injected sinusoidal components.

Console Output

IMSL Power Spectral Density Example
Signal: 2·sin(2π·10t) + 1.5·sin(2π·50t) + sin(2π·120t) + noise
  Sample rate: 1024 Hz   N=2048
-------------------------------------------------------
Rank     Frequency (Hz)    PSD (V┬▓/Hz)
-------------------------------------------------------
1                 10.00       4.060980
2                 50.00       2.276353
3                120.00       0.987153
-------------------------------------------------------