An Improved Empirical Mode Decomposition Based on Local Integral Mean and Its Application in Signal Processing

Empirical mode decomposition (EMD) is an effective method to deal with nonlinear nonstationary data, but the lack of orthogonal decomposition theory and mode-mixing are the main problems that limit the application of EMD. In order to solve these two problems, we propose an improved method of EMD. The most important part of this improved method is to change the mean value by envelopes of signal in EMD to the mean value by the definite integral, which enables the mean value to be mathematically expressed strictly. Firstly, we prove that the signal is orthogonally decomposed by the improved method. Secondly, the Monte Carlo method of white noise is used to explain that the improved method can effectively alleviate mode-mixing. In addition, the improved method is adaptive and does not need any input parameters, and the intrinsic mode functions (IMFs) generated from it is robust to sifting. We have carried out experiments on a series of artificial and real data, the results show that the improved method is the orthogonal decomposition method and can effectively alleviate mode-mixing, and it has better decomposition performance and physical meaning than EMD, ensemble EMD (EEMD), and complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN). In addition, the improved method is generally more time-consuming than EMD, but far less than EEMD and CEEMDAN.

However, because EMD is an experience-based method, there are many problems [11]. ese problems include the mixing of multiscale modes, and the sifting process for intrinsic mode functions (IMFs) requires an appropriate stopping criterion, the end effect, and the orthogonality of IMFs. Among them, EMD has been criticized by signal processing experts because of the lack of orthogonal decomposition theory, which has become a major problem restricting the application of EMD [12,13]. Various research studies based on the orthogonal or approximate orthogonal decomposition methods of EMD [12,[14][15][16][17] are not convincing. References [18,19] use Schmidtʼs formula to make IMFs strictly orthogonal, but it cannot guarantee that the IMF is obtained from the orthogonal decomposition of the signal. Mode-mixing between IMFs is another major problem that limits EMD applications. Although a large number of improved methods of EMD are devoted to solving this problem, they often need to determine the input parameters a priori. Inappropriate parameters will reduce the accuracy of decomposition [20] and even lead to decomposition failure. For example, ensemble EMD (EEMD) [21], complementary EEMD (CEEMD) [22], and complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [23] effectively alleviate mode-mixing, but need to set parameters of the amplitude and quantity of auxiliary noise; moreover, their related procedures are time-consuming, which is often intolerable for large-length signal; empirical wavelet transform (EWT) [24], variational mode decomposition (VMD) [25], and the methods based on Hankel matrix eigenvalue decomposition (EVD) [26] build EMD on additional classical theoretical methods to alleviate mode-mixing, but they also need to carefully select appropriate parameters, especially the number of mode components.
In fact, the orthogonality and mode-mixing between IMFs are related. If the signal will be orthogonally decomposed and the mode components will orthogonal, mode-mixing will be effectively alleviated. For this reason, an improved EMD method is proposed in this study. e most important change in the improved method is to change the envelope mean method (in EMD) [1] to the integral mean value theorem (IMVT) to find the data's mean, which makes the mean be strictly expressed mathematically. Although some EMD improvement methods have tried to change the mean method, none of these methods can be expressed in a strict analytical formula [27][28][29][30]. Local mean decomposition (LMD) constructs a mean curve with the mean of adjacent extreme points of signal, but this will lose the detailed fluctuation information between these two adjacent extreme points [31]. We call the improved method the integral mean mode decomposition (IMMD). In the IMMD method, we try to prove the orthogonality of IMFs and show that mode-mixing is effectively alleviated. In addition, IMMD uses spline interpolation to directly predict the mean value at both ends of the data and uses a stop criterion with a fixed number of sifting [32] to reduce complexity. e rest of the manuscript is organized as follows: in Section 2, we give in detail the steps of the IMMD method and the corresponding flowchart (Figure 1), the stopping criterion of IMF, and the method of restraining the end effect. And in the part of the stopping criterion, we will take Gaussian white noise as an example to analyze the robustness of IMFs from IMMD to sifting. In Section 3, a widely used nonlinear nonstationary signal is used as an example to prove the orthogonality of IMMD decomposition in three steps. In Section 4, the Monte Carlo method is used to show that IMMD can alleviate mode-mixing than EMD. In Section 5, a series of artificial and real signals are used to verify that the IMMD method is the orthogonal decomposition method and can alleviate mode-mixing. Finally, discussion and conclusion are given in Section 6 ( Figure 1).

Steps
(1) e original signal (or protomode function (PMF) [33]) is divided into many local parts by all extrema. (2) For each local part, the value of the local mean is calculated by IMVT, and its position is fixed at the midpoint of the local time series. All the local means are interpolated by cubic splines to get the mean curve (in EMD, the mean curve of the signal is obtained by the difference between the upper and lower envelopes, and the upper (lower) envelope is obtained by smoothing the maxima (minima) of the signal. e remaining steps of EMD are the same as those of IMMD (see Figure 1)). (3) PMF is obtained by subtracting the mean curve from the signal. (4) Repeat steps (1)-(3) on PMF (as a signal) iteratively, until PMF satisfies the stoppage criterion. is process of subtracting the mean iteratively from the signal is called sifting. When the sifting is over, PMF is IMF 1 . (5) New signal is obtained by subtracting IMF 1 from the original signal, and IMF 2 is obtained by repeating steps (1)-(4) on the new signal. (6) Repeat steps (1)- (5) to get other IMFs of the original signal.
e flowchart of IMMD based on EMD is shown in Figure 1. e method to construct the mean curve of signal in [34] is the same as that in IMMD. However, in [34], the definition of "centroid" is used to determine the time series position of the local mean t, which makes the time series position of t be restricted by the signal amplitude. erefore, IMMD uses the simplest time series midpoint as the time series position of t.

2.2.
e Stoppage Criterion. ere are two purposes of sifting: to eliminate riding waves and to make the wave profiles more symmetric [1]. e commonly used sifting stoppage criterion was called the Cauchy type of stoppage criterion [11], and the equation is as follows: . (1) Among them, PMF k is the PMF obtained by k-th sifting. If SD < 0.2∼0.3, sifting stops. In fact, although 0.2∼0.3 as the SD threshold is accepted by most users of EMD, it is controversial [11]. It is generally accepted that, under a limited sifting number, the smaller the SD value, the better. In order to reduce the complexity, EMD uses the stoppage criterion of fixing sifting number to 10, which also satisfies the Cauchy type of stoppage criterion [32].
Let IMMD adopt the stoppage criterion of fixed sifting number and decompose a randomly generated Gaussian white noise. e SDs of all IMF vary with the sifting number, as shown in Figure 2(b). It can be seen that SD 1 ∼SD 14 shows fast and stable strong convergence. When the sifting number >6, all SD values <0.01; when the sifting number >9, all SD values <0.001; and when 100 < the sifting number <5000, almost all SD values <1 × 10 −30 .
Let EMD decompose the same Gaussian white noise, and the SDs of IMF vary with the sifting number, as shown in Figure 2(a). Almost all SDs show weak convergence of slow oscillation. Except for SD 8 If PMF satisfies the stopping criterion ?
If r is a monotonic curve ? End Using integral to find the local mean between two successive extremes All the local mean points are smoothed by the spline to get the mean curve m(t) EMD IMMD    values < 0.001; and when the sifting number < 5000, almost all SD values are still more than 1 × 10 −10 .
e results of the two methods show that the sifting process of IMFs in IMMD is stable and can meet the Cauchy stop criterion very well. In addition, we have tested all the signals decomposed by IMMD in this study about sifting, and the SD curve of them is almost exactly similar to that of Gaussian white noise, and there is no qualitative difference. erefore, we recommend that the IMMD method adopts the stop criterion of fixed sifting number.

Method of Alleviating End Effects.
End-effect processing methods are essentially predicting the end mean of data [21]. e error caused by prediction will spread from the end to the interior of the data with the sifting, so it will affect the stability of IMF. We propose a method of predicting the end mean which is suitable for IMMD: using the spline interpolation of the data's mean to directly predict the end mean of the data.
Compared with the other prediction methods, the proposed method has the following characteristics: (1) it directly predicts the end mean; (2) it does not need to introduce additional end-effect processing methods; and (3) it ensures that IMFs have good robustness to the sifting. We present the SD curves (see Figure 3) obtained when Gaussian white noise is decomposed by IMMD with enlargement of the ends (one common method to alleviate end effects in EMD) [21], and make a simple visual comparison between Figures 3 and 2(b). e details of the comparison are similar to that in Section 2.2, and there is no essential difference, so we will not repeat it here. e comparative analysis between the proposed method and other methods is beyond the scope of this study.

The Proof of IMMD Orthogonal Decomposition
Let IMMD decompose the signal x(t) into a total of n components of IMF 1 , . . ., IMF j , . . ., IMF (n − 1), r n , where IMF j is the j-th mode component and its corresponding residue is r j . e proof of the orthogonal decomposition of IMMD is divided into three steps: first, IMF j is orthogonal to r j ; second, IMF j is orthogonal to r k (j < k ≤ n); and third, IMF j is orthogonal to IMF k .

e Proof of Orthogonality of IMF j and r j
3.1.1. e Proof. Take IMF 1 as an example. Suppose that PMF k (k � 0, 1, . . ., PMF 0 � x(t)) is obtained after the k-th sifting and the number of its extrema is n, let PMF k � h k (t), h k (t) can be described as follows: where h ki (t) is the local part of h k (t) from t i to t i+1 , t i is the position (on the time series) of the i-th extremum of h k (t) (see Figure 4). Let m k (t) is the mean curve of h k (t); m k (t) can also be described as follows: According to the IMVT, we have 10 100 1000 5000 e number of si ings SD 1 SD 2 SD 3 SD 4 SD 5 SD 6 SD 7 SD 8 SD 9 SD 10 SD 11 SD 12 SD 13 SD 14 SD 15   where m ki (t) is a constant, described as m ki (see Figure 4). After one shifting on h k (t), we get h k+1 (t) � h k (t) − m k (t), and h k+1 (t) can also be described as follows: h l (t) is obtained after the l-th (l ≤ k) sifting, and its mean curve is m l (t). m l (t) can also be described as m l (t) � [m l1 (t), . . . , m li (t), . . . , m l(n−1) (t)], where m li (t) is also a constant, described as m li (see Figure 4). en, we have erefore, we have i.e., h k+1 (t)⊥m l (t), (l � 0, 1, . . ., k; m 0 (t) is the mean curve of PMF 0 ). Here, it should be noted that the sifting will produce new extrema (the result of overshoot and undershoot), but the new extrema have no effect on the above conclusion. e reason is that the mean curve m li (t) is a constant in any local part of h k+1 (t) (take the annual global surface temperature anomaly (GSTA) [35] as an example ( Figure 5): the sifting produced a new extremum at 1912 in PMF 3 , but in the local parts [1911,1912] and [1912,1913], m 0 (t), m 1 (t), andm 2 (t) are all constant).

(9)
In the same way, it is easy to get erefore, any IMF is orthogonal to its corresponding residue.

e Effect of Sifting on the Orthogonality of IMF j and r j .
Although the new extremum produced by sifting has no effect on the above conclusion, the spline interpolation of all local means in the sifting process will affect the orthogonality of IMF and its corresponding residue. In the above proof, the spline interpolation of all local means is ignored. e explanation is as follows: suppose mean 2 (t) is the mean curve of the signal obtained by spline interpolating all local means of the signal. e error of the mean curve is the difference between mean 2 (t) and the stair step curve of mean, mean 1 (t) (the default mean curve in the above proof). Similar to the SD of two continuous PMFs, the standard deviation of mean is used to evaluate this error, which is defined as follows: Take the IMF 1 of GSTA as an example (see Figure 6). With the increase of the sifting number, the SD value increases linearly and reaches the maximum 0.028 after 20 siftings. For GSTA with low sampling frequency, 0.028 is enough small (in practice, the error can be reduced by increasing the sampling frequency). erefore, the error of the mean curve has a small effect on the orthogonality, and it can be ignored in theoretical proof.
3.2. e Proof of Orthogonality of IMF j and r k (j < k ≤ n) e average period of the IMF j (j � 1, 2, ...) from EMD and its improved methods will increase with the increase of order j. which means In particular, for Gaussian white noise, T IMF(j+1) � 2T IMF j [32,36,37]. Let L j be any local time region [t i , t i+1 ] of IMF j (t i is the time series point of the i-th extremum of IMF j ) and L k be any local time region [t i′ , t i′+1 ] of IMF k (t i′ is the time series point of the i′-th extremum of IMF j ). Generally, if L j ∩ L k ≠ 0, then according to equation (12), t i′ ≤ t i and Mathematical Problems in Engineering t i+1 ≤ t i′+1 , that is, L j ⊂ L k . Taking GSTA as an example, when the fixed sifting number is 1, IMMD decomposes into IMF 1 ∼IMF 7 of GSTA, as shown in Figure 7. In most local parts, L j ≤ L k (a typical part is marked in the black dotted box: L 1 < L 2 <· · ·<L 7 ), but there are a few local parts (marked by red box), L 1 > L 2 or L 1 > L 3 . Section 3.1 has proved that ignoring the smoothing of the mean, the residue r k of IMF k is a constant within L k . Because L j ⊂ L k , r k should also be a constant in L j . Due to the importance of this conclusion, we take GSTA as an example to explain and verify. First of all, in order to reduce the influence of the spline smoothing mean curve, and to better illustrate the above conclusion, the sifting number is fixed at 1. Figure 8 gives all the IMF j curves and r k (k > j; j, k � 1∼6) step curves of GSTA. It can be seen in Figure 8 that r k is a constant in most locals except for a few locals (marked by color dotted line in Figure 8) of IMF j .
Data sampling and spline interpolation produce a few locals of IMF j where r k is not constant. Discrete sampling of data and local definite integral to obtain the mean value makes the mean curve of the data discontinuous and needs to be smoothed. Smoothing may cause a small number of extrema to move or even generate new extrema. For example, in Figures 8(e) and 8(f ), r 5 becomes rs 5 after smoothing. As a new data, rs 5 is decomposed by IMMD to get IMF 6 . Let r 6 be the step mean curve of rs 5 . If r 6 is not smooth, r 6 is constant in the two parts of rs 5 -r 6 (should be equal to IMF 6 ); however, r 6 must be smoothed, so the extremum of real IMF 6 has shifted slightly to the right, which causes r 6 to be no longer a constant in the first local of IMF 6 .
After the above analysis, it can be determined that if the smoothing of the mean curve is ignored, not only r j but also r k (k > j) are constant in any local of IMF j (j � 1, 2, ..., n). In addition, according to equation (6), in any local of IMF j , the definite integral of IMF j is 0, so the inner product of IMF j and r k is 0. erefore, globally, IMF j and r k are orthogonal, which is

Index of Orthogonality of IMF j and r k .
Because of the smoothing of the mean curve, the finite length of the data will cause leakage, and equations (10) and (13) are not strictly equal to 0. erefore, define an index for checking the orthogonality of IMF j and r k (k ≥ j) as When the sifting number is fixed to 1, the IRO jk of IMF 1 ∼IMF 6 and r 1 ′∼r 6 ′ (sum of step mean curves) and the IRO jk of IMF 1 ∼IMF 6 and r 1 ∼r 6 (sum of smooth mean curves obtained by smoothing the step mean curves) of GSTA are shown in Table 1. Overall, these values verify IMF j and r k is orthogonal, and the smoothing of the mean curve reduces the orthogonality of IMF j and r k .

e Effect of Sifting on the Orthogonality of IMF j and r k .
In the process of sifting, the influence on the orthogonality of IMF j and r k includes the smoothing of the mean curve of PMF and the symmetry of PMF with respect to the time axis. It has been explained and verified in Sections 3.2.1 and 3.2.2 that smoothing the mean of PMF will reduce the orthogonality of IMF j and r k . However, when the PMF is more symmetric about the time axis, its local integral tends to 0, so the orthogonality is stronger. As the sifting number increases, the mean value of PMF and the change of PMF both tend to 0, and the orthogonality no longer changes. erefore, sifting makes the orthogonality index oscillate and eventually tend to be constant. Figure 9 shows the curve of the IRO 1k value of IMF 1 and r 1 ∼r 6 of GSTA with the change of fixed sifting number (1∼5000). It can be seen visually that when the fixed sifting number is less than 20, the IRO 1k value oscillates; when it is greater than 20, it quickly converges to a fixed value. Table 2 shows the IRO jk value of IMF j and r k (k ≥ j; j, k � 1∼9) when the fixed sifting number � 10. Compared with Table 1, it can be seen that sifting enhances orthogonality.    (10) and (13), it can be obtained that . . . Kelvin 1 (f ) Figure 8: IMF j (j � 1∼6) and their corresponding residual components r j of GSTA. In the locality of IMF j , the variable r k (k � j, j + 1, . . ., 6) is marked by the corresponding color dotted line.
Mathematical Problems in Engineering 9 And it is equivalent to: So far, we have proved that any two IMF components are orthogonal.

3.3.2.
Index of Orthogonality of IMF j and IMF k . Similar to equation (13), equation (16) is not strictly equal to zero. Reference [1] has defined an index to check the orthogonality of any two IMFs: When the fixed sifting number is 1, the IO jk values of IMF 1 ∼IMF 7 of GSTA are shown in Table 3.
ere are 5 values with the order of magnitude of 10 −1 in Table 3, mainly due to the poor symmetry of the IMFs obtained by sifting once (especially IMF 4 and IMF 5 in Figures 8(d) and 8(e)). erefore, the integral of IMF in each local is 0 and cannot be satisfied very well.

e Effect of Sifting on the Orthogonality of IMF j and IMF k .
With the increase of sifting number, the symmetry of IMF and the orthogonality of IMFs increased. Figure 10 shows that the IO 1k of GSTA varies with the fixed sifting number (1∼5000). It can be seen that, with the increase of sifting number, the values of IO 1k oscillate but converge rapidly to the constant. In fact, Figure 10 is consistent with Figure 9, and it shows that IMF from IMMD has good robustness to sifting.
When the fixed sifting number is 10, the IO jk values of IMF 1 ∼IMF 10 of GSTA are shown in Table 4. In Table 4, the five IO jk values with an order of magnitude of 10 −1 (in Table 3) do decrease. In addition, the increase in the number of IMFs also shows the performance of orthogonal decomposition. erefore, sifting can enhance the orthogonality, which is consistent with the conclusion of Section 3.2.3.
In addition, reference [1] also has defined an overall index of orthogonality (IO): .
(18) Figure 11(a) shows the IO value of GSTA, which varies with the fixed sifting number (1∼5000). When sifting number <20, IO value oscillates; when sifting number ≥20, IO � 5.1 × 10 −2 . Figure 11(b) shows the number of IMFs of GSTA, which varies with the fixed sifting number (1∼5000). Compared Figure 11(a) with Figure 11(b), the IO value is approximately inversely proportional to the number of IMFs, which indicates that the number of IMFs increases, the IO value is smaller, and the orthogonality is stronger.

Alleviation of Mode-Mixing
IMFs from IMMD are orthogonal to each other, so it can effectively alleviate the mode-mixing. Using the Monte Carlo method of Gaussian white noise, references [32,36,37] indirectly explain that EMD, EEMD, B-spline interpolationbased EMD (B-EMD) [38], and trigonometric cardinal spline interpolation-based EMD (C-EMD) [39] are all equivalent to a set of binary filter banks with a constant quality factor Q. is section will use the same method to illustrate that IMMD can effectively alleviate the modemixing.
Let IMMD and EMD decompose 5000 samples of Gaussian white noise with length � 512, mean � 0, and variance � 1, respectively. IMMD (EMD) decomposes each Gaussian white noise into at least 10 (8) IMFs. e average Fourier power spectrum of the corresponding order IMF of all samples is shown in Figure 12. e IMF 1 from IMMD or EMD is equivalent to a high-pass filter, and the other IMF is equivalent to a set of overlapping band-pass filters. e center frequency of the latter band-pass filter from IMMD (EMD) is approximately 2/3 (1/2) of the center frequency of the previous band-pass filter. erefore, the decomposition of Gaussian white noise by the IMMD method is also equivalent to filter banks and has the characteristic of trisection.
Except for IMF 1 , the bandwidth of the equivalent bandpass filter of the latter IMF from IMMD (EMD) is roughly 2/ 3 (1/2) times that of the previous IMF, so the equivalent band-pass filter bank has a constant-Q property. It means that the power spectra of IMFs have self-similarity [36]:  Figure 9: IRO 1k (k � 1∼6) of IMF 1 and r k of GSTA which is decomposed by IMMD with the fixing sifting number (1∼5000).      Among them, S k′ (f) (S k (f)) is the power spectrum of the k ′ -th (k-th) IMF; ρ is a constant, for IMMD, ρ � (3/2); for EMD, ρ � 2. Based on equation (5), the power spectra of all IMFs are standardized and collapsed and coincident into a single curve [36]. For EMD, the mode-mixing causes the band-pass filter to mix more low-frequency components. erefore, the filter has poor symmetry about the center frequency and poor narrow-band characteristics, and the ρ value needs to be adjusted to 2.01 [40]. IMMD does not have these problems, and the power spectra normalization result of IMFs is shown in Figure 13.
When IMMD decomposes Gaussian white noise, its equivalent filter bank has the characteristic of trisection. e number of equivalent filters from IMMD is more than that from EMD, and the narrow-band characteristics and constant-Q properties of the filter bank are better than EMD.
erefore, it has better multiresolution analysis capabilities and can effectively alleviate mode-mixing.

Experiments and Results
As mentioned in the Introduction, the orthogonality of signal decomposition will lead to the alleviation of mode-mixing, and we believe that they are unified. ey are not only jointly displayed in the orthogonal index, but more importantly reflected in the decomposition capabilities of good decorrelation, energy concentration, and better physical meaning. IMMD has a good orthogonal decomposition theory, which can effectively alleviate mode-mixing. erefore, IMMD has good orthogonal decomposition capabilities.
In this section, we will conduct experiments on the decomposition of a series of test signals. Firstly, we start with multicomponent nonstationary signals. Secondly, in order to highlight the capabilities and shortcomings of EMD, the authors in [41] successfully used a two-tone signal, and it is often used in many literature studies [15,25,27,28,41]. erefore, the two-tone signal is used to evaluate the IMMD method. In addition, EMD has a critical frequency ratio of 0.67 that cannot successfully decompose two-tone signal [41], so we specially construct a three-sine superimposed signal with a frequency ratio of 1 : 0.67 : 0.5 and compare IMMD with EMD to evaluate the IMMD method. Finally, two real complex signals are used to test and evaluate IMMD. One is the annual GSTA data, which is one of the most widely studied nonlinear nonstationary climate time series. It can explain the general method well and is often used by Huang to illustrate the performance of EMD and EEMD methods [42][43][44]. e other is the electrocardiogram (ECG) signal provided by Moody and Mark [45], which is often used by researchers to illustrate the performance of signal processing methods [46][47][48][49][50][51]. For these two realistic nonstationary and nonlinear signals, we compare the   performance of EMD, EEMD, CEEMDAN, and IMMD methods at the same time and demonstrate the ability of the IMMD method to orthogonally decompose and alleviate mode-mixing.

Multicomponent Nonstationary Signals
5.1.1. Example 1. Signal 1 (Sig 1 ) is the superposition of two different frequency-modulation (FM) signals: Let IMMD and EMD decompose Sig 1 into Sig 3 mode components. e waveforms and time-frequency spectra of these components are shown in Figure 14. Figure 14(a) shows that IMF 1 and IMF 2 from EMD have severe waveform distortion. In the frequency domain (see Figure 14(c)), IMF 1 and IMF 2 from EMD are mixed with more frequency components; at t � 0.7 s, mode-mixing occurs between these two IMFs. IMMD can separate mode components well, and only the frequency mixing occurs at the end of IMF 2 (see Figures 14(b) and 14(d)). Compared with EMD, IMMD alleviates mode-mixing, separates two FM signals well, and shows better orthogonal decomposition capability than EMD. Table 5 shows the values of orthogonal indexes of Sig 1 decomposed by two methods. ese values verified the theory that IMFs from IMMD are orthogonal, and IMFs from EMD are posteriori orthogonal (the orthogonality of the EMD method cannot be proved, and it can only be verified by the IO value of data [1]). In addition, the three IO jk values of IMMD are all smaller than those of EMD, but the IO value is larger than that of EMD. Overall, IMMD has better orthogonality than EMD.

Example 2.
Signal 2 (Sig 2 ) is the superposition of three monocomponent signals including a FM signal, an amplitude-modulated (AM) signal, and a linear signal: Let IMMD and EMD decompose Sig 2 into Sig 3 IMFs, respectively, and waveforms and time-frequency spectra of the components are shown in Figure 15. e characteristic of Sig 2 is that the frequency of the FM signal is equal to the frequency of amplitude change of the AM signal when t ⟶ 0. When 0 < t < 2, two mode components from EMD have serious waveform distortion (see Figure 15(a)). In the frequency domain, when t < 2, IMF 1 and IMF 2 from EMD are mixed with more frequency components; when t � 0.17 s and 0.6 s, mode-mixing occurs between IMF 1 and IMF 2 (see Figure 15(c)). IMF 1 from IMMD restores the FM signal well; except for the end, IMF 2 also restores the AM signal well, but more frequency components are mixed (see Figure 15(d)). In addition, the residual components from the two methods  Figure 13: Collapse and coincidence of the average power spectrum of IMFs from EMD (a) and IMMD (b) based on equation (19).

Mathematical Problems in Engineering 13
restore the linear signal well. erefore, IMMD alleviates mode-mixing. All IO jk and IO values of Sig 2 decomposed by IMMD are smaller than that by EMD (see Table 6). It is verified that IMFs from IMMD are orthogonal, and IMMD has better orthogonality than EMD.

Two-Tone Signal.
Compared with EMD, this section uses a commonly used discrete-time two-tone signal to evaluate the decomposition capability of IMMD. e twotone signal is as follows: where a ∈ [0.01, 100] is the amplitude ratio and f ∈ (0, 1) is the frequency ratio of low-frequency (LF) component, a cos(2πft + φ), and high-frequency (HF) component, cos(2πt). For the convenience of discussion, let φ � 0 [28]. e criterion for judging the correct separation of two components is as follows: IMF 1 is the first IMF extracted from Sig 3 with n-th siftings and ‖ · ‖ L 2 (T) stands for the Euclidean norm on functions defined over [0, T]. e value of equation (23) equals 0 when the two components are correctly separated; it is close to 1 when the two components are badly separated.
However, this criterion has one obvious drawback: there may be some discrete data points with large errors at the ends of IMF 1 , which makes the value of equation (23) far greater than 1, especially when a is small (for example, in Figure 14(b), the error between IMF and the FM signal is  great only at the ends, so that the value of equation (23) is much greater than 1, but the high-and low-frequency components are well separated). erefore, we change the criterion to the correlation coefficients described in the following equation: If the HF and LF components are separated correctly, the difference between IMF 1 and HF does not contain any information of LF, and the value of ρ is 0. When they are not well separated, the difference between IMF 1 and HF is close to LF and the value of ρ is close to 1. reshold 0.5 is used to determine whether HF and LF are separated.
Based on equation (24) When EMD and IMMD decompose Sig 3 , there is a critical cutoff frequency ratio (see Figure 16). If f exceeds the cutoff frequency ratio, no matter what the amplitude ratio is, the two components cannot be separated. e cutoff frequency ratio from EMD is about 0.65 (0.67 [41] and 0.65 [28]), and that from IMMD is about 0.95. e transition areas (af � 1; af 2 � 1) [41] of the two methods are almost identical. Obviously, in the area where the critical frequency ratio is 0.65∼0.95, IMMD can successfully separate two-tone signals against mode-mixing. erefore, IMMD shows superior orthogonal decomposition capability than EMD.

16
Mathematical Problems in Engineering

ree-Tone Signal.
e cutoff frequency ratio for twotone signals decomposed by IMMD is greater than that by EMD. erefore, we particularly construct the two-tone signal with a frequency ratio of 1 : 0.67, to evaluate the decomposition performance of the IMMD method in detail. Furthermore, as there is a dyadic filter bank of EMD, a signal whose frequency is 0.5 times of frequency of the HF signal is superposed on the signal, which increases the difficulty of separation. e three-tone signal model is as follows: � sin(20πt) + sin(0.67 × 20πt) + sin(0.5 × 20πt).

(25)
Let two methods decompose Sig 4 into Sig 3 IMFs and a residual component r. e Fourier spectra of three IMFs are shown in Figure 17.
In the results of Sig 4 decomposed by EMD (see Figures 17(a) and 17(c)), tone with a frequency of 10 Hz only exists in IMF 1 . Undoubtedly, EMD breaks the tone with a frequency of 6.7 Hz into two parts, one in IMF 1 and the other in IMF 2 . To our surprise, the tone with a frequency of 5 Hz does not exist only in IMF 2 , but is broken in both IMF 1 and IMF 2 . So there is severe mode-mixing between IMF 1 and IMF 2 from EMD. IMF 3 is a pseudocomponent from EMD. On the contrary, although the three sinusoidal signals in the frequency domain have mode-mixing, the amount of mixing is very small. IMMD separates three-tone signal well (see Figures 17(b) and 17(d)).
Finally, Table 7 shows the orthogonal values of Sig 4 decomposed by two methods. It is verified that IMFs from IMMD are orthogonal. Although some IO jk and IO values from EMD are better than those from IMMD, the decomposition of Sig 4 by the EMD method is a failure. erefore, there is no doubt that the orthogonal decomposition capability of IMMD is better than that of EMD about Sig 4 .

Two Commonly Used Realistic Nonstationary and Nonlinear Signals.
In this section, we will conduct experiments on two commonly used real-world signals to compare the decomposition performance of the three common EMD domain methods of EMD, EEMD, and CEEMDAN with the IMMD method, including orthogonal performance and mode-mixing.
We carry out a statistical significance test on IMFs of annual GSTA with the posteriori test method proposed in [42][43][44]. e premise of applying the posteriori test method is that IMF 1 of any well-sampled data is almost always the noise, so the noise contained in the data was estimated based on the IMF 1 [42]. In Figure 19, the solid line is the expectation of variance of IMFs of the white noise, and IMF 1 of white noise contains the same variance as the IMF 1 of GSTA (the expected case); the upper (lower) dotted line is the expectation of variance of IMFs of white noise of three (one-third) times that of the expected case [42]; the upper (lower) dash-dotted line is the expectation of variance of IMFs of white noise of two (one-second) times that of the expected case. e results show that IMMD and EMD each have 4 IMFs and EEMD and CEEMDAN each have 3 IMFs which are useful information of the data because they are beyond the three-time variance of the noise. EEMDʼs IMF 4 and IMMDʼs IMF 5 are also useful information for data, but it is beyond the twice variance of noise. In addition, there is a pseudocomponent IMF 6 (with minimal energy) from EMD.
For IMFs with insignificant noise in the statistical sense, it is indicated that these IMFs contain physically meaningful information, that is, various trends of GSTA [42] (see Figure 20). EMD obtains the overall trend and the change trend about 60 years of GSTA [42]; in addition to the above two trends, EEMD also obtains a trend of about 20 years of GSTA [42]; CEEMDAN obtains a suspected noised 60-year trend of GSTA; in addition to the above three different trends, IMMD also obtains a linear trend of GSTA [42]. Because the decomposition of GSTA by IMMD is adaptive and orthogonal, the number of mode components obtained is more, so it can effectively alleviate mode-mixing (remarkably, the mixing of IMF 9 and r of IMMD is approximately the r component of the other three methods), making GSTA have more physical meaning. Tables 8 and 9 give the IO jk values of the first 7 IMFs of GSTA, which verifies that GSTA is orthogonally decomposed by IMMD and posteriori orthogonally decomposed by the other three methods. IO j6 and IO 67 from EMD are very small, which are caused by the pseudocomponent IMF 6 , and have no practical significance. e number of IMFs from IMMD is more than that from the other three methods, resulting in more leakage, so most of the IO jk values from IMMD are greater than that from the other three methods. Table 10 shows the IO value of GSTA. Among them, EEMD has the best IO value, followed by EMD, IMMD, and CEEMDAN. Except for IMMD, the order of the magnitude of the IO value corresponds to the ability of the other three methods to decompose GSTA (the number of trends). We believe that if we evaluate the decomposition orthogonality of the signal, then the decomposition result with better practical physical meaning should be the most important. erefore, IMMD has better orthogonal decomposition capability than the other three methods.
In addition, Table 10 also shows the time-consuming of GSTA decomposed by four methods related procedures; EMD has the least time-consuming, followed by IMMD, EEMD, and CEEMDAN.

Electrocardiogram (ECG)
Signal. ECG signal is a typical nonlinear and nonstationary weak signal. Muscle artifact (MA) and baseline wander (BW) are the main noises of ECG signal. Figure 21 shows the original ECG100 (numbered as 100 in the database), MA, BW waveforms [52], and the noised ECG100 waveform formed by superimposing these three signals, and Figure 22 shows all IMFs of the noised ECG100 decomposed by EMD, EEMD, CEEMDAN, and IMMD. In ECG reconstruction, IMF 1 contains a lot of high-frequency noise, which is generally recognized as MA, so it is often removed [53][54][55][56]. For the remaining IMFs, the QRS characteristic wave is used to identify the mode components of ECG. IMF 2 ∼IMF 5 from EMD, IMF 2 ∼IMF 5 from EEMD, IMF 2 ∼IMF 7 from CEEMDAN, and IMF 2 ∼IMF 11 from IMMD contain obvious QRS waves, which are identified as mode components of ECG100 and used to reconstruct ECG100. Figure 23 shows the waveform and spectrum of the ECG100 reconstructed by four methods. In the frequency domain, the four methods have removed BW well, and for low-frequency BW with f < 1 Hz, IMMD has the best elimination effect; the best method to eliminate MA is EMD, and the worst is CEEMDAN, which may be caused by the largest amount of auxiliary noise in CEEMDAN. In the time domain, except for the IMMD, the first, second, third, and

f(Hz)
Amplitude IMF 1 from EMD IMF 2 from EMD IMF 3 from EMD Amplitude IMF 1 from IMMD IMF 2 from IMMD IMF 3 from IMMD  sixth characteristic T waves of ECG100 reconstructed by the other three methods have all moved forward, so pseudo-T waves are generated, which may cause doctors to misdiagnose.
In order to quantitatively evaluate the performance of the four methods, the correlation coefficient R is used to quantitatively describe the reconstruction accuracy, the signal-to-noise ratio (SNR) and the mean square error   Sixty-year or longer timescale trend (IMF 5 + IMF 7 ) Twenty-year or longer timescale trend (IMF 4 + IMF 5 + IMF 7 ) 1860 1880 1900 1920 1940 1960 1980 2000 (b) GSTA Overall trend (IMF 7 ) Year Sixty-year or longer timescale trend (IMF 6 + IMF 7 ) Sixty-year or longer timescale trend with suspected noise (IMF 2 + IMF 6  Year GSTA Linear trend (IMF 10 ) Overall trend (IMF 9 + IMF 10 ) Sixty-year or longer timescale trend (IMF 7 + IMF 8 + IMF 9 + IMF 10 ) Twenty-year or longer timescale trend (IMF 5 + IMF 7 + IMF 8 + IMF 9 + IMF 10    Mathematical Problems in Engineering Table  9: IO jk of GSTA decomposed by CEEMDAN and IMMD.      , h). In the green dashed box, pseudo-T waves are obviously caused by the forward movement of the T wave from reconstructed ECG100 by EMD, EEMD, and CEEMDAN.       (MSE) are used to quantitatively describe the ability of denoising, and the corresponding formulas are as follows:

Mathematical Problems in Engineering
where x [n] is ECG100 and y [n] is ECG100 reconstructed by methods. e three characteristic quantities of each of the four methods are shown in Table 11. e ECG100 reconstructed by IMMD has the largest ρ and SNR values and the smallest MSE value. erefore, IMMD causes the smallest ECG100 distortion, its denoising performance is the best, and it shows the best reconstruction performance, followed by CEEMDAN, EEMD, and EMD. IMMD has better orthogonal decomposition capability than the other three methods, noise can be well separated from ECG100, and mode-mixing can be alleviated. In addition, the related procedure of noised ECG100 decomposed by EMD has the least time-consuming, followed by IMMD, EEMD, and CEEMDAN. Tables 12-15 show the IO jk values of the first 11 IMF components of noised ECG100, which verifies that the IMMD method decomposing complex ECG is orthogonal and that the other three methods decomposing ECG have posterior orthogonality. In addition, most of the IO jk values of IMMD are larger than those of the other three methods.
is is because the number of IMFs from IMMD is nearly 1.7 times more than that from the other three methods, which causes more leakage. Table 11 shows the IO value of ECG100 decomposed by the four methods. IMMD has the best IO value, followed by EMD, CEEMDAN, and EEMD.

Discussions and Conclusions
6.1. Discussions. In theory, the signal is orthogonally decomposed by IMMD. e mean value by the envelope used in EMD and its existing improved methods is difficult to be expressed mathematically, so the decomposition orthogonality cannot be proved. IMMD overcomes this problem by using local integral to calculate the mean. eoretical proof and experiments involving signals in manuscripts verify that the decomposition of the signal by IMMD is orthogonal, but this needs to be verified by more signals, especially real signals.
Sifting affects the orthogonal decomposition of the signal by IMMD method. With the increase of the sifting number, (1) the smoothing of the mean curve leads to the decrease of orthogonality and (2) IMF symmetry leads to the enhancement of orthogonality. For a certain two mode components, the orthogonality may decrease with the increase of sifting number, or even to the extent that the two mode components are no longer considered to be orthogonal, but it is precisely this way that the IMMD (or EMD) method can separate two mode components that are not strictly orthogonal (for example, signal Sig 2 ). erefore, the role of sifting needs to be studied more deeply.
Because the signal is orthogonally decomposed by IMMD, mode-mixing can be alleviated. If the signal is orthogonally decomposed, the mode components of the signal are orthogonal, and there is no mode-mixing between them, so IMMD can alleviate mode-mixing.
e Monte Carlo method of Gaussian white noise and the signal experiments given in this study show that IMMD effectively alleviates mode-mixing and has better multiresolution analysis characteristics. In particular, the characteristic of IMMD equivalent constant-Q filter bank is trisection, which is different from EMD, EEMD, etc.
PMF from IMMD has good robustness to the sifting, but it also needs to be verified by more signals, especially real signals. Due to the envelope mean method, the PMF from EMD, EEMD, and CEEMDAN is weakly convergent to the sifting, and the PMF will continue to change with the increase of sifting number and eventually become uniform in amplitude (the problem of oversifting), and IMF loses the substantial significance of amplitude modulation [57]. In addition, due to the influence of auxiliary noise, any two decomposition experiments of the same signal by EEMD or CEEMDAN have small differences in the two decomposition results [20]. If the signal needs to be decomposed very accurately, this small difference is very likely to cause wrong conclusions.  In short, compared with the EMD, EEMD, and CEEMDAN methods, the decomposition of signal by the IMMD method is orthogonal, and the number of IMFs from IMMD is more, so it has better multiresolution analysis capabilities, and the decomposition has better physical meaning; IMF from IMMD is robust to the sifting and will not cause oversifting problem; in addition, IMMD can effectively alleviate mode-mixing; compared to EEMD and CEEMDAN, IMMD does not require any parameter settings, and the decomposition is data-driven and adaptable; due to the large number of IMFs produced by IMMD, the time-consuming of IMMD-related procedures is greater than that of EMD, but it is much less than that of EEMD and CEEMDAN procedures, which is caused by the hundreds of auxiliary noises of them.

6.2.
Conclusions. An improved method of EMD, IMMD, is proposed in this study. Compared with EMD, EEMD, and CEEMDAN, IMMD has the following advantages: smoothing of the mean curve is ignored, IMMD method decomposes the signal orthogonally; without any input parameters, IMMD can effectively alleviate mode-mixing; and IMMD method has good robustness to sifting.
We suggest that IMMD is suitable for nonlinear and nonstationary signal research fields that require multiresolution and high stability, such as geophysics and bioelectric signal processing.

Data Availability
e data used to support the findings of this study are publicly available on the Internet, and the detailed website address or information is given in the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.