Getting started#

This notebook shows the basic usage of the PyRASA package.

[1]:

from neurodsp.sim import sim_combined
from neurodsp.utils import create_times
import numpy as np
import scipy.signal as dsp
import matplotlib.pyplot as plt
import pandas as pd

We first simulate a signal with a single oscillation and a spectral slope of -1

[2]:

fs = 1000
n_seconds = 60
duration=4
overlap=0.5

sim_components = {'sim_powerlaw': {'exponent' : -1},
                  'sim_oscillation': {'freq' : 10}}


sig = sim_combined(n_seconds=n_seconds, fs=fs, components=sim_components)
times = create_times(n_seconds=n_seconds, fs=fs)

max_times = times < 1
f, axes = plt.subplots(ncols=2, figsize=(8, 4))
axes[0].plot(times[max_times], sig[max_times])
axes[0].set_ylabel('Amplitude (a.u.)')
axes[0].set_xlabel('Time (s)')
freq, psd = dsp.welch(sig, fs=fs, nperseg=duration*fs, noverlap=duration*fs*overlap)
axes[1].loglog(freq, psd)
axes[1].set_ylabel('Power (a.u.)')
axes[1].set_xlabel('Frequency (Hz)')

plt.tight_layout()

../_images/examples_basic_functionality_3_0.png

Now we can use IRASA to seperate the signal in its periodic and aperiodic components

[3]:

from pyrasa.irasa import irasa

irasa_out = irasa(sig,
                    fs=fs,
                    band=(1, 100),
                    psd_kwargs={'nperseg': duration*fs,
                                'noverlap': duration*fs*overlap
                            },
                    hset_info=(1, 2, 0.05))

f, axes = plt.subplots(ncols=2, figsize=(8, 4))
axes[0].set_title('Periodic')
axes[0].plot(irasa_out.freqs, irasa_out.periodic[0,:])
axes[0].set_ylabel('Power (a.u.)')
axes[0].set_xlabel('Frequency (Hz)')
axes[1].set_title('Aperiodic')
axes[1].loglog(irasa_out.freqs, irasa_out.aperiodic[0,:])
axes[1].set_ylabel('Power (a.u.)')
axes[1].set_xlabel('Frequency (Hz)')

f.tight_layout()

../_images/examples_basic_functionality_5_0.png

Now we can further analyse the periodic and aperiodic components using the get_peak_params and compute_slope functions, which will return pandas dataframes containing specific information about the slope or the oscillatory parameters.

[4]:

# %% get periodic stuff
irasa_out.get_peaks()

[4]:

	ch_name	cf	bw	pw
0	0	10.0	1.18842	0.490839

[5]:

# %% get aperiodic stuff
aperiodic_fit = irasa_out.fit_aperiodic_model()
aperiodic_fit.aperiodic_params

[5]:

	Offset	Exponent	fit_type	ch_name
0	-1.328678	0.984115	fixed	0

[6]:

aperiodic_fit.gof

[6]:

	mse	R2	R2_adj.	BIC	BIC_adj.	AIC	fit_type	ch_name
0	0.000363	0.997522	0.997509	-35.430795	-41.776853	-43.398668	fixed	0

But how does all of this work in practice? The beauty of IRASA lies in its simplicity. Essentially, its just up/downsampling and averaging. I will deconstruct the algorithm below to show you its inner workings and highlight potential pitfalls.

We start by simply computing a psd

[7]:

kwargs_psd = {'nperseg': duration*fs,
              'noverlap': duration*fs*overlap}

freq, psd = dsp.welch(sig, fs=fs, **kwargs_psd)

f, ax = plt.subplots(figsize=(4,4))
ax.loglog(freq, psd)
ax.set_ylabel('Power')
ax.set_xlabel('Frequency (Hz)')

[7]:

Text(0.5, 0, 'Frequency (Hz)')

../_images/examples_basic_functionality_11_1.png

Now we need to create two other psds from an up-/downsampled version of the data. Note that the data is up-/downsampled by the same factor

[8]:

import fractions

resampling_factor = 1.5

def simple_rasa(resampling_factor):
    rat = fractions.Fraction(str(resampling_factor))
    up, down = rat.numerator, rat.denominator

    # Much faster than FFT-based resampling
    data_up = dsp.resample_poly(sig, up, down, axis=-1)
    data_down = dsp.resample_poly(sig, down, up, axis=-1)

    # Calculate an up/downsampled version of the PSD using same params as original
    _, psd_up = dsp.welch(data_up, fs * resampling_factor,  **kwargs_psd)
    _, psd_dw = dsp.welch(data_down, fs / resampling_factor, **kwargs_psd)

    return psd_up, psd_dw


psd_up, psd_dw = simple_rasa(resampling_factor)

f, ax = plt.subplots(figsize=(4,4))
f_max = freq < 100
ax.loglog(freq[f_max], psd[f_max], color='k', label='original')
ax.loglog(freq[f_max], psd_up[f_max], color='b', label='upsampled (factor = 1.5)')
ax.loglog(freq[f_max], psd_dw[f_max], color='r', label='downsampled (factor = 1.5)')
ax.set_ylabel('Power')
ax.set_xlabel('Frequency (Hz)')
plt.legend()

[8]:

<matplotlib.legend.Legend at 0x111a35400>

../_images/examples_basic_functionality_13_1.png

We can see that up-/downsampling shifted the peak (oscillation) in the spectrum relative to the original data. Now we compute the geometric mean of the up-/downsampled version of the data and repeat the procedure for a different resampling factor. We can see that this creates a version of the original data with 2 peaks that are shifted around the original peak by a factor of x.

[9]:

psd_up_19, psd_dw_19 = simple_rasa(1.9)
psd_up_11, psd_dw_11 = simple_rasa(1.1)

gmean = lambda x, y : np.sqrt(x * y)

f, ax = plt.subplots(figsize=(4,4))
f_max = freq < 100
ax.loglog(freq[f_max], psd[f_max], label='original')
ax.loglog(freq[f_max], gmean(psd_up_11, psd_dw_11)[f_max], label='factor 1.1')
ax.loglog(freq[f_max], gmean(psd_up, psd_dw)[f_max], label='factor 1.5')
ax.loglog(freq[f_max], gmean(psd_up_19, psd_dw_19)[f_max], label='factor 1.9')
ax.set_ylabel('Power')
ax.set_xlabel('Frequency (Hz)')

plt.legend()

[9]:

<matplotlib.legend.Legend at 0x111e85820>

../_images/examples_basic_functionality_15_1.png

Now we can compute the median, between our 3 geometric means to obtain our aperiodic spectrum. In reality we use a bit more than 3 up-/downsampling factors, but on this data its enough to get a decent spectrum

[10]:

aperiodic_spectrum = np.median([gmean(psd_up_11, psd_dw_11),
                                gmean(psd_up, psd_dw),
                                gmean(psd_up_19, psd_dw_19)], axis=0)

f, ax = plt.subplots(figsize=(4,4))
f_max = freq < 100
ax.loglog(freq[f_max], psd[f_max], label='original')
ax.loglog(freq[f_max], aperiodic_spectrum[f_max], label='aperiodic spectrum')
ax.set_ylabel('Power')
ax.set_xlabel('Frequency (Hz)')

plt.legend()

[10]:

<matplotlib.legend.Legend at 0x111dbd5e0>

../_images/examples_basic_functionality_17_1.png

But how do we get the periodic spectrum? Its actually quite simple. We just need to subtract our aperiodic spectrum from the original spectrum.

[11]:

f, ax = plt.subplots(figsize=(4,4))
f_max = freq < 100
ax.plot(freq[f_max], psd[f_max] - aperiodic_spectrum[f_max], label='periodic spectrum')
ax.set_ylabel('Power')
ax.set_xlabel('Frequency (Hz)')

plt.legend()

[11]:

<matplotlib.legend.Legend at 0x112ae2ff0>

../_images/examples_basic_functionality_19_1.png

Now we can do slope fits on the aperiodic spectrum or do peak detection on the periodic spectrum. However, this was quite a simple spectrum and reality is usually much messsier and noisier. For instance we might get spectra that dont linearly decrease with frequency by the same value but have a deflection point (spectral knee) at which the slope starts changing. We can also deal with those using IRASA. See below for how we would do this.

[12]:

# %% Lets check the knee
knee_freq = 15
exp = 1.5
knee = knee_freq ** exp
duration=4
overlap=0.99

sim_components = {'sim_knee': {'exponent1' : -.0, 'exponent2': -1*exp, 'knee':knee },
                  'sim_oscillation': {'freq' : 10}}

sig = sim_combined(n_seconds=n_seconds, fs=fs, components=sim_components)

# %%
max_times = times < 1
f, axes = plt.subplots(ncols=2, figsize=(8, 4))
axes[0].plot(times[max_times], sig[max_times])
axes[0].set_ylabel('Amplitude (a.u.)')
axes[0].set_xlabel('Time (s)')
freq, psd = dsp.welch(sig, fs=fs, nperseg=duration*fs, noverlap=duration*fs*overlap)
axes[1].loglog(freq, psd)
axes[1].set_ylabel('Power (a.u.)')
axes[1].set_xlabel('Frequency (Hz)')

plt.tight_layout()

../_images/examples_basic_functionality_21_0.png

[13]:

irasa_out = irasa(sig,
                    fs=fs,
                    band=(.1, 100),
                    psd_kwargs={'nperseg': duration*fs,
                                'noverlap': duration*fs*overlap
                            },
                    hset_info=(1, 2.5, 0.05))

[14]:

f, axes = plt.subplots(ncols=2, figsize=(8, 4))
axes[0].plot(irasa_out.freqs, irasa_out.periodic[0,:])
axes[0].set_ylabel('Power (a.u.)')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_title('Periodic Spectrum')


axes[1].loglog(freq[1:], psd[1:], label='psd')
axes[1].loglog(irasa_out.freqs, irasa_out.aperiodic[0,:])
axes[1].set_ylabel('Power (a.u.)')
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_title('Original + \n Aperiodic Spectrum')

f.tight_layout()

../_images/examples_basic_functionality_23_0.png

[15]:

# %% get periodic stuff
irasa_out.get_peaks()

[15]:

	ch_name	cf	bw	pw
0	0	9.75	1.188926	0.501843

[16]:

aps = irasa_out.fit_aperiodic_model(fit_func='knee').aperiodic_params
aps

[16]:

	Offset	Knee	Exponent_1	Exponent_2	fit_type	Knee Frequency (Hz)	tau	ch_name
0	3.610455e-15	62.715188	0.047433	1.469384	knee	14.093899	0.011292	0

[17]:

# %% get aperiodic stuff
ap_f = irasa_out.fit_aperiodic_model(fit_func='fixed')
ap_k = irasa_out.fit_aperiodic_model(fit_func='knee')

[18]:

pd.concat([ap_f.gof, ap_k.gof])

[18]:

	mse	R2	R2_adj.	BIC	BIC_adj.	AIC	fit_type	ch_name
0	0.014386	0.887723	0.887157	-13.429702	-19.775835	-21.412631	fixed	0
0	0.000120	0.999062	0.999053	-30.116565	-42.808830	-46.082423	knee	0

[20]:

import seaborn as sns
f, ax = plt.subplots(figsize=(4,4))
sns.lineplot(data=ap_k.model, x='Frequency (Hz)', y='aperiodic_model', ax=ax)
ax.set_xscale('log')
ax.set_yscale('log')

../_images/examples_basic_functionality_28_0.png

[ ]: