# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Aperiodic and oscillatory components of the power spectrum
#
# This notebook demonstrates how to use YASA to separate the aperiodic (= fractal = $1/f$) components of the EEG power spectra using the IRASA method.
#
# Please make sure to install the latest version of YASA first by typing the following line in your terminal or command prompt:
#
# `pip install --upgrade yasa`

import mne
import yasa
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
sns.set(style='white', font_scale=1.2)

# ## Data loading
#
# Let's load 6 min of eyes open EEG data, sampled at 200 Hz. We use the [MNE library](https://mne.tools/stable/index.html), which includes various functions and utilities for reading EEG data and electrode locations. For a description of all the file formats supported by MNE, please refer to https://mne.tools/stable/auto_tutorials/io/plot_20_reading_eeg_data.html
#
# Note that during loading, MNE automatically scale the EEG data from microVolts (uV) to Volts. However, YASA requires the data to be in microVolts (uV), and therefore we need to reverse this scaling by multiplying our data by 1e6 (1 V = 1,000,000 uV).

# +
# Load data as a MNE Raw file
raw = mne.io.read_raw_fif('data_resting_EO_200Hz_raw.fif', preload=True, verbose=0)

# Extract data, sf, and chan
data = raw.get_data() * 1e6  # Convert data from V to uV
sf = raw.info['sfreq']
chan = raw.ch_names

# Let's have a look at the data
print('Chan =', chan)
print('Sampling frequency =', sf, 'Hz')
print('Data shape =', data.shape)
print('Duration =', data.shape[1] / sf, 'seconds')
# -

# **************
#
# ## Calculate the original power spectrum
#
# For a primer on EEG spectral power please refer to https://raphaelvallat.com/bandpower.html or [08_bandpower.ipynb](08_bandpower.ipynb)

# +
from scipy.signal import welch

win = int(4 * sf)  # Window size is set to 4 seconds
freqs, psd = welch(data, sf, nperseg=win)  # Works with single or multi-channel data

print(freqs.shape, psd.shape)  # psd has shape (n_channels, n_frequencies)

# Plot
plt.plot(freqs, psd[1, :], 'k', lw=2)
plt.fill_between(freqs, psd[1, :], cmap='Spectral')
plt.xlim(1, 30)
plt.yscale('log')
sns.despine()
plt.title(chan[1])
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD log($uV^2$/Hz)');
# -

# As can be seen on the plot above, there is a strong alpha oscillatory peak around 10-11 Hz.
#
# ## Extract the aperiodic and periodic component of the spectrum
#
# The Irregular-Resampling Auto-Spectral Analysis (IRASA) method ([Wen & Liu 2016](https://doi.org/10.1007/s10548-015-0448-0)) allows to separate the
# fractal and oscillatory components in the power spectrum of EEG signals.
#
# The steps are:
#
# 1. Compute the original power spectral density (PSD).
# 2. Resample the EEG data by multiple non-integer factors and their
# reciprocals ($h$ and $1/h$).
# 3. For every pair of resampled signals, calculate the PSD and take the
# geometric mean of both. In the resulting PSD, the power associated with the
# oscillatory component is redistributed away from its original
# (fundamental and harmonic) frequencies by a frequency offset that varies
# with the resampling factor, whereas the power solely attributed to the
# fractal component remains the same power-law statistical distribution
# independent of the resampling factor.
# 4. It follows that taking the median of the PSD of the variously
# resampled signals can extract the power spectrum of the fractal component,
# and the difference between the original power spectrum and the extracted
# fractal spectrum offers an approximate estimate of the power spectrum of
# the oscillatory component.
#
# **For more details on the IRASA method, please refer to**
#
# - <NAME>., & <NAME>. (2016). Separating Fractal and Oscillatory Components in the Power Spectrum of Neurophysiological Signal. Brain Topography, 29(1), 13–26. https://doi.org/10.1007/s10548-015-0448-0
# - https://github.com/fieldtrip/fieldtrip/blob/master/specest/
#
# **For another method to extract fractal noise, please refer to**
# - https://github.com/fooof-tools/fooof
# - https://www.biorxiv.org/content/10.1101/299859v1

# Apply the IRASA technique
freqs, psd_aperiodic, psd_osc = yasa.irasa(data, sf, ch_names=chan, band=(1, 30), win_sec=4, return_fit=False)

# Plot the aperiodic component on a linear-log scale
plt.plot(freqs, psd_aperiodic[1, :], 'k', lw=2)
plt.xlim(1, 30)
plt.yscale('log')
sns.despine()
plt.title('Aperiodic component, chan = ' + chan[1])
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD log($uV^2$/Hz)');

# Plot the oscillatory component on a linear-linear scale
plt.plot(freqs, psd_osc[1, :], 'k', lw=2)
plt.xlim(1, 30)
sns.despine()
plt.title('Oscillatory component, chan = ' + chan[1])
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD log($uV^2$/Hz)');

# +
# Plot the oscillatory + aperiodic component on a linear-log scale
# psd_combined = psd_aperiodic[1, :] + psd_osc[1, :]
# plt.plot(freqs, psd_combined, 'k', lw=2)
# plt.fill_between(freqs, psd_combined, cmap='Spectral')
# plt.xlim(1, 30)
# plt.yscale('log')
# sns.despine()
# plt.title('Aperiodic + oscillatory component, chan = ' + chan[1])
# plt.xlabel('Frequency [Hz]')
# plt.ylabel('PSD log($uV^2$/Hz)');
# -

# ## Fit the fractal component (1/f)
#
# If `return_fit=True` (default), `yasa.irasa` [fits an exponential function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) to the aperiodic power spectrum and return the fit parameters (intercept, slope), the $R^2$ of the fit, and the standard deviation of the oscillatory component.
#
# The aperiodic signal, $L$, is modeled using an exponential function in semilog-power space (linear frequencies and log PSD) as:
#
# $L = a + log(F^b)$
#
# where $a$ is the intercept, $b$ is the slope, and $F$ the vector of input frequencies.

freqs, psd_aperiodic, psd_osc, fit_params = yasa.irasa(data, sf, ch_names=chan)
fit_params

StarcoderdataJupyterScriptsDedupFiltered