[DOC] - Add tutorials for aperiodic sub-module

neurodsp/aperiodic/autocorr.py

@@ -31,9 +31,9 @@ def compute_autocorr(sig, max_lag=1000, lag_step=1, demean=True):
         sig = sig - sig.mean()
 
     autocorrs = np.correlate(sig, sig, "full")[len(sig)-1:]
-    autocorrs = autocorrs[:max_lag] / autocorrs[0]
+    autocorrs = autocorrs[:max_lag+1] / autocorrs[0]
     autocorrs = autocorrs[::lag_step]
 
-    timepoints = np.arange(0, max_lag, lag_step)
+    timepoints = np.arange(0, max_lag+1, lag_step)
 
     return timepoints, autocorrs

neurodsp/aperiodic/irasa.py

@@ -1,4 +1,4 @@
-"""The IRASA method for separating periodic and aperiodic activity."""
+"""The IRASA algorithm for separating periodic and aperiodic activity."""
 
 import fractions
 
@@ -12,20 +12,23 @@
 ###################################################################################################
 ###################################################################################################
 
-def compute_irasa(sig, fs=None, f_range=(1, 30), hset=None, **spectrum_kwargs):
-    """Separate the aperiodic and periodic components using the IRASA method.
+def compute_irasa(sig, fs, f_range=None, hset=None, thresh=None, **spectrum_kwargs):
+    """Separate aperiodic and periodic components using IRASA.
 
     Parameters
     ----------
     sig : 1d array
         Time series.
     fs : float
         The sampling frequency of sig.
-    f_range : tuple or None
-        Frequency range.
-    hset : 1d array
+    f_range : tuple, optional
+        Frequency range to restrict the analysis to.
+    hset : 1d array, optional
         Resampling factors used in IRASA calculation.
         If not provided, defaults to values from 1.1 to 1.9 with an increment of 0.05.
+    thresh : float, optional
+        A relative threshold to apply when separating out periodic components.
+        The threshold is defined in terms of standard deviations of the original spectrum.
     spectrum_kwargs : dict
         Optional keywords arguments that are passed to `compute_spectrum`.
 
@@ -40,9 +43,9 @@ def compute_irasa(sig, fs=None, f_range=(1, 30), hset=None, **spectrum_kwargs):
 
     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.
+    Irregular-Resampling Auto-Spectral Analysis (IRASA) is an algorithm that aims to
+    separate 1/f and periodic components by resampling time series and computing power spectra,
+    averaging away any activity that is frequency specific to isolate the aperiodic component.
 
     References
     ----------
@@ -73,19 +76,25 @@ def compute_irasa(sig, fs=None, f_range=(1, 30), hset=None, **spectrum_kwargs):
         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
+        # 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
+        # Calculate the 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
+    # 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
 
+    # Apply a relative threshold for tuning which activity is labelled as periodic
+    if thresh is not None:
+        sub_thresh = np.where(psd_periodic - psd_aperiodic < thresh * np.std(psd))[0]
+        psd_periodic[sub_thresh] = 0
+        psd_aperiodic[sub_thresh] = psd[sub_thresh]
+
     # Restrict spectrum to requested range
     if f_range:
         psds = np.array([psd_aperiodic, psd_periodic])

neurodsp/tests/aperiodic/test_autocorr.py

@@ -7,5 +7,6 @@
 
 def test_compute_autocorr(tsig):
 
-    timepoints, autocorrs = compute_autocorr(tsig)
-    assert len(timepoints) == len(autocorrs)
+    max_lag = 500
+    timepoints, autocorrs = compute_autocorr(tsig, max_lag=max_lag)
+    assert len(timepoints) == len(autocorrs) == max_lag + 1

neurodsp/tests/aperiodic/test_irasa.py

@@ -17,7 +17,6 @@ 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
@@ -28,7 +27,7 @@ def test_compute_irasa(tsig_comb):
 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))
+    freqs, psd_ap, _ = compute_irasa(tsig_comb, FS, f_range=[1, 30], noverlap=int(2*FS))
     b0, b1 = fit_irasa(freqs, psd_ap)
 
     assert round(b1) == EXP1

tutorials/aperiodic/README.txt

@@ -0,0 +1,4 @@
+Analyzing Aperiodic Signal Properties
+-------------------------------------
+
+Tutorials for the aperiodic module.

tutorials/aperiodic/plot_Autocorr.py

@@ -0,0 +1,168 @@
+"""
+Autocorrelation Measures
+========================
+
+Apply autocorrelation measures to neural signals.
+
+Autocorrelation is the correlation of a signal with delayed copies of itself.
+Autocorrelation measures can be useful to investigate properties of neural signals.
+
+This tutorial covers ``neurodsp.aperiodic.autocorr``.
+"""
+
+###################################################################################################
+
+# sphinx_gallery_thumbnail_number = 1
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from neurodsp.sim import sim_powerlaw, sim_oscillation
+
+# Import the function for computing autocorrelation
+from neurodsp.aperiodic import compute_autocorr
+
+###################################################################################################
+# Autocorrelation Measures
+# ------------------------
+#
+# Autocorrelation is computed as the correlation between the original signal and delayed
+# copies, across different lags.
+#
+# The result is a measure of how correlated a signal is to itself, across time,
+# and the timescale of autocorrelation.
+#
+# Algorithm Settings
+# ~~~~~~~~~~~~~~~~~~
+#
+# Settings for computing autocorrelation are:
+#
+# - `max_lag` : the maximum lag to compute autocorrelation for
+# - `lag_step` : the step size to advance across when computing autocorrelation
+#
+# Both parameters are defined in samples, with defaults of using a step size of 1 sample,
+# stepping up to a maximum lag of 1000 samples.
+#
+# Autocorrelation can be computed with :func:`~.compute_autocorr` function, which
+# returns `timepoints` at which autocorrelation was calculated, and `autocorrs`, which are
+# the resulting correlation coefficients.
+#
+
+###################################################################################################
+# Autocorrelation of Periodic Signals
+# -----------------------------------
+#
+# First, let's examine periodic signals, specifically sine waves.
+#
+# Note that periodic signals are by definition rhythmic. This means that we should expect
+# the autocorrelation to also have a rhythmic pattern across time.
+#
+# Imagine, for example, moving one sine wave across another. At some points, when they are
+# in phase, they will line up exactly (have a high correlation), while at others, when they
+# are out of phase, they will have either low or anti-correlation.
+#
+
+###################################################################################################
+
+# Simulation settings
+n_seconds = 10
+fs = 1000
+
+# Define the frequencies for the sinusoids
+freq1 = 10
+freq2 = 20
+
+# Simulate sinusoids
+sig_osc1 = sim_oscillation(n_seconds, fs, freq1)
+sig_osc2 = sim_oscillation(n_seconds, fs, freq2)
+
+###################################################################################################
+
+# Compute autocorrelation on the periodic time series
+timepoints_osc1, autocorrs_osc1 = compute_autocorr(sig_osc1)
+timepoints_osc2, autocorrs_osc2 = compute_autocorr(sig_osc2)
+
+###################################################################################################
+# Autocorrelation of a sinusoidal signal
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+# The autocorrelation of the first sine wave is plotted below.
+#
+
+###################################################################################################
+
+# Plot autocorrelations
+_, ax = plt.subplots(figsize=(6, 4))
+ax.plot(timepoints_osc1, autocorrs_osc1)
+ax.set(xlabel='lag (samples)', ylabel='autocorrelation');
+
+###################################################################################################
+#
+# As we can see, the autocorrelation of a sinusoid is itself a sinusoid!
+#
+# This reflects that a sinusoid related to itself will oscillate between being
+# positively and negatively correlated with itself.
+#
+# Next, let's compare the autocorrelation of different sinusoids.
+#
+
+###################################################################################################
+
+# Plot autocorrelations for two different sinusoids
+_, ax = plt.subplots(figsize=(6, 4))
+ax.plot(timepoints_osc1, autocorrs_osc1, alpha=0.75, label='10 Hz')
+ax.plot(timepoints_osc2, autocorrs_osc2, alpha=0.75, label='20 Hz')
+ax.set(xlabel='lag (samples)', ylabel='autocorrelation')
+plt.legend(loc='upper right')
+
+###################################################################################################
+#
+# In the above, we can see that the autocorrelation of sinusoids with different frequencies
+# leads to autocorrelation results with different timescales.
+#
+# If you compare to the number of samples on the x-axis, keeping in mind the sampling
+# rate (1000 Hz), you can check that the autocorrelation of a sinusoidal signal is
+# a sinusoid of the same frequency.
+#
+
+###################################################################################################
+# Autocorrelation of Aperiodic Signals
+# ------------------------------------
+#
+# Next, lets consider the autocorrelation of aperiodic signals.
+#
+# Here we will use white noise, as an example of a signal without autocorrelation, and
+# pink noise, which does, by definition, have temporal auto-correlations.
+#
+
+###################################################################################################
+
+# Simulate a white noise signal
+sig_wn = sim_powerlaw(n_seconds, fs, exponent=0)
+
+# Simulate a pink noise signal
+sig_pn = sim_powerlaw(n_seconds, fs, exponent=-1)
+
+###################################################################################################
+
+# Compute autocorrelation on the aperiodic time series
+timepoints_wn, autocorrs_wn = compute_autocorr(sig_wn)
+timepoints_pn, autocorrs_pn = compute_autocorr(sig_pn)
+
+###################################################################################################
+
+# Plot the autocorrelations of the aperiodic signals
+_, ax = plt.subplots(figsize=(5, 4))
+ax.plot(timepoints_wn, autocorrs_wn, label='White Noise')
+ax.plot(timepoints_pn, autocorrs_pn, label='Pink Noise')
+ax.set(xlabel="lag (samples)", ylabel="autocorrelation")
+plt.legend()
+
+###################################################################################################
+#
+# In the above, we can see that for white noise, the autocorrelation is only high at
+# a lag of 0 samples, with all other lags being uncorrelated.
+#
+# By comparison, the pink noise signal has a pattern of decreasing autocorrelation
+# across increasing lags. This is characteristic of powerlaw data.
+#