Spectral Analysis

This tutorial covers cogpy’s spectral analysis stack: PSD estimation, time-frequency spectrograms, and derived features.

PSD estimation

cogpy provides two PSD backends:

from cogpy.spectral.psd import psd_welch, psd_multitaper

# Welch's method (good for long signals)
psd_w, freqs = psd_welch(signal, fs=1000.0)

# Multitaper (better frequency resolution, recommended for short windows)
psd_m, freqs = psd_multitaper(signal, fs=1000.0)

Both return (psd, freqs) tuples. The PSD shape matches the input with the time axis replaced by a frequency axis: (..., freq).

xarray interface

For xarray DataArrays, use the specx module:

from cogpy.spectral.specx import psdx, spectrogramx

# PSD with xarray coordinates preserved
psd_da = psdx(sig, NW=4.0)
# <xarray.DataArray (AP: 16, ML: 16, freq: 513)>

# Time-frequency spectrogram
spec = spectrogramx(sig, window_size=512, window_step=128, NW=4.0)
# <xarray.DataArray (AP: 16, ML: 16, time_win: ..., freq: 257)>

Spectral features

All spectral features follow the PSD-first convention: they accept pre-computed (psd, freqs) and reduce over the frequency axis.

from cogpy.spectral.features import (
    band_power,
    relative_band_power,
    spectral_entropy,
    spectral_flatness,
    spectral_edge,
    broadband_snr,
    line_noise_ratio,
    narrowband_ratio,
    am_artifact_score,
    aperiodic_exponent,
)

# Band power — integrate PSD over a frequency band
gamma = band_power(psd, freqs, band=(30, 100))      # (...,) scalar
rel_gamma = relative_band_power(psd, freqs, band=(30, 100))  # fraction

# Spectral shape
entropy = spectral_entropy(psd, freqs)     # high = broadband, low = peaked
flatness = spectral_flatness(psd, freqs)   # 1.0 = white noise, 0.0 = pure tone
edge_95 = spectral_edge(psd, freqs, p=0.95)  # frequency below which 95% of power

# Signal quality
snr = broadband_snr(psd, freqs)
lnr = line_noise_ratio(psd, freqs, line_freq=60.0)

# Narrowband artifact detection
nb_ratio = narrowband_ratio(psd, freqs, flank_hz=5.0)  # (..., freq)
# Values >> 1 indicate narrowband peaks (line noise, artifact)

# Amplitude modulation
am_score = am_artifact_score(psd, freqs)

# Aperiodic fit (requires specparam/FOOOF)
exponent = aperiodic_exponent(psd, freqs)

Batch dimensions

All features broadcast over leading dimensions. For a grid PSD of shape (16, 16, 513) (AP, ML, freq), each feature returns (16, 16).

Line noise detection

For detecting narrowband lines at unknown frequencies, use the F-test scanner:

from cogpy.spectral.features import ftest_line_scan

# Scans all frequency bins for sinusoidal components
fstat, freqs, sig_mask = ftest_line_scan(signal_1d, fs=1000.0, NW=4.0)

# sig_mask[i] == True means a significant line at freqs[i]
line_freqs = freqs[sig_mask]

This implements Thomson’s (1982) eigenvalue-weighted F-test across all frequency bins simultaneously.

Next steps