Comparing pyrasa with neurodsp on their IRASA implementations.

irasa_test.py

import neurodsp
import neurodsp.aperiodic
import neurodsp.plts
import neurodsp.sim
import neurodsp.spectral
import numpy as np
from pyrasa.irasa import irasa

n_seconds = 10
fs = 500  # sampling rate in Hz

# Parameters of the simulated components
cf = 10  # frequency of the oscillatory component
exp = -2  # exponent of the power-law
print(f"simulated exponent: {exp:1.2f}")

# Components for the simulated signal
components = dict(sim_oscillation=dict(freq=cf), sim_powerlaw=dict(exponent=exp))

# Frequency range of interest
f_range = (1, 40)

# Create the stimulated time series
sig = neurodsp.sim.sim_combined(n_seconds, fs, components)
time = np.arange(len(sig)) / fs

# Power spectrum
freqs, psd = neurodsp.spectral.compute_spectrum(sig, fs, nperseg=4 * fs)
freqs, psd = neurodsp.spectral.trim_spectrum(freqs, psd, f_range)
neurodsp.plts.plot_power_spectra(freqs, psd, title="Original spectrum")

# Compute the IRASA decomposition of the data
freqs, psd_aperiodic, psd_periodic = neurodsp.aperiodic.compute_irasa(
    sig, fs, f_range=f_range, thresh=1
)
neurodsp.plts.plot_power_spectra(
    freqs,
    [psd_aperiodic, psd_periodic],
    labels=["aperiodic", "periodic"],
    title="neurodsp IRASA components",
)
# Fit the aperiodic component of the IRASA results
intercept, fit_sl = neurodsp.aperiodic.fit_irasa(freqs, psd_aperiodic)
print(f"neurodsp computed exponent: {fit_sl:1.2f}")

irasa_out = irasa(
    sig,
    fs=fs,
    band=f_range,
    psd_kwargs=dict(nperseg=4 * fs, noverlap=4 * fs // 8),
    hset_info=(1.1, 1.9, 0.05),
)
neurodsp.plts.plot_power_spectra(
    irasa_out.freqs,
    [irasa_out.aperiodic[0], irasa_out.periodic[0]],
    labels=["aperiodic", "periodic"],
    title="pyrasa IRASA components",
)
aperiodic_fit = irasa_out.fit_aperiodic_model()
fit_sl_new = aperiodic_fit.aperiodic_params.at[0, "Exponent"]
print(f"pyrasa computed exponent: {fit_sl:1.2f}")