From 7fb9f338d5b3b0336f8da4475edde38a9ec30e7e Mon Sep 17 00:00:00 2001
From: rdgao <r.dg.gao@gmail.com>
Date: Thu, 5 Sep 2019 11:10:57 +0200
Subject: [PATCH 01/11] added dfa in aperiodic module

---
 neurodsp/aperiodic/dfa.py | 109 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 109 insertions(+)
 create mode 100644 neurodsp/aperiodic/dfa.py

diff --git a/neurodsp/aperiodic/dfa.py b/neurodsp/aperiodic/dfa.py
new file mode 100644
index 00000000..4140e973
--- /dev/null
+++ b/neurodsp/aperiodic/dfa.py
@@ -0,0 +1,109 @@
+import numpy as np
+import scipy as sp
+
+def dfa(data, fs, n_scales=10, min_scale=0.01, max_scale=1.0, deg=1, method='dfa'):
+    """ Perform Detrended Fluctuation Analysis (DFA) on a given time series by
+    dividing the data into non-overlapping windows at log-spaced scales and
+    computing the mean RMS of fits over all window at each scale.
+
+    Parameters
+    ----------
+    data : array, 1-D
+        Data to compute DFA over.
+    fs : float, Hz
+        Sampling frequency of data.
+    n_scales : int (default=10)
+        Number of scales to estimate fluctuations over.
+    min_scale : float, seconds (default=0.01s)
+        Shortest scale to compute DFA over.
+    max_scale : float, seconds (default=1.0s)
+        Longest scale to compute DFA over.
+    deg : int (default=1)
+        Polynomial degree for detrending. 1 for regular DFA, 2 or higher for
+        generalized DFA.
+    method : str (default='dfa')
+        Method to compute per-window fluctuation.
+        'dfa' : detrended fluctuation (OLS fit of data window). Empirically, DFA
+            seems to work for a larger range of exponent values.
+        'rs' : rescaled range (range divided by std)
+
+    Returns
+    -------
+    t_scales : array, 1-D
+        Time-scales over which detrended fluctuations was performedself.
+    DFs : array, 1-D
+        Mean RMS (fluctuations) at each scale.
+    alpha : float
+        Slope of line in loglog when plotting t_scales against DFs (Hurst exp.).
+    """
+    # 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')
+    # get culmulative sum
+    data_walk = sp.cumsum(data - np.mean(data))
+    DFs = np.zeros_like(t_scales)
+    # step through each scale and get RMS
+    for idx, win_len in enumerate(win_lens):
+        if method is 'dfa':
+            DFs[idx] = compute_DF(data_walk, win_len=win_len, deg=deg)
+        elif method is 'rs':
+            DFs[idx] = compute_RS(data_walk, win_len=win_len)
+
+    # get DFA exponent
+    alpha = np.polyfit(np.log10(t_scales), np.log10(DFs), deg=1)[0]
+    return t_scales, DFs, alpha
+
+def compute_RS(data, win_len):
+    """Compute rescaled range of a given time-series at a given scale.
+
+    Parameters
+    ----------
+    data : array, 1-D
+        Data to compute R/S over.
+    win_len : int
+        Window length for each R/S computation.
+
+    Returns
+    -------
+    x : float
+        Mean R/S over all sub-windows of data.
+
+    """
+    # gather all windows
+    n_win = int(np.floor(len(data) / win_len))
+    # vectorize the data so we can call math functions in one go
+    Z_rect = np.reshape(data[:n_win * win_len], (n_win, win_len)).T
+    # get back the data by taking the derivative of data_walk
+    X = np.concatenate((data[:1],np.diff(data)))
+    X_rect = np.reshape(X[:n_win * win_len], (n_win, win_len)).T
+
+    RS_seg = np.ptp(Z_rect, axis=0)/np.std(X_rect, axis=0)
+    return np.mean(RS_seg)
+
+
+def compute_DF(data, win_len, deg=1):
+    """ Compute detrended fluctuation of the data at the given window length.
+
+    Parameters
+    ----------
+    data : array, 1-D
+        Data to compute DFA over.
+    win_len : int
+        Window length for each DF fit.
+    deg : int (default=1)
+        Polynomial degree for detrending.
+
+    Returns
+    -------
+    x : float
+        mean RMS across fits.
+    """
+    # gather all windows
+    n_win = int(np.floor(len(data) / win_len))
+    # vectorize the data so we can call np.polyfit in one go
+    data_rect = np.reshape(data[:n_win * win_len], (n_win, win_len)).T
+    # fitting
+    coef, fluc, _, _, _ = np.polyfit(
+        np.arange(win_len), data_rect, deg=deg, full=True)
+    # fluc cntains sum of squared error, we want mean RMS across fits
+    return np.mean((fluc / win_len)**0.5)

From 4943acbd15639fbafede3fd28ea62d0555bcbe45 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Thu, 5 Sep 2019 10:57:10 -0700
Subject: [PATCH 02/11] Some refactoring, to fit NDSP

---
 neurodsp/aperiodic/dfa.py | 172 +++++++++++++++++++++++---------------
 1 file changed, 104 insertions(+), 68 deletions(-)

diff --git a/neurodsp/aperiodic/dfa.py b/neurodsp/aperiodic/dfa.py
index 4140e973..b2429331 100644
--- a/neurodsp/aperiodic/dfa.py
+++ b/neurodsp/aperiodic/dfa.py
@@ -1,109 +1,145 @@
+"""Fluctuation analyses to measure fractal properties of time series."""
+
 import numpy as np
 import scipy as sp
 
-def dfa(data, fs, n_scales=10, min_scale=0.01, max_scale=1.0, deg=1, method='dfa'):
-    """ Perform Detrended Fluctuation Analysis (DFA) on a given time series by
-    dividing the data into non-overlapping windows at log-spaced scales and
-    computing the mean RMS of fits over all window at each scale.
+###################################################################################################
+###################################################################################################
+
+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
     ----------
-    data : array, 1-D
-        Data to compute DFA over.
+    sig : 1d array
+        Time series.
     fs : float, Hz
-        Sampling frequency of data.
-    n_scales : int (default=10)
+        Sampling rate, in Hz.
+    n_scales : int, optional, default=10
         Number of scales to estimate fluctuations over.
-    min_scale : float, seconds (default=0.01s)
-        Shortest scale to compute DFA over.
-    max_scale : float, seconds (default=1.0s)
-        Longest scale to compute DFA over.
-    deg : int (default=1)
-        Polynomial degree for detrending. 1 for regular DFA, 2 or higher for
-        generalized DFA.
-    method : str (default='dfa')
-        Method to compute per-window fluctuation.
-        'dfa' : detrended fluctuation (OLS fit of data window). Empirically, DFA
-            seems to work for a larger range of exponent values.
-        'rs' : rescaled range (range divided by std)
+    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 : array, 1-D
-        Time-scales over which detrended fluctuations was performedself.
-    DFs : array, 1-D
-        Mean RMS (fluctuations) at each scale.
-    alpha : float
-        Slope of line in loglog when plotting t_scales against DFs (Hurst exp.).
+    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 loglog when plotting t_scales against fluctuations.
+        This is the alpha value for DFA, or the Hurst exponent for rescaled range.
+
+    Notes
+    -----
+    These analysis 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.
+
+    Measures available 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
+
+    Empirically, DFA seems to work for a larger range of exponent values.
     """
-    # get log10 equi-spaced scales and translate that into window lengths
+
+    # 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')
-    # get culmulative sum
-    data_walk = sp.cumsum(data - np.mean(data))
-    DFs = np.zeros_like(t_scales)
-    # step through each scale and get RMS
+
+    # Calculate culmulative sum of the signal
+    sig_walk = sp.cumsum(sig - np.mean(sig))
+
+    # Step through each scale and get RMS
+    fluctuations = np.zeros_like(t_scales)
     for idx, win_len in enumerate(win_lens):
         if method is 'dfa':
-            DFs[idx] = compute_DF(data_walk, win_len=win_len, deg=deg)
+            fluctuations[idx] = compute_detrended_fluctuation(sig_walk, win_len=win_len, deg=deg)
         elif method is 'rs':
-            DFs[idx] = compute_RS(data_walk, win_len=win_len)
+            fluctuations[idx] = compute_rescaled_range(sig_walk, win_len=win_len)
+        else:
+            raise ValueError('Fluctuation method not understood.')
 
-    # get DFA exponent
-    alpha = np.polyfit(np.log10(t_scales), np.log10(DFs), deg=1)[0]
-    return t_scales, DFs, alpha
+    # Calculate the relationship between between fluctuations & time scales
+    exp = np.polyfit(np.log10(t_scales), np.log10(fluctuations), deg=1)[0]
 
-def compute_RS(data, win_len):
+    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
     ----------
-    data : array, 1-D
-        Data to compute R/S over.
+    sig : 1d array
+        Time series.
     win_len : int
-        Window length for each R/S computation.
+        Window length for each rescaled range computation.
 
     Returns
     -------
-    x : float
-        Mean R/S over all sub-windows of data.
-
+    rs : float
+        Average rescaled range over windows.
     """
-    # gather all windows
-    n_win = int(np.floor(len(data) / win_len))
-    # vectorize the data so we can call math functions in one go
-    Z_rect = np.reshape(data[:n_win * win_len], (n_win, win_len)).T
-    # get back the data by taking the derivative of data_walk
-    X = np.concatenate((data[:1],np.diff(data)))
+
+    # Gather windows & vectorize, so we can call math functions in one go
+    n_win = int(np.floor(len(sig) / win_len))
+    sig_rect = np.reshape(sig[:n_win * win_len], (n_win, win_len)).T
+
+    # get back the sig by taking the derivative of sig_walk
+    X = np.concatenate((sig[:1], np.diff(sig)))
     X_rect = np.reshape(X[:n_win * win_len], (n_win, win_len)).T
 
-    RS_seg = np.ptp(Z_rect, axis=0)/np.std(X_rect, axis=0)
-    return np.mean(RS_seg)
+    # Calculate rescaled range as range divided by std, and take mean across windows
+    rs_win = np.ptp(sig_rect, axis=0) / np.std(X_rect, axis=0)
+    rs = np.mean(rs_win)
 
+    return rs
 
-def compute_DF(data, win_len, deg=1):
-    """ Compute detrended fluctuation of the data at the given window length.
+
+def compute_detrended_fluctuation(sig, win_len, deg=1):
+    """Compute detrended fluctuation of a time-series at the given window length.
 
     Parameters
     ----------
-    data : array, 1-D
-        Data to compute DFA over.
+    sig : 1d array
+        Time series.
     win_len : int
-        Window length for each DF fit.
-    deg : int (default=1)
+        Window length for each detrended fluctuation fit.
+    deg : int, optional, default=1
         Polynomial degree for detrending.
 
     Returns
     -------
-    x : float
-        mean RMS across fits.
+    det_fluc : float
+        Measured detrended fluctuaiton, as the average error fits of the window.
     """
-    # gather all windows
-    n_win = int(np.floor(len(data) / win_len))
-    # vectorize the data so we can call np.polyfit in one go
-    data_rect = np.reshape(data[:n_win * win_len], (n_win, win_len)).T
-    # fitting
-    coef, fluc, _, _, _ = np.polyfit(
-        np.arange(win_len), data_rect, deg=deg, full=True)
-    # fluc cntains sum of squared error, we want mean RMS across fits
-    return np.mean((fluc / win_len)**0.5)
+
+    # Gather windows & vectorize, so we can call math functions in one go
+    n_win = int(np.floor(len(sig) / win_len))
+    sig_rect = np.reshape(sig[:n_win * win_len], (n_win, win_len)).T
+
+    # Calculate local trend, as the line of best fit within the time window
+    _, fluc, _, _, _ = np.polyfit(np.arange(win_len), sig_rect, 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

From 5d38ab06fd014241cba7e8cdcfed9887a92460b9 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Thu, 5 Sep 2019 10:57:28 -0700
Subject: [PATCH 03/11] Add smoke tests aperiodic

---
 neurodsp/tests/test_aperiodic_dfa.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)
 create mode 100644 neurodsp/tests/test_aperiodic_dfa.py

diff --git a/neurodsp/tests/test_aperiodic_dfa.py b/neurodsp/tests/test_aperiodic_dfa.py
new file mode 100644
index 00000000..c87fe81c
--- /dev/null
+++ b/neurodsp/tests/test_aperiodic_dfa.py
@@ -0,0 +1,20 @@
+"""Tests for fractal analysis using fluctuation measures."""
+
+from neurodsp.tests.settings import FS
+
+from neurodsp.aperiodic.dfa import *
+
+###################################################################################################
+###################################################################################################
+
+def test_compute_fluctuations(tsig):
+
+    t_scales, flucs, exp = compute_fluctuations(tsig, FS)
+
+def test_compute_rescaled_range(tsig):
+
+    rs = compute_rescaled_range(tsig, 10)
+
+def test_compute_detrended_fluctuation(tsig):
+
+    out = compute_detrended_fluctuation(tsig, 10)

From f61f339dbc1782f7c5cd6cee6e3a5cd4758622bb Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Wed, 9 Oct 2019 15:41:23 -0700
Subject: [PATCH 04/11] update DFA with review comments & associated updates

---
 neurodsp/aperiodic/__init__.py |  3 ++
 neurodsp/aperiodic/dfa.py      | 52 ++++++++++++++++------------------
 neurodsp/burst/__init__.py     |  1 +
 neurodsp/utils/data.py         | 16 +++++------
 4 files changed, 37 insertions(+), 35 deletions(-)
 create mode 100644 neurodsp/aperiodic/__init__.py

diff --git a/neurodsp/aperiodic/__init__.py b/neurodsp/aperiodic/__init__.py
new file mode 100644
index 00000000..62af65c0
--- /dev/null
+++ b/neurodsp/aperiodic/__init__.py
@@ -0,0 +1,3 @@
+""""Functions for calculating aperiodic properties of signals."""
+
+from .dfa import compute_fluctuations
\ No newline at end of file
diff --git a/neurodsp/aperiodic/dfa.py b/neurodsp/aperiodic/dfa.py
index b2429331..4dc04e56 100644
--- a/neurodsp/aperiodic/dfa.py
+++ b/neurodsp/aperiodic/dfa.py
@@ -3,6 +3,8 @@
 import numpy as np
 import scipy as sp
 
+from neurodsp.utils.data import split_signal
+
 ###################################################################################################
 ###################################################################################################
 
@@ -39,12 +41,12 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     fluctuations : 1d array
         Average fluctuation at each scale.
     exp : float
-        Slope of line in loglog when plotting t_scales against fluctuations.
+        Slope of line in loglog when plotting time scales against fluctuations.
         This is the alpha value for DFA, or the Hurst exponent for rescaled range.
 
     Notes
     -----
-    These analysis compute fractal properties by analyzing fluctuations across windows.
+    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.
@@ -52,28 +54,25 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     information on the fractal properties of a signal.
 
     Measures available 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
-
-    Empirically, DFA seems to work for a larger range of exponent values.
     """
 
     # 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')
 
-    # Calculate culmulative sum of the signal
-    sig_walk = sp.cumsum(sig - np.mean(sig))
-
-    # Step through each scale and get RMS
+    # Step through each scale and measure fluctuations
     fluctuations = np.zeros_like(t_scales)
     for idx, win_len in enumerate(win_lens):
+
         if method is 'dfa':
-            fluctuations[idx] = compute_detrended_fluctuation(sig_walk, win_len=win_len, deg=deg)
+            fluctuations[idx] = compute_detrended_fluctuation(sig, win_len=win_len, deg=deg)
         elif method is 'rs':
-            fluctuations[idx] = compute_rescaled_range(sig_walk, win_len=win_len)
+            fluctuations[idx] = compute_rescaled_range(sig, win_len=win_len)
         else:
             raise ValueError('Fluctuation method not understood.')
 
@@ -84,14 +83,14 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
 
 
 def compute_rescaled_range(sig, win_len):
-    """Compute rescaled range of a given time-series at a given scale.
+    """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.
+        Window length for each rescaled range computation, in samples.
 
     Returns
     -------
@@ -99,45 +98,44 @@ def compute_rescaled_range(sig, win_len):
         Average rescaled range over windows.
     """
 
-    # Gather windows & vectorize, so we can call math functions in one go
-    n_win = int(np.floor(len(sig) / win_len))
-    sig_rect = np.reshape(sig[:n_win * win_len], (n_win, win_len)).T
+    # 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
 
-    # get back the sig by taking the derivative of sig_walk
-    X = np.concatenate((sig[:1], np.diff(sig)))
-    X_rect = np.reshape(X[:n_win * win_len], (n_win, 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)
 
-    # Calculate rescaled range as range divided by std, and take mean across windows
-    rs_win = np.ptp(sig_rect, axis=0) / np.std(X_rect, 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.
+    """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.
+        Window length for each detrended fluctuation fit, in samples.
     deg : int, optional, default=1
         Polynomial degree for detrending.
 
     Returns
     -------
     det_fluc : float
-        Measured detrended fluctuaiton, as the average error fits of the window.
+        Measured detrended fluctuation, as the average error fits of the window.
     """
 
-    # Gather windows & vectorize, so we can call math functions in one go
-    n_win = int(np.floor(len(sig) / win_len))
-    sig_rect = np.reshape(sig[:n_win * win_len], (n_win, win_len)).T
+    # 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), sig_rect, deg=deg, full=True)
+    _, 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)
diff --git a/neurodsp/burst/__init__.py b/neurodsp/burst/__init__.py
index ec014d5a..6b27cad4 100644
--- a/neurodsp/burst/__init__.py
+++ b/neurodsp/burst/__init__.py
@@ -1,3 +1,4 @@
+"""Burst detection functions."""
 
 from .dualthresh import detect_bursts_dual_threshold
 from .utils import compute_burst_stats
diff --git a/neurodsp/utils/data.py b/neurodsp/utils/data.py
index 3a6af3aa..6bb4b040 100644
--- a/neurodsp/utils/data.py
+++ b/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

From 01b2120dc4e3ef32794f8bf424ed72e4c74b40c8 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Wed, 9 Oct 2019 15:45:42 -0700
Subject: [PATCH 05/11] lints

---
 neurodsp/aperiodic/__init__.py | 2 +-
 neurodsp/aperiodic/dfa.py      | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/neurodsp/aperiodic/__init__.py b/neurodsp/aperiodic/__init__.py
index 62af65c0..f1c5c6aa 100644
--- a/neurodsp/aperiodic/__init__.py
+++ b/neurodsp/aperiodic/__init__.py
@@ -1,3 +1,3 @@
 """"Functions for calculating aperiodic properties of signals."""
 
-from .dfa import compute_fluctuations
\ No newline at end of file
+from .dfa import compute_fluctuations
diff --git a/neurodsp/aperiodic/dfa.py b/neurodsp/aperiodic/dfa.py
index 4dc04e56..cd65e311 100644
--- a/neurodsp/aperiodic/dfa.py
+++ b/neurodsp/aperiodic/dfa.py
@@ -69,9 +69,9 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     fluctuations = np.zeros_like(t_scales)
     for idx, win_len in enumerate(win_lens):
 
-        if method is 'dfa':
+        if method == 'dfa':
             fluctuations[idx] = compute_detrended_fluctuation(sig, win_len=win_len, deg=deg)
-        elif method is 'rs':
+        elif method == 'rs':
             fluctuations[idx] = compute_rescaled_range(sig, win_len=win_len)
         else:
             raise ValueError('Fluctuation method not understood.')

From 7bf7002641318eca582265932a34daccc67468ae Mon Sep 17 00:00:00 2001
From: ryanhammonds <ryan.hammonds@utexas.edu>
Date: Tue, 28 Jul 2020 17:42:36 -0700
Subject: [PATCH 06/11] nan when win_len <= 1

---
 neurodsp/aperiodic/dfa.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/neurodsp/aperiodic/dfa.py b/neurodsp/aperiodic/dfa.py
index cd65e311..964342a7 100644
--- a/neurodsp/aperiodic/dfa.py
+++ b/neurodsp/aperiodic/dfa.py
@@ -69,7 +69,9 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     fluctuations = np.zeros_like(t_scales)
     for idx, win_len in enumerate(win_lens):
 
-        if method == 'dfa':
+        if win_len <= 1:
+            fluctuations[idx] = np.nan
+        elif 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)
@@ -77,7 +79,8 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
             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]
+    exp = np.polyfit(np.log10(t_scales[~np.isnan(fluctuations)]),
+                     np.log10(fluctuations[~np.isnan(fluctuations)]), deg=1)[0]
 
     return t_scales, fluctuations, exp
 

From 6cf0971f0907485513126c00d61baff673100e14 Mon Sep 17 00:00:00 2001
From: ryanhammonds <ryan.hammonds@utexas.edu>
Date: Tue, 28 Jul 2020 18:05:38 -0700
Subject: [PATCH 07/11] tests updated

---
 neurodsp/tests/aperiodic/test_dfa.py | 35 ++++++++++++++++++++++++++++
 neurodsp/tests/test_aperiodic_dfa.py | 20 ----------------
 2 files changed, 35 insertions(+), 20 deletions(-)
 create mode 100644 neurodsp/tests/aperiodic/test_dfa.py
 delete mode 100644 neurodsp/tests/test_aperiodic_dfa.py

diff --git a/neurodsp/tests/aperiodic/test_dfa.py b/neurodsp/tests/aperiodic/test_dfa.py
new file mode 100644
index 00000000..0caa2236
--- /dev/null
+++ b/neurodsp/tests/aperiodic/test_dfa.py
@@ -0,0 +1,35 @@
+"""Tests for fractal analysis using fluctuation measures."""
+
+import pytest
+
+from neurodsp.tests.settings import FS
+
+from neurodsp.aperiodic.dfa import (compute_fluctuations, compute_rescaled_range,
+                                    compute_detrended_fluctuation)
+
+###################################################################################################
+###################################################################################################
+
+@pytest.mark.parametrize("method",
+    ['dfa', 'rs', pytest.param(None, marks=pytest.mark.xfail(raises=ValueError))]
+)
+def test_compute_fluctuations(tsig, method):
+
+    t_scales, flucs, exp = compute_fluctuations(tsig, FS, method=method)
+
+    assert len(t_scales) == len(flucs)
+    assert exp > 0
+
+
+def test_compute_rescaled_range(tsig):
+
+    rs = compute_rescaled_range(tsig, 10)
+
+    assert isinstance(rs, float)
+
+
+def test_compute_detrended_fluctuation(tsig):
+
+    out = compute_detrended_fluctuation(tsig, 10)
+
+    assert out >= 0
diff --git a/neurodsp/tests/test_aperiodic_dfa.py b/neurodsp/tests/test_aperiodic_dfa.py
deleted file mode 100644
index c87fe81c..00000000
--- a/neurodsp/tests/test_aperiodic_dfa.py
+++ /dev/null
@@ -1,20 +0,0 @@
-"""Tests for fractal analysis using fluctuation measures."""
-
-from neurodsp.tests.settings import FS
-
-from neurodsp.aperiodic.dfa import *
-
-###################################################################################################
-###################################################################################################
-
-def test_compute_fluctuations(tsig):
-
-    t_scales, flucs, exp = compute_fluctuations(tsig, FS)
-
-def test_compute_rescaled_range(tsig):
-
-    rs = compute_rescaled_range(tsig, 10)
-
-def test_compute_detrended_fluctuation(tsig):
-
-    out = compute_detrended_fluctuation(tsig, 10)

From 10102bb8b5a76b029abe7b70ca65c993bf50d02d Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Wed, 29 Jul 2020 21:45:05 -0700
Subject: [PATCH 08/11] update DFA from review

---
 neurodsp/aperiodic/dfa.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/neurodsp/aperiodic/dfa.py b/neurodsp/aperiodic/dfa.py
index 964342a7..545c70b2 100644
--- a/neurodsp/aperiodic/dfa.py
+++ b/neurodsp/aperiodic/dfa.py
@@ -15,7 +15,7 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     ----------
     sig : 1d array
         Time series.
-    fs : float, Hz
+    fs : float
         Sampling rate, in Hz.
     n_scales : int, optional, default=10
         Number of scales to estimate fluctuations over.
@@ -41,7 +41,7 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     fluctuations : 1d array
         Average fluctuation at each scale.
     exp : float
-        Slope of line in loglog when plotting time scales against fluctuations.
+        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
@@ -53,7 +53,7 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     The relationship of this fluctuation measure across window sizes provides
     information on the fractal properties of a signal.
 
-    Measures available are:
+    Available measures are:
 
     - DFA: detrended fluctuation analysis
         - computes ordinary least squares fits across signal windows
@@ -65,13 +65,17 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
     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 win_len <= 1:
-            fluctuations[idx] = np.nan
-        elif method == 'dfa':
+        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)
@@ -79,8 +83,7 @@ def compute_fluctuations(sig, fs, n_scales=10, min_scale=0.01, max_scale=1.0, de
             raise ValueError('Fluctuation method not understood.')
 
     # Calculate the relationship between between fluctuations & time scales
-    exp = np.polyfit(np.log10(t_scales[~np.isnan(fluctuations)]),
-                     np.log10(fluctuations[~np.isnan(fluctuations)]), deg=1)[0]
+    exp = np.polyfit(np.log10(t_scales), np.log10(fluctuations), deg=1)[0]
 
     return t_scales, fluctuations, exp
 
@@ -141,6 +144,6 @@ def compute_detrended_fluctuation(sig, win_len, deg=1):
     _, 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)
+    det_fluc = np.mean((fluc / win_len))**0.5
 
     return det_fluc

From 3fcf009653646fb39b0e0de0a2a40fcaae7585b1 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Wed, 29 Jul 2020 21:45:31 -0700
Subject: [PATCH 09/11] add coloured noise test signals

---
 neurodsp/tests/conftest.py | 15 ++++++++++++---
 neurodsp/tests/settings.py |  1 +
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py
index ebc2dfec..36821d7e 100644
--- a/neurodsp/tests/conftest.py
+++ b/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."""
diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py
index 07d2a146..652f89c8 100644
--- a/neurodsp/tests/settings.py
+++ b/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
 

From 7d35d04cf0db08d6af28e4f753050b4f081c4522 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Wed, 29 Jul 2020 21:45:46 -0700
Subject: [PATCH 10/11] update tests for fluctuations

---
 neurodsp/tests/aperiodic/test_dfa.py | 48 ++++++++++++++++++++--------
 1 file changed, 35 insertions(+), 13 deletions(-)

diff --git a/neurodsp/tests/aperiodic/test_dfa.py b/neurodsp/tests/aperiodic/test_dfa.py
index 0caa2236..cc8d8132 100644
--- a/neurodsp/tests/aperiodic/test_dfa.py
+++ b/neurodsp/tests/aperiodic/test_dfa.py
@@ -1,8 +1,11 @@
 """Tests for fractal analysis using fluctuation measures."""
 
-import pytest
+from pytest import raises
 
-from neurodsp.tests.settings import FS
+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)
@@ -10,26 +13,45 @@
 ###################################################################################################
 ###################################################################################################
 
-@pytest.mark.parametrize("method",
-    ['dfa', 'rs', pytest.param(None, marks=pytest.mark.xfail(raises=ValueError))]
-)
-def test_compute_fluctuations(tsig, method):
-
-    t_scales, flucs, exp = compute_fluctuations(tsig, FS, method=method)
+def test_compute_fluctuations(tsig):
 
+    t_scales, flucs, exp = compute_fluctuations(tsig, 500)
     assert len(t_scales) == len(flucs)
-    assert exp > 0
 
+    # Check error for if the settings create window lengths that are too short
+    with raises(ValueError):
+        t_scales, flucs, exp = compute_fluctuations(tsig, 100)
 
-def test_compute_rescaled_range(tsig):
+    # 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)
 
-    rs = compute_rescaled_range(tsig, 10)
+    # 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)
 
-    assert isinstance(rs, float)
+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

From b7d9eaba328d24759e451b89bfc12e589f9c4641 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Wed, 29 Jul 2020 21:45:59 -0700
Subject: [PATCH 11/11] add fluctuations to API listing

---
 doc/api.rst | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/doc/api.rst b/doc/api.rst
index 7624c8d7..d91ed24a 100644
--- a/doc/api.rst
+++ b/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
 -----------