Fast algorithm for scaling analysis with higher-order detrending moving average method

Among scaling analysis methods based on the root-mean-square deviation from the estimated trend, it has been demonstrated that centered detrending moving average (DMA) analysis with a simple moving average has good performance when characterizing long-range correlation or fractal scaling behavior. Furthermore, higher-order DMA has also been proposed; it is shown to have better detrending capabilities, removing higher-order polynomial trends than original DMA. However, a straightforward implementation of higher-order DMA requires a very high computational cost, which would prevent practical use of this method. To solve this issue, in this study, we introduce a fast algorithm for higher-order DMA, which consists of two techniques: (1) parallel translation of moving averaging windows by a fixed interval; (2) recurrence formulas for the calculation of summations. Our algorithm can significantly reduce computational cost. Monte Carlo experiments show that the computational time of our algorithm is approximately proportional to the data length, although that of the conventional algorithm is proportional to the square of the data length. The efficiency of our algorithm is also shown by a systematic study of the performance of higher-order DMA, such as the range of detectable scaling exponents and detrending capability for removing polynomial trends. In addition, through the analysis of heart-rate variability time series, we discuss possible applications of higher-order DMA.

Figure 1

Moving average polynomials $\left\{{\tilde{y}}_n\left(i\right)\right\}$ of degree $m$ at scale $n$ (dashed lines) in centered DMA, where $\left\{y\right(i\left)\right\}$ is a discrete sample path of Brownian motion.

Figure 2

Illustration of parallel translation of a moving averaging window. In centered DMA, each moving average window over $[i,i+n-1]$ is shifted by a fixed interval $[-n^{\prime },n^{\prime }]$, where $n$ is the window size and $n=2n^{\prime }+1$. In backward DMA, each moving average window over $[i,i+n-1]$ is shifted by a fixed interval $[-n+1,0]$.

Figure 3

Illustration of recurrence calculations of $Y^{\left(k\right)}(i,n)$. Every $r$ point in the smallest scale $n_1, Y^{\left(k\right)}(lr+1,n_1)$ ($l=0,1,\cdots$) is calculated using ${\sum }_j{j}^{k}y\left(j\right)$ (gray shaded terms). Otherwise, subsequent $Y^{\left(k\right)}(i,n)$ values are calculated using recurrence formulas.

Figure 4

Relative error induced by iterative calculations of recurrence formulas [Eq. (12)]. The relative error is defined in the text. The mean values of the relative error when $n=11$ were calculated using 100 samples of white Gaussian noise.

Figure 5

Comparison of computation times. The CPU time of the process calculating ${\sigma }_{\mathrm{DMA}}\left(n\right)$ were measured on a Windows PC with an Intel Core i7-5930K (3.5 GHz). The computation times were averaged over 100 runs.

Figure 6

Estimated scaling exponent $\alpha$ for $m\mathrm{th}$ order DMA versus the scaling exponent $\beta$ of the analyzed time series. Numerically generated time series displaying the power-law scaling of the power spectral density, $S\left(f\right)\sim f^{-\beta }$, were analyzed. The length of the time series is $N={10}^5$. The mean value of $\alpha$, estimated using linear regression of ${log}_{10}{\sigma }_{\mathrm{DMA}}\left(n\right)$ versus ${log}_{10}n$ in the range ${10}^2\le n\le {10}^4$, was calculated based on 10 samples. Dashed lines represent the relation $\alpha =(\beta +1)/2$.

Figure 7

Detrending capability for removing a polynomial trend. The analyzed time series were generated by superposition of white Gaussian noise and a $q\mathrm{th}$ order polynomial function, and analyzed by centered DMA (center) and backward DMA (right). The estimated values of ${\sigma }_{\mathrm{DMA}}\left(n\right)$ for each value of $n$ were averaged over 10 samples. From top to bottom, linear, quadratic, cubic, quartic, and quintic trends were analyzed. The left panels show sample time series, where polynomial functions are described using dashed lines. The dashed lines in the center and right panels indicate lines with slope 0.5.

Figure 8

$1/f$ scaling behavior of heart-rate variability (HRV). (a) HRV time series of a normal subject (a 35-year-old male) [37]. (b) Log-log plots of ${\sigma }_{\mathrm{DMA}}$ vs. $n$ for centered DMA combined with DFA (solid lines) [25]. In the case of DFA, the fluctuation functions $F\left(n\right)$ were plotted instead of ${\sigma }_{\mathrm{DMA}}$; In solid lines from top to bottom, the orders of DFA are 1, 2, 3, 4, and 5. The curves have been shifted vertically for better visibility.