The red traces in Figure 1 show an ordinary random walk (lower panel), a monofractal random walk (middle panel) and a multifractal random walk (upper panel). The fractal property of these random walks is reflected by their picture-in-picture similarity as illustrated in the upper panel of Figure 1. Small “hills” and “valleys” with similar structure appear when you zoom on the large “hills” and “valleys” of the random walk. The DFA is employed to time series with a random walk like structure (Peng et al., 1995). However, most biomedical time series have fluctuations that are more similar to the increments of the random walks (see the blue traces in Figure 1). If the biomedical time series has the noise like structure of the blue traces in Figure 1, it should be converted to a random walk like time series before employing DFA. Noises can be converted to random walks by subtracting the mean value and integrate the time series. Time series whitenoise, monofractal, and multifractal are all noise like time series and are converted to random walk like time series by Matlab code 1 below:

A conventional analysis of variation in biomedical time series is to compute the average variation as a RMS. The reader can use Matlab code 2 below to compute RMS for the time series whitenoise, monofractal, and multifractal:

Figure 2 illustrates that the average amplitude of variation (i.e., RMS) is equal for all three time series, whitenoise, monofractal, and multifractal, even though they have quite different structures. MFDFA will be able to distinguish between these structures as we will see in the sections below.

The multifractal time series in the upper panel have local fluctuations with both large and small magnitudes. RMS in Matlab code 2 can be computed for segments of the time series to differentiate between the magnitudes of the local fluctuations. A simple procedure is to cut the time series into equal-sized non-overlapping segments and compute a local RMS for each segment. This can be done by Matlab code 3 below and is the core procedure of MFDFA:

The first line of Matlab code 3 converts the noise like time series, multifractal, to a random walk like time series X (i.e., Matlab code 1). The third line of Matlab code 3 set the parameter scale that defines the sample size of the non-overlapping segments in which the local RMS, RMS{1}, are computed. The fifth line is the number of segments that the time series X can be divided into where length(X) is the sample size of time series X. Thus, segments = 8 when length(X) = 8000 and scale = 1000. The sixth to fourteenth line is a loop that computes the local RMS around a trend fit{v} for each segment by updating the time Index (see red v arguments in Matlab code 3). In the first loop v = 1, the Index{1} goes from sample 1 to sample 1000. In the second loop v = 2, the Index{2} goes from sample 1001 to sample 2000. In the final loop v = 8, the Index{8} goes from sample 7001 to 8000.

Slow varying trends are present in biomedical time series and detrending of the signal is therefore necessary to quantify the scale invariant structure of the variation around these trends. In Matlab code 3, a polynomial trend fit{v} is fitted to X within each segment v (see blue command lines in Matlab code 3). The first blue command line is the parameter m that defines the order of the polynomial. The polynomial trend are linear when m = 1, quadratic when m = 2, and cubic when m = 3 (see Figures 3A–C). The first blue command line within the loop defines the polynomial coefficients C used to create the polynomial trend fit{v} of each segment (see dashed red lines in Figure 3). The local fluctuation, RMS{1}(v), is then computed for the residual variation, X(Index{v})-fit{v}, within each segment v. The local fluctuation, RMS{1}(v), is illustrated in Figure 3 as the distance between the red dashed trends and the red solid lines.

In the DFA the variation of the local RMS{1} are quantified by an overall RMS (F) in Matlab code 4 below:

The fast changing fluctuations in the time series X will influence the overall RMS, F, for segments with small sample sizes (i.e., small scale) whereas slow changing fluctuations will influence F for segments with large sample sizes (i.e., large scale). The scaling function, F, should therefore be computed for multiple segments sizes (i.e., multiple scales) to emphasize both fast and slow evolving fluctuations that influence the structure of the time series. The scaling function, F(ns), can be computed for multiple scales by including Matlab code 3 and 4 within a new loop marked as red command lines and arguments ns below:

In the first red command line, a vector with multiple segment sizes (i.e., scales) is set by the reader. In the second red command line, a loop is initiated where Matlab code 3 is computed from the smallest to the largest scale. The segment sample size, scale(ns), are updated by the red index ns. The local fluctuation, RMS{ns}, is a set of vectors where each vector have a length equal to the number of segments [i.e., segments(ns)]. In the first loop for ns = 1, the local fluctuation RMS{1} is a vector of local RMS values computed for 500 segments [i.e., segments(1)] each containing 16 samples [i.e., scale(1)]. In the last loop for ns = 7, the local fluctuation RMS{7} is a vector with local RMS values computed for seven segments [i.e., segments(7)] each containing 1024 samples [i.e., scale(7)]. In the last command line, the scaling function (i.e., overall RMS), F(ns), are computed for multiple scales by Matlab code 4. Figure 4 illustrates the local fluctuations, RMS{ns}, and the overall RMS, F(ns), for multiple scales.

DFA indentifies the monofractal structure of a time series as the power law relation between the overall RMS (i.e., F in Matlab code 4) computed for multiple scales (i.e., scale in Matlab code 5). The power law relation between the overall RMS is indicated by the slope (H) of the regression line (RegLine) computed by Matlab code 6 below:

The slope, H, of the regression line, RegLine, is called the Hurst exponent (Hurst, 1951). The Hurst exponent defines the monofractal structure of the time series by how fast the overall RMS, F, of local fluctuations, RMS, grows with increasing segment sample size (i.e., scale). Figure 5 shows that the overall RMS, F, of local fluctuations, RMS, is growing faster with the segment sample size for the monofractal and multifractal time series compared with whitenoise time series. The larger Hurst exponent, H, is visually seen as more slow evolving variations (i.e., more persistent structure) in monfractal and multifractal time series compared with whitenoise. Figure 6 illustrates that the Hurst exponents defines a continuum between a noise like time series and a random walk like time series. The Hurst exponent will be in the interval between 0 and 1 for noise like time series whereas above 1 for a random walk like time series. A time series has a long-range dependent (i.e., correlated) structure when the Hurst exponent is in the interval 0.5–1 and an anti-correlated structure when the Hurst exponent is in the interval 0–0.5. The time series has an independent or short-range dependent structure in the special case when the Hurst exponent is equal to 0.5. According to Figure 5, time series whitenoise has a time independent structure with Hurst exponent close to 0.5 whereas monofractal, and multifractal has a long-range dependent structure with Hurst exponent between 0.7 and 0.8. The reader should notice that short-range dependent processes can mimic the scale invariance in Figure 5 for certain scaling ranges (cf. Gao et al., 2006).

The structure of the monofractal and multifractal time series are different even though they have similar overall RMS and slopes H in Figure 5. The multifractal time series have local fluctuations with both extreme small and large magnitudes that is absent in the monofractal time series. The absence of fluctuations with extreme large and small magnitudes results in a normal distribution for the monofractal time series where the variation is described by the second order statistical moment (i.e., variance) alone. The monofractal DFA are therefore based on the second order statistics of the overall RMS (i.e., F in Matlab code 4). In the multifractal time series, local fluctuation, RMS{ns}(v), will be extreme large magnitude for segments v within the time periods of large fluctuations and extreme small magnitude for segments v within the time periods of small fluctuations. Consequently, the multifractal time series are not normal distributed and all q-order statistical moments should to be considered. Thus, it’s necessary to extend the overall RMS in the monofractal DFA (i.e., F in Matlab code 4) to the following q-order RMS of the multifractal DFA (Fq in Matlab code 7 below):

The first command line in Matlab code 7 defines a set of q-orders from −5 to 5. The second line initiates a loop that computes the overall q-order RMS, Fq(nq), from negative to positive q’s (see the red arguments nq). The q-order weights the influence of segments with large and small fluctuations, RMS{1}, as illustrated in Figure 7. Fq(nq) for negative q’s (i.e., nq = 1–3) are influenced by segments v with small RMS{1}(v). In contrast, Fq(nq) for positive q’s (i.e., nq = 4–6) are influenced by segments v with large RMS{1}(v). The local fluctuations RMS{1} with large and small magnitudes is graded by the magnitude of the negative or positive q-order, respectively. Fq for q = −3 and 3 is more influenced by the segments v with the smallest and largest RMS{1}(v), respectively, compared to Fq for q = -1 and 1 (see Figure 7). The midpoint q = 0 are neutral to influence of segments with small and large RMS{1}. Notice that the last line of Matlab code 7 redefines the special case q(nq) = 0 because 1/0 goes to infinity [i.e., 1/q(q = = 0) = inf in Matlab]. The reader should also notice that Fq(q = = 2) are equal to second order statistics F in Matlab code 4 because sqrt(x) = x^(1/2) in Matlab. The monofractal DFA in Matlab code 5 can now be extended to a MFDFA by simply changing Matlab code 4 to Matlab code 7 in the last line of Matlab code 5. This change is highlighted with red command lines in Matlab code 8 below:

The q-order Hurst exponent can now be defined as the slopes (Hq) of regression lines (qRegLine) for each q-order RMS (Fq). Both Hq(nq) and qRegLine{nq} is computed by looping Matlab code 6 for each q-order (see red command lines and arguments nq in Matlab code 9 below):

Figure 8 shows that the slopes Hq of the regression lines are q-dependent for the multifractal time series (see Figure 8A). The difference between the q-order RMS for positive and negative q’s are more visual apparent at the small segment sizes compared to the large segment sizes (see Figure 8A). The small segments are able to distinguish between the local periods with large and small fluctuations (i.e., positive and negative q’s, respectively) because the small segments are embedded within these periods. In contrast, the large segments cross several local periods with both small and large fluctuations and will therefore average out their differences in magnitude. Thus, the relationship between the q-order RMS of the multifractal time series becomes similar to the monofractal time series at the largest segment sizes (compare Figures 8A,B). Both the monofractal and whitenoise time series have no periods with small and large fluctuations and, consequently, the same difference between the q-order RMS irrespective of the segment sample sizes (see Figures 8B,C). The growing similarity between the q-order RMS of multifractal and monofractal time series with increasing segment sample size leads to a decreasing Hq for multifractal time series (see Figure 8D). The decreasing Hq indicates that the segments with small fluctuations have a random walk like structure whereas segments with large fluctuations have a noise like structure (see the continuum of Hurst exponents in Figure 6). The similarity between the scaling function F of monofractal and multifractal time series in Figure 5 is indicated by the intercept of Hq around q = 2 (compare blue and red traces in Figure 8D). Thus, the monofractal DFA in Matlab code 5 and 6 will not distinguish between the structure of the monofractal and multifractal time series.

The q-order Hurst exponent Hq is only one of several types of scaling exponents used to parameterize the multifractal structure of time series. The typical procedure in the literature of MFDFA is to first convert Hq to the q-order mass exponent (tq) and thereafter convert tq to the q-order singularity exponent (hq) and q-order singularity dimension (Dq; Kantelhardt et al., 2002). The plot of hq versus Dq is referred to as the multifractal spectrum (see Figure 9C). The q-order mass exponent tq can be computed from Hq by the Matlab code 10 below (see Figure 9B):

The mass exponent tq is used to compute the q-order singularity exponent (hq) and the q-order singularity dimension (Dq) by Matlab code 11 below (see upper right Figure 9):

The monofractal and whitenoise time series has a mass exponent tq with a linear q-dependency. The linear q-dependency of tq leads to a constant hq of these time series because hq is the tangent slope of tq (see the first command line in Matlab code 11). The constant hq reduces the multifractal spectrum to a small arc for the monofractal and whitenoise time series (see Figure 9C). In contrast, the multifractal time series has mass exponents tq with a curved q-dependency and, consequently, a decreasing singularity exponent hq. The resulting multifractal spectrum is a large arc where the difference between the maximum and minimum hq are called the multifractal spectrum width (see arrow in Figure 9C).

The reader should notice that the q-order singularity exponent hq and corresponding dimension Dq computed by Matlab code 11 are referred to as α and f(α) in Kantelhardt et al. (2002), but as h and D(h) in other literature (e.g. Ihlen and Vereijken, 2010). Furthermore, the singularity dimension can be confused with the generalized dimension and the box counting dimension that is other ways to parameterize the multifractal structure of time series (see Equation 14 in Kantelhardt et al., 2002).

The Hurst exponent defined by the monofractal DFA represents the average fractal structure of the time series as illustrated in Figure 6 and is closely related to the central tendency of multifractal spectrum. The deviation from average fractal structure for segments with large and small fluctuations is represented by the multifractal spectrum width. Thus, each average fractal structure in the continuum of Hurst exponents (see Figure 6) has a new continuum of multifractal spectrum widths that represents the deviations from the average fractal structure (see Figure 10). Furthermore, the shape of the multifractal spectrum in Figure 10 does not have to be symmetric. The multifractal spectrum can also have either a left or a right truncation that originate from a leveling of the q-order Hurst exponent for negative or positive q’s, respectively (see Figure 11). The leveling of q-order Hurst exponent reflects that the q-order RMS is insensitive to the magnitude of the local fluctuations. The multifractal spectrum will have a long left tail when the time series have a multifractal structure that is insensitive to the local fluctuations with small magnitudes (see upper Figure 11). In contrast, the multifractal spectrum will have a long right tail when the time series have a multifractal structure that is insensitive to the local fluctuations with large magnitudes (see lower Figure 11). Another continuum of right and left truncations exists for each multifractal spectrum width in the continuum illustrated in Figure 10. Thus, the width and shape of the multifractal spectrum is able to classify a wide range of different scale invariant structures of biomedical time series.

The transformation of the q-order Hurst exponent Hq to the mass exponent tq and, finally, to the multifractal spectrum Dq and hq is stated as a “just-the-way-it-is” argument in the above section without mathematical details. The reader might ask at this point why one should define and interpret the multifractal spectrum Dq and hq and not only Hq that are directly estimated by Matlab code 8 and 9. Estimating the multifractal spectrum directly from the local fluctuation, will answer this question and give a less abstract definition of the multifractal spectrum. A local Hurst exponent can be defined directly from, RMS{ns}(v), for each time instant v. The local Hurst exponent estimated for a multifractal time series will fluctuate in time in contrast to the time independent Hurst exponent estimated by the monofractal DFA (see Matlab code 5 and 6; Figure 5). The temporal variation of the local Hurst exponent can be summarized in a probability distribution and the multifractal spectrum is just the normalized probability distribution in log-coordinates. Thus, the width and shape of the multifractal spectrum reflect the temporal variation of the local Hurst exponent or, in other words, the temporal variation in the local scale invariant structure of the time series. In order to estimate the local Hurst exponent, the local fluctuation, RMS{ns}(v), has to be defined within a translating segment centered at sample v instead of within non-overlapping segments. Thus, Matlab code 8 has to be modified to Matlab code 12 below (see red command lines):

The first red command line defines a vector of small segment sizes (i.e., scale_small) where the segment sizes increases with two samples. This increase of two samples is necessary to align the center of segments according to the Time_index. The third red line set the Time_index for the local fluctuation, RMS{ns}, that exclude the halfmax number of samples at the start and the end of the time series (see second red command line). Then a loop are initiated for each segment size where the local fluctuations, RMS{ns}(v), are computed for a translating segment centered at sample v = Time_index. The translating segment includes the local samples that are updated by T_index (see last red command line). The local Hurst exponent (Ht) can now be computed from the local fluctuation, RMS{ns}, by the Matlab code 13 below:

The first two command line defines the regression line Regfit equal to the regression line qRegLine{q }= = 0 computed by Matlab code 8. Regfit represents the center of the spread of local RMS and are the regression line of the overall q-order RMS, Fq(q = = 0,:). A loop for each scale ns computes the residual fluctuation resRMS{ns} of log2(RMS{ns}) around the regression line Regfit for each sample v of the time series. In Figure 12B, the differences between the overall q-order RMS, Fq, converge toward each other with increasing scale. This convergence is inevitable for multifractal variation by the linear relationship between Fq for all q-order and the assumption of monotonical decreasing q-order Hurst exponent, Hq (see Figure 8D). The same convergence is seen for the local RMS in Figure 12A and is used to estimate the local Hurst exponents, Ht(ns,:). Ht(ns,:) is estimated as the slope of the line from local RMS in log-coordinates to the endpoint of the regression line, Regfit, at the largest scale, maxL (see Figure 12). Consequently, Ht(ns,:) are obtained by dividing the residuals resRMS{ns} by logscale (i.e., the difference between maximal scale maxL and the scale(ns) in log-coordinates) and adding the slope Hq(q = = 0) of regression line, Regfit (see Figure 12A). Figure 13A illustrates the local Hurst exponent Ht(ns,:) for ns = 5 (i.e., scale(ns) = 15) for the multifractal, monofractal, and whitenoise time series. The local Hurst exponent Ht(ns,:) has larger variation for the multifractal time series compared to the monofractal and whitenoise time series. The small Ht(ns,:) in the periods of the multifractal time series with local fluctuations of large magnitudes (i.e., large RMS{ns}) reflects the noise like structure of the local fluctuations (see red dashed lines in Figure 13A). In contrast, the larger Ht(ns,:) in the periods with local fluctuations of small magnitudes (i.e., small RMS{ns}) reflects the random walk like structure of the local fluctuations (see black dashed lines in Figure 13A). The local Hurst exponent Ht in periods with fluctuations of small and large magnitudes is therefore consistent with the q-order Hurst exponent Hq for negative and positive q’s, respectively. The advantage of local Hurst exponent Ht compared with q-order Hurst exponent Hq is the ability of Ht to identify the time instant of structural changes within the time series. In studies where the physiological phenomenon is perturbed at some time instant v, the local Hurst exponent Ht(ns,v) can identify how this perturbation affects the local scale invariant structure of the biomedical time series. The temporal variation of local Hurst exponent Ht can be summarized in a histogram representing the probability distribution (Ph) of Ht (see Figure 13B). The multifractal spectrum (Dh) is defined simply by the log-transformation of the normalized probability distribution (Ph_norm). The probability distribution (Ph) and multifractal spectrum (Dh) are computed by Matlab code 14 below:

The first line in Matlab code 14 convert the matrix Ht to the vector Ht_row that are the input argument in hist function used to compute the histogram for Ht_row. The second input argument in hist function are BinNumb that set the number of bins in the histogram. A sufficient choice for BinNumb is the square root of the sample size of Ht_row (see the second command line). The output variables of hist function are the center of each bin Htbin and the number freq of local Hurst exponents that fall into each bin. The probability distribution Ph are computed by dividing the number freq of local Hurst exponents in each bin by the total number of local Hurst exponents, sum(freq) (see Figure 13B). The multifractal spectrum Dh are computed by first define Ph_norm by normalizing Ph to the maximum probability max(Ph) and then divide log(Ph_norm) by –log(mean(scale)) (cf. Struzik, 2000; Scafetta et al., 2003). The multifractal spectrum Dh are therefore directly related to the distribution Ph of the local fractal structure of the time series. The distribution Ph is the same for the local scale invariant structure of the time series as the conventional probability distribution are for the local amplitudes of the time series. The present state of the physiological system is connected to both past and future states that influence the local scale invariant structure of time series. Thus, distribution Ph and the multifractal spectrum Dh of biomedical time series might reflect important properties of the self-regulation of physiological processes.

MFDFA1 and MFDFA2 have the best performance when signal are a noise like time series. However, it can be difficult according to Figure 6 to visually differentiate between random walk and noise like time series. A possible solution suggested by Eke et al. (2002) is to run a monofractal DFA (i.e., Matlab code 5 and 6) before running MFDFA1 and MFDFA2. The time series are noise like if Hurst exponent H is between 0.2 and 0.8. In this case, MFDFA1 and MFDFA2 can be employed directly without transformation of the time series. However, the time series are random walk like when H is between 1.2 and 1.8. In these cases, the time series should either be differentiated before entering the MFDFA1 or MFDFA2 or the conversion to random walk in the first line of Matlab code 8 and 12 should be eliminated. If the time series are random walk like + 1 should be added to the output variables Hq, hq, and tq from MFDFA1 and Ht and Htbin from MFDFA2. Table 1 summarize the categories of the Hurst exponent estimated by a monofractal DFA with corresponding conversion of the biomedical time series that should be performed before entering it into MFDFA1 and MFDFA2.

The local fluctuation in the time series is defined as a local RMS within both MFDFA1 and MFDFA2. Large error appears in the multifractal spectrum when RMS is close to zero because both log2(Fq) for negative q’s in Matlab code 8 and log2(RMSt) in Matlab code 12 becomes infinitely small (i.e., -inf in Matlab). Extreme large Hq will be present for negative q’s as output from MFDFA1. Equivalently, extreme large outliers in Ht will be present as output from MFDFA2. Consequently, local RMS close to zero will lead to large right tails for the multifractal spectrum. The problem of segments with RMS close to zero can be solved by eliminating RMS below a certain threshold (eps). The threshold eps can be set to the precision of the measurement device that is recording the biomedical time series. As an example, the measurement of the inter-beat intervals of the human heart is measured as the time interval between R-peaks in ECG and has a typical precision of 1 millisecond. Thus, RMS below 1 millisecond can be eliminated from further analysis when MFDFA1 and MFDFA2 are employed to series of inter-beat intervals. Elimination of local fluctuations below the measurement error is possible in MFDFA1 by setting eps = 1 and RMS{ns}(RMS < eps) = [] in Matlab code 8.

There are two main reasons why the local fluctuation RMS becomes zero in segments with small sample sizes. First, the polynomial trend fit of the time series can be overfitted in segments with small sample sizes (i.e., small scale). An overfitted trend will be similar to the time series and cause the residual fluctuations, RMS, to be close to zero. The sample size of the smallest segment (i.e., scale) should therefore be much larger than the polynomial order m to prevent an overfitted trend. Secondly, the biomedical time series might be smooth with little apparent variation and therefore similar to the polynomial trend even for low order m. In these cases, the value of the smallest scales should be raised and the scale invariance checked carefully (see “Is the Time Series Scale Invariant?” below).

The application of both Matlab function MFDFA1 and MFDFA2 assumes that the biomedical time series are scale invariant. This means that plot(log2(scale),log2(Fq)) yield a linear relationship between log2(scale) and log2(Fq) (see Figure 8)., The q-order Hurst exponent Hq should not be estimated by a linear regression if the relationship between log2(scale) and log2(Fq) is curved or S-shaped. Consequently, the first output from MFDFA1 to be visually checked should be plot(log2(scale),log2(Fq)). Non-linear relation in this plot might arise from several reasons. First, an insufficient order m for the polynomial detrending will yield a non-linear relationship between log2(scale) and log2(Fq) for scale invariant time series with a trend. The solution is to run MFDFA1 or MFDFA2 multiple times with different m and compare their plot(log2(scale),log2(Fq). Secondly, local fluctuations RMS close to zero for small scales would yield a non-linear dip in lower end of plot(log2(scale),log2(Fq)). This dip can be prevented by elimination of RMS below the measurement error suggested in Section “Local Fluctuations Close to Zero?” or by choosing a larger minimum input scale. Finally and most importantly, a non-linear relationship in plot(log2(scale),log2(Fq)) might originate from the phenomenon recorded in the biomedical time series. As an example, respiratory frequency creates distinct oscillations in the fast fluctuations of the heart rate variability (Stein and Kleiger, 1999) and cause the scale invariance to break down at the smallest scales. Another example is postural sway in humans where the variation of the center of pressure has two distinct scaling regions thought to represent two distinct modes for human balance control (Collins and De Luca, 1993). One way to detect the sub-regions with scale invariance is to look for periods with approximately constant log2(Fq(q = = 1,:)./scale) in plot(log2(scale),log2(Fq(q = = 1,:)./scale)) within the entire scaling range (cf. Gao et al., 2006). The scales where log2(Fq(q = = 1,:)./scale) are no longer constant indicates the segment sizes above and below which the local fluctuations (i.e., RMS) are no longer scale invariant. These points will in many cases have phenomenological explanations and should not be ignored.

The Matlab functions MFDFA1 and MFDFA2 have input parameters scale, q and m. The estimation of the multifractal spectra is dependent on these parameter settings. The rest of this section gives guidelines to the parameter settings in MFDFA1 and MFDFA2:

The basic component of both MFDFA1 and MFDFA2 are the local fluctuation, RMS. Statistical parameters other than RMS can be used to define the local fluctuation in a time series. In multifractal analyses based on wavelet transformations, the local fluctuation is defined as the convolution product between the time series and a waveform fitted within local segments of the time series (cf. Daubechies, 1992; Mallat, 1999). The results from analyses called wavelet transformation modulus maxima (Muzy et al., 1991), multifractal analysis with wavelet leaders (Jaffard et al., 2006; Wendt, 2008), and gradient modulus wavelet projection (Turiel et al., 2006) can therefore be directly compared with the results from MFDFA1 and MFDFA2. In an entropy-based estimation of the multifractal spectrum, the local fluctuation is defined as the sum of the time series within the local segment relative to the total sum of the entire time series (Chhabra and Jensen, 1989). This method uses a q-order entropy function instead of a q-order RMS and estimates hq and Dq, directly, as the regression slope of the q-order entropy functions. The MFDFA has been shown to perform as well as or better than these multifractal analyses (Kantelhardt et al., 2002; Oświęcimka et al., 2006; Serrano and Figliola, 2009; Huang et al., 2011). However, extensions of detrending procedure in MFDFA1 and MFDFA2 should be considered when the biomedical time series contains strong oscillatory or ramp-like trends (Hu et al., 2009; Huang et al., 2011).