As described previously (Eke et al., 2000), statistically self-affine time series of the fGn type with exact long-range-dependent structuring were generated by the Davies and Harte method (DHM; Davies and Harte, 1987) at length, N, and pre-set Hurst exponent, ^fGnH_true. fBm-type time series were produced by cumulatively summing (Eke et al., 2000) their fGn counterparts, yielding a true extended Hurst exponent of H_true = ^fBmH_true = ^fGnH_true + 1. fBm-type time series were also produced directly at N and pre-set H_true = ^fBmH_true by the spectral synthesis method (SSM; Saupe, 1988).

Cantor sets and—by their cumulative summation—Cantor functions as examples of exact multifractal structuring were generated at pre-set weight factors. Statistical self-similar multifractals with known H(q) were synthesized for testing purposes using the generalized binomial multifractal model (Oświȩcimka et al., 2006; Schumann and Kantelhardt, 2011). The generating algorithm is an iterative process in which a stochastic binomial cascade is created at a given weight factor and at a given length. In the generalized method, the degree of correlation and the strength of multifractality can be separately tuned, the former being set with the Fourier Filtering Method; for further details see Schumann and Kantelhardt (2011) and Mukli et al. (2015).

After considering various numerical testing frameworks, we chose the DHM algorithm and Cantor functions as offering the best control over the cardinal parameters in our testing. The above-listed signals represent cases of unimodality with a single SR. Multimodal synthetic (or mock) signals with multiple SRs were created by adding these unimodal fractal time series of known attributes (N, H_true). Positioning of crossovers was controlled by setting the standard deviation (SD) ratios (i.e., focus ratio) of the signal components in addition to the differences in H_true.

Human NIRS measurements using a NIRO 500 Cerebral Oxygen Monitor (Hamamatsu Photonics, Hersching, Germany) at a rate of 2 Hz were carried out to record the relative change in total hemoglobin concentration with a length of N = 16,384 data points (for details, see Eke et al., 2006). Human EEG signals were sampled with a length of N = 16,384 data points with eyes closed during random hand movements at 500 Hz using a Neurofax EEG System (see Acknowledgment; Pattnaik and Sarraf, in press). These empirical records were acquired from healthy volunteers above the pre-frontal area. Rat fMRI-BOLD data with a length of N = 4,096 data points were obtained by using a modified 11.7 T Bruker horizontal-bore spectrometer (Bruker AVANCE, Billerica, MA, USA) using a 1 H surface coil (1.4 cm diameter) with sequential sampling gradient echo planar imaging (EPI) sequence (Hyder et al., 1995): field of view of 2.56 × 2.56 cm²; image matrix of 64 × 64; slice thickness of 2 mm; repetition time of 200 ms (i.e., 5 Hz of sampling frequency) and echo time of 13 ms; and voxel size of 400 × 400 × 2000 μm³ (for further details see Herman et al., 2011).

The FMF (Mukli et al., 2015) exploits the fact that the family of q-wise scaling functions yields a convergent structure merging at the coarsest scale, s = N. At this scale, the value defined in Equation (1) becomes q-independent yielding a “focus,” S[X_i](N). Enforcing S[X_i](N) as an analytical constraint improves the extraction of q-wise regression slopes for finite-length (by excluding inherently conflicting estimates due to the merely explicit emergence of this restriction) and thus offers a robust method for multifractal analysis of empirical data (Mukli et al., 2015; Ali et al., 2016; Delignières et al., 2016), the very approach adopted in this work. Instead of independently performing repetitive monofractal analyses on a set of empirical scaling functions (Equation 2), FMF methods—through minimizing the residual sum of squared errors, SSE—iteratively find the best-fitting, true multifractal and assign it to the family of the evaluated scaling functions (Equation 3).

In Equation (2), x has two specific values (at x = 0 and in the case of the “focus” x = N); otherwise, it represents the scale, where the exact scale-dependent statistic is being evaluated. The case of x = N and represents the enforced constraint in Equation (3). Thus, according to FMF, a set of model (i.e., exact) scaling functions with iterated parameters—Ĥ[X_i](q) and log Ŝ[X_i](N)—are fitted all at once to the actual data set of the scaling functions. In order to obtain an overall measure of the goodness-of-fit of the FMF regression procedure, its mean squared error (MSE) was calculated according to Equation (21) of Mukli et al. (2015).

The standard segmented line regression method is capable of finding breakpoints, s_b, and also in the case of a superimposed signal (Figure 3) approximating crossovers, s_x. Equation (4) is an adaptation of the segmented line regression method for a bimodal scaling function, where s′ could be any particular temporal scale. To capture q-dependent breakpoints, we introduced a q-wise regression algorithm, broken down into three steps of Equations (4a–c)

where indices ^f and ⁿ stand for different fractal processes: in our particular case uncorrelated (noise) and correlated (fractal) signals within a co-sampled arrangement, respectively. We chose noise and fractal signals as the constituents of a bimodal signal in describing our methods because this was the case for bimodal cerebral hemodynamic data reported earlier (Eke et al., 2006; Herman et al., 2011) and is used in this study as an exemplary dataset. The breakpoint for a given moment is obtained at the minimum value of the SSE(q,s′) function as the estimates, Ĥ[X_i](q) and log Ŝ[X_i](q,x) for the fractal and noise components, respectively, are being refined during the iteration process. Further away from the breakpoint, in the low range of scales, the underlying noise with low H will dominate, while in the high range of scales, the fractal component with high H will dominate. However, in between the two, where the respective function values are commensurable, their fluctuations blend into a non-fractal segment, a so-called exclusion range (ER; Figure 3). Consequently, by excluding the ER from the analysis, the error in assessing the low and high Hs will be smaller, while for the FMF analysis yielding H(q), the error will increase due to SR being contracted by the exclusion range. To find the ER, we used Equation (5) (a similar approach as Equation (20) of Ge and Leung, 2013) at a specific q and selected a range of scales with SSE(q,s′) lower than SSE_tolerance calculated at a preset level of tolerance (0 ≤ tolerance ≤ 1)

Any moment-to-moment inconsistencies in the regression analysis will upset the expected structural aspect of multifractal scaling functions known as the “H(q)-monotonicity” (Mukli et al., 2015), the monotonous drop of regression slopes of H(q) from q = −∞ to q = +∞, because in subsequent steps of multifractal formalism this automatically results in “inversed” or “corrupted” multifractal spectra. To eliminate this eventuality, the analytical constraints of H[X_i](q) > H[X_i](q-k) with k > 0 were enforced on the nested iterative processes of minimization. MSE was obtained from the sum of SSEs from respective FMF analyses see Equation (3) on the noise and fractal components separated by q-wise breakpoints.

As fluctuations from the two underlying signal components mutually contribute to each other's scaling functions near the breakpoint, they hold estimates deviating from the power-law relationship (Figure 4). Taking the exemplary case of the SSC algorithm—where the statistical measure is the standard deviation, SD—this relationship is readily seen as a realization of the Bienaymé formula stating that in the case of uncorrelated variables, the variance (SD²) of their sum equals the sum of the respective variances (Bienaymé, 1853). A generalization for the case of (anti)correlated signals is given in the Appendix (See Supplementary Material). Thus, the emerging scaling function (the power mean of SD) can be well estimated (Figure A1 in Supplementary Material) as the root sum square of the composing scaling functions as

where the signal (^cX_i) used in the calculation of S is in square brackets with c being a positive integer referring to each and every of the N_c constituent signals.

Earlier—for the cases of resting-state cerebral hemodynamic fluctuations—we showed that a fractally correlated signal is typically interwoven by uncorrelated noise (Eke et al., 2006; Herman et al., 2011). In this simplest case of bimodality, the above derivation can be reduced to the following direct approximation

Accordingly, instead of fitting the two constituting fractally correlated components of a bimodal scaling function separately in two distinct processes, an exact bimodal model scaling function is reconstructed from two properly fitting power-law sets, based on the rule of addition Equation (7). Performing this one-pass regression on a log-log scale, the minimization of SSE—with the generalized Hurst exponent and the focus of the scaling function being iterated—results in the best fit of the exact bimodal model as follows

A special application of this procedure is when one of the constituting components of the composite signal is known, which obviously reduces the number of tuning factors in the minimization process. Specifically, this component could be uncorrelated noise [i.e., instrument and/or biological noise (Peng et al., 1995; Blesic et al., 2003; Eke et al., 2006; Herman et al., 2011)] with a fractal estimate of Ĥ ≅ H_true = 0.5 [as will be seen later in Section Performance of the SFD Method on an Empirical Bimodal Signal with Limited Definition (fMRI-BOLD)].

SFD is not at all limited to q-wise applications, but can also be performed along with FMF. In this case, the process of minimization of the FMF analysis needs to be modified by raising the number of tuning parameters Equation (9). Thus, both the two sets (ⁿX_i and ^fX_i) of H(q) and their associated two foci, S(N), see Equations (9a and 9b) are being simultaneously adjusted in the same iterative process Equation (9c)

Similarly to the qSRA method, H(q)-monotonicity was granted by applying the same analytical constraints. The calculation of MSE from SSE Equation (9) was as explained in Section Focus-Based Multifractal Method.

The crossover scale, s_x, of the decomposed scaling functions—where the respective statistical values are in principle the same—can be determined as the common value of the equations of the two underlying regression lines

Thus, the crossover scale can be calculated as

When enforcing the respective foci of the underlying components, the best value of the crossover scale is obtained as

Additive random or correlated noise readily upsets multifractal analysis as demonstrated by Ludescher et al. (2011). In particular, when random noise with increasing amplitudes is added to a multifractal signal, the q-wise scaling function obtained by standard tools such as, MF-DFA will exhibit crossovers gradually shifting to larger scales. Accordingly, H(q) evaluated within the range of chosen scales becomes dominated by that of the added noise component. While the aim of these authors was not to provide a solution for handling noise contamination, they made an important contribution demonstrating that added noise propagates across scales ultimately leading to spurious results in multifractal analysis.

As for the impact of the component focus ratio (Figure 6), from the above-mentioned geometrical properties of multifractal scaling functions and the relationships shown in Figure 5A2 it follows that the crossover scale is low when both H and the component focus ratio are low (Figure 6, vignettes in lower left corner). Conversely, it is high when both H and the component focus ratio are high (Figure 6, vignettes in upper right corner). In between these extremes, a diagonal band of low bias due to the impact of mid-range crossover scales in the data is seen (Figure 6, bottom row) where the presence of merging scale-free patterns can be statistically confirmed (Clauset et al., 2009). When—due to the actual representation of scaling in the empirical data—the analysis is not performed within this optimized range, the signal definition should be improved: crossover scales that are much too low require an increased sampling rate, those that are much too high call for a longer signal to be collected (Eke et al., 2002).

There are cases when the superposition of two fractal components yields a composite signal with a crossover falling outside the observed range of scales (Figure 10). The multifractal spectrum in this case is typically asymmetric (Drożdż and Oświęcimka, 2015). While under these conditions the crossover is not directly accessible to our SFD-based analysis, our additive model still allows for its characterization and offers an explanation for the asymmetry in D(h). This way, a composite process yielding asymmetric D(h) can be modeled too. Asymmetric D(h) can also be interpreted as a phase transition (Grassberger et al., 1988; Muzy et al., 1993; Arneodo et al., 1995; Radons and Stoop, 1996) based on the deep analogy that exists between the multifractal formalism and equilibrium statistical thermodynamics (Stanley and Meakin, 1988; Tel, 1988; Arneodo et al., 1995). According to this model, the superimposed partition functions under and above a critical q and the signals mutually perturb each other to a slight extent (Arneodo et al., 1995); thus, the partition function, τ(q) and D(h) are always dominated by a single component under which condition the Bienaymé formula yields similar results in describing the superposition (Ssuperimposed≈Scomponent2+0=Scomponent).

In a scaling function representation of empirical temporal multifractality, the crossover scale for the chosen smallest negative moment is the largest and it becomes the smallest at the largest positive moment (see for example Figure 5B2). This moment-wise distribution of crossover scales emerges from the geometrical underpinnings of FMF (see Figure 2C) and the way H(q)-dependence is formulated in Equation (12) yielding the crossover scale itself. As crossover scales and breakpoints are similar manifestations of scaling, breakpoints should also be captured in a moment-wise manner.

The significance of the breakpoint in the analysis of bi- or multimodal signals has already been recognized in the literature (Peng et al., 1995; Kantelhardt et al., 2001; Eke et al., 2006; Herman et al., 2011; Ge and Leung, 2013). However, disregarding their moment-wise nature (Ge and Leung, 2013) distorts the acquired H(q)s with breakpoints falling outside the proper SR. This generates uneven error propagation across the moments in the case of standard multifractal methods (Ludescher et al., 2011), or introduces a slight but global error when FMF methods are applied. To the best of our knowledge, the qSRA- and SFD-FMF methods should be regarded as first attempts to carry out genuinely adaptive multifractal analyses on bi- or multimodal signals in a properly designed, moment-wise manner. They can readily be combined with various fractal and multifractal tools as their adaptive step, thus opening new possibilities for future applications.

The crossover between the EEG signal components was found at the boundary between the δ and θ bands (Figure 8) of EEG classification. An independent δ and θ rhythm has already been proposed due to the significant interregional gap in synchrony (Mormann et al., 2008; Fetterhoff et al., 2015). Our current findings also prove the presence of an independent multifractal δ rhythm. The NIRS signal was also found to be the sum of temporally correlated, scale-free fluctuations dominated by low frequencies (i.e., low-frequency fluctuations) and uncorrelated [but not instrumental (Eke et al., 2006)] noise.

Resting-state brain dynamics as captured in fMRI-BOLD fluctuations is powered by ongoing neurodynamics spreading across the functional connections of a fractally organized anatomical network of an immense neuronal pool (Bullmore et al., 2009; Werner, 2010). Previously, we demonstrated that the resting-state fMRI-BOLD signal in an animal (Herman et al., 2011) and the fNIRS signal in a human model (Eke et al., 2006), alike, were bimodal temporal fractals. Thus, subject to future but likely advances in fMRI-BOLD technology, an improved signal definition may likely reveal bimodality in human fMRI-BOLD signals too. However, due to limitations inherent to the current fMRI-BOLD technology (Eke et al., 2012), their scale-free character—especially in cases of bimodality—is hard to capture (see Figure 12). Accordingly, while the SFD analysis yielded valid estimates of crossovers in a sizable fraction of the scanned voxels, still—depending on the correlation level and component focus ratio—numerous crossover estimates—due to sub-optimal data acquisition—proved invalid at extreme scales. This severely upset the topology of the multifractal parametric maps and prompted us to look for justifiable grounds to reduce the degree of freedom in finding breakpoints and/or crossovers in an attempt to improve the performance of the analysis. We did indeed find that the uncorrelated components of the scaling functions could be taken as default throughout the voxels of the scan. This uncorrelated component was determined at a site (Figure 9) where the precision of the method has been shown to be maximal in tests on synthetic signals (Figure 6). This procedure effectively compensated for the limited definition in the BOLD signals as seen in the enhanced topology of the parametric maps (Figure 9D).