[ENH] - Add DFA to a new aperiodic module

doc/api.rst

@@ -183,6 +183,20 @@ Lagged Coherence
 
     compute_lagged_coherence
 
+Aperiodic Analyses
+------------------
+
+Functions and utilities in the ``aperiodic`` module, for analyzing aperiodic activity in time series.
+
+Fluctuation Analyses
+~~~~~~~~~~~~~~~~~~~~
+
+.. currentmodule:: neurodsp.aperiodic
+.. autosummary::
+    :toctree: generated/
+
+    compute_fluctuations
+
 Simulations
 -----------
 

neurodsp/aperiodic/__init__.py

@@ -0,0 +1,3 @@
+""""Functions for calculating aperiodic properties of signals."""
+
+from .dfa import compute_fluctuations

neurodsp/aperiodic/dfa.py

@@ -0,0 +1,149 @@
+"""Fluctuation analyses to measure fractal properties of time series."""
+
+import numpy as np
+import scipy as sp
+
+from neurodsp.utils.data import split_signal
+
+###################################################################################################
+###################################################################################################
+
+def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, deg=1, method='dfa'):
+    """Compute a fluctuation analysis on a signal.
+
+    Parameters
+    ----------
+    sig : 1d array
+        Time series.
+    fs : float
+        Sampling rate, in Hz.
+    n_scales : int, optional, default=10
+        Number of scales to estimate fluctuations over.
+    min_scale : float, optional, default=0.01
+        Shortest scale to compute over, in seconds.
+    max_scale : float, optional, default=1.0
+        Longest scale to compute over, in seconds.
+    deg : int, optional, default=1
+        Polynomial degree for detrending. Only used for DFA.
+
+        - 1 for regular DFA
+        - 2 or higher for generalized DFA
+    method : {'dfa', 'rs'}
+        Method to use to compute fluctuations:
+
+        - 'dfa' : detrended fluctuation
+        - 'rs' : rescaled range
+
+    Returns
+    -------
+    t_scales : 1d array
+        Time-scales over which fluctuation measures were computed.
+    fluctuations : 1d array
+        Average fluctuation at each scale.
+    exp : float
+        Slope of line in log-log when plotting time scales against fluctuations.
+        This is the alpha value for DFA, or the Hurst exponent for rescaled range.
+
+    Notes
+    -----
+    These analyses compute fractal properties by analyzing fluctuations across windows.
+
+    Overall, these approaches involve dividing the time-series into non-overlapping
+    windows at log-spaced scales and computing a fluctuation measure across windows.
+    The relationship of this fluctuation measure across window sizes provides
+    information on the fractal properties of a signal.
+
+    Available measures are:
+
+    - DFA: detrended fluctuation analysis
+        - computes ordinary least squares fits across signal windows
+    - RS: rescaled range
+        - computes the range of signal windows, divided by the standard deviation
+    """
+
+    # Get log10 equi-spaced scales and translate that into window lengths
+    t_scales = np.logspace(np.log10(min_scale), np.log10(max_scale), n_scales)
+    win_lens = np.round(t_scales * fs).astype('int')
+
+    # Check that all window sizes are fit-able
+    if np.any(win_lens <= 1):
+        raise ValueError("Some of window sizes are too small to run. "
+                         "Try updating `min_scale` to a value that works "
+                         "better for the current sampling rate.")
+
+    # Step through each scale and measure fluctuations
+    fluctuations = np.zeros_like(t_scales)
+    for idx, win_len in enumerate(win_lens):
+
+        if method == 'dfa':
+            fluctuations[idx] = compute_detrended_fluctuation(sig, win_len=win_len, deg=deg)
+        elif method == 'rs':
+            fluctuations[idx] = compute_rescaled_range(sig, win_len=win_len)
+        else:
+            raise ValueError('Fluctuation method not understood.')
+
+    # Calculate the relationship between between fluctuations & time scales
+    exp = np.polyfit(np.log10(t_scales), np.log10(fluctuations), deg=1)[0]
+
+    return t_scales, fluctuations, exp
+
+
+def compute_rescaled_range(sig, win_len):
+    """Compute rescaled range of a given time series at a given scale.
+
+    Parameters
+    ----------
+    sig : 1d array
+        Time series.
+    win_len : int
+        Window length for each rescaled range computation, in samples.
+
+    Returns
+    -------
+    rs : float
+        Average rescaled range over windows.
+    """
+
+    # Demean signal
+    sig = sig - np.mean(sig)
+
+    # Calculate cumulative sum of the signal & split the signal into segments
+    segments = split_signal(sp.cumsum(sig), win_len).T
+
+    # Calculate rescaled range as range divided by standard deviation (of non-cumulative signal)
+    rs_win = np.ptp(segments, axis=0) / np.std(split_signal(sig, win_len).T, axis=0)
+
+    # Take the mean across windows
+    rs = np.mean(rs_win)
+
+    return rs
+
+
+def compute_detrended_fluctuation(sig, win_len, deg=1):
+    """Compute detrended fluctuation of a time series at the given window length.
+
+    Parameters
+    ----------
+    sig : 1d array
+        Time series.
+    win_len : int
+        Window length for each detrended fluctuation fit, in samples.
+    deg : int, optional, default=1
+        Polynomial degree for detrending.
+
+    Returns
+    -------
+    det_fluc : float
+        Measured detrended fluctuation, as the average error fits of the window.
+    """
+
+    # Calculate cumulative sum of the signal & split the signal into segments
+    segments = split_signal(sp.cumsum(sig - np.mean(sig)), win_len).T
+
+    # Calculate local trend, as the line of best fit within the time window
+    _, fluc, _, _, _ = np.polyfit(np.arange(win_len), segments, deg=deg, full=True)
+
+    # Convert to root-mean squared error, from squared error
+    det_fluc = np.mean((fluc / win_len))**0.5
+
+    return det_fluc

neurodsp/tests/aperiodic/test_dfa.py

@@ -0,0 +1,57 @@
+"""Tests for fractal analysis using fluctuation measures."""
+
+from pytest import raises
+
+import numpy as np
+
+from neurodsp.sim import sim_powerlaw
+from neurodsp.tests.settings import FS, FS_HIGH
+
+from neurodsp.aperiodic.dfa import (compute_fluctuations, compute_rescaled_range,
+                                    compute_detrended_fluctuation)
+
+###################################################################################################
+###################################################################################################
+
+def test_compute_fluctuations(tsig):
+
+    t_scales, flucs, exp = compute_fluctuations(tsig, 500)
+    assert len(t_scales) == len(flucs)
+
+    # Check error for if the settings create window lengths that are too short
+    with raises(ValueError):
+        t_scales, flucs, exp = compute_fluctuations(tsig, 100)
+
+    # Check error for nonsense method input
+    with raises(ValueError):
+        t_scales, flucs, exp = compute_fluctuations(tsig, 500, method='nope.')
+
+def test_compute_fluctuations_dfa(tsig_white, tsig_brown):
+    """Tests fluctuation analysis for the DFA method."""
+
+    # Test white noise: expected DFA of 0.5
+    t_scales, flucs, exp = compute_fluctuations(tsig_white, FS_HIGH, method='dfa')
+    assert np.isclose(exp, 0.5, atol=0.1)
+
+    # Test brown noise: expected DFA of 1.5
+    t_scales, flucs, exp = compute_fluctuations(tsig_brown, FS_HIGH, method='dfa')
+    assert np.isclose(exp, 1.5, atol=0.1)
+
+def test_compute_fluctuations_rs(tsig_white):
+    """Tests fluctuation analysis for the RS method."""
+
+    # Test white noise: expected RS of 0.5
+    t_scales, flucs, exp = compute_fluctuations(tsig_white, FS_HIGH, method='rs')
+    assert np.isclose(exp, 0.5, atol=0.1)
+
+def test_compute_rescaled_range(tsig):
+
+    out = compute_rescaled_range(tsig, 10)
+    assert isinstance(out, float)
+    assert out >= 0
+
+def test_compute_detrended_fluctuation(tsig):
+
+    out = compute_detrended_fluctuation(tsig, 10)
+    assert isinstance(out, float)
+    assert out >= 0

neurodsp/tests/conftest.py

@@ -6,12 +6,11 @@
 
 import numpy as np
 
+from neurodsp.sim import sim_oscillation, sim_powerlaw
 from neurodsp.utils.sim import set_random_seed
-from neurodsp.tests.settings import (FS, N_SECONDS, N_SECONDS_LONG, FREQ_SINE,
+from neurodsp.tests.settings import (FS, FS_HIGH, N_SECONDS, N_SECONDS_LONG, FREQ_SINE,
                                      BASE_TEST_FILE_PATH, TEST_PLOTS_PATH)
 
-from neurodsp.sim import sim_oscillation
-
 ###################################################################################################
 ###################################################################################################
 
@@ -39,6 +38,16 @@ def tsig_sine_long():
 
 	yield sim_oscillation(N_SECONDS_LONG, FS, freq=FREQ_SINE, variance=None, mean=None)
 
+@pytest.fixture(scope='session')
+def tsig_white():
+
+    yield sim_powerlaw(N_SECONDS_LONG, FS_HIGH, exponent=0)
+
+@pytest.fixture(scope='session')
+def tsig_brown():
+
+    yield sim_powerlaw(N_SECONDS_LONG, FS_HIGH, exponent=-2)
+
 @pytest.fixture(scope='session', autouse=True)
 def check_dir():
     """Once, prior to session, this will clear and re-initialize the test file directories."""

neurodsp/tests/settings.py

@@ -10,6 +10,7 @@
 
 # Define general settings for simulations & tests
 FS = 100
+FS_HIGH = 1000
 N_SECONDS = 1.0
 N_SECONDS_LONG = 10.0
 

neurodsp/utils/data.py

@@ -67,7 +67,7 @@ def create_samples(n_samples, start_val=0):
 
 
 def split_signal(sig, n_samples):
-    """Split a signal into non-overlapping chunks.
+    """Split a signal into non-overlapping segments.
 
     Parameters
     ----------
@@ -78,17 +78,17 @@ def split_signal(sig, n_samples):
 
     Returns
     -------
-    chunks : 2d array
-        The signal, split into chunks, with shape [n_chunks, chunk_size].
+    segs : 2d array
+        The signal, split into segments, with shape [n_segment, segment_size].
 
     Notes
     -----
-    If the signal does not divide evenly into the number of chunks, this approach
-    will truncate the signal, returning the maximum number of chunks, and dropping
+    If the signal does not divide evenly into the number of segments, this approach
+    will truncate the signal, returning the maximum number of segments, and dropping
     any leftover samples.
     """
 
-    n_chunks = int(np.floor(len(sig) / float(n_samples)))
-    chunks = np.reshape(sig[:int(n_chunks * n_samples)], (n_chunks, int(n_samples)))
+    n_segments = int(np.floor(len(sig) / float(n_samples)))
+    segments = np.reshape(sig[:int(n_segments * n_samples)], (n_segments, int(n_samples)))
 
-    return chunks
+    return segments