Merge pull request #212 from neurodsp-tools/irasa

[ENH] - Add IRASA

doc/api.rst

@@ -196,6 +196,8 @@ Fluctuation Analyses
     :toctree: generated/
 
     compute_fluctuations
+    compute_irasa
+    fit_irasa
 
 Simulations
 -----------

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

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

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()

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():
 

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)