speckit 1.0.0

A Python toolkit for Spectral Analysis.

SpecKit: Python library for dynamic-resolution spectral analysis

Miguel Dovale
Gerhard Heinzel

SpecKit is a high-performance Python toolkit for advanced spectral analysis, designed for scientists and engineers who require high-fidelity, high-dynamic-range spectral estimates.

The library's core feature is its logarithmically-spaced frequency scheduling algorithm. Unlike traditional methods such as Welch’s, which use a fixed segment length and resolution across all frequencies, SpecKit dynamically adjusts the frequency resolution. The result is that, for a given dataset, this method provides many more high-resolution spectral estimates at low frequencies, while minimizing statistical uncertainty (via more averages) at high frequencies. This method is ideally suited for ana...

Built on rigorous, textbook-backed theory, SpecKit also provides comprehensive error analysis for all its estimates, making it a reliable tool for scientific research.

---

Key Features

• Logarithmic Frequency Scheduling
  Dynamically adapts segment length and averaging to deliver optimal resolution and stability across the entire frequency spectrum.

• High-Performance Engine
  Leverages numba for JIT-compiling and parallelizing performance-critical code. Provides pure NumPy fallbacks for full portability.

• Comprehensive Spectral Quantities
  Calculates a full suite of auto-spectral and cross-spectral estimates, including:

  • Power Spectral Density (PSD) and Amplitude Spectral Density (ASD)
  • Cross Spectral Density (CSD)
  • Transfer Functions and Coherence
  • Conditioned spectra for MISO/SISO systems

• Rigorous Error Analysis
  Provides analytically derived standard deviations and normalized random errors for all key estimates, following Bendat & Piersol.
  Additionally, SpecKit computes empirical errors directly from segment scatter using Welford-style accumulation (see below).

• Advanced DSP Utilities
  Includes tools for:

  • Multi-file data loading/resampling
  • High-order Lagrange interpolation for fractional time shifting
  • Optimal linear combination for noise subtraction
  • Specialized flat-top windows with known Equivalent Noise Bandwidth (ENBW)

• High-Quality Noise Generation
  Ultra-fast, numerically stable generators for white, red (Brownian), pink, and general 1/f^α colored noise.

• Publication-Ready Plotting
  Built-in, user-friendly plotting interface for creating clear and professional figures.

---

Installation

Install SpecKit directly from PyPI:

pip install speckit

For development:

git clone https://github.com/mdovale/speckit.git
cd speckit
pip install -e .[dev]

---

Quickstart

The main user-facing function is speckit.compute_spectrum.
Here’s a simple example of computing the Amplitude Spectral Density (ASD) of a noisy signal and comparing it to SciPy’s Welch method:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
from speckit import compute_spectrum

N = int(1e6)  # samples
fs = 2.0      # Hz
x = np.random.randn(N)

# Welch (fixed resolution)
f_welch, psd_welch = welch(x, fs, window="hann", nperseg=N // 8)

# SpecKit (dynamic resolution)
result = compute_spectrum(x, fs=fs, Jdes=1000, Kdes=100, win="Kaiser", psll=200)

fig, ax = plt.subplots(figsize=(8, 5), dpi=150)
ax.loglog(f_welch, np.sqrt(psd_welch), label="Welch", color="gray", alpha=0.7)
result.plot(which="asd", errors=True, sigma=1, ax=ax, label="SpecKit", color="crimson")
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel(r"ASD (units / $\sqrt{\rm Hz}$)")
ax.legend()
ax.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()

---

API Overview

Top-level entry point

• speckit.compute_spectrum(data, fs, **kwargs)
  Configures, plans, and executes a spectral analysis in one line.

Core classes

• SpectrumAnalyzer
  Configuration + computation object.
  Methods:

  • .plan() → returns frequency/segmentation plan.
  • .compute() → full spectrum analysis.
  • .compute_single_bin(freq, fres or L) → analyze a single frequency bin.

• SpectrumResult
  Immutable container with all results and helper methods.
  Attributes are lazily computed and cached on demand:

  • psd, asd, csd, tf, coh, ENBW, …
  • Errors: Gxx_dev, Hxy_dev, coh_error, etc.
  • Empirical errors: XY_emp_dev, Gxy_emp_dev, …
  • Utilities: .plot(), .to_dataframe(), .get_rms(), .get_measurement().

---

Performance Strategies

SpecKit is designed for datasets with millions of samples. Performance comes from:

1. Numba JIT Kernels

   • Core bin-by-bin statistics are computed with hand-optimized Goertzel recurrences.
   • Parallelized across segments using numba.prange.

2. Selective Detrending

   • Order = -1: window only
   • Order = 0: mean removal
   • Order = 1,2: polynomial detrend via orthonormal basis (_build_Q)
   • The detrend mode changes which JIT kernel is dispatched.

3. Window/Polynomial Caching

   • Expensive objects like windows and polynomial bases are cached per segment length.

4. Fallback to NumPy

   • If Numba is unavailable, pure-NumPy helpers vectorize segment processing.
   • NumPy path uses highly optimized polyfit/polyval for detrending.

5. Single-Bin Evaluation

   • .compute_single_bin() avoids planning full spectra and is highly efficient for probing specific frequencies.

---

Error Analysis and Welford’s Algorithm

Textbook formulas

SpecKit implements Bendat & Piersol-style analytic formulas:

• Variances of PSD, CSD, and coherence estimates scale as 1/√K, where K is the number of averages.
• Coherence errors scale with (1 - γ²) / √K.

Empirical scatter (M2)

In addition, the JIT kernels compute:

• MXX, MYY: mean segment powers
• M2: mean squared deviation of per-segment cross-spectra about the mean

This is essentially Welford’s algorithm for streaming variance, adapted to complex quantities.

From M2, SpecKit derives:

• XY_emp_var: variance of the averaged cross-spectrum
• XY_emp_dev: empirical standard deviation
• Gxx_emp_dev, Gxy_emp_dev: empirical deviations mapped into spectral density units

These provide non-parametric uncertainty estimates, complementing the analytic formulas. They are particularly useful when:

• The Gaussian assumptions of textbook theory may not hold
• Segments are correlated due to overlap
• One wants to validate that statistical convergence is as expected

Why both?

• Analytic errors are fast, interpretable, and match the literature.
• Empirical errors reflect actual scatter in the data.
  Together, they provide a complete picture of uncertainty.

---

References and Acknowledgements

1. General Spectral Analysis and Uncertainty

   • Bendat, Piersol. Random Data: Analysis and Measurement Procedures. John Wiley & Sons.
   • Bendat, Piersol. Engineering Applications of Correlation and Spectral Analysis. John Wiley & Sons.

2. Logarithmic Frequency Scheduling

   • Tröbs, M., Heinzel, G. Improved spectrum estimation from digitized time series on a logarithmic frequency axis. doi:10.1016/j.measurement.2005.10.010
   • Heinzel, G., Rüdiger, A., & Schilling, R. (2002). Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows. Max-Planck-Institut für Gravitationsphysik.
   • Basis for SpecKit’s lpsd_plan and ltf_plan scheduler, inspired by algorithms developed for LISA Pathfinder data analysis.

3. Colored Noise Generation

   • Plaszczynski, S. (2007). Generating long streams of 1/f^α noise. Fluctuation and Noise Letters. doi:10.1142/S0219477507003635
   • Provides the stable cascade-of-filters method used in SpecKit.
   • Numba-optimized version of original code by Jan Waldmann

4. High-order Lagrange Interpolation

   • 
   • dsp.timeshift and dsp.lagrange_taps adopted from the PyTDI project. Credit to: Staab, Martin, Bayle, Jean-Baptiste, Hartwig, Olaf

---

License

BSD 3-Clause License. See the LICENSE file.

---

Contributing

Contributions are welcome!
Please open issues or PRs for bug reports, feature requests, or ideas.