The first step is to estimate the power spectral density of the observation vector. This estimation can be carried out by calculating the periodogram of the signal corresponding to the squared modulus of the discrete Fourier transform of the signal. The periodogram is formally written as P ω j = 2 ∗ 1 2 π N ∑ j = 1 m x j e − i ω j j 2 , 2 π N < ω j < 2 π m N , (2) where the set of variables used in Equation 2 is the same as that described in Equation 1. To meet the requirements of Equation 1, we note that the frequency range ω is limited. The frequency ω = 0 is excluded, and the maximum frequency is ω = 2 π m / N , where m is the integer part of N − 1 / 2 . When the length of the signal is equal to the power of 2, this procedure excludes the Nyquist frequency; otherwise, it limits the length of the spectrum to half the length of the signal minus 1. Finally, note that the value of the periodogram is doubled. This is a method described at the end of the help paragraph of the periodogram() function on the MathWorks website ² , which conserves the total power in one-sided periodograms. The following MATLAB codes are used to estimate the periodogram of the signal:

…………………………………………………………

The first line standardizes the observation vector x by setting its mean to 0 and its standard deviation to 1. This operation is relatively common in the field and is essential for the rest of this tutorial because Equations 3, 4 in the following sections are derived from the autocorrelation function (and not from the autocovariance function) and therefore assume a variance equal to 1. This operation also normalizes the total power of the spectrum P , but it does not change the shape of the spectrum, so it does not alter the information of the fractal exponents.

The second line returns the length N of the input signal x, and the third computes the upper bound m. In the fourth line, we use the Signal Processing Toolbox command periodogram() to estimate the power spectral density of the input signal x, and an alternative code using the fft() command is given if the periodogram() command is not included in your MATLAB version. Finally, the fifth and the sixth lines bound P and w within the interval presented in Equation 2. Figure 1 represents the plot of the estimated periodograms of the series choleskyfgn, arfima0d0, whitenoise, and empirical.

Type Fig1_PSD in the command window to access Figure 1.

…………………………………………………………

In this sub-step, we present the equations and their computation for the two theoretical spectral densities derived from the fGn/fBm and ARFIMA (0,d,0) models, respectively.

The theoretical fGn spectral density was given by Li and Lim (2006) as follows: T f G n ′ ω j ; H = sin H π Γ 2 H + 1 ω j 1 − 2 H , (3) where Γ is the gamma function, and the other variables used in Equation 3 are the same as those described in Equation 1. The following MATLAB code is the computation of Equation 3.

…………………………………………………………

The theoretical ARFIMA (0,d,0) spectral density was given by Taqqu et al. (1995) as follows: T A R F I M A ′ ω j ; d = 1 2 π 2 ⁡ sin ω j 2 − 2 d , (4) where d is the integration parameter belonging to − 0.5,0.5 . One can easily convert the exponent H to the exponent d via the equation d = H − 0.5 . The following MATLAB code converts H to d and computes Equation 4.

…………………………………………………………

Figure 2 repeats Figure 1 by adding the theoretical spectral densities calculated using MATLAB codes 2 and 3. To better illustrate the global functioning of the algorithm, in Figure 2, we present the theoretical power spectral densities computed with the estimated values of H and d obtained via whittle.m. These values are reported in Table 1.

Type Fig2_TPSD in the command window to access Figure 2.

As advised by Jennane et al. (2001), in this sub-step, we calculate the constant of proportionality c to adjust the power of the model process to that of the signal as follows: c = ∑ ω P ω j ∑ ω T ′ ω j , H . (5)

…………………………………………………………

In the first line, the constant c is calculated, and in the second line, the theoretical spectrum Tp is adjusted to the empirical periodogram P. Figure 3 shows the plot of Figure 2 with adjusted theoretical spectrum. The total power of the theoretical spectrum Tp is determined by the value of H in Equation 3 or d in Equation 4, but the function of fractal exponents is to shape the spectrum, not to modulate its power, hence the need for adjustment between theoretical power and signal power.

Type Fig3_TPSD_adjusted in the command window to access Figure 3.

In this second step, we construct Whittle’s log-likelihood function by injecting the two previous sub-steps into Equation 1, which is the function that we want to maximize. However, the optimization toolbox of MATLAB only allows minimizing functions, so we will have to minimize the inverse of Equation 1. This inverse function is written in the MATLAB code as follows:

…………………………………………………………

where lwH is the inverse of Whittle’s likelihood function, N is the length of the observation vector, T is the scaled theoretical periodogram computed via MATLAB codes 2 or 3 and then MATLAB code 4, and P is the estimated periodogram of the observation vector computed via MATLAB code 1. MATLAB codes 6 and 7 are given in the following sections, which are the functions declared in MATLAB and which we will have to minimize. These codes correspond to the combination of MATLAB codes 2–5. MATLAB code 6 gives the function based on the theoretical spectrum of fGn, while MATLAB code 7 gives the function based on the theoretical spectrum of ARFIMA (0,d,0). These functions will have to be placed either at the end of the mother code whittle.m (as is the case with the provided code) or in two separate files that can be named WLLFfgn.m and WLLFarf.m.

…………………………………………………………

In the first line, the function is used to declare the functions WLLFfgn and WLLFarf. The outputs are lwHfgn for code 6 and lwHarf for code 7 and correspond to l W H in Equation 1. The inputs are as follows: H is the Hurst exponent, w is the Fourier frequencies, P is the estimated periodogram of the observation vector, and N is the length of the observation vector.

In Figure 4, we present the plots corresponding to these two functions for our four test signals: choleskyfgn, arfima0d0, whitenoise, and empirical, with H values ranging from 0.05 to 0.95 by steps of 0.05.

Type Fig4_lwH in the command window to access Figure 4.

In this third step, we describe the method to solve the optimization problem. In other words, we must find the value of H for which we reach the minimum of the inverse of Whittle’s likelihood: α ^ = argmin 0 < H < 1 2 N ∑ j = 1 m ln ⁡ c T ′ ω j ; H + P ω j c T ′ ω j ; H , (6) where α ^ is the estimate of the fractal exponent of the observation vector x j . As stated in the introduction of this first part, the present section 2.5 marks the transition from the exponents H and d characterizing the stationary processes fGn and ARFIMA (0,d,0) to the exponent α that can characterize the full set of stationary and non-stationary processes on a continuum from α = 0 to α = 2 . The α metric allows the whittle.m algorithm to work without making the a priori distinction between stationary and non-stationary signals, as DFA does.

In order to implement the operation of Equation 6, we use the MATLAB function fminbnd(), which is a minimizer based on golden section search and parabolic interpolation. The MATLAB codes to find the minimum of the WLfgn and Wlarf functions are as follows:

…………………………………………………………

Afgn and Aarf are the two values of α ^ ; the first one is estimated by fGn-based Whittle’s likelihood, and the second is estimated by ARFIMA-based Whittle’s likelihood. The syntax @(H) intimates the function fminbnd() that H is the variable to be optimized, while WLfgn (H,w,P,N) and WLarf(H,w,P,N) are the functions that are optimized (corresponding to MATLAB codes 6 and 7). Finally, 0 is the lower bound and 1 is the upper bound in the optimization problem. This algorithm never optimizes on the bounds, so using 0 as a lower bound and 1 as an upper bound satisfies both theoretical conditions 0 < H < 1 and − 0.5 < d < 0.5 .

By definition, Whittle’s log-likelihood function only holds for stationary signals such as fGn or ARFIMA (0,d,0). As a result, when analyzing non-stationary signals, the minimum of this function occurs when α ^ is almost equal to 1. More precisely, when we refer to the help provided with the fminbnd() function ³ , this translates into values of α ^ greater than 1 − 2 . 10 − 4 − 6 2,2204 × 10 − 16 ) ∼ 0.9998 . Thus, if the algorithm returns a value greater than 0.9998, we can classify the signal as non-stationary. To calculate the fractal exponent for the non-stationary signals, we apply the method proposed by Diebolt and Guiraud (2005), which consists in analyzing a differentiated version of the signal, and then add 1 to α ^ . This translates in MATLAB code as

…………………………………………………………

The first line of both codes is the logical test. If the answer is true, then the four consecutive lines are computed. In the second line, the periodogram is estimated on the differentiated observation vector diff(x). In the third line, mdiff is calculated because diff(x) is one point shorter than x. In the fourth and fifth lines, the zero frequency and those greater than mdiff are excluded. Finally, in the sixth line, Afgn and Aarf are estimated by adding 1 to the value returned by fminbnd().

We would like to warn the reader that in order to make this tutorial progressive, the steps are not ordered in the same way as in the whittle.m algorithm. Moreover, some algorithmic intricacies were unfolded in sub-steps, giving more lines of code than necessary to build the all-in-one whittle.m algorithm. To construct whittle.m, the reader should proceed as follows:

1. Begin with the header: function [Afgn, Aarf] = whittle(x), where x is the observation vector, Aarf is the exponent estimated using ARFIMA (0,d,0) theoretical power spectral density, and Afgn is the exponent estimated using fGn theoretical power spectral density.

To test Whittle’s maximum likelihood estimator, we generated two sets of signals, one for each theoretical model (fGn/fBm and ARFIMA (0,d,0)). The first one, based on the fGn properties, is the Cholesky decomposition, whose algorithm is named cholfgn.m. The second one consists of ARFIMA (0,d,0) filtering, whose algorithm is named ARFIMA0d0.m. We thus generated, for each of these generators, 42 subsets of 120 signals of length N = 1,024 for a total of 5,040 signals for each generator. These 42 subsets correspond to 42 different α -values; the first value is 0.01 because the value 0 is excluded from the theoretical models, the following 19 values range from 0.05 to 0.95 by steps of 0.05, and the 21st value is 0.99 because the value 1 is also excluded from the models. The last 21 values are ranged in the same way from 1.01 to 1.99.

We then estimated α ^ using three analysis methods: fGn-based Whittle’s likelihood and ARFIMA-based Whittle’s likelihood, presented in this tutorial and implemented in the whittle.m function, and DFA, which is implemented in the DFA.m function. Note that the DFA algorithm we used is the improved version presented by Almurad and Delignières (2016). The algorithm of generation and analysis is presented in signal_generation_and_analysis.m, and the results are saved in generatedsignals.mat. The first striking result is the analysis time. It took 57.24 s to compute α ^ on the set of signals using the two Whittle’s functions, while it took 2052.30 s to compute α ^ on the set of signals using DFA. In sum, the analysis for Whittle’s function is approximately 70 times faster than DFA. We conducted the analysis on a Windows laptop equipped with an Intel i77700HQ processor, 16 GB of RAM, an Nvidia GTX 1050 graphics card, and a 512 GB SSD hard drive using MATLAB version 2019a.

Figure 5 shows the whole set of computed α ^ . This figure is arranged in two columns corresponding to the cholfgn and ARFIMA0d0 generators and in three rows corresponding to the three analysis methods: fGn-based Whittle’s maximum likelihood estimator, ARFIMA-based Whittle’s maximum likelihood estimator, and DFA. The figure shows that the cholfgn generator is strongly biased around the frontier between noise and motion (for α = 1 ). We had already encountered this phenomenon in the study by Roume et al. (2019) using the Davies and Harte algorithm and can easily conclude that for generators based on the autocorrelation function of fGn, the integration of fractional Gaussian noises with an α close to 0 to obtain fractional Brownian motions with an α close to 1 is a technique that does not work properly. However, we can observe that this technique holds relatively well when the signals have been generated via ARFIMA filtering, as shown by the continuity between noise and motion observed in the bottom row. In addition, we will rely on the ARFIMA0d0 generator to estimate the efficiency of the analysis methods in the following section.

Type Fig5_Genbiases in the command window to access Figure 5.

In this section, we compare the efficiency of fGn-based Whittle’s likelihood, ARFIMA-based Whittle’s likelihood, and DFA. To assess the efficiency of these analysis methods, one must account for both their bias and variance. We therefore calculated Mean Squared Error (MSE), which is defined as the average of the squared error values. These analysis methods, being estimators in the sense of statistics, MSE can also be written as the sum of the variance and the squared bias of the values of α ^ , allowing us to directly compare their quality. Thus, when comparing two estimators, the estimator characterized by the smallest MSE value can be considered better.

First, for each estimator, we calculated the 5,040 (42 α -values × 120 signals) squared error values, which were then averaged. We obtained an MSE of 0.0014±0.0019 for fGn-based Whittle’s likelihood ⁴ , 0.0007±0.0011 for ARFIMA-based Whittle’s likelihood, and 0.0078±0.0127 for DFA. We performed a one-way ANOVA that confirmed the significant differences between these groups (F (2,15,117) = 1,385.8; p < 0.001; η ² = .15). Thus, ARFIMA-based Whittle’s likelihood is better than fGn-based Whittle’s likelihood, which itself is a better estimator than DFA. To go beyond the reductionist nature of this elementary comparison, a box plot of all the squared error values is constructed, as shown in Figure 6.

Type Fig6_SquarredError in the command window to access Figure 6.

In line with the results of the comparison of the MSE values, the first observation that can be made is that regardless of the value of α , the ARFIMA method is characterized by a lower median error than the other two analysis methods, as well as by a lower dispersion of these errors. We can also observe biases that could already be predicted from Figure 5; for example, for the fGn-based Whittle’s likelihood, the overestimation bias for α is less than 0.3, and the high variance around α is equal to 1. For DFA, the overestimation bias for α is less than 0.3, the underestimation bias for α is greater than 1.8, and there is a simultaneous increase between the variance of α ^ and the true α value.

In this section, we discuss the performance of ARFIMA-based Whittle’s likelihood depending on signal lengths. Fractal analysis, like DFA, often requires a large series (N > 500) to provide reliable alpha estimates (Damouras et al., 2010; Phinyomark et al., 2020). This can lead to experimental issues when working with physiological signals, such as the difficulty of older adults to walk for more than 5 min, which is between 100 and 250 strides (Moon et al., 2016; Phinyomark et al., 2020).

To assess the performance, we first made four sets of signals by reducing the length of tsarf from 1 to N, with the following four N values: 32, 64, 128, and 256. We then estimated α ^ for these four sets and computed MSE. We obtained an MSE value of 0.7709±0.7954 for N = 32, 0.0884±0.1540 for N = 64, 0.0124±0.0207 for N = 128, and 0.0036±0.0053 for N = 256. Similar to Figure 6, the box plot of α ^ squared error values for the four reduced sets is constructed, as shown in Figure 7.

Type Fig7_SignalLength in the command window to access Figure 7.

The first observation that can be made from this figure is that the analysis method provides unreliable estimates for lengths N = 32 and 64. The second is that for N = 256, the ARFIMA-based Whittle’s likelihood gives more reliable estimates than DFA with N = 1,024. Finally, for N = 128, this method gives slightly less reliable estimates than DFA for 1,024 points.

In addition, it can be observed that, particularly for N = 128 and N = 256, the errors are maximized to approximately α = 1. These errors are linked to misclassification of signals as stationary or non-stationary. This is even more predominant for signals with an α exponent just above 1 for N = 128, although it should be noted that these signals are truncated and the variation in mean and standard deviation was distributed over the entire original signal. It would be interesting to add a size variable to the generated signals in future studies.

For further comparison, we calculated the value of MSE for α ^ obtained with DFA for N = 512 points, which gives 0.0203 ± 0.355, i.e., approximately twice the value of MSE for ARFIMA-based Whittle’s likelihood with N = 128 (0.0124 ± 0.0207). Considering that the variability of the version of DFA we use is optimized, for N = 512, the variability is smaller than that of the original DFA with N =1,024 (Almurad and Delignières, 2016), and given that several studies using the original DFA with N lengths of approximately 1,000 data points have given elegant results, it is reasonable to predict that an α ^ value estimated via ARFIMA-based Whittle’s likelihood on a series of 128 points could be acceptable even if not optimal.

During the development of this tutorial, we tested Whittle’s maximum likelihood estimator on several physiological signals, including those made available by J. Hausdorff on the PhysioNet platform (Goldberger et al., 2000; Hausdorff, 2001). The results obtained are summarized in Table 2.

Under the conditions without metronome (slow, normal, and fast), ARFIMA-based Whittle’s likelihood and DFA provide similar results for the mean α ^ , but a slightly higher standard deviation was observed for DFA. This result becomes very interesting under metronome conditions, especially when the ARFIMA-based Whittle’s likelihood and DFA provide diametrically opposed results, the first detecting persistence ( α ^ > 0.5) and the second detecting anti-persistence ( α ^ < 0.5).

To understand why we obtained these results, we grouped the different series into three groups: the first group, “persistent behavior,” assembles the series for which DFA and Whittle’s likelihood analysis provide an α ^ greater than 0.5. In the second group, “anti-persistent behavior,” we group the time series for which both DFA and Whittle likelihood provide an α ^ below 0.5. Finally, in the third group, “mixed behavior,” we group the series with α ^ lower than 0.5 from DFA and α ^ higher than 0.5 from Whittle’s likelihood. Then, we estimated the power spectra using the periodogram method, and finally, we averaged these spectra for each group. These averaged spectra are presented in the first three columns in Figure 8, and the data are saved in spectralbehavior. mat.

Type Fig8_SpectralBehavior in the command window to access Figure 8.

In the first two columns, when the Whittle’s and DFA methods provided similar α ^ (i.e., both greater than 0.5 for the first column and both lower than 0.5 for the second column), we observe that the power spectrum is dominated by a trend that is persistent in the first column and anti-persistent in the second column. In the third column, when Whittle’s and DFA provided opposite results, we can observe a mixed behavior where the low frequency part of the spectrum is dominated by anti-persistence, while the high frequencies are dominated by persistence.

This non-monotonicity of the spectrum in the logarithmic space, which was already observed by Torre and Delignières (2008), is not predicted in the two models, fGn/fBm and ARFIMA (0,d,0); therefore, we cannot apply these models, from a theoretical point of view, on a signal being characterized by such a spectrum. However, the ARFIMA (p,d,q) model allows the creation of non-monotonic spectra, as shown in the fourth column in Figure 8. It will therefore be interesting to improve the method that we have presented so far by incorporating in Equation 4 the AR and MA components known as short-memory components, all of which will give a theoretical spectrum with three parameters: T A R F I M A ′ ω j ; φ , d , θ = 1 2 π 2 ⁡ sin ω j 2 − 2 d θ e − i ω j 2 φ e − i ω j 2 , (7) where

• φ e − i ω j = 1 − ∑ j = 1 p φ j e − i ω j j is the autoregressive (AR) component.

Finally, it will be interesting to compare the proposed method in this tutorial with Higuchi’s fractal dimension method that Phinyomark et al. (2020) has shown to be better than DFA, which is on our agenda for future work.