diff --git a/doc/api.rst b/doc/api.rst index d91ed24a..b9ec41ef 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -196,6 +196,8 @@ Fluctuation Analyses :toctree: generated/ compute_fluctuations + compute_irasa + fit_irasa Simulations ----------- diff --git a/neurodsp/aperiodic/__init__.py b/neurodsp/aperiodic/__init__.py index f1c5c6aa..867a5dd7 100644 --- a/neurodsp/aperiodic/__init__.py +++ b/neurodsp/aperiodic/__init__.py @@ -1,3 +1,4 @@ """"Functions for calculating aperiodic properties of signals.""" from .dfa import compute_fluctuations +from .irasa import compute_irasa, fit_irasa diff --git a/neurodsp/aperiodic/irasa.py b/neurodsp/aperiodic/irasa.py new file mode 100644 index 00000000..2170b11f --- /dev/null +++ b/neurodsp/aperiodic/irasa.py @@ -0,0 +1,128 @@ +"""The IRASA method for separating periodic and aperiodic activity.""" + +import fractions + +import numpy as np + +from scipy import signal +from scipy.optimize import curve_fit + +from neurodsp.spectral import compute_spectrum, trim_spectrum + +################################################################################################### +################################################################################################### + +def compute_irasa(sig, fs=None, f_range=(1, 30), hset=None, **spectrum_kwargs): + """Separate the aperiodic and periodic components using the IRASA method. + + Parameters + ---------- + sig : 1d array + Time series. + fs : float + The sampling frequency of sig. + f_range : tuple or None + Frequency range. + hset : 1d array + Resampling factors used in IRASA calculation. + If not provided, defaults to values from 1.1 to 1.9 with an increment of 0.05. + spectrum_kwargs : dict + Optional keywords arguments that are passed to `compute_spectrum`. + + Returns + ------- + freqs : 1d array + Frequency vector. + psd_aperiodic : 1d array + The aperiodic component of the power spectrum. + psd_periodic : 1d array + The periodic component of the power spectrum. + + Notes + ----- + Irregular-Resampling Auto-Spectral Analysis (IRASA) is described in Wen & Liu (2016). + Briefly, it aims to separate 1/f and periodic components by resampling time series, and + computing power spectra, effectively averaging away any activity that is frequency specific. + + References + ---------- + Wen, H., & Liu, Z. (2016). Separating Fractal and Oscillatory Components in the Power Spectrum + of Neurophysiological Signal. Brain Topography, 29(1), 13–26. DOI: 10.1007/s10548-015-0448-0 + """ + + # Check & get the resampling factors, with rounding to avoid floating point precision errors + hset = np.arange(1.1, 1.95, 0.05) if not hset else hset + hset = np.round(hset, 4) + + # The `nperseg` input needs to be set to lock in the size of the FFT's + if 'nperseg' not in spectrum_kwargs: + spectrum_kwargs['nperseg'] = int(4 * fs) + + # Calculate the original spectrum across the whole signal + freqs, psd = compute_spectrum(sig, fs, **spectrum_kwargs) + + # Do the IRASA resampling procedure + psds = np.zeros((len(hset), *psd.shape)) + for ind, h_val in enumerate(hset): + + # Get the up-sampling / down-sampling (h, 1/h) factors as integers + rat = fractions.Fraction(str(h_val)) + up, dn = rat.numerator, rat.denominator + + # Resample signal + sig_up = signal.resample_poly(sig, up, dn, axis=-1) + sig_dn = signal.resample_poly(sig, dn, up, axis=-1) + + # Calculate the power spectrum using the same params as original + freqs_up, psd_up = compute_spectrum(sig_up, h_val * fs, **spectrum_kwargs) + freqs_dn, psd_dn = compute_spectrum(sig_dn, fs / h_val, **spectrum_kwargs) + + # Geometric mean of h and 1/h + psds[ind, :] = np.sqrt(psd_up * psd_dn) + + # Now we take the median resampled spectra, as an estimate of the aperiodic component + psd_aperiodic = np.median(psds, axis=0) + + # Subtract aperiodic from original, to get the periodic component + psd_periodic = psd - psd_aperiodic + + # Restrict spectrum to requested range + if f_range: + psds = np.array([psd_aperiodic, psd_periodic]) + freqs, (psd_aperiodic, psd_periodic) = trim_spectrum(freqs, psds, f_range) + + return freqs, psd_aperiodic, psd_periodic + + +def fit_irasa(freqs, psd_aperiodic): + """Fit the IRASA derived aperiodic component of the spectrum. + + Parameters + ---------- + freqs : 1d array + Frequency vector, in linear space. + psd_aperidic : 1d array + Power values, in linear space. + + Returns + ------- + intercept : float + Fit intercept value. + slope : float + Fit slope value. + + Notes + ----- + This fits a linear function of the form `y = ax + b` to the log-log aperiodic power spectrum. + """ + + popt, _ = curve_fit(fit_func, np.log(freqs), np.log(psd_aperiodic)) + intercept, slope = popt + + return intercept, slope + + +def fit_func(freqs, intercept, slope): + """A fit function to use for fitting IRASA separated 1/f power spectra components.""" + + return slope * freqs + intercept diff --git a/neurodsp/tests/aperiodic/test_irasa.py b/neurodsp/tests/aperiodic/test_irasa.py new file mode 100644 index 00000000..46af4664 --- /dev/null +++ b/neurodsp/tests/aperiodic/test_irasa.py @@ -0,0 +1,44 @@ +"""Tests for IRASA functions.""" + +import numpy as np + +from neurodsp.tests.settings import FS, N_SECONDS_LONG, EXP1 + +from neurodsp.sim import sim_combined +from neurodsp.spectral import compute_spectrum, trim_spectrum + +from neurodsp.aperiodic.irasa import * + +################################################################################################### +################################################################################################### + +def test_compute_irasa(tsig_comb): + + # Estimate periodic and aperiodic components with IRASA + f_range = [1, 30] + freqs, psd_ap, psd_pe = compute_irasa(tsig_comb, FS, f_range, noverlap=int(2*FS)) + + assert len(freqs) == len(psd_ap) == len(psd_pe) + + # Compute r-squared for the full model, comparing to a standard power spectrum + _, powers = trim_spectrum(*compute_spectrum(tsig_comb, FS, nperseg=int(4*FS)), f_range) + r_sq = np.corrcoef(np.array([powers, psd_ap+psd_pe]))[0][1] + assert r_sq > .95 + +def test_fit_irasa(tsig_comb): + + # Estimate periodic and aperiodic components with IRASA & fit aperiodic + freqs, psd_ap, _ = compute_irasa(tsig_comb, FS, noverlap=int(2*FS)) + b0, b1 = fit_irasa(freqs, psd_ap) + + assert round(b1) == EXP1 + assert np.abs(b0 - np.log((psd_ap)[0])) < 1 + +def test_fit_func(): + + freqs = np.arange(30) + intercept = -2 + slope = -2 + + fit = fit_func(freqs, intercept, slope) + assert (fit == slope * freqs + intercept).all() diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py index 36821d7e..eadd8c7b 100644 --- a/neurodsp/tests/conftest.py +++ b/neurodsp/tests/conftest.py @@ -6,9 +6,10 @@ import numpy as np -from neurodsp.sim import sim_oscillation, sim_powerlaw +from neurodsp.sim import sim_oscillation, sim_powerlaw, sim_combined from neurodsp.utils.sim import set_random_seed -from neurodsp.tests.settings import (FS, FS_HIGH, N_SECONDS, N_SECONDS_LONG, FREQ_SINE, +from neurodsp.tests.settings import (N_SECONDS, N_SECONDS_LONG, + FS, FS_HIGH, FREQ_SINE, FREQ1, EXP1, BASE_TEST_FILE_PATH, TEST_PLOTS_PATH) ################################################################################################### @@ -38,6 +39,13 @@ def tsig_sine_long(): yield sim_oscillation(N_SECONDS_LONG, FS, freq=FREQ_SINE, variance=None, mean=None) +@pytest.fixture(scope='session') +def tsig_comb(): + + components = {'sim_powerlaw': {'exponent' : EXP1}, + 'sim_oscillation': {'freq' : FREQ1}} + yield sim_combined(n_seconds=N_SECONDS_LONG, fs=FS, components=components) + @pytest.fixture(scope='session') def tsig_white(): diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index 652f89c8..dceb615c 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -14,12 +14,13 @@ N_SECONDS = 1.0 N_SECONDS_LONG = 10.0 -# Define frequency options for simulations +# Define parameter options for simulations FREQ1 = 10 FREQ2 = 25 FREQ_SINE = 1 FREQS_LST = [8, 12, 1] FREQS_ARR = np.array([5, 10, 15]) +EXP1 = -1 # Define error tolerance levels for test comparisons EPS = 10**(-10)