Using Custom functions to model aperiodic signals

Aperiodic changes in a power spectrum can be modeled using a variety of different linear models. This can be easily accomplished using scipy’s curve_fit function. PyRASA aims to give users both some preset models that are similar to what has been implemented in specparam, as well as ability to fit custom models. The preset models are called fixed and knee and can be entered to the compute_slope method in the IrasaSpectrum returned by the pyrasa.irasa function. However, users can also import the class AbstractFitFun from pyrasa.utils.fit_funcs and inherit from it. Below we will illustrate how this is done and how the resultant child class can be used as a Custom fit function to model aperiodic signals.

First we simulated a signal with a spectral exponent of 1.5

import scipy.signal as dsp
from pyrasa.utils.aperiodic_utils import compute_aperiodic_model
from pyrasa.utils.fit_funcs import AbstractFitFun
from pyrasa import irasa
import numpy as np
from neurodsp.sim import sim_powerlaw
import pandas as pd

n_secs = 60
fs=1000
f_range = [1.5, 150]
exponent = -1.5

sig = sim_powerlaw(n_seconds=n_secs, fs=fs, exponent=exponent)
freqs, psd = dsp.welch(sig, fs, nperseg=int(4 * fs))
freq_logical = np.logical_and(freqs >= f_range[0], freqs <= f_range[1])
psd, freqs = psd[freq_logical], freqs[freq_logical]

Now we need to simply overwrite the func method and enter it as a fit_func to either the compute_slope function imported fromfrom pyrasa.utils.aperiodic_utils. Or we can use the method of the IrasaSpectrum returned by the pyrasa.irasa function.

class CustomFitFun(AbstractFitFun):
    log10_aperiodic = True
    log10_freq = True

    def func(self, x: np.ndarray, a: float, b: float) -> np.ndarray:
        """
        Fixed fitting function.
        Use this to model aperiodic activity without a spectral knee
        """
        y_hat = a - b * x

        return y_hat


slope_fit = compute_aperiodic_model(psd, freqs, fit_func=CustomFitFun)

slope_fit.aperiodic_params

	a	b	fit_type	ch_name
0	-1.246405	1.480257	custom	0

slope_fit_2 = irasa(sig, fs=fs, band=f_range, psd_kwargs={'nperseg': 4 * fs}).fit_aperiodic_model(fit_func=CustomFitFun)

slope_fit_2.aperiodic_params

	a	b	fit_type	ch_name
0	-1.237573	1.483932	custom	0

pd.concat([slope_fit.gof,
           slope_fit_2.gof])

	mse	R2	R2_adj.	BIC	BIC_adj.	AIC	fit_type	ch_name
0	0.006865	0.979519	0.979450	-19.046324	-25.395720	-27.823447	custom	0
0	0.000361	0.998907	0.998903	-37.861877	-44.211273	-46.638999	custom	0