Source code for neurokit2.complexity.fractal_dfa

# -*- coding: utf-8 -*-
from warnings import warn

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from ..misc import find_knee

[docs] def fractal_dfa( signal, scale="default", overlap=True, integrate=True, order=1, multifractal=False, q="default", maxdfa=False, show=False, **kwargs, ): """**(Multifractal) Detrended Fluctuation Analysis (DFA or MFDFA)** Detrended fluctuation analysis (DFA) is used to find long-term statistical dependencies in time series. For monofractal DFA, the output *alpha* :math:`\\alpha` corresponds to the slope of the linear trend between the scale factors and the fluctuations. For multifractal DFA, the slope values under different *q* values are actually generalised Hurst exponents *h*. Monofractal DFA corresponds to MFDFA with *q = 2*, and its output is actually an estimation of the **Hurst exponent** (:math:`h_{(2)}`). The Hurst exponent is the measure of long range autocorrelation of a signal, and :math:`h_{(2)} > 0.5` suggests the presence of long range correlation, while :math:`h_{(2)} < 0.5`suggests short range correlations. If :math:`h_{(2)} = 0.5`, it indicates uncorrelated indiscriminative fluctuations, i.e. a Brownian motion. .. figure:: ../img/douglas2022a.png :alt: Illustration of DFA (Douglas et al., 2022). Multifractal DFA returns the generalised Hurst exponents *h* for different values of *q*. It is converted to the multifractal **scaling exponent** *Tau* :math:`\\tau_{(q)}`, which non-linear relationship with *q* can indicate multifractility. From there, we derive the singularity exponent *H* (or :math:`\\alpha`) (also known as Hölder's exponents) and the singularity dimension *D* (or :math:`f(\\alpha)`). The variation of *D* with *H* is known as multifractal singularity spectrum (MSP), and usually has shape of an inverted parabola. It measures the long range correlation property of a signal. From these elements, different features are extracted: * **Width**: The width of the singularity spectrum, which quantifies the degree of the multifractality. In the case of monofractal signals, the MSP width is zero, since *h*\\(q) is independent of *q*. * **Peak**: The value of the singularity exponent *H* corresponding to the peak of singularity dimension *D*. It is a measure of the self-affinity of the signal, and a high value is an indicator of high degree of correlation between the data points. In the other words, the process is recurrent and repetitive. * **Mean**: The mean of the maximum and minimum values of singularity exponent *H*, which quantifies the average fluctuations of the signal. * **Max**: The value of singularity spectrum *D* corresponding to the maximum value of singularity exponent *H*, which indicates the maximum fluctuation of the signal. * **Delta**: the vertical distance between the singularity spectrum *D* where the singularity exponents are at their minimum and maximum. Corresponds to the range of fluctuations of the signal. * **Asymmetry**: The Asymmetric Ratio (AR) corresponds to the centrality of the peak of the spectrum. AR = 0.5 indicates that the multifractal spectrum is symmetric (Orozco-Duque et al., 2015). * **Fluctuation**: The *h*-fluctuation index (hFI) is defined as the power of the second derivative of *h*\\(q). See Orozco-Duque et al. (2015). * **Increment**: The cumulative function of the squared increments (:math:`\\alpha CF`) of the generalized Hurst's exponents between consecutive moment orders is a more robust index of the distribution of the generalized Hurst's exponents (Faini et al., 2021). This function can be called either via ``fractal_dfa()`` or ``complexity_dfa()``, and its multifractal variant can be directly accessed via ``fractal_mfdfa()`` or ``complexity_mfdfa()``. .. note :: Help is needed to implement the modified formula to compute the slope when *q* = 0. See for instance Faini et al. (2021). See Parameters ---------- signal : Union[list, np.array, pd.Series] The signal (i.e., a time series) in the form of a vector of values. scale : list A list containing the lengths of the windows (number of data points in each subseries) that the signal is divided into. Also referred to as Tau :math:`\\tau`. If ``"default"``, will set it to a logarithmic scale (so that each window scale has the same weight) with a minimum of 4 and maximum of a tenth of the length (to have more than 10 windows to calculate the average fluctuation). overlap : bool Defaults to ``True``, where the windows will have a 50% overlap with each other, otherwise non-overlapping windows will be used. integrate : bool It is common practice to convert the signal to a random walk (i.e., detrend and integrate, which corresponds to the signal 'profile') in order to avoid having too small exponent values. Note that it leads to the flattening of the signal, which can lead to the loss of some details (see Ihlen, 2012 for an explanation). Note that for strongly anticorrelated signals, this transformation should be applied two times (i.e., provide ``np.cumsum(signal - np.mean(signal))`` instead of ``signal``). order : int The order of the polynomial trend for detrending. 1 corresponds to a linear detrending. multifractal : bool If ``True``, compute Multifractal Detrended Fluctuation Analysis (MFDFA), in which case the argument ``q`` is taken into account. q : Union[int, list, np.array] The sequence of fractal exponents when ``multifractal=True``. Must be a sequence between -10 and 10 (note that zero will be removed, since the code does not converge there). Setting ``q = 2`` (default for DFA) gives a result of a standard DFA. For instance, Ihlen (2012) uses ``q = [-5, -3, -1, 0, 1, 3, 5]`` (default when for multifractal). In general, positive q moments amplify the contribution of fractal components with larger amplitude and negative q moments amplify the contribution of fractal with smaller amplitude (Kantelhardt et al., 2002). maxdfa : bool If ``True``, it will locate the knee of the fluctuations (using :func:`.find_knee`) and use that as a maximum scale value. It computes max. DFA (a similar method exists in :func:`entropy_rate`). show : bool Visualise the trend between the window size and the fluctuations. **kwargs : optional Currently not used. Returns ---------- dfa : float or pd.DataFrame If ``multifractal`` is ``False``, one DFA value is returned for a single time series. parameters : dict A dictionary containing additional information regarding the parameters used to compute DFA. If ``multifractal`` is ``False``, the dictionary contains q value, a series of windows, fluctuations of each window and the slopes value of the log2(windows) versus log2(fluctuations) plot. If ``multifractal`` is ``True``, the dictionary additionally contains the parameters of the singularity spectrum. See Also -------- fractal_hurst, fractal_tmf, entropy_rate Examples ---------- **Example 1:** Monofractal DFA .. ipython:: python import neurokit2 as nk signal = nk.signal_simulate(duration=10, frequency=[5, 7, 10, 14], noise=0.05) @savefig p_fractal_dfa1.png scale=100% dfa, info = nk.fractal_dfa(signal, show=True) @suppress plt.close() .. ipython:: python dfa As we can see from the plot, the final value, corresponding to the slope of the red line, doesn't capture properly the relationship. We can adjust the *scale factors* to capture the fractality of short-term fluctuations. .. ipython:: python scale = nk.expspace(10, 100, 20, base=2) @savefig p_fractal_dfa2.png scale=100% dfa, info = nk.fractal_dfa(signal, scale=scale, show=True) @suppress plt.close() **Example 2:** Multifractal DFA (MFDFA) .. ipython:: python @savefig p_fractal_dfa3.png scale=100% mfdfa, info = nk.fractal_mfdfa(signal, q=[-5, -3, -1, 0, 1, 3, 5], show=True) @suppress plt.close() .. ipython:: python mfdfa References ----------- * Faini, A., Parati, G., & Castiglioni, P. (2021). Multiscale assessment of the degree of multifractality for physiological time series. Philosophical Transactions of the Royal Society A, 379(2212), 20200254. * Orozco-Duque, A., Novak, D., Kremen, V., & Bustamante, J. (2015). Multifractal analysis for grading complex fractionated electrograms in atrial fibrillation. Physiological Measurement, 36(11), 2269-2284. * Ihlen, E. A. F. E. (2012). Introduction to multifractal detrended fluctuation analysis in Matlab. Frontiers in physiology, 3, 141. * Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Havlin, S., Bunde, A., & Stanley, H. E. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. Physica A: Statistical Mechanics and its Applications, 316(1-4), 87-114. * Hardstone, R., Poil, S. S., Schiavone, G., Jansen, R., Nikulin, V. V., Mansvelder, H. D., & Linkenkaer-Hansen, K. (2012). Detrended fluctuation analysis: a scale-free view on neuronal oscillations. Frontiers in physiology, 3, 450. """ # Sanity checks if isinstance(signal, (np.ndarray, pd.DataFrame)) and signal.ndim > 1: raise ValueError( "Multidimensional inputs (e.g., matrices or multichannel data) are not supported yet." ) n = len(signal) scale = _fractal_dfa_findscales(n, scale) # Sanitize fractal power (cannot be close to 0) q = _sanitize_q(q, multifractal=multifractal) # Store parameters info = {"scale": scale, "q": q} # Preprocessing if integrate is True: # Get signal profile signal = np.cumsum(signal - np.mean(signal)) # Function to store fluctuations. For DFA this is an array. For MFDFA, this # is a matrix of size (len(scale),len(q)) fluctuations = np.zeros((len(scale), len(q))) # Start looping over scale for i, window in enumerate(scale): # Get window segments = _fractal_dfa_getwindow(signal, n, window, overlap=overlap) # Get polynomial trends trends = _fractal_dfa_trends(segments, window, order=order) # Get local fluctuation fluctuations[i] = _fractal_dfa_fluctuation(segments, trends, q) if len(fluctuations) == 0: return np.nan, info # Max. DFA --------------------- if maxdfa is True: # Find knees of fluctuations knee = np.repeat(len(scale), fluctuations.shape[1]) for i in range(fluctuations.shape[1]): knee[i] = find_knee( y=np.log2(fluctuations[:, i]), x=np.log2(scale), show=False, verbose=False ) knee = np.exp2(np.nanmax(knee)) # Cut fluctuations fluctuations = fluctuations[scale <= knee, :] scale = scale[scale <= knee] # ------------------------------ # Get slopes slopes = _slopes(scale, fluctuations, q) if len(slopes) == 1: slopes = slopes[0] # Prepare output info["Fluctuations"] = fluctuations info["Alpha"] = slopes # Extract features if multifractal is True: info.update(_singularity_spectrum(q=q, slopes=slopes)) out = pd.DataFrame( { k: v for k, v in info.items() if k in [ "Peak", "Width", "Mean", "Max", "Delta", "Asymmetry", "Fluctuation", "Increment", ] }, index=[0], ) else: out = slopes # Plot if show is True. if show is True: if multifractal is True: _fractal_mdfa_plot( info, scale=scale, fluctuations=fluctuations, q=q, ) else: _fractal_dfa_plot(info=info, scale=scale, fluctuations=fluctuations) return out, info
# ============================================================================= # Utils # ============================================================================= def _fractal_dfa_findscales(n, scale="default"): # Convert to array if isinstance(scale, list): scale = np.asarray(scale) # Default scale number if scale is None or isinstance(scale, str): scale = int(n / 10) # See if isinstance(scale, int): scale = np.exp(np.linspace(np.log(10), np.log(int(n / 10)), scale)).astype(int) scale = np.unique(scale) # keep only unique # Sanity checks (return warning for too short scale) if len(scale) < 2: raise ValueError("NeuroKit error: more than one window is needed. Increase 'scale'.") if np.min(scale) < 2: raise ValueError( "NeuroKit error: there must be at least 2 data points in each window. Decrease 'scale'." ) if np.max(scale) >= n: raise ValueError( "NeuroKit error: the window cannot contain more data points than the time series. Decrease 'scale'." ) return scale def _sanitize_q(q=2, multifractal=False): # Turn to list if isinstance(q, str): if multifractal is False: q = [2] else: q = [-5, -3, -1, 0, 1, 3, 5] elif isinstance(q, (int, float)): q = [q] # Fractal powers as floats q = np.asarray_chkfinite(q, dtype=float) return q def _slopes(scale, fluctuations, q): # Extract the slopes of each `q` power obtained with MFDFA to later produce # either the singularity spectrum or the multifractal exponents. # Note: added by Leonardo Rydin (see # Ensure mfdfa has the same q-power entries as q if fluctuations.shape[1] != q.shape[0]: raise ValueError("Fluctuation function and q powers don't match in dimension.") # Allocated array for slopes slopes = np.zeros(len(q)) # Find slopes of each q-power for i in range(len(q)): # if fluctiations is zero, log2 wil encounter zero division old_setting = np.seterr(divide="ignore", invalid="ignore") slopes[i] = np.polyfit(np.log2(scale), np.log2(fluctuations[:, i]), 1)[0] np.seterr(**old_setting) return slopes def _fractal_dfa_getwindow(signal, n, window, overlap=True): # This function reshapes the segments from a one-dimensional array to a # matrix for faster polynomail fittings. 'Overlap' reshapes into several # overlapping partitions of the data # TODO: see whether this step could be integrated within complexity_coarsegraining if overlap: segments = np.array([signal[i : i + window] for i in np.arange(0, n - window, window // 2)]) else: segments = signal[: n - (n % window)] segments = segments.reshape((signal.shape[0] // window, window)) return segments def _fractal_dfa_trends(segments, window, order=1): # TODO: can we rely on signal_detrend? x = np.arange(window) coefs = np.polyfit(x[:window], segments.T, order).T # TODO: Could this be optimized? Something like np.polyval(x[:window], coefs) trends = np.array([np.polyval(coefs[j], x) for j in np.arange(len(segments))]) return trends def _fractal_dfa_fluctuation(segments, trends, q=2): # Detrend detrended = segments - trends # Compute variance var = np.var(detrended, axis=1) # Remove where var is zero var = var[var > 1e-08] if len(var) == 0: warn("All detrended segments have no variance. Retuning NaN.") return np.nan # Find where q is close to zero. Limit set at |q| < 0.1 # See is0 = np.abs(q) < 0.1 # Reshape q to perform np.float_power q_non0 = q[~is0].reshape(-1, 1) # Get the fluctuation function, which is a function of the windows and of q # When q = 2 (i.e., multifractal = False) # The formula is equivalent to np.sqrt(np.mean(var)) # And corresponds to the Root Mean Square (RMS) fluctuation = np.float_power(np.mean(np.float_power(var, q_non0 / 2), axis=1), 1 / q_non0.T) if np.sum(is0) > 0: fluc0 = np.exp(0.5 * np.mean(np.log(var))) fluctuation = np.insert(fluctuation, np.where(is0)[0], [fluc0]) return fluctuation # ============================================================================= # Utils MFDFA # ============================================================================= def _singularity_spectrum(q, slopes): """Extract the slopes of the fluctuation function to futher obtain the singularity strength `α` (or Hölder exponents) and singularity spectrum `f(α)` (or fractal dimension). This is iconically shaped as an inverse parabola, but most often it is difficult to obtain the negative `q` terms, and one can focus on the left side of the parabola (`q>0`). Note that these measures rarely match the theoretical expectation, thus a variation of ± 0.25 is absolutely reasonable. The parameters are mostly identical to `fractal_mfdfa()`, as one needs to perform MFDFA to obtain the singularity spectrum. Calculating only the DFA is insufficient, as it only has `q=2`, and a set of `q` values are needed. Here defaulted to `q = list(range(-5,5))`, where the `0` element is removed by `_cleanse_q()`. This was first designed and implemented by Leonardo Rydin in `MFDFA <>`_ and ported here by the same. """ # Components of the singularity spectrum # --------------------------------------- # The generalised Hurst exponents `h(q)` from MFDFA, which are simply the slopes of each DFA # for various `q` values out = {"h": slopes} # The generalised Hurst exponent h(q) is related to the Scaling Exponent Tau t(q): out["Tau"] = q * slopes - 1 # Calculate Singularity Exponent H or α, which needs tau out["H"] = np.gradient(out["Tau"]) / np.gradient(q) # Calculate Singularity Dimension Dq or f(α), which needs tau and q # The relation between α and f(α) (H and D) is called the Multifractal (MF) spectrum or # singularity spectrum, which resembles the shape of an inverted parabola. out["D"] = q * out["H"] - out["Tau"] # Features (Orozco-Duque et al., 2015) # --------------------------------------- # The width of the MSP quantifies the degree of the multifractality # the spectrum width delta quantifies the degree of the multifractality. out["Width"] = np.nanmax(out["H"]) - np.nanmin(out["H"]) # the singularity exponent, for which the spectrum D takes its maximum value (α0) out["Peak"] = out["H"][np.nanargmax(out["D"])] # The mean of the maximum and minimum values of singularity exponent H out["Mean"] = (np.nanmax(out["H"]) + np.nanmin(out["H"])) / 2 # The value of singularity spectrum D, corresponding to the maximum value of # singularity exponent H, indicates the maximum fluctuation of the signal. out["Max"] = out["D"][np.nanargmax(out["H"])] # The vertical distance between the singularity spectrum *D* where the singularity # exponents are at their minimum and maximum. out["Delta"] = out["D"][np.nanargmax(out["H"])] - out["D"][np.nanargmin(out["H"])] # the asymmetric ratio (AR) defined as the ratio between h calculated with negative q and the # total width of the spectrum. If the multifractal spectrum is symmetric, AR should be equal to # 0.5 out["Asymmetry"] = (np.nanmin(out["H"]) - out["Peak"]) / out["Width"] # h-fluctuation index (hFI), which is defined as the power of the second derivative of h(q) # hFI tends to zero in high fractionation signals. if len(slopes) > 3: # Help needed to double check that! out["Fluctuation"] = np.sum(np.gradient(np.gradient(out["h"])) ** 2) / ( 2 * np.max(np.abs(q)) + 2 ) else: out["Fluctuation"] = np.nan # hFI tends to zero in high fractionation signals. hFI has no reference point when a set of # signals is evaluated, so hFI must be normalisedd, so that hFIn = 1 is the most organised and # the most regular signal in the set # BUT the formula in Orozco-Duque (2015) is weird, as HFI is a single value so cannot be # normalized by its range... # Faini (2021): new index that describes the distribution of the generalized Hurst's exponents # without requiring the Legendre transform. This index, αCF, is the cumulative function of the # squared increments of the generalized Hurst's exponents between consecutive moment orders. out["Increment"] = np.sum(np.gradient(slopes) ** 2 / np.gradient(q)) return out # ============================================================================= # Plots # ============================================================================= def _fractal_dfa_plot(info, scale, fluctuations): polyfit = np.polyfit(np.log2(scale), np.log2(fluctuations), 1) fluctfit = 2 ** np.polyval(polyfit, np.log2(scale)) plt.loglog(scale, fluctuations, "o", c="#90A4AE") plt.xlabel(r"$\log_{2}$(Scale)") plt.ylabel(r"$\log_{2}$(Fluctuations)") plt.loglog(scale, fluctfit, c="#E91E63", label=r"$\alpha$ = {:.3f}".format(info["Alpha"])) plt.legend(loc="lower right") plt.title("Detrended Fluctuation Analysis (DFA)") return None def _fractal_mdfa_plot(info, scale, fluctuations, q): # Prepare figure fig = plt.figure(constrained_layout=False) spec = matplotlib.gridspec.GridSpec(ncols=2, nrows=2) ax_fluctuation = fig.add_subplot(spec[0, 0]) ax_spectrum = fig.add_subplot(spec[0, 1]) ax_tau = fig.add_subplot(spec[1, 0]) ax_hq = fig.add_subplot(spec[1, 1]) n = len(q) colors =, 1, n)) for i in range(n): # Plot the points ax_fluctuation.loglog( scale, fluctuations[:, i], "o", fillstyle="full", markeredgewidth=0.0, c=colors[i], alpha=0.2, markersize=6, base=2, zorder=1, ) # Plot the polyfit line polyfit = np.polyfit(np.log2(scale), np.log2(fluctuations[:, i]), 1) fluctfit = 2 ** np.polyval(polyfit, np.log2(scale)) ax_fluctuation.loglog(scale, fluctfit, c=colors[i], base=2, label="_no_legend_", zorder=2) # Add labels for max and min if i == 0: ax_fluctuation.plot( [], label=f"$h$ = {polyfit[0]:.3f}, $q$ = {q[0]:.1f}", c=colors[0], ) elif i == (n - 1): ax_fluctuation.plot( [], label=f"$h$ = {polyfit[0]:.3f}, $q$ = {q[-1]:.1f}", c=colors[-1], ) ax_fluctuation.set_xlabel(r"$\log_{2}$(Scale)") ax_fluctuation.set_ylabel(r"$\log_{2}$(Fluctuations)") ax_fluctuation.legend(loc="lower right") # Plot the singularity spectrum # ax.set_title("Singularity Spectrum") ax_spectrum.set_ylabel(r"Singularity Dimension ($D$)") ax_spectrum.set_xlabel(r"Singularity Exponent ($H$)") ax_spectrum.axvline( x=info["Peak"], color="black", linestyle="dashed", label=r"Peak = {:.3f}".format(info["Peak"]), ) ax_spectrum.plot( [np.nanmin(info["H"]), np.nanmax(info["H"])], [np.nanmin(info["D"])] * 2, color="red", linestyle="solid", label=r"Width = {:.3f}".format(info["Width"]), ) ax_spectrum.plot( [np.nanmin(info["H"]), np.nanmax(info["H"])], [info["D"][-1], info["D"][0]], color="#B71C1C", linestyle="dotted", label=r"Delta = {:.3f}".format(info["Delta"]), ) ax_spectrum.plot(info["H"], info["D"], "o-", c="#FFC107") ax_spectrum.legend(loc="lower right") # Plot tau against q # ax.set_title("Scaling Exponents") ax_tau.set_ylabel(r"Scaling Exponent ($τ$)") ax_tau.set_xlabel(r"Fractal Exponent ($q$)") ax_tau.plot(q, info["Tau"], "o-", c="#E91E63") # Plot H against q # ax.set_title("Generalised Hurst Exponents") ax_hq.set_ylabel(r"Singularity Exponent ($H$)") ax_hq.set_xlabel(r"Fractal Exponent ($q$)") ax_hq.plot(q, info["H"], "o-", c="#2196F3") fig.suptitle("Multifractal Detrended Fluctuation Analysis (MFDFA)") return None