Abstract
Chaotic data analysis is important in many areas of science and engineering. However, the chaotic signals are inevitably contaminated by complicated noise in the collection process which greatly interferes with the analysis of chaos identification. The chaotic vibration is extremely nonlinear and has a broad range of frequencies; linear filtering methods are not effective for chaotic signal noise reduction. Then an improved ensemble empirical mode decomposition (EEMD) based on singular value decomposition (SVD) and Savitzky-Golay (SG) filtering method was proposed. Firstly, the noise energy of first level intrinsic mode function (IMF) was estimated by “” criterion, and then SVD was used to extract the signal details from first IMF, and the singular value was selected to reconstruct the IMF according to noise energy of the first IMF. Secondly, the remaining IMFs are divided into high frequency and low frequency components based on consecutive mean square error (CMSE), and the useful signals of high frequency components and low frequency components are extracted based on SVD and SG filtering method, respectively. The superiority of the proposed method is demonstrated with simulated signal, two-degree-of-freedom chaotic vibration signals, and the experimental signals based on double potential well theory.
1. Introduction
Chaotic data analysis is very important in many areas of science and engineering such as parameter estimation, system identification, computation of correlation dimension, and prediction [1]. However, the chaotic signals are inevitably contaminated by complicated noise in the collection process which greatly interferes with the analysis of chaos identification and prediction. Due to the complexity, nonstationarity, and nonlinearity of vibration signals, the linear low-pass filtering is not suitable for the noise reduction. In fact, linear filtering based on frequency separations often causes distortion because some of the suppressed frequency components are part of the dynamics, while chaos-based approaches often are computationally expensive and may not be very effective either, especially when noise is very serious [2–4]. Therefore, the effective denoising method plays a decisive role in the analysis of chaotic signals [5].
Lots of research works on denoising method of nonlinear and nonstationary chaotic signals have been done in different fields [3, 6, 7]. Davies applied the gradient descent algorithm to reduce noise in a time series from nonlinear dynamical system [8]. Holzfuss and Kadtke researched the noise reduction method with the radial basis functions. The interpolation of the global flow of the nonlinear dynamical systems is conducted with the global radial basis function, which provides a considerable increase of the signal-to-noise ratio [9]. The shadowing method is applied by Doyne Farmer and Sidorowich to find the deterministic orbit as close as possible to a given noisy orbit, and the optimal solution in the sense of least-mean-squares is presented. It provides an effective and convenient numerical method for noise reduction for data generated by a nonlinear dynamical system [10]. In order to increase the signal noise ratio of the measured signals, Lee and Williams presented the preprocessing approach to improve noise reduction for chaotic signals [11]. Han and Liu proposed an effective noise reduction method to remove Gaussian noise embedded in chaotic signal. The noise residual ration based on dual-wavelet is applied to choose wavelet decomposition scales and the singular spectrum analysis and spatial correlation theory is used to estimate wavelet coefficients [12]. Kang and Harlim presented the autoregressive linear stochastic models to filter the nonlinear spatiotemporal chaos noise [13]. Shin et al. studied the iterative SVD method for noise reduction for low dimensional chaotic time series. The technique is particularly concerned with the improvement of the phase space reconstruction diagram. The method can be used blindly in the case of very noisy signal [14].
In recent years empirical mode decomposition (EMD) has been widely applied to noise reduction which is an adaptive time-frequency signal processing method [15]. The original signal is decomposed into a series of intrinsic mode functions (IMFs) according to the signal characteristics with EMD, and EMD can well reveal the nonstationary and nonlinearity information of the signal. EMD has been applied to denoising for vibration signal successfully [16–19]. However, EMD cannot extract true information of vibration signal accurately for the problem of mode mixing. The phenomenon that the different or similar scales appear in the same mode is regarded as mode mixing. As the frequency of noise and useful signal is too close to cause energy mutual penetration, then the noise and useful signal cannot be separated accurately based on EMD [20]. To alleviate mode mixing, Huang et al. [21] developed EEMD to improve EMD. By adding white noise to the original signal and calculating the means of IMFs repeatedly, EEMD is more accurate and effective for noise reduction of vibration signal.
However, in the traditional EEMD noise reduction methods, the selection method of IMF components is not mentioned in the majority of literatures. An automatic algorithm of choosing IMF components to reconstruct signal was designed according to the characteristic when a white noise signal is decomposed by EMD; the product of the energy density of IMF component and its corresponding averaged period is a constant in [22]. There are two shortcomings in this method. The useful information in the high frequency IMF components is eliminated but the noise in the low frequency IMF components is remained. Then the conventional EEMD noise reduction method cannot achieve the desired denoising target effectively. In order to overcome these shortcomings, an improved EEMD denoising method was proposed in this paper. The main denoising idea of the proposed method is to suppress the noise and extract the useful information for each IMF component, which can remove noise and retain useful signal as much as possible at the same time. The implementation processes of the improved EEMD denoising method are as follows. In Section 3.1, the combination of “3σ” criterion and SVD was used for extracting the useful signal from the first IMF. Then, SVD refines the high frequency IMFs in Section 3.2, and the Savitzky-Golay (SG) is used for making the low frequency IMFs smoother in Section 3.3. Finally, the process of the proposed method is discussed in detail with a simulated signal analysis and two-degree-of-freedom chaotic vibration signals. In order to generate reliable and repetitive chaotic vibration to verify the theoretical results about extraction of the chaos characteristics and evaluate the performance of the proposed method, the experiment of leaf spring based on the double potential well theory was carried out in this paper, and the results show that the contaminated noise can be filtered effectively.
2. EEMD Principle
The empirical mode decomposition as a new signal processing method is used for analyzing nonstationary and nonlinear signals which are decomposed into several stationary data and the trend superposition based on EMD, and the IMFs are the stationary data and represent the different scale characteristic of the original signal. Compared with the Fourier transform and wavelet transform, the EMD decomposes any given data into intrinsic mode functions that are not set analytically and are instead determined by an analyzed sequence alone. The basis functions are in this case derived adaptively and directly from input data. For an IMF resulting from the EMD, these requirements should be met; that is, the number of IMF extrema (the sum of the maxima and minima) and the number of zero-crossings must either be equal or differ at most by one; at any point of an IMF the mean value of the envelope defined by the local maxima and the envelope defined by the local minima should be zero.
Generally, the following steps are included in the process of EMD. The maxima and minima of the sequence are identified, respectively, at first, and the envelopes are obtained. The mean is calculated based on the two envelopes, and the mean value so calculated is further subtracted from the initial sequence, and then the above steps result in the extraction of the required intrinsic mode function in the first approximation. To obtain the final IMF, new maxima and minima shall again be identified and all the above steps repeated. This repeated process is called sifting. The first IMF is obtained when the sifting process is successfully completed, and the next IMF can be obtained by subtracting the previously extracted IMF from the original signal and repeating the above described procedure once again. This continues until all IMFs are extracted.
EEMD is a substantial improvement of EMD. In this method, the added white noise is distributed in the whole time-frequency space uniformly, which facilitates a separation of the frequency scales and reduces the occurrence of mode mixing. The procedures are presented as follows [23].(1)Initialize the number of ensemble .(2)Given the amplitude of the added white noise, .(3)Given is an original signal, add a random white noise signal to , where denotes the th added white noise series and represents the th noise-added signal, ~.(4)Decompose the noise-added signal into IMFs by using the EMD, where denotes the th IMF of the th noise-added signal and is the number of IMFs.(a)Localize the maximal points of the signal and connect them with a cubic spline. The upper envelope of the signal is represented by .(b)Localize the minimal points of the signal and connect them with a cubic spline. The lower envelope of the signal is represented by .(c)Subtract the mean envelop from the original signal to obtain a pro-IMF , and is denoted as .(d)Regard as the new signal when meets these two properties: the mean of the upper and lower envelope is equal to zero, and the number of extrema and the number of zero-crossings are equal or different at most by one.(e)Regard substeps (a)–(d) of step (4) until the resulting signal is a proper IMF .(f)Subtract IMF from the signal . The residual is considered as new data and is fed back to step (1): (g)Complete the algorithm when the residual of substep (f) is a monotonic function. The last residual can be considered as the trend: If then go to step (3) with .(5)Calculate the ensemble mean of the trials for each IMF: (6)Regard the mean of each of the IMFs as the final IMFs.
Evidently, the added white noise series cancel each other in the final mean of the corresponding IMFs and the mean IMFs stay within the natural dyadic filter windows and thus significantly reduce the chance of mode mixing and preserve the dyadic property.
3. Improved EEMD Method
3.1. “3σ” Criterion and Singular Value Decomposition
In the process of signal denoising by conventional EEMD, all of the first IMF (IMF1) is usually considered as noise. However, through the in-depth study, it is found that there is a certain amount of useful signals in IMF1 [24], and the refine of IMF1 is still a difficult problem because of lack of a prior knowledge. Therefore, extracting the useful signal components of IMF1 is very important for reducing the signal distortion. Consider that the noise interfered in the most part of IMF1 approximately obeys the normal distribution with zero mean value; the noise energy of first level intrinsic mode function (IMF) was estimated by using the “3σ” criterion. Then singular value decomposition (SVD) was used to extract the signal details from first IMF, and the singular value was selected to reconstruct the IMF according to noise energy of the IMF1 in this paper.
According to “3σ” criterion, the distribution of noise meets , where is the mean value of . Suppose ; that means the probability of is 0.9973. Therefore we can think the signals falling outside must contain the useful information to be retained. The detailed expression is as follows:where expresses the extracted signals from , the noise variance is , and HH is the high frequency coefficient of IMF1.
According to the theorem of the SVD, for a real matrix with the size of , it always satisfies where the matrix is a diagonal matrix. All its diagonal entries are nonnegative and sorted in descending order; that is, . These elements are called the singular values of the matrix . The key problem is how to judge which of diagonal elements in the matrix should be regarded as those with small values that are replaced by zero.
Since the original signal is decomposed by EEMD and the trend component is extracted from the original signal, SVD can be used to remove the detrended signal, and the difference spectrum of singular value is used to select the singular values adaptively [25].
The difference spectrum of singular value is defined aswhere .
The difference spectrum series is . Suppose the first singular values correspond to the useful signal, and then the difference spectrum of th element is the maximum peakwhere is the peak value of . The useful signal can be achieved by reconstructing the first singular values.
“3σ” criterion has been used to extract the useful signal details of first IMF directly in [24]. As can be seen from formula (5), the extracted signal is a high amplitude impulse signal. As an example, the details extracted from IMF1 of 5 dB noisy Lorenz time series after EEMD by using “3σ” criterion are shown in Figure 1. Obviously, it is untrue to take the extracted details as useful signal of IMF1.
In order to obtain accurate details from IMF, the energy of extracted signal based on “3σ” criterion is estimated as a priori knowledge. Then SVD is used to extract the signal details from IMF1, and the singular value was selected to reconstruct the IMF1 according to noise energy in the IMF1. The steps are given as follows.(1)The details are extracted based on formula (5).(2)The energy of extracted details is estimated based on the following formula: where is the number of IMFs which meet “3σ” criterion.(3)Suppose the singular value of IMF1 after SVD is , and the reconstructed signal of the first singular value is . Then the value of can be calculated according to
3.2. The Extraction of High Frequency IMFs by Using SVD
After processing signals based on EEMD, the IMF1 whose spectrum is half of the original signal is equivalent of a high-pass filtering, and the noise energy density of each IMF is decreasing 2 times of original one; that is, the noise strength is substantial reduced with the increasing of the order of IMF [26]. Therefore, the singular values of IMF1 are evenly distributed. For the remaining IMFs, the singular values of signals are significantly greater than the singular values of noise. According to this feature, the singular values of useful signal can be selected based on difference spectrum method and then they are used to reconstruct the original signal after denoising.
For the chaotic vibration signals, the noise is mainly distributed in the high frequency IMFs. On the contrary, the useful signals are mainly distributed in the low frequency IMFs. Therefore, the IMFs dominated by noise and useful signals can be distinguished by using consecutive mean square error (CMSE) [26]. The main principle is to find an index value and make the reconstruction error of IMFs when is the minimumwhere is the length of signal. The index value is given as follows:
After high and low frequency IMFs being distinguished, the high frequency IMFs can be processed by SVD.
3.3. The Smooth of Low Frequency IMFs by SG
The IMFs of low frequency still have few noise components through EEMD processing. If the signals can be smoothed further, the signal-to-noise ratio (SNR) will be improved effectively. Reference [27] indicates that the Savitzky-Golay method has a great denoising function for the low frequency IMFs. Therefore, the useful signals are extracted from low frequency IMFs by SG method.
Suppose is one of the IMF coefficients, and a polynomial with th order is fitted by points in the vicinity of at the sense of least square. The smooth value shown in the following is ’s value in :where is the number of the left points of and is the number of the right points of . is the coefficient of . Suppose the measurement data is ; in order to fit the testing data with , we need to optimize
The flow chart of improved EEMD is given in Figure 2.
4. Simulation Analysis of Proposed Method
4.1. Denoising for Noisy Lorenz
In order to test the denoising performance of improved EEMD, the denoising for the Lorenz signal with noise is taken as an example and the process of improved EEMD method is described in detail.
The SNR of Lorenz time series is set to 5 dB. The Lorenz time history is shown in Figure 3. The EEMD results of the noisy Lorenz signal are shown in Figure 4.
(a)
(b)
The EEMD results comprise 7 IMFs (IMF1~7) and a residue component. The index value based on (11) and (12); that is, IMF1~3 are dominated by noise and IMF4~7 are dominated by signal. Therefore, the improved EEMD process of noisy Lorenz time series can be described as follows.
(1) The useful signal details of IMF1 are extracted based on “3σ” criterion and SVD. The singular values of IMF1 are shown in Figure 5, where represents the number of singular values and indicates the amplitude of the singular value. We can see that the singular values are evenly distributed, and the number of singular values cannot be determined by using spectrum difference effectively. Therefore, the details information of IMF1 is extracted based on “3σ” criterion; then the noise energy which provides a priori knowledge for selecting the singular values can be estimated.
The extracted useful signal from 5 dB IMF1 based on “3σ” criterion and SVD is shown in Figure 6. Compared with Figure 1, the extracted signal is more uniform and reasonable.
(2) The useful details of high frequency are extracted by SVD. The singular values of IMF2 and IMF3 are shown in Figure 7; it is shown that the first several singular values are significantly larger than the other ones and the trend component is extracted from the original signal as the EEMD. Then the singular values of useful details are selected based on spectrum difference and they are used to reconstruct , which are shown in Figure 8.
(a)
(b)
(a)
(b)
(3) The low frequency IMFs are smoothed based on Savitzky-Golay filter algorithm. The fitting polynomial order is 3, and data window is 15. Then ~ are obtained.
(4) The denoised result is by reconstructing signals.
In order to compare the proposed method with the traditional noise reduction, dB Lorenz time series are denoised by using wavelet weighted threshold method [28], the traditional EEMD method, and the proposed method, respectively. In the process of wavelet weighted threshold method, “sym8” is chosen as wavelet basis function and the wavelet decomposition level is taken as 5. In the process of traditional EEMD method, the hard threshold method is used in denoising. The denoising results for noisy Lorenz time series with three methods are shown in Figure 9. Since the model is known for Lorenz time series, SNR and MSE (mean square error) are used to evaluated the denoising effectiveness. It is obvious that the greater SNR and the smaller MSE, the better the denoising effectiveness. The denoised results for Lorenz time series with dB, 5 dB, 10 dB, and 15 dB are shown in Table 1.
(a)
(b)
(c)
As seen in Figure 9, there are a lot of jitters in (a) and (b), which means the nonsmooth, while Figure 9(c) is smooth and in a good agreement with the original signal. As shown in Table 1, the proposed method is superior to the wavelet weighted threshold method and the traditional EEMD method, and the advantage is particularly prominent for low SNR Lorenz time series.
4.2. Denoising for Noisy Two-Degree-of-Freedom Chaotic Vibration Signal
In order to validate the effectiveness of the proposed method in chaotic vibration analysis, the two-degree-of-freedom chaotic signal is denoised as an essential example for further confirmations. The power system physical model of two-degree-of-freedom chaotic vibration is shown in Figure 10, where is the mass of vibration isolated equipment, is mass of basis, represents the displacement of the isolated equipment, and represents the basis displacement. The elastic supporting of basis is assumed to be linear and stiffness and damping coefficient are denoted as , , respectively. The vibration isolation elements have characteristics of cubic stiffness and linear damping. The linear and nonlinear stiffness coefficients are and , respectively, and is damping coefficient. According to Newton’s law, the differential equation of motion is derived:
Let , , , , and then
Therefore (15) can be rewritten as dimensionless equation form,where
When parameter , , the system works in a chaotic state [29]. The phase diagram of the system is shown in Figure 10. The white Gaussian noise is superposed on the system. Then the noisy signal is denoised by wavelet weighted threshold method, the traditional EEMD method, and the proposed method, respectively. The denoised results are shown in Figures 11(b), 11(c), and 11(d). As shown in Figure 11, the signal curve is coarse in Figures 11(b) and 11(c), and the signal curve in Figure 11(d) is more close to the original one in Figure 11(a). The SNR and MSE of denoised signal based on the three methods above are listed in Table 2. It is shown that the SNR of denoised signal based on the proposed method is the largest, and the MSE is least. It can be concluded that a wonderful noise reduction effect is obtained based on the proposed method.
(a)
(b)
(c)
(d)
4.3. Denoising for Noisy Two-Well Potential Magnetic Leaf Spring Chaotic Vibration Signal
In order to verify the effectiveness of the proposed method in practical chaotic vibration analysis, experiment of leaf spring based on the double potential well theory was carried out in this paper.
The cantilever structure of mechanical chaotic platform based on the double potential well theory is designed in Figure 12, which includes eight components. (1) Exciter is applied to provide exciting forces. (2) Supporter is used to keep the exciter in horizontal state. (3) Leaf spring is designed to generate chaotic signal. (4) Rectangular magnet is chosen to produce the double potential well. (5) Mass block is installed at the end of the leaf spring and it will be attracted by the magnets which acts as two potential wells. (6) External frame is given to install magnets and connect with the exciter. (7) Clamp is applied to fix the leaf spring tensely. (8) Laser displacement sensor model is CD33-30NV. The maximal displacement is ±4 mm (the maximum permissible error is ±0.1% FS). The distance between the sensor and the leaf spring is 30 mm.
(a)
(b)
The frequency sweeping experiment is conducted in the range of 5 Hz~27 Hz containing the natural frequency. Sampling rate is set as 2 kHz; the measuring time length is 5 seconds for each measurement. The procedures of the experiment are as follows: (1) LabVIEW is functioned as signal acquisition, and the time history and power spectrum are monitored online; (2) the sinusoid signal from the generator is amplified, and then the exciter is driven; (3) Design a recording rule of the measuring data, such as F××L××S×××.lvm, where is the gain of the amplifier, is the distance between the magnet and the leaf spring central line, and is the height of spring; (4) the frequency of the excitation is changed in 0.1 Hz steps; (5) the nonlinear signal analysis is carried out, such as visual inspection; extraction of the frequency of the measured signal; calculating the reconstruction of the phase space and reconstructing attractor; and computation of the characteristic exponent; (6) chaos identification software based on the nonlinear time series analysis is programmed, and the improved Poincaré section and LE exponent are obtained; (7) determine the chaotic parameters range of the system.
The system enters into stable chaotic state automatically when the gain is 0.98. As shown in Figure 13(a), the time history is random-like motion and the power spectrum is broadband one. The reconstructed phase space parameters are calculated with false nearest neighbors (FNN) and first minimal mutual information (FMMI) methods. The reconstructed embedding dimension and delay time are 4 and 16, respectively. The chaotic attractor is obtained in Figure 13(e), and the maximal Lyapunov exponent (LE) is 0.313. The measured signal is polluted by ambient noise apparently.
(a)
(b)
(c)
(d)
(e)
The measured signal is denoised by using wavelet weighted threshold method, the traditional EEMD method, and the proposed method, respectively. The denoised results are shown in Figures 14(a), 14(b), and 14(c). By comparison, it is obvious that the denoised phase diagram Figure 14(c) based on the proposed is smoother than the others and exhibits the geometry of the original signal attractor more clearly. Since the original signal derives from a mechanical vibration signal and the signal and the noise are unknown, the SNR and MSE cannot be used to quantitatively evaluate the effectiveness of denoised for unknown system. However, the autocorrelation coefficient value of chaotic signal is much greater than that of noise signal [30]. Therefore, the autocorrelation coefficient value of denoised signal can be applied to verify the effectiveness of the proposed method in practical chaotic vibration analysis.
(a)
(b)
(c)
The autocorrelation coefficient is defined aswhere is denoised signal, is time delay, and is the length of signal.
The autocorrelation coefficient values of noisy and denoised signal based on the three methods above are shown in Table 3. As can be seen that the autocorrelation coefficient values of noisy signal are much smaller than that of denoised signal, and the autocorrelation coefficient values of the denoised signal based on the proposed method are the greatest, which further demonstrates its superiority for denoising.
The parameter regime of the system is obtained through the experiment, as shown in Figure 15, where the horizontal axis is excitation frequency and the vertical axis is the gain. The reconstruction parameters and characteristic exponents of measured signals are shown in Table 4. As all we know, if the maximum LE of the system is positive and the correlation dimension is fractal, the system is a chaotic motion. When Hz and 15 Hz, the correlation dimension noisy signal is close to integer; it illustrates that the unexpected noise seriously degrades the calculation of the characteristic exponents for measured signals and the proposed method is beneficial to the further analysis of chaotic vibration signals.
5. Conclusions
In order to develop a more reasonable and practical method for chaotic signal denoising, an improved EEMD method is proposed in this paper. The detail components are extracted from the first IMF based on “3σ” criterion and SVD. The signals of high and low frequency IMF are extracted based on SVD and SG method, respectively, and thus it can be used not only to reduce the loss of useful signal in high frequency IMF but also to remove part noise of low frequency IMF effectively. The effectiveness of the proposed method is verified by a series of simulated and practical experiments. From the calculation results and discussions, the following conclusions can be drawn:(1)The useful details extracted from the first IMF based on combination of “3σ” criterion and SVD is more reasonable than that based on “3σ” criterion.(2)The original signal is decomposed by EEMD and the trend component is extracted; then the SVD can be used to denoise the high frequency IMFs effectively and the difference spectrum can be used to select the singular values adaptively.(3)The proposed method is effective in many cases no matter whether the system equation of motion is known or not. It can extract the useful signal from noisy signal and even the SNR of the signal is very low as well.(4)The chaotic experiment setup based on the double potential well theory is designed, which is applied to produce repetitive, stable chaos, and the chaotic parameter regime of the system is determined, in which the vibration signal is denoised by using the proposed method.(5)The calculation of characteristic exponents of measured signals is degraded by the noise and the proposed method provides an important preprocessing way for further analysis of chaotic vibration signals.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work is supported by National Nature Science of China (nos. 51179197, 51579242, and 51509253).