SL_specdecomp

SL_specdecomp is a PyMC-based package for Bayesian decomposition of neural power spectra into broadband and rhythmic components.

The package was developed for the manuscript:

Patrick F. Bloniasz and Emily P. Stephen. Spectral decompositions of neural voltage recordings are susceptible to model misspecifications that cause meaningful estimation error. DOI: https://doi.org/10.64898/2026.06.12.718232

The package is intended for model-based analysis of neural power spectra, such as spectra estimated from ECoG, EEG, LFP, or related field-potential recordings. It treats spectral decomposition as an estimation problem: the user specifies a latent spectral model, and the observed power spectrum is fit with a Gamma observation model.

Overview

Neural power spectra often contain both:

• broadband structure, often summarized by aperiodic height, slope, or knee parameters; and
• rhythmic structure, often summarized by center frequency, bandwidth, and amplitude.

SL_specdecomp fits these components using Bayesian inference. The default aperiodic predictor is a Lorentzian / 1/f-like function,

P_ap(f) = 10**b / (knee + f**chi)

and the default rhythmic predictor is a mirrored Gaussian bump. The model uses a Gamma likelihood for the observed linear-power spectrum.

Public API

The package exports:

from SL_specdecomp import Decompose, Decomposition, from_multitaper

Decompose

Fit a decomposition to a frequency grid and an already-computed power spectrum.

from_multitaper

Compute a multitaper spectrum from a time-domain signal using spectral_connectivity, then call Decompose.

Decomposition

Dataclass returned by fitted models. It stores the observed spectrum, posterior mean estimates, PyMC InferenceData, and plotting helpers.

Model modes

SL_specdecomp currently supports two model modes.

Mode	Latent spectrum	Returned components
additive	Broadband power plus one or more rhythmic components in linear power: mu = P_ap + P_rh.	estimated_spectrum, broadband, rhythms, broadband_components when present, and idata.
FOOOF_spectrum	FOOF-like multiplicative model, equivalent to adding rhythmic bumps in log10-power space.	Posterior estimate of the total spectrum and aperiodic baseline. Additional internal diagnostic attributes may be present, but this mode is mainly intended for comparison with multiplicative/log-power parametrizations.

Inputs

At minimum, Decompose expects:

freqs

One-dimensional frequency grid in Hz. Frequencies must be non-negative. The first value may be DC/0 Hz, but otherwise frequencies should be positive.

psd

One-dimensional power spectrum in linear power, not log power. It must have the same shape as freqs and should contain positive finite values.

mode

Either "additive" or "FOOOF_spectrum".

n_rhythms

Number of rhythmic components to fit.

rhythm_bands

List of (low, high) frequency bands constraining the rhythm centers. The length should match n_rhythms. If rhythms are included, it is best to provide these bands explicitly.

k_tapers

Number, or effective number, of approximately uncorrelated direct spectral estimates being averaged. This becomes the Gamma shape parameter.

sample_kwargs

Optional keyword arguments forwarded to pymc.sample.

The default sampler settings are:

dict(
    draws=1000,
    tune=1000,
    chains=2,
    target_accept=0.9,
    progressbar=False,
    random_seed=42,
    cores=1,
)

If blackjax is installed and the user does not explicitly set nuts_sampler, the package uses nuts_sampler="blackjax". Otherwise it uses the default PyMC NUTS sampler.

Quick start

Fit an additive decomposition to an existing power spectrum:

from SL_specdecomp import Decompose

model = Decompose(
    freqs,
    psd,
    mode="additive",
    n_rhythms=2,
    rhythm_bands=[(0.1, 4), (8, 12)],
    k_tapers=3,
    sample_kwargs={
        "draws": 800,
        "tune": 800,
        "chains": 2,
        "target_accept": 0.9,
    },
    plot=True,
)

total = model.estimated_spectrum
broadband = model.broadband
rhythms = model.rhythms
idata = model.idata

model.plot()

The returned object also provides model.total as an alias for model.estimated_spectrum.

FOOOF-like mode

Fit the multiplicative / FOOOF-like model:

from SL_specdecomp import Decompose

model = Decompose(
    freqs,
    psd,
    mode="FOOOF_spectrum",
    n_rhythms=2,
    rhythm_bands=[(0.1, 4), (8, 12)],
    k_tapers=3,
    sample_kwargs={
        "draws": 800,
        "tune": 800,
        "chains": 2,
        "target_accept": 0.9,
    },
    plot=True,
)

total = model.estimated_spectrum

When plot=True is passed to Decompose, component plotting is shown by default for additive fits. For FOOOF_spectrum fits, call model.plot() afterward if you want to inspect baseline and rhythmic contribution terms.

Multitaper convenience function

If you have a time-domain signal and want the package to compute the multitaper spectrum before fitting, use from_multitaper:

from SL_specdecomp import from_multitaper

model = from_multitaper(
    signal,
    fs,
    mode="additive",
    n_rhythms=2,
    rhythm_bands=[(0.1, 4), (8, 12)],
    time_halfbandwidth_product=3,
    n_tapers=5,
    freq_min=1.0,
    freq_max=200.0,
    downsample_independent=True,
    sample_kwargs={
        "draws": 800,
        "tune": 800,
        "chains": 2,
        "target_accept": 0.9,
    },
    plot=True,
)

from_multitaper flattens the input signal, estimates a multitaper spectrum with spectral_connectivity.Multitaper, averages non-frequency dimensions of the returned power array if necessary, optionally thins the frequency grid, and then calls Decompose.

The effective number of tapers passed to the Gamma likelihood is determined as:

n_tapers if n_tapers is not None else k_tapers

If downsample_independent=True, the frequency grid is thinned by an approximate multitaper-resolution stride before fitting. This is useful because the likelihood treats frequency bins as approximately independent.

Outputs

A fitted Decomposition object contains:

Attribute	Meaning
freqs	Frequency grid used for the fit.
observed	Observed linear-power PSD used for the fit.
estimated_spectrum	Posterior mean of the model-implied total spectrum.
total	Alias for estimated_spectrum.
broadband	Posterior mean aperiodic / baseline spectrum.
rhythms	In additive mode, posterior mean summed rhythmic power. In FOOOF_spectrum mode, posterior mean additive contribution above baseline.
broadband_components	Posterior mean individual aperiodic components with shape (M, F) when available.
rhythms	In additive mode, posterior mean summed rhythmic power. In FOOOF_spectrum mode, posterior mean additive contribution above baseline.
idata	Full PyMC / ArviZ InferenceData object.
mode	Fit mode: "additive" or "FOOOF_spectrum".

Plotting

Use the fitted object's plotting method:

ax = model.plot()

or, to suppress components:

ax = model.plot(plot_components=False)

The plot is drawn on log-log axes and includes the observed PSD and posterior mean spectrum. When components are requested and available, it also shows the broadband estimate, rhythmic contribution, and individual broadband components.

Installation

Recommended development installation

git clone https://github.com/Stephen-Lab-BU/SL_specdecomp.git
cd SL_specdecomp
conda env create -f environment.yml
conda activate SL_specdecomp
python -m pip install -e .

If Conda dependency resolution is slow, use Mamba:

conda install mamba -n base -c conda-forge
mamba env create -f environment.yml
mamba activate SL_specdecomp
python -m pip install -e .

Install directly from GitHub

python -m pip install git+https://github.com/Stephen-Lab-BU/SL_specdecomp.git

Requirements

See environment.yml for the recommended environment.

The package uses:

• numpy
• pymc
• pytensor
• arviz
• matplotlib
• spectral-connectivity for the from_multitaper workflow
• blackjax optionally, when available, for external NUTS sampling

Repository layout

.
├── SL_specdecomp/
│   ├── __init__.py
│   ├── api.py
│   ├── plotting.py
│   ├── predictors.py
│   ├── pymc_models.py
│   ├── types.py
│   └── version.py
├── Notebooks/
├── environment.yml
├── environment_jax.yml
├── setup.py
├── LICENSE
└── README.rst

Advanced usage

Multiple aperiodic components

Decompose accepts n_aperiodics. Values greater than 1 are supported by the code, but they trigger a warning because component identifiability and sampler convergence are experimental. Inspect posterior diagnostics carefully when using multiple aperiodic components.

Custom predictors

The default predictors are:

• APERIODIC_1OVERF with parameters b, knee, and chi;
• MIRRORED_GAUSSIAN with parameters center, sigma, and nyquist.

The package includes a lightweight predictor registry. Advanced users can pass custom predictor objects or registry keys through aperiodic_predictor and rhythm_predictor.

Custom parameter specifications

Advanced users can override parameter behavior with aperiodic_param_specs and rhythm_param_specs. A parameter specification may be a fixed numeric constant or a dictionary containing a PyMC distribution/factory.

Methodological notes

Gamma observation model

The likelihood is:

y(f) ~ Gamma(alpha=K, beta=K / mu(f))

where mu(f) is the latent mean spectrum and K is k_tapers. This parameterization gives:

E[y(f)] = mu(f)
Var[y(f)] = mu(f)**2 / K

Thus, the model treats power spectral estimates as heteroscedastic: frequencies with larger power have larger sampling variance.

Additive model

The additive model writes the latent spectrum as:

mu(f) = P_ap(f) + P_rh(f)

where P_ap is the aperiodic / broadband component and P_rh is the sum of linear-power rhythmic components.

FOOOF-like multiplicative model

The FOOOF_spectrum model writes the latent spectrum as a multiplicative model in linear power, equivalently an additive model in log10-power space:

log10(mu(f)) = log10(P_ap(f)) + sum_j a_j * G_j(f)

or equivalently:

mu(f) = P_ap(f) * 10**(sum_j a_j * G_j(f))

This is useful for comparing against log-power, FOOOF-like spectral parametrizations while still using the Gamma likelihood.

Choosing rhythm bands

Rhythm bands should reflect the experimental context and the scientific question. For example:

rhythm_bands = [
    (0.1, 4),    # slow-wave / delta
    (8, 20),     # alpha
    (20, 30),    # beta
]

If n_rhythms is greater than zero and rhythm_bands is omitted, the code falls back to repeated (1.0, 20.0) bands. This fallback is mainly a safety default; explicit bands are recommended for interpretable analyses.

Diagnostics

Because the package uses Bayesian inference, inspect diagnostics before interpreting results. At minimum, check:

• sampler warnings and divergences;
• effective sample sizes;
• r_hat values;
• posterior uncertainty for broadband and rhythm parameters;
• sensitivity to rhythm-band choices and priors;
• whether the posterior mean spectrum captures the observed PSD.

For expensive analyses, first run a short test fit, then increase draws, tune, and chains for final inference.

Version

The uploaded code reports:

__version__ = "0.0.3"

Citation

If you use SL_specdecomp, please cite the manuscript and this repository.

@misc{bloniasz_stephen_2026_spectral_decomposition,
  title  = {Spectral decompositions of neural voltage recordings are susceptible to model misspecifications that cause meaningful estimation error},
  author = {Bloniasz, Patrick F. and Stephen, Emily P.},
  year   = {2026},
  doi    = {10.64898/2026.06.12.718232},
  url    = {https://doi.org/10.64898/2026.06.12.718232},
  note   = {Preprint}
}

Related repositories

• SL_GPsim: Gaussian-process simulation tools for generating time-domain signals from analytically specified power spectra.
• Bloniasz_Stephen_Estimation: analysis and figure-generation code for the manuscript.

License

This repository is released under the MIT License. See LICENSE for details.

Contact

For questions, bug reports, or feature requests, please open an issue on the GitHub repository.