A Novel MFDFA Algorithm and Its Application to Analysis of Harmonic Multifractal Features

Abstract

A power grid harmonic signal is characterized as having both nonlinear and nonstationary features. A novel multifractal detrended fluctuation analysis (MFDFA) algorithm combined with the empirical mode decomposition (EMD) theory and template movement is proposed to overcome some shortcomings in the traditional MFDFA algorithm. The novel algorithm is used to study the multifractal feature of harmonic signals at different frequencies. Firstly, the signal is decomposed and the characteristics of wavelet transform multiresolution analysis are employed to obtain the components at different frequency bands. After this, the local fractal characteristic of the components is studied by utilizing the novel MFDFA algorithm. The experimental results show that the harmonic signals exhibit obvious multifractal characteristics and that the multifractal intensity is related to the signal frequency. Compared with the traditional MFDFA algorithm, the proposed method is more stable in curve fitting and can extract the multifractal features more accurately.

1. Introduction

With a large number of nonlinear electrical equipment and electronic devices being applied in the power system, more and more harmonics are appearing in the power grid. Investigations show that the pollutions from harmonics have become increasingly severe in the power grid [1,2,3,4,5,6]. Therefore, the study and analysis of harmonics has become an important research topic both at home and abroad. Traditional harmonic analysis methods mainly employ artificial neural network [7,8,9,10,11,12,13], Hilbert Huang Transform (HHT) [14,15], instantaneous reactive power [16], wavelet transform [17,18] and Fourier transform (FT) [19], etc. However, these methods ignore the self-similarity that exists in the time-domain waveform and the trend of harmonics, i.e., the so-called fractal characteristics [20]. Therefore, bringing fractal characteristics into harmonic analysis is of great significance.

As one of the emerging theories for studying nonlinear systems, monofractal analysis has made great progress in the nonlinear signal analysis. However, monofractal is mainly a description of the overall average of the object of study, with the local characteristics of the signal being insufficiently characterized [21]. Multifractal analysis [22,23] provides the ability to describe the local characteristics of the signal in detail. By combining multifractal and detrended fluctuation analysis, a multifractal detrended fluctuation analysis (MFDFA) algorithm [24] was proposed by Kantelhardt and Zschiegner et al. Since then, MFDFA has been quickly applied in a number of research areas, including market [25], temperature time series [26], seismic wave signal [27], vibration fault diagnosis [28], image processing [29], etc. For instance, in 2016, Zhao and He [30] analyzed the characteristics of speech signals based on MFDFA with moving, overlapping windows. To solve the shortcomings of the periodic trend affecting the Hurst index estimation in MFDFA, Fourier transform and MFDFA were combined to analyze the fluctuations in a high-frequency power load [31]. Aiming to solve the subjectivity, lack of dynamics and theoretical basis for setting the risk threshold, the authors of Reference [32] proposed a method that combined surrogate data method and MFDFA to determine the prewarning threshold of electric load risk threshold based on historical load data.

The above-mentioned various methods can effectively describe the nonlinear system, especially the multifractal features of time series, but the analysis of time-series signals requires the process of detrending. New pseudofluctuation errors occur during the process. Such errors are mainly caused by two reasons. One is that the sequence is overfitting or underfitting due to the uncertainty of the order in the fitting polynomial function, while the other is that MFDFA uses a uniform sequence to segment data, resulting in discontinuity at sequence segmentation points.

To solve the above-mentioned problem, this paper proposes a new MFDFA algorithm based on empirical mode decomposition (EMD) and template movement. The multifractal characteristics of harmonic signals are analyzed by the new algorithm and a new method for determining the harmonic signals’ characteristics is obtained. Since harmonic signals are susceptible to factors such as randomness, distribution factors and nonstationary factors [33], this paper first uses the wavelet multiresolution method to decompose the harmonic signals and the harmonic components, obtained at different frequencies. After this, the multifractal characteristics of the harmonic components are determined by the novel MFDFA algorithm. Finally, the proposed method and the traditional method are analyzed and compared. The analysis results verify the effectiveness of the proposed method.

2. MFDFA Algorithm Based on Empirical Modality and Template Movement

2.1. Empirical Mode Decomposition Algorithm

The EMD algorithm is an adaptive signal time-frequency processing algorithm [34], which is mainly used for the processing of nonlinear and nonstationary signals [35,36]. The procedure of EMD analysis is described as follows:

• Find all local extreme points in signal $x\left(t\right)$ and connect the extreme points smoothly through the cubic spline function to get the upper envelope $e_{-}\left(t\right)$, lower envelope $e_{+}\left(t\right)$ and the average envelope $m\left(t\right)=(e_{+}\left(t\right)+e_{-}\left(t\right))/2$ of the original signal;
• Find the difference function $z\left(t\right)=x\left(t\right)-m\left(t\right)$;
• Determine whether $z\left(t\right)$ meets intrinsic mode function’s (IMF) conditions or not. If it is not satisfied, repeat the above steps. If it is satisfied, $z\left(t\right)$ is the first IMF and is recorded as $im{f}_i\left(t\right)$;
• Finally, $x\left(t\right)$ is decomposed by EMD into n frequencies from high to low $im{f}_i\left(t\right)$ and a remainder $r_n\left(t\right)$. Essentially, $x\left(t\right)={{\displaystyle \sum }}_{i=1}^{n}im{f}_i\left(t\right)+r_n\left(t\right)$.

2.2. Proposed Novel MFDFA Algorithm

The EMD algorithm successively decomposes n $im{f}_i\left(t\right)$ and one residual $r_n\left(t\right)$ from the signal and the residual part is used to replace the detrend fluctuation item in MFDFA. To be specific, the residual $r_n\left(t\right)$ replaces the detrend polynomial $y_v\left(t\right)$ and the result, after the least squares polynomial fitting, eliminates the residual sequences to obtain the fluctuations polynomial. After this, the fluctuation polynomial is divided into several segments. After calculating the fluctuation function $F_q\left(s\right)$ of the signal according to different $q$ values, the generalized Hurst index of the signal can be obtained by finding the logarithmic least squares fitting slope of $F_q\left(s\right)$ and $q$. Finally, the multifractal spectrum of the signal can be obtained through the Legendre transform. The proposed algorithm is illustrated in the flowchart in Figure 1.

Step-1: Set the time series as $x\left(k\right), k=1,2,\dots , N$ and subsequently, construct a new sequence:

$$Y\left(i\right)={\displaystyle \sum }_{k=1}^i\left(x_k-\overline{x}\right)$$

(1)

where $\overline{x}$ is the mean value of $x_k$, with $\overline{x}=\frac{1}{N}{{\displaystyle \sum }}_{k=1}^N{x}_k$.

Step-2: The traditional MFDFA algorithm divides the new sequences $Y_i$ into $N_s=int\left(N/s\right)$ nonoverlapping intervals, each of which contains $s$ data. Since $N$ may not be divisible by $s$, there will be a residual value. To fit all the data in the sequence in the calculation, another segmentation process is performed from the end of the sequence and thereafter, $2N_s$ equal length segments are obtained.

In this paper, the method of sequence segmenting MFDFA is improved by introducing the template movement and the number of sequences is increased from $N_s$ or $2N_s$ to $N-s+1$.

Step-3: For the $s$ points in each $v \left(v=1,2,\dots ,2N_s\right)$ interval, the least square method is used to fit the polynomial of the order k and the results are obtained:

$$Y_v\left(i\right)=a_1{i}^k+a_2{i}^{k-1}+\dots +a_{k}i+a_{k+1},(i=1,2,\dots ,s;k=1,2,\dots )$$

(2)

This paper introduces the EMD algorithm to replace the least squares fitting, to detrend sequence $Y\left(i\right)$. The remainder of the EMD decomposition is the trend polynomial that reflects the general trend of the signal, with $Y_v\left(i\right)=r_v\left(i\right)$.

Step-4: Calculate the mean square error of $F^2\left(s,v\right)$ by the following:

$$F^2\left(s,v\right)=\frac{1}{s}{{\displaystyle \sum }}_{i=1}^s{\left\{Y\left[\left(v-1\right)s+i\right]-r_v\left(i\right)\right\}}^2$$

(3)

Step-5: For $N-s+1$ segments, find the mean value of $F^2\left(s,v\right)$ and get the $q-order$ fluctuation function $F_q\left(s\right)$ by the following:

$$F_q\left(s\right)=\{\begin{array}{c}{\left\{\frac{1}{N-s+1}{{\displaystyle \sum }}_{v=1}^{N-s+1}{[F^2\left(s,v\right)]}^{\frac{q}{2}}\right\}}^{\frac{1}{q}}, q\ne 0\\ exp\left\{\frac{1}{2\left(N-s+1\right)}{{\displaystyle \sum }}_{v=1}^{N-s+1}ln\left[F^2\left(s,v\right)\right]\right\}, q=0\end{array}$$

(4)

$F_q\left(s\right)$ is a function of data length $s$ and $q-order$. With an increase in $s$, $F_q\left(s\right)$ increases in the power–law relationship

$$s^{h\left(q\right)}\propto F_q\left(s\right)$$

(5)

where $h\left(q\right)$ is the generalized Hurst index.

$h\left(q\right)$ is related to the multifractal quality index $\tau \left(q\right)$ and their relation is described as follows:

$$\tau \left(q\right)=qh\left(q\right)-1$$

(6)

According to the Legendre transform, the relationship between the singular exponent $\alpha$ and the multifractal spectrum $f\left(\alpha \right)$ and $h\left(q\right),$ is obtained in Reference [37]

$$\alpha =h\left(q\right)+q{h}^{\prime }\left(q\right)=h\left(q\right)+q\frac{dh\left(q\right)}{dq}$$

(7)

$$f\left(\alpha \right)=q\left[\alpha -h\left(q\right)\right]+1$$

(8)

From the above equations, three characteristics in judging multifractal features based on the MFDFA method can be obtained [38]:

• Judging by $q$ and $h\left(q\right)$:
  • If $q$ is not related to $h\left(q\right)$, the signal is monofractal.
  • If $q$ is related to $h\left(q\right)$, the signal is multifractal.
• Judging by $q$ and $\tau \left(q\right)$:
  • If $\tau \left(q\right)$ is a straight line, the signal is monofractal.
  • If $q$ and $\tau \left(q\right)$ are nonlinear, the signal is multifractal.
• Judging by $\alpha$ and $f\left(\alpha \right)$:
  • If $f\left(\alpha \right)$ is a constant, the signal is monofractal.
  • If the curve of $\alpha$ and $f\left(\alpha \right)$ has a single-peak bell shape, the signal is multifractal.

2.3. MFDFA Feature Extraction Parameters

The singular exponent $\alpha$ reflects the degree of unevenness of the fractal sequence in the local probability measure distribution [39,40] while the multifractal spectrum width $\Delta \alpha ={\alpha }_{max}-{\alpha }_{min}$ reflects the signal’s intensity of the multifractal property. A larger $\Delta \alpha$ highlights stronger multifractal characteristics in the signal and more severe fluctuations in the signal. The singular exponent ${\alpha }_0,$ corresponding to the maximum value of multifractal spectrum $f\left(\alpha \right)$, reflects the randomness of the signal. ${\alpha }_{min}$ is the smallest singular exponent, which reflects the intensity of local changes in the signal. A smaller ${\alpha }_{min}$ indicates a stronger local singularity in the signal and a more intense local variation in the signal.

Therefore, based on the above analysis, $\Delta \alpha$, ${\alpha }_0$ and ${\alpha }_{min}$ can be used as the feature quantities after extracting the multifractal property of harmonic signals.

3. Signal Acquisition and Analysis

3.1. Signal Acquisition

In this paper, the electromagnetic flow meter contaminated by power harmonics was taken as the study object and the fractal characteristics of power system harmonics were analyzed. The advantages of the new MFDFA algorithm in the harmonic analysis of the power system were studied.

The electromagnetic flow meter adopted Faraday’s law of electromagnetic induction and the conductor moving in the magnetic field generated an induced voltage. As shown in Figure 2, in the acquisition system, an alternating magnetic field was generated on the coil by controlling Switch 1, Switch 2, Switch 3 and Switch 4. Furthermore, the voltage induced by the flowing liquid was connected to the signal amplifying circuit through ELECTRODE+ and ELECTRODE–. Finally, the amplified signal was sent to the computer through the signal acquisition system after the analysis. The system adopted a low-frequency rectangular wave excitation mode, with the excitation frequency at 6.25 Hz and the excitation current at 25 mA. It is important to note that the induced voltage between ELECTRODE+ and ELECTRODE– is extremely susceptible to interference from power system harmonics as its amplitude is generally in the range of tens of uV and several mV.

Figure 3 shows a waveform of the voltage after removing the DC offset signal between ELECTRODE+ and ELECTRODE–. The signal acquisition system adopted a sampling rate of 500 Hz, and the spectrum of the acquired signal was calculated and is presented in Figure 4. The spectrogram indicated that the signal contained obvious harmonic components. In addition to the fundamental component, the harmonic components were dominated by second-order (100 Hz) and fourth-order (200 Hz) harmonics.

The original signal in Figure 3 was decomposed by the Mallat algorithm [41,42,43] and the actual frequency decomposition layer $p$ was:

$$p=lo{g}_2\left(\frac{f_s}{2f_0}\right)-1,$$

(9)

$f_0$ was selected as 6.25 Hz and $f_s$ was the sampling frequency of 500 Hz. The db24 wavelet was used to decompose the signal into 5 layers to obtain the high-frequency part of the original signal.

The total bandwidth of the signal was 250 Hz and the frequency band range is shown in

Table 1. The frequency band of each layer by wavelet decomposition.

Layer	Frequency Band
a5	0–7.8125 Hz
a4	0–15.625 Hz
a3	0–31.25 Hz
a2	0–62.5 Hz
a1	0–125 Hz
d5	7.8125–15.625 Hz
d4	15.625–31.25 Hz
d3	31.25–62.5 Hz
d2	62.5–125 Hz
d1	125–250 Hz

. The approximate signal (a5–a1) c­onsisted of the low-frequency components of the signal and the detail signal (d5–d1) consisted of the high-frequency components of the signal.

The decomposition coefficient d3 contained the fundamental wave, as shown in Table 1, the decomposition coefficient d2 contained second-order harmonics, the decomposition coefficient d1 contained fourth-order harmonics and the waveform of d3–d1 is shown in Figure 5.

3.2. Signal Analysis

Three components, d1, d2 and d3, were selected to calculate and analyze the local fractal characteristics by using the traditional MFDFA algorithm [44] and the proposed MFDFA algorithm.

Figure 6a,b shows the multifractal spectrum of the three signal components after wavelet decomposition, using the traditional MFDFA algorithm and the new MFDFA algorithm. $\alpha$ and $f\left(\alpha \right)$ are single-peak bell-shaped graphs. The slope of the curve is obvious and the value of the singular index $\alpha$ is not unique. Therefore, all three signal components have typical multifractal characteristics.

Figure 7a,b is the $q-h\left(q\right)$ curves of the three signal components, derived from the traditional MFDFA algorithm and new MFDFA algorithm, respectively. $h\left(q\right)$ of the three signal components gradually decreased with an increase in $q,$ as seen in in Figure 7, and the value of $q$ varied with $h\left(q\right)$. In other words, $q$ was related to $h\left(q\right)$. The local structure of the fundamental signal, both the second-order harmonic and fourth-order harmonic, were not uniform and therefore indicated the existence of obvious multifractal characteristics. From the comparison displayed in Figure 7a,b, the curve of the traditional MFDFA algorithm had obvious bumps and inflection points, and the curve using the new algorithm was smoother and more stable than the traditional algorithm. Moreover, the density of the three curves of the new algorithm was also much lower than that of the traditional algorithm. This indicated that the proposed algorithm reduced the pseudofluctuation error, caused by the discontinuity of the data segmentation and obtained a more stable fluctuation function.

Figure 8a,b exhibits the $q-\tau \left(q\right)$ curves of the three signal components, obtained using the conventional MFDFA algorithm and the new MFDFA algorithm, respectively. $\tau \left(q\right)$ is a curve where $q$ and $\tau \left(q\right)$ have a nonlinear relationship. Furthermore, this was a convex-increasing function, which proved that the fundamental signal, the second-order harmonic and the fourth-order harmonic signals had multifractal characteristics. When $q$ was greater than 0, the $\tau \left(q\right)$ values of the three components were more densely distributed. When $q$ was less than 0, the $\tau \left(q\right)$ value distribution of the three signals was relatively sparse, and the $\tau \left(q\right)$ value of the high-frequency component d1 was lower than the $\tau \left(q\right)$ value of the high-frequency component d2. Furthermore, the $\tau \left(q\right)$ value of the low-frequency component d3 signal was the largest. From comparing Figure 8a,b, we can see that the $q-\tau \left(q\right)$ curve of the new algorithm was smoother than the traditional algorithm, which verified the superiority of the new algorithm.

The fluctuation function was calculated according to the $q$ values of different fluctuation orders, where $q$ ranged from −5 to 5 and different $q$ values were calculated to obtain different $log{F}_q\left(s\right)-log\left(s\right)$ scatter plots. All $log{F}_q\left(s\right)-log\left(s\right)$ calculated at all q values are plotted on Figure 9, Figure 10 and Figure 11. Figure 9a with Figure 9b, Figure 10a with Figure 10b and Figure 11a with Figure 11b are $log{F}_q\left(s\right)-log\left(s\right)$, double logarithmic graphs corresponding to different $q$-valued harmonic signal components obtained by the traditional MFDFA algorithm and the novel MFDFA algorithm, respectively.

In Figure 9, Figure 10 and Figure 11, the fluctuation function $F_q\left(s\right)$ that was obtained by the two algorithms, changes with the scale $s$, which indicated that the signal had multifractal characteristics. Using the traditional MFDFA algorithm and the new MFDFA algorithm to analyze each harmonic component, the extracted multifractal characteristics parameters $\Delta \alpha$, ${\alpha }_0$ and ${\alpha }_{min}$ are listed in

Table 2. $\Delta \alpha$, ${\alpha }_0$ and ${\alpha }_{min}$ values of the harmonic component computed analytically through traditional MFDFA and new MFDFA.

Harmonic Component	Traditional MFDFA			New MFDFA
	$\Delta \alpha$	${\alpha }_0$	${\alpha }_{min}$	$\Delta \alpha$	${\alpha }_0$	${\alpha }_{min}$
d1	0.957	0.113	−0.126	0.397	0.066	−0.122
d2	0.782	0.128	−0.039	0.342	0.051	−0.071
d3	0.969	0.324	0.103	0.800	0.061	−0.074

.

In Table 2, the $\Delta \alpha$ value of d2 is the smallest and the $\Delta \alpha$ value of d3 is the largest. This represents that the singular exponential distribution range of the second-order harmonic component was the narrowest and the singular exponential distribution range of d1 was wider than that of d3. It shows that the fractal distribution of the second-order harmonic component signal was better than that of the fundamental signal and the fourth-order harmonic. The fourth-order harmonic showed better uniformity than the fundamental signal. Compared with the second-order harmonic, the fundamental wave of the fourth-order harmonic had more obvious multifractal features and the fluctuations were more obvious with scale changes. The ${\alpha }_{min}$ of d1 is the smallest, indicating that the harmonic signal with the highest frequency had the strongest local fluctuation. It was known that the multifractal intensity of the fundamental signal was the largest and the signal distribution is relatively uniform. The multifractal intensity of the fourth-order harmonic signal was larger than that of the second-order harmonic and the signal distribution was the most uneven. These features can be well reflected in Figure 5.

Compared with the traditional MFDFA method, the characteristic extraction parameter values obtained by the new MFDFA method showed a significant decreasing trend. The multifractal spectrum curve distribution was narrowed in the new method, indicating that the new algorithm can reduce the error in the multifractal-characteristics analysis of harmonic signals and thus, make the curve fitting more stable.

The $h\left(q\right)$ values obtained by analyzing the harmonic components using traditional MFDFA and the new MFDFA are presented in

Table 3. h(q) values of the harmonic signal computed analytically through MFDFA and the new MFDFA.

$\mathit{q}$	Traditional MFDFA			New MFDFA
	d1	d2	d3	d1	d2	d3
−5	0.583	0.495	0.747	0.171	0.151	0.395
−4	0.521	0.433	0.666	0.142	0.121	0.313
−3	0.426	0.336	0.553	0.116	0.094	0.202
−2	0.258	0.192	0.418	0.091	0.072	0.098
−1	0.113	0.128	0.324	0.066	0.051	0.061
0	0.067	0.097	0.272	0.042	0.033	0.041
1	0.034	0.071	0.238	0.02	0.016	0.022
2	0.006	0.049	0.212	0	0.001	0.007
3	−0.017	0.032	0.192	−0.019	−0.012	−0.006
4	−0.037	0.018	0.175	−0.036	−0.023	−0.018
5	−0.055	0.006	0.160	−0.051	−0.032	−0.029
$\Delta \mathrm{h}$	0.638	0.489	0.587	0.222	0.183	0.4244

.

The $h\left(q\right)$ values corresponding to all $q$ values in the Table 3 were less than 0.5 for the new MFDFA algorithm, which indicated that the low-frequency harmonic signal component and the high-frequency harmonic signal component were anticorrelated. The $\Delta h$ of d2 was the smallest and the $\Delta h$ value of d3 was the largest. Thus, if $\Delta h$ is larger, the multifractal is stronger, which resulted in the multifractal degree of the d2 component being the smallest and the multifractal degree of the fundamental signal being the strongest. Comparing the $\Delta h$ values obtained by the traditional MFDFA algorithm and the new algorithm, the results obtained by the new MFDFA algorithm were significantly smaller than the traditional algorithms, which indicated that the new algorithm could reduce the pseudofluctuation error and had higher precision than the traditional algorithm. The new MFDFA could more accurately characterize the multifractal features of the signal.

By using the characteristics of multifractal, the h(q) value distribution of different signals can be calculated under the condition of a fixed range of q values to measure the irregularity of multifractal in different frequency harmonic signals.

4. Conclusions

By using the new MFDFA algorithm proposed in this paper to analyze the multifractal characteristics of harmonic signals, the following conclusions are drawn:

• The power grid harmonic signals in the flow meter signal exhibit multifractal characteristics. Furthermore, the multifractal intensity of the fundamental signal is the largest and the multifractal intensity of the higher-order harmonic is larger than that of the lower harmonic.
• Compared with the traditional MFDFA algorithm, the new algorithm can effectively reduce the pseudofluctuation error caused by the discontinuity of the traditional algorithm, making the fitting curve more stable and more accurately revealing the multifractal characteristics of harmonic signals.
• $\Delta \mathsf{\alpha }$, ${\mathsf{\alpha }}_0$, ${\mathsf{\alpha }}_{\min }$ and $\mathrm{h}\left(\mathrm{q}\right)$ can provide theoretical and algorithmic support for grid harmonic management.
• Although the algorithm shows good performance in the multifractal-characteristics analysis of harmonic signals, the new algorithm is mainly for integer subharmonics. Using the algorithm for analyzing noninteger harmonics still remains an interesting issue that needs to be studied.