# Chaotic Signal Denoising Based on Adaptive Smoothing Multiscale Morphological Filtering

Authors: Guiji Tang; Xiaoli Yan; Xiaolong Wang
Journal: Complexity
 (2020)

---

## Abstract

Nonlinear time series denoising is the prerequisite for extracting effective information from observation sequence. An effective chaotic signal denoising method not only has a good signal-to-noise ratio (SNR) enhancement performance, but also can remain as a good unpredictable denoised signal. However, the inherent characteristics of chaos, such as extreme sensitivity to initial values and broadband spectrum, pose challenges for noise reduction of polluted chaotic signals. To address these issues, an adaptive smoothing multiscale morphological filtering (ASMMF) is proposed to reconstruct chaotic signals. In the process of noise reduction for contaminated chaotic signals, firstly, a multiscale morphological filter is constructed adaptively according to the multiscale permutation entropy (MPE) of morphological filter residuals, and the contaminated signals are filtered. Secondly, the weight coefficients of scale structural element are calculated by the residual sum of squares operation, and the chaotic signals are reconstructed. Finally, the resampled filter signals are smoothed by cubic B-spline interpolation operation. In the experiment, the Lorenz signal with white Gaussian noise, the measured sunspot, and the chaotic vibration signal are reconstructed by four comparison methods. The test results show that the proposed ASMMF method has obvious advantages in noise suppression and topological trajectory restoration.


---

## Body

## 1. Introduction

Chaos is a seemingly random irregular motion that occurs in a deterministic system [1]. It is vital to preprocess observation data to effectively extract desired chaotic information. In the process of analyzing the existing chaotic phenomena, the measured chaotic behavior is susceptible to noise interference. Therefore, it is desirable to preprocess the noise of observation data without distorting the dynamic characteristics of underlying signals [2, 3]. Owing to chaotic signals usually having a broadband spectrum that overlaps with the spectrum of noise [4], the conventional linear low-pass filtering method is ineffective for chaotic signal denoising. If these methods are applied in noise reduction of chaos signals, the clean chaotic signals may be oversmoothed severely [5]. Therefore, it is of great significance to study the corresponding noise reduction methods according to the characteristics of chaotic signals.Over the past decades, a number of approaches for chaotic signal noise reduction have been proposed [3, 6–13]. The singular spectrum analysis method [6] and local projection method [7, 8] are denoising algorithms based on phase space reconstruction; due to their different dynamic characteristics after reconstructing, the chaotic signals are separated from noise. But it is difficult to accurately determine the noise boundary points of noise reduction by singular spectrum analysis, and the parameter selection method of local projection has limitations. Therefore, the noise reduction ability will decrease if the noises of denoising chaotic signals are severe. The threshold denoising method based on time-frequency analysis, such as the wavelet threshold (WT) denoising method [3], and empirical mode decomposition (EMD) [9, 10] are noise suppression algorithms which decompose the noisy signal in time-frequency domain and separate the signal from the noise by threshold processing. However, the denoising performance of the WT method is greatly impacted by the wavelet basis and the number of layers, which reduces the adaptability of the method. Although the EMD threshold denoising overcomes the issue that WT must select the appropriate wavelet basis according to the signal characteristics, it is still hard to determine the threshold properly. The local curve fitting method [11], which smooths the noisy signal by segments, uses the least square polynomial fitting to realize the signal denoising, and it is arduous to reconstruct precisely by local linear approximation due to the high nonlinearity of chaos. Cooperative filtering denoising [12] uses the self-similarity of chaotic signals in time domain to group similar blocks, and the noise suppression effect of the reconstructed signal after filtering is obvious. Nevertheless, for the reason that the selection of filtering parameters of this method is not adaptive, the application of this method is still limited. Wang et al. [13] further improved the cooperative filtering denoising method; the complexity of chaotic signals is analyzed and the optimal filtering parameters are adaptively selected according to the smallest permutation entropy (PE) [14]. Most of the existing noise suppression methods are proposed on the basis of low-pass filter [15, 16]. Analyzing the above research results, the study trend of chaotic signal noise suppression is developing from global noise reduction to local noise reduction, and the noise reduction methods proposed for dynamic characteristics of chaotic signal have better effect.Morphological filtering (MF) is a nonlinear filtering technique [17]. It can be employed to make translation matching and local correction of the original signal from front to back by constructing a certain structuring element (SE), and the morphological characteristics of the original signal are preserved while the noises are suppressed [18]. Currently, morphological filtering has been widely used in image processing [19, 20], signal analysis [21], feature extraction [22], and other fields. In this paper, the MF is first introduced to reduce the noise of contaminated chaotic signal. Combined with the self-similar characteristics of chaotic signals, an adaptive smooth multiscale morphological filtering (ASMMF) method for chaotic signal denoising is proposed on the basis of existing algorithms. Within this scheme, the scale range of structural element is determined in accordance with the self-similarity of temporal waveform of chaotic signal. The multiscale permutation entropy (MPE) [23] of the residual signal is calculated to select the optimal combination of morphological filtering scales, and the residual variance is calculated to determine the weight coefficients of different scale filtering results. The chaotic signal is reconstructed by weighted optimal scale combination of morphological filtering results. To solve the problem of topping distortion in morphological filtering, the reconstructed signal is resampled and smoothed by cubic B-spline interpolation. In order to distinguish ASMMF from the adaptive multiscale morphological filtering without smoothing, the latter is called multiscale morphological filtering for short. The ASMMF is an improved algorithm on the basis of AMMF.The compendium of this paper is organized as follows. In Section2, the process of ASMMF is described in detail. Section 3 analyses the parameters of the algorithm and proposes the corresponding optimization scheme. Moreover, a noisy Lorenz signal is taken to describe the optimizing process. Section 4 applies the proposed method and three comparison methods to denoise different noisy chaotic signals. The ultimate conclusions are drawn in Section 5.


## 2. Adaptive Smoothing Multiscale Morphological Filtering

In order to denoise the contaminated chaotic signals, a novel method called ASMMF is presented in this paper. And this algorithm is divided into three steps: multiscale morphological filtering, signal reconstruction, and signal smoothing. Figure1 shows the noise reduction process.Figure 1
Flowchart of the denoising algorithm for chaotic signals based on ASMMF.


### 2.1. Multiscale Morphological Filtering

Setfn as a one-dimensional (1D) signal, which is a discrete function over the domain F=0,…,N−1, set gm as the SE, which is a function over the domain G=0,…,M−1, and the stipulated value of N should be larger than that of M. The dilation operation and the erosion operation are two kinds of basic operations of mathematical morphology (MM) [24], which can be denoted as symbols ⊕ and Θ, respectively, and the signal process can be represented as follows [17]:(1)f⊕gn=maxfn−m+gmm=0,1,…,N−1,fΘgn=minfn+m−gmm=0,1,…,N−1.A number of basic operators are used in MF. The opening operation and closing operation of a setfn by SE gm are the most widely used; they can be defined as(2)f∘gn=fΘg⊕gn,f•gn=f⊕gΘgn,where symbols ∘ and • stand for the opening operation and closing operation, respectively. The opening operation flattens the positive impulses and matches the negative ones, and the closing operation can flatten the negative impulses and matches the positive ones [20]. On the basis of the sequential operation of opening and closing, the open-closing operator and the close-opening operator are defined as [25, 26](3)FOCfn=f∘g•gn,FCOfn=f•g∘gn.The output signals of the open-closing operator and the close-opening operator have the problem of median bias; thus, the combined morphological filter (CMF) is used in signal processing. The combined morphological filter is able to smooth the signal by eliminating the positive and the negative impacts in the signal simultaneously. The combination morphological filter is further defined as(4)CMFfn=FOCfn+FCOfn2.For complex noisy signals, morphological filtering of single-scale structural element has limited accuracy in reconstructing. It is necessary to construct multiscale morphological filters to denoise target signals at different levels and scales. Letλg be the SE at scale λ, which is obtained by λ−1 times dilation operation; λg is described as [27](5)λg=g⊕g⊕⋯⊕g︸λ−1times.The multiscale CMF can be redefined as(6)yλn=f∘λg•λgn+f•λg∘λgn2.CMF noise reduction results can effectively restore the original signal waveform. Also, the weighted average operation can restrain the statistical bias of denoising results. Indeed, compared with a single-scale CMF, the multiscale CMF is more effective in noise reduction.


### 2.2. Signal Reconstruction

These results of morphological filtering at different scales can be obtained by multiscale morphological analysis of noisy signal. Each filtering result at single scale reflects a specific component contained in the signal. The different components of chaotic signals should appear in multiple scales as a processed signal. The reconstructed signal is obtained by aggregating the filtering results of different scales as follows:(7)Y=∑j=1jωjyj∑j=1jωj=1,where yj is a result of morphological filtering on scale j, ωj is the weight coefficient of yj, and Y is the reconstruction result.The weight coefficient directly affects the effect of multiscale morphological filtering. For the purpose of obtaining a better denoising result, the square operation of residual amplitude difference (RAD) before and after filtering processing is used to determine the weight coefficients adaptively. The RAD can be defined as(8)Rjn=fn−yjn.The weight coefficient can be expressed as(9)ωj=∑nRjn2∑j∑nRjn2.


### 2.3. Signal Smoothing

It is found that the waveform of the aggregated signal is not smooth enough, and the local waveform is slightly rough. Owing to the morphological filter with large-scale SE blurring the detail information of signal and bringing the top-weakening distortion in the process of reducing noise, the signal directly gathered by the results of morphological filtering at different scales is not smooth enough. However, the occurrence of the top-weakening distortion is regular all the way. The distortion only made some local waveforms a little rough, but did not change the general shape of the whole waveform. In order to reduce the top-weakening distortion, a new method of resampling the filtered signal is proposed according to the regularity of the distortion. The cubic B-spline interpolation method is proposed to smooth the filtered signal after resampling.Letdii=0,1,…,n be the de Boor point, where i is the serial number and n is the number of the de Boor point. The k-th B-spline function is defined as follows [28, 29]:(10)Pt=∑i=0ndiBi,kt,t∈0,1,Bi,0t=1,ti≤t≤ti+1,0,others,Bi,kt=t−titi+k−1−tiBi,k−1t+ti+k−tti+k−ti+1Bi+1,k−1t,where t0,t1,…,tn+k+1 are node vectors and Bi,kt is a base function of the node vector. When the denominator in equation (10) is 0, the value is specified as 0. When k is 3, the cubic uniform B-spline function defined in [0, 1] is described as(11)pt=1−t36d0+3t3−6t2+46d1+−3t3+3t2+3t+16d2+t36d3.Equation (11) is a cubic polynomial divided into uniform segments; the signal processed by cubic B-spline interpolation has the characteristics of detail smoothness and linear independence.


## 2.1. Multiscale Morphological Filtering

Setfn as a one-dimensional (1D) signal, which is a discrete function over the domain F=0,…,N−1, set gm as the SE, which is a function over the domain G=0,…,M−1, and the stipulated value of N should be larger than that of M. The dilation operation and the erosion operation are two kinds of basic operations of mathematical morphology (MM) [24], which can be denoted as symbols ⊕ and Θ, respectively, and the signal process can be represented as follows [17]:(1)f⊕gn=maxfn−m+gmm=0,1,…,N−1,fΘgn=minfn+m−gmm=0,1,…,N−1.A number of basic operators are used in MF. The opening operation and closing operation of a setfn by SE gm are the most widely used; they can be defined as(2)f∘gn=fΘg⊕gn,f•gn=f⊕gΘgn,where symbols ∘ and • stand for the opening operation and closing operation, respectively. The opening operation flattens the positive impulses and matches the negative ones, and the closing operation can flatten the negative impulses and matches the positive ones [20]. On the basis of the sequential operation of opening and closing, the open-closing operator and the close-opening operator are defined as [25, 26](3)FOCfn=f∘g•gn,FCOfn=f•g∘gn.The output signals of the open-closing operator and the close-opening operator have the problem of median bias; thus, the combined morphological filter (CMF) is used in signal processing. The combined morphological filter is able to smooth the signal by eliminating the positive and the negative impacts in the signal simultaneously. The combination morphological filter is further defined as(4)CMFfn=FOCfn+FCOfn2.For complex noisy signals, morphological filtering of single-scale structural element has limited accuracy in reconstructing. It is necessary to construct multiscale morphological filters to denoise target signals at different levels and scales. Letλg be the SE at scale λ, which is obtained by λ−1 times dilation operation; λg is described as [27](5)λg=g⊕g⊕⋯⊕g︸λ−1times.The multiscale CMF can be redefined as(6)yλn=f∘λg•λgn+f•λg∘λgn2.CMF noise reduction results can effectively restore the original signal waveform. Also, the weighted average operation can restrain the statistical bias of denoising results. Indeed, compared with a single-scale CMF, the multiscale CMF is more effective in noise reduction.


## 2.2. Signal Reconstruction

These results of morphological filtering at different scales can be obtained by multiscale morphological analysis of noisy signal. Each filtering result at single scale reflects a specific component contained in the signal. The different components of chaotic signals should appear in multiple scales as a processed signal. The reconstructed signal is obtained by aggregating the filtering results of different scales as follows:(7)Y=∑j=1jωjyj∑j=1jωj=1,where yj is a result of morphological filtering on scale j, ωj is the weight coefficient of yj, and Y is the reconstruction result.The weight coefficient directly affects the effect of multiscale morphological filtering. For the purpose of obtaining a better denoising result, the square operation of residual amplitude difference (RAD) before and after filtering processing is used to determine the weight coefficients adaptively. The RAD can be defined as(8)Rjn=fn−yjn.The weight coefficient can be expressed as(9)ωj=∑nRjn2∑j∑nRjn2.


## 2.3. Signal Smoothing

It is found that the waveform of the aggregated signal is not smooth enough, and the local waveform is slightly rough. Owing to the morphological filter with large-scale SE blurring the detail information of signal and bringing the top-weakening distortion in the process of reducing noise, the signal directly gathered by the results of morphological filtering at different scales is not smooth enough. However, the occurrence of the top-weakening distortion is regular all the way. The distortion only made some local waveforms a little rough, but did not change the general shape of the whole waveform. In order to reduce the top-weakening distortion, a new method of resampling the filtered signal is proposed according to the regularity of the distortion. The cubic B-spline interpolation method is proposed to smooth the filtered signal after resampling.Letdii=0,1,…,n be the de Boor point, where i is the serial number and n is the number of the de Boor point. The k-th B-spline function is defined as follows [28, 29]:(10)Pt=∑i=0ndiBi,kt,t∈0,1,Bi,0t=1,ti≤t≤ti+1,0,others,Bi,kt=t−titi+k−1−tiBi,k−1t+ti+k−tti+k−ti+1Bi+1,k−1t,where t0,t1,…,tn+k+1 are node vectors and Bi,kt is a base function of the node vector. When the denominator in equation (10) is 0, the value is specified as 0. When k is 3, the cubic uniform B-spline function defined in [0, 1] is described as(11)pt=1−t36d0+3t3−6t2+46d1+−3t3+3t2+3t+16d2+t36d3.Equation (11) is a cubic polynomial divided into uniform segments; the signal processed by cubic B-spline interpolation has the characteristics of detail smoothness and linear independence.


## 3. Algorithmic Parameter Analysis and Adaptive Optimization

In this section, the Lorenz signal with superimposed Gaussian white noise is taken as an example to discuss the influence of algorithm parameters on noise reduction performance. The equation of Lorenz system is expressed as [30](12)x˙=−ax−y,y˙=−xz+cx−y,z˙=xy−bz,where a=10, b=8/3, and c=28; the system is in chaos. The fourth-order Runge–Kutta method is used to solve the equation, and the sampling time is set to 0.01. Chaotic signals are generated by the state variable x. Figure 2(a) shows the waveform of state variables changing with time. The chaotic signal is superimposed with white Gaussian noise of different intensities. Figure 2(b) shows a Lorenz signal with a superimposed signal to noise ratio (SNR) of 15 dB noise. The SNR is defined as [31](13)SNR=10log10∑ts2t∑txt−st2,where st is the original chaotic time series. When calculating the signal-to-noise ratio of the input signal (SNRin), the xt is regarded as a noisy chaotic signal. Nevertheless, if the calculated object is the output signal to noise ratio (SNRout), the xt is regarded as the denoised signal.Figure 2
Time domain waveform of (a) noise-free signal and (b) noisy signals.
(a)
(b)The ASMMF method makes full use of the self-similar structure characteristics of chaotic signals and has good chaotic attractor reduction performance. Just like all the MF based algorithms, ASMMF algorithm also needs to specify the shape and select the optimal scale combination of SE.


### 3.1. The Shape Selection of Structural Element

The SE can be used as a filter window in the process of MF, and its geometric feature matches the signal much better, and then main waveform of this part of the signal can be retained. The shape and size of structural element play an important role in the application of mathematical morphology. The triangular SE, the sinusoidal SE, and the linear SE are widely used SEs for signal analysis [32–36].Morphological operation is a numerical operation of adding or subtracting extreme values. Flat SE is defined as the linear SE with height of 0, which can retain the shape feature of the impulse as much as possible [34]. And it has the simplest algorithm and higher operation efficiency [35]. However, if the raw signal is processed by a scale far away from the theoretical central scale, the filtered result is often heavily polluted. Fortunately, multiscale morphological filtering with adaptive weighting can solve this problem [36]. Therefore, the flat SE is therefore chosen in filtering chaotic signal.


### 3.2. The Selection of Optimal Scale Combination

Set the initial SE to {0, 0, 0}. Scale of flat SE is used to represent the length of SE, and define the scale of initial SE to beλ=1. The difference between flat structural element at different scales is the number of zeros. The selection of scale coefficients of structural element of morphological filters has a significant effect on the denoising results. If the scale of structural element is relatively small, as shown in Figure 3(a), the effect of morphological filter on noise suppression is not obvious. On the contrary, when the scale of structural element is too large, the effect of noise reduction is obvious as shown in Figure 3(b), but the waveform has serious distortion.Figure 3
Time domain waveform of (a) morphological filtering result at scale 1 and (b) morphological filtering result at scale 25.
(a)
(b)The experimental analysis shows that the optimal scale combination is affected by the signal characteristics and noise level, and structural element directly affects the effect of morphological filtering. In order to achieve the best filtering result, it is necessary to select the appropriate scale combination of structural element. This optimal filtering parameter is constantly changing in practical application. How to realize the automatic optimization of filter parameters determines the adaptability of the algorithm. In order to increase the adaptability of the algorithm, an adaptive method for determining the optimal scale combination parameter is proposed in this paper. The flowchart to determine the optimal size combination is shown in Figure4.Figure 4
Flowchart of the optimal scale combination selection.As shown in Figure4, the adaptive selection process of optimal scale combination is divided into two steps.Step 1.
The selection of the maximum scale structural element. An appropriate maximum scale SE is able to reduce the computational complexity and operating time. The maximum scale of structural element is denoted asλmax. The experiment shows that the optimal length of the SE is closely related to the sampling frequency of the signal and the number of sampling points in the fault period. λmax should be smaller than the sampling point in the fault period, which is corresponding to a period of Lorenz signal. The number of sampling points corresponding to a period is denoted by fs/f0, where fs is sampling frequency (fs=100Hz) and f0 is the characteristic frequency of signal. According to the calculation method of Reference [37], the characteristic frequency is defined as(14)f0=∑t2πNtt,where NT is the number of turns performed in T.(15)λmax=fsf0,where symbol • represents round down the number.
The characteristic frequency of Lorenz signal isf≈7.04 Hz; therefore, fs/f≈1.40Hz Hz. λmax should be less than 15, so set λmax to 14 adaptively.Step 2.
The multiscale permutation entropy of residual signal is used to select the optimal scale combination of SE adaptively. MPE is an average entropy parameter developed on the basis of PE [14] to describe the irregularity degree and self-similarity of time series at different scales. The calculation process of MPE is described as follows.
Setxt,t=1,2,…,N to be a time series, and the scale τ represents the window length of coarsening. In the first place, the time series is coarsely granulated under the condition of scale factor τ. And the resampling of scale τ is carried out; the sequence after coarse granulation is denoted as(16)zjτ=1τ∑t=j−1τ+1jτxt,where j=1,2,…,N/τ, and • denotes downward rounding.
Then, the phase space of time serieszτj,j=1,2,…,N is reconstructed as(17)Zτj=zτj,zτj+ε,…,zτj+m−1ε,where m indicates the embedded dimension and ε indicates the delay time. The m data points of zτj are rearranged in ascending order as(18)Zτj=zτj+j1−1ε≤zτj+j2−1ε≤⋯≤zτj+jm−1ε,where jd,d=1,2,…,m is the column index of each element in the reconstructed component. In case there are m! symbol sequences at most, k≤m!. The frequency of occurrence of each permutation is recorded as fkτ,k=1,2,…k, ∑kfk=1; afterwards, the probability of each sequence in ε scale is described as(19)pkτ=fkτN−m−1ε.
According to Reference [14, 15], the MPE is ultimately defined as(20)Pτ=−∑k=1m!pkτlnpkτ.
In order to facilitate the analysis, the MPE is normalized by the maximum value oflnm!, that is,(21)Pτ=Pτlnm!0<Pτ<1.
The MPE vector can be described asP=P1,P2,…. The MPE consists of four parameters: embedding dimension m, time delay ε, length of time series N, and scale factor τ. The PE value highly depends on the selected embedding dimension m and the scale factor τ. The PE method works with the embedding dimension 3≤m≤7 [38], and the embedding dimension m and time series N should satisfy the condition 5m!≤N [39]; then, the embedding dimension is set to 4 in this paper. According to Reference [40], the delay time is set to 1. The selection criterion of scale factor is given as τ≤N/m+1! in Reference [41], and the maximum value of scale factor is 10 in this paper. Figure 5 shows the MPE analysis of morphological filtering results at different scales.Figure 5
The MPE of morphological filtering results: (a)τ=1,2,…,5; (b) τ=6,7,…,10.
(a)
(b)According to the analysis, with the increase of structural element scale, we find that the MPE values decrease gradually. And the complexity of the sequence decreases with the increase of the scale. The increase in scale makes the filtered sequence more orderly; however, it does not mean that the signal reduction effect is better. As shown in Figure3(b), the filtering results with too large scale of structural element will also have smaller PE values.In this paper, the optimal filtering scale is determined by the residual multiscale permutation entropy of the filtering result. If the residual MPE is the largest, it means that the randomness of residual is the strongest and the filtering effect is the best. The MPE of multiscale morphological filtering residuals is shown in Figure6.Figure 6
The MPE of morphological filter residuals: (a)τ=1,2,…,5; (b) τ=6,7,…,10.
(a)
(b)As shown, PE of the residuals corresponding to differentτ changes with the increase of λ, and all go through the process of gradually increasing from small to a maximum and then gradually decreasing. λ corresponding to the maximum of τ scale can be used as the optimal filtering scale coefficient of chaotic signal at this scale λbestτ. Λ=λbest1,λbest2,…,λbestτis initially determined to be the optimal combination of scale coefficients. Those repeat scales should be removed to ensure that each scale appears only one time. Finally, the optimal combination of scale coefficients for noisy signals as shown in Figure 2(b) turns out to be Λ=1,4,6,7,8,9,13.After determining the optimal combination of filtering scales, the weighting coefficients of filtering results are calculated according to equation (9), and the chaotic signal is aggregated according to equation (7). The aggregation result is the noise reduction result by AMMF, which is shown in Figure 7(a). Compared with the noisy chaotic signal and the original signal, the waveform of the aggregated signal is similar to the original signal as a whole, but the local details are slightly rough. However, after smoothing, as shown in Figure 7(b), the reconstructed signal waveform is basically the same as the original signal.Figure 7
Results of denoising: (a) unsmoothed morphological filtering and (b) smooth morphological filtering.
(a)
(b)


## 3.1. The Shape Selection of Structural Element

The SE can be used as a filter window in the process of MF, and its geometric feature matches the signal much better, and then main waveform of this part of the signal can be retained. The shape and size of structural element play an important role in the application of mathematical morphology. The triangular SE, the sinusoidal SE, and the linear SE are widely used SEs for signal analysis [32–36].Morphological operation is a numerical operation of adding or subtracting extreme values. Flat SE is defined as the linear SE with height of 0, which can retain the shape feature of the impulse as much as possible [34]. And it has the simplest algorithm and higher operation efficiency [35]. However, if the raw signal is processed by a scale far away from the theoretical central scale, the filtered result is often heavily polluted. Fortunately, multiscale morphological filtering with adaptive weighting can solve this problem [36]. Therefore, the flat SE is therefore chosen in filtering chaotic signal.


## 3.2. The Selection of Optimal Scale Combination

Set the initial SE to {0, 0, 0}. Scale of flat SE is used to represent the length of SE, and define the scale of initial SE to beλ=1. The difference between flat structural element at different scales is the number of zeros. The selection of scale coefficients of structural element of morphological filters has a significant effect on the denoising results. If the scale of structural element is relatively small, as shown in Figure 3(a), the effect of morphological filter on noise suppression is not obvious. On the contrary, when the scale of structural element is too large, the effect of noise reduction is obvious as shown in Figure 3(b), but the waveform has serious distortion.Figure 3
Time domain waveform of (a) morphological filtering result at scale 1 and (b) morphological filtering result at scale 25.
(a)
(b)The experimental analysis shows that the optimal scale combination is affected by the signal characteristics and noise level, and structural element directly affects the effect of morphological filtering. In order to achieve the best filtering result, it is necessary to select the appropriate scale combination of structural element. This optimal filtering parameter is constantly changing in practical application. How to realize the automatic optimization of filter parameters determines the adaptability of the algorithm. In order to increase the adaptability of the algorithm, an adaptive method for determining the optimal scale combination parameter is proposed in this paper. The flowchart to determine the optimal size combination is shown in Figure4.Figure 4
Flowchart of the optimal scale combination selection.As shown in Figure4, the adaptive selection process of optimal scale combination is divided into two steps.Step 1.
The selection of the maximum scale structural element. An appropriate maximum scale SE is able to reduce the computational complexity and operating time. The maximum scale of structural element is denoted asλmax. The experiment shows that the optimal length of the SE is closely related to the sampling frequency of the signal and the number of sampling points in the fault period. λmax should be smaller than the sampling point in the fault period, which is corresponding to a period of Lorenz signal. The number of sampling points corresponding to a period is denoted by fs/f0, where fs is sampling frequency (fs=100Hz) and f0 is the characteristic frequency of signal. According to the calculation method of Reference [37], the characteristic frequency is defined as(14)f0=∑t2πNtt,where NT is the number of turns performed in T.(15)λmax=fsf0,where symbol • represents round down the number.
The characteristic frequency of Lorenz signal isf≈7.04 Hz; therefore, fs/f≈1.40Hz Hz. λmax should be less than 15, so set λmax to 14 adaptively.Step 2.
The multiscale permutation entropy of residual signal is used to select the optimal scale combination of SE adaptively. MPE is an average entropy parameter developed on the basis of PE [14] to describe the irregularity degree and self-similarity of time series at different scales. The calculation process of MPE is described as follows.
Setxt,t=1,2,…,N to be a time series, and the scale τ represents the window length of coarsening. In the first place, the time series is coarsely granulated under the condition of scale factor τ. And the resampling of scale τ is carried out; the sequence after coarse granulation is denoted as(16)zjτ=1τ∑t=j−1τ+1jτxt,where j=1,2,…,N/τ, and • denotes downward rounding.
Then, the phase space of time serieszτj,j=1,2,…,N is reconstructed as(17)Zτj=zτj,zτj+ε,…,zτj+m−1ε,where m indicates the embedded dimension and ε indicates the delay time. The m data points of zτj are rearranged in ascending order as(18)Zτj=zτj+j1−1ε≤zτj+j2−1ε≤⋯≤zτj+jm−1ε,where jd,d=1,2,…,m is the column index of each element in the reconstructed component. In case there are m! symbol sequences at most, k≤m!. The frequency of occurrence of each permutation is recorded as fkτ,k=1,2,…k, ∑kfk=1; afterwards, the probability of each sequence in ε scale is described as(19)pkτ=fkτN−m−1ε.
According to Reference [14, 15], the MPE is ultimately defined as(20)Pτ=−∑k=1m!pkτlnpkτ.
In order to facilitate the analysis, the MPE is normalized by the maximum value oflnm!, that is,(21)Pτ=Pτlnm!0<Pτ<1.
The MPE vector can be described asP=P1,P2,…. The MPE consists of four parameters: embedding dimension m, time delay ε, length of time series N, and scale factor τ. The PE value highly depends on the selected embedding dimension m and the scale factor τ. The PE method works with the embedding dimension 3≤m≤7 [38], and the embedding dimension m and time series N should satisfy the condition 5m!≤N [39]; then, the embedding dimension is set to 4 in this paper. According to Reference [40], the delay time is set to 1. The selection criterion of scale factor is given as τ≤N/m+1! in Reference [41], and the maximum value of scale factor is 10 in this paper. Figure 5 shows the MPE analysis of morphological filtering results at different scales.Figure 5
The MPE of morphological filtering results: (a)τ=1,2,…,5; (b) τ=6,7,…,10.
(a)
(b)According to the analysis, with the increase of structural element scale, we find that the MPE values decrease gradually. And the complexity of the sequence decreases with the increase of the scale. The increase in scale makes the filtered sequence more orderly; however, it does not mean that the signal reduction effect is better. As shown in Figure3(b), the filtering results with too large scale of structural element will also have smaller PE values.In this paper, the optimal filtering scale is determined by the residual multiscale permutation entropy of the filtering result. If the residual MPE is the largest, it means that the randomness of residual is the strongest and the filtering effect is the best. The MPE of multiscale morphological filtering residuals is shown in Figure6.Figure 6
The MPE of morphological filter residuals: (a)τ=1,2,…,5; (b) τ=6,7,…,10.
(a)
(b)As shown, PE of the residuals corresponding to differentτ changes with the increase of λ, and all go through the process of gradually increasing from small to a maximum and then gradually decreasing. λ corresponding to the maximum of τ scale can be used as the optimal filtering scale coefficient of chaotic signal at this scale λbestτ. Λ=λbest1,λbest2,…,λbestτis initially determined to be the optimal combination of scale coefficients. Those repeat scales should be removed to ensure that each scale appears only one time. Finally, the optimal combination of scale coefficients for noisy signals as shown in Figure 2(b) turns out to be Λ=1,4,6,7,8,9,13.After determining the optimal combination of filtering scales, the weighting coefficients of filtering results are calculated according to equation (9), and the chaotic signal is aggregated according to equation (7). The aggregation result is the noise reduction result by AMMF, which is shown in Figure 7(a). Compared with the noisy chaotic signal and the original signal, the waveform of the aggregated signal is similar to the original signal as a whole, but the local details are slightly rough. However, after smoothing, as shown in Figure 7(b), the reconstructed signal waveform is basically the same as the original signal.Figure 7
Results of denoising: (a) unsmoothed morphological filtering and (b) smooth morphological filtering.
(a)
(b)


## 4. Experimental Analysis

In order to verify the effectiveness of proposed ASMMF for different chaotic signals, the existing EMD-CIT [9], WT [3], and AMMF methods are applied for comparison, and the noise suppression results of the Lorenz signals with different intensities of Gaussian white noises, the sunspot chaotic signals, and chaotic vibration signals are, respectively, shown in this section.


### 4.1. Lorenz Chaotic Signal Denoising

The total length of the original chaotic Lorenz signal is set at 8000.SNRout and the root mean square error (RMSE) [42] indexes are both used to evaluate the noise reduction performance. RMSE can be defined as follows:(22)RMSE=1N∑tRt2.Figure8 shows the denoising results of different noisy signals processed by four different algorithms. In the results, when the signal-to-noise ratio of input signal is low, SNRout of the unsmoothed adaptive multiscale morphological filtering is slightly lower than that of the EMD threshold method and the wavelet threshold method. As the noise intensity increases gradually, its performance index gradually exceeds the wavelet threshold method, but the overall denoising performance is poor. The morphological filtering algorithm after smoothing is superior to the other three methods in both the aspects.Figure 8
Analysis of signal denoised by different algorithms: (a)SNRout; (b) RMSE.
(a)
(b)In order to intuitively reflect the denoising results of different algorithms, the phase diagrams based on the chaotic characteristic are drawn ulteriorly. Figure9 shows the denoising output results by four methods for noisy chaotic signal with input signal-to-noise ratio of 10 dB. The two-dimensional phase diagrams of the clean signal and the noisy signal are plotted in Figures 9(a) and 9(b). The noise reduction results of WT and EMD-CIIT are shown in Figures 9(c) and 9(d). According the phase diagrams, both of them can restore the attractor shape of chaotic signals smoothly, but the self-similarity structure of the system cannot be clearly expressed due to the orbital deformation. Unsmoothed morphological filtering result is shown in Figure 9(e) The shape of the attractor restored is not smooth enough, and the effect is not good enough. Furthermore, the cancellation distortion of local waveform in the process of noise reduction is obvious. The ASMMF is an improvement algorithm of AMMF. After smoothing the rough waveform of morphological filtering, the shape and self-similar structure of chaotic attractors can be restored more clearly.Figure 9
Phase diagram of Lorenz attractor: (a) original signal; (b) noisy signals; (c) signal denoised by EMD thresholding; (d) signal denoised by wavelet thresholding; (e) signal denoised by AMMF; (f) signal denoised by ASMMF.
(a)
(b)
(c)
(d)
(e)
(f)


### 4.2. Sunspot Chaotic Signal Denoising

Sunspots are a phenomenon of temperature change on the surface of the solar photosphere, and their activity will directly affect the changes of the Earth's magnetic field and climate. Recently, chaos theory has become an important means to study the rule of sunspot activity. The observed sunspot signal is always chaotic and noisy [10, 31], so it is necessary to effectively reduce the noise of the observed data. The sunspot sequence from January 1749 to September 2016 provided by NASA is selected as experimental data (http://solarscience.msfc.nasa.gov/green-wch.shtml). The Lyapunov exponent (LE) of the attractor represents the predictability and the sensitivity of dynamical system. Therefore, the LE is an important feature for the study of chaotic behavior. For a chaotic system, the largest Lyapunov exponent (LLE) is positive [43]. Referring to the calculation method in Reference [36], the LLE is 0.7615.To investigate the method for different kinds of noise, EMD-CIT, WT, AMMF, and ASMMF are used to denoise the sunspot data. Compared with the original sunspot observed data shown in Figure10(a), there is obvious measurement noise in the observed data. Figures 10(b)–10(e) show the denoising results. As shown in the figures, we can notice that all the techniques significantly reduce the noise obviously. But the waveform of EMD-CIT, WT, and AMMF noise reduction results is not smooth enough. And the waveform of ASMMF after noise reduction is the smoothest.Figure 10
(a) Observation data of monthly mean number of sunspot plane. (b) The noise reduction result by EMD thresholding. (c) The noise reduction result by wavelet thresholding. (d) The noise reduction result by AMMF. (e) The noise reduction result by ASMMF.
(a)
(b)
(c)
(d)
(e)The projection of the attractor reconstructed from the data sequence on the two-dimensional plane before and after noise reduction is shown in Figure11. It is clear that the chaotic attractor phase diagram of the measured lunar sunspot signal is very disordered due to the noise interference, and the attractor trajectory in the phase diagram space is not smooth. No regular geometric structure is shown in Figure 11(a). Noise reduction results of four mentioned methods are shown in Figures 11(b)–11(e); all of them obviously reduce the noise and obtain the geometric structure of attractor with strong regularity. Moreover, the attractor orbits obtained by the ASMMF method are smoother and clearer than those of the other three methods. It shows that this method has certain advantages in reducing noise and enhancing the inherent deterministic components of chaotic sequences.Figure 11
Phase diagrams for the sunspot signal: (a) observation data; (b) the noise reduction result by EMD thresholding; (c) the noise reduction result by wavelet thresholding; (d) the noise reduction result by AMMF; (e) the noise reduction result by ASMMF.
(a)
(b)
(c)
(d)
(e)In order to intuitively verify the advantages and disadvantages of noise reduction methods, MPE is used to describe the complexity and regularity of time series on multiple time scales. Figure12 shows the analysis results by MPE.Figure 12
The MPE of denoising results.By direct comparison, the four methods all perform obvious effect in chaotic signal denoising. The result of noise reduction by these methods is that the MPE value decreases compared with the signal before noise reduction, and the signal after noise reduction is more regular. If the time scale is small, the PE of the unsmoothed adaptive multiscale morphological filtering is lower than that of the EMD threshold method and the wavelet threshold method. As the time scale increases gradually, its performance index gradually exceeds the wavelet threshold method, but the overall denoising performance is poor. The denoising result using ASMMF algorithm has the lowest complexity in different time scales. It indicates that the intrinsic deterministic components of chaotic sequences are enhanced after noise reduction. The simulation results show that compared with WT, EMD-CIIT, and AMMF, ASMMF has a better denoising effect on the observed sunspot sequence, with a smoother phase diagram, and higher self similarity.


### 4.3. Chaotic Vibration Signal Denoising

Detection of faults of rotating machinery by the complex noisy vibration signals is very important and difficult. In recent years, features of chaotic behavior are very good for fault diagnosis and identification. But the features of the LLE, approximate entropy, and correlation dimension are easily affected by noise. Clear attractor orbits and signals that do not damage nonlinear dynamical structures are helpful to obtain chaotic behavior characteristics accurately. So, detecting chaotic vibration signals in noise environment is a significant task.In the experiment, the Bently RK-4 rotor system test bed was used to simulate the fault; it is equipped with a signal preadapter, speed control device, and bearing and oil pump system for oil film instability, which can realize the oil film whirl and oil film oscillation test. Eddy current sensors were installed on both sides of the single-disk rotor to measure radial vibration displacement of rotating shaft. ZonicBook/618E produced by Iotech Company in the United States collects vibration signals of rotors during rotation with sampling frequency of 1280 Hz. When the rotation speed is maintained at 4560 r/min, the rotor system shows obvious resonance oil film oscillation phenomenon, and the time-domain waveform is shown in Figure13.Figure 13
The time-domain waveform of vibration signals.As mentioned inSection 4.2, the LE is an important feature for the study of chaotic behavior. For a chaotic system, the LLE is positive. Figure 14(a) shows the largest Lyapunov exponents, and phase diagram for the vibration signal is positive. So, the collected vibration signal is chaotic. Figure 14(b) shows the phase diagram for the vibration signal. Due to the influence of noise, the chaotic attractors cannot be clearly displayed.Figure 14
The largest Lyapunov exponents and phase diagram for the vibration signal.
(a)
(b)Due to the influence of noise, the chaotic attractors cannot be clearly displayed. It is necessary to denoise the collected vibration signal. In the experiment, the contaminated chaotic vibration signal is denoised by the 4 noise reduction methods mentioned. Figure15 shows the denoising results.Figure 15
Phase diagrams for the vibration signal: (a) the noise reduction result by EMD thresholding; (b) the noise reduction result by wavelet thresholding; (c) the noise reduction result by AMMF; (d) the noise reduction result by ASMMF.
(a)
(b)
(c)
(d)The phase diagram of the ASMMF method after noise reduction is smoother, clearer, and more regular than other methods. As shown in Figure16, as the time scale increases gradually, with the increase of time scale, the PE values of AMMF, WT, and EMD CIIT are not significantly lower than before noise reduction. The PE value of noise reduction result by WT algorithm on some time scales is even greater than that before noise reduction, and the overall denoising performance is poor. The MPE of the ASMMF method’s denoising result is the minimum, which represents that the result of denoising by this algorithm is the most regular. The experimental results show that this method can also achieve good noise reduction effect on the vibration signal, and the denoising performance is better than the proposed three other methods.Figure 16
The MPE of denoising results for the chaotic vibration signal.


## 4.1. Lorenz Chaotic Signal Denoising

The total length of the original chaotic Lorenz signal is set at 8000.SNRout and the root mean square error (RMSE) [42] indexes are both used to evaluate the noise reduction performance. RMSE can be defined as follows:(22)RMSE=1N∑tRt2.Figure8 shows the denoising results of different noisy signals processed by four different algorithms. In the results, when the signal-to-noise ratio of input signal is low, SNRout of the unsmoothed adaptive multiscale morphological filtering is slightly lower than that of the EMD threshold method and the wavelet threshold method. As the noise intensity increases gradually, its performance index gradually exceeds the wavelet threshold method, but the overall denoising performance is poor. The morphological filtering algorithm after smoothing is superior to the other three methods in both the aspects.Figure 8
Analysis of signal denoised by different algorithms: (a)SNRout; (b) RMSE.
(a)
(b)In order to intuitively reflect the denoising results of different algorithms, the phase diagrams based on the chaotic characteristic are drawn ulteriorly. Figure9 shows the denoising output results by four methods for noisy chaotic signal with input signal-to-noise ratio of 10 dB. The two-dimensional phase diagrams of the clean signal and the noisy signal are plotted in Figures 9(a) and 9(b). The noise reduction results of WT and EMD-CIIT are shown in Figures 9(c) and 9(d). According the phase diagrams, both of them can restore the attractor shape of chaotic signals smoothly, but the self-similarity structure of the system cannot be clearly expressed due to the orbital deformation. Unsmoothed morphological filtering result is shown in Figure 9(e) The shape of the attractor restored is not smooth enough, and the effect is not good enough. Furthermore, the cancellation distortion of local waveform in the process of noise reduction is obvious. The ASMMF is an improvement algorithm of AMMF. After smoothing the rough waveform of morphological filtering, the shape and self-similar structure of chaotic attractors can be restored more clearly.Figure 9
Phase diagram of Lorenz attractor: (a) original signal; (b) noisy signals; (c) signal denoised by EMD thresholding; (d) signal denoised by wavelet thresholding; (e) signal denoised by AMMF; (f) signal denoised by ASMMF.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.2. Sunspot Chaotic Signal Denoising

Sunspots are a phenomenon of temperature change on the surface of the solar photosphere, and their activity will directly affect the changes of the Earth's magnetic field and climate. Recently, chaos theory has become an important means to study the rule of sunspot activity. The observed sunspot signal is always chaotic and noisy [10, 31], so it is necessary to effectively reduce the noise of the observed data. The sunspot sequence from January 1749 to September 2016 provided by NASA is selected as experimental data (http://solarscience.msfc.nasa.gov/green-wch.shtml). The Lyapunov exponent (LE) of the attractor represents the predictability and the sensitivity of dynamical system. Therefore, the LE is an important feature for the study of chaotic behavior. For a chaotic system, the largest Lyapunov exponent (LLE) is positive [43]. Referring to the calculation method in Reference [36], the LLE is 0.7615.To investigate the method for different kinds of noise, EMD-CIT, WT, AMMF, and ASMMF are used to denoise the sunspot data. Compared with the original sunspot observed data shown in Figure10(a), there is obvious measurement noise in the observed data. Figures 10(b)–10(e) show the denoising results. As shown in the figures, we can notice that all the techniques significantly reduce the noise obviously. But the waveform of EMD-CIT, WT, and AMMF noise reduction results is not smooth enough. And the waveform of ASMMF after noise reduction is the smoothest.Figure 10
(a) Observation data of monthly mean number of sunspot plane. (b) The noise reduction result by EMD thresholding. (c) The noise reduction result by wavelet thresholding. (d) The noise reduction result by AMMF. (e) The noise reduction result by ASMMF.
(a)
(b)
(c)
(d)
(e)The projection of the attractor reconstructed from the data sequence on the two-dimensional plane before and after noise reduction is shown in Figure11. It is clear that the chaotic attractor phase diagram of the measured lunar sunspot signal is very disordered due to the noise interference, and the attractor trajectory in the phase diagram space is not smooth. No regular geometric structure is shown in Figure 11(a). Noise reduction results of four mentioned methods are shown in Figures 11(b)–11(e); all of them obviously reduce the noise and obtain the geometric structure of attractor with strong regularity. Moreover, the attractor orbits obtained by the ASMMF method are smoother and clearer than those of the other three methods. It shows that this method has certain advantages in reducing noise and enhancing the inherent deterministic components of chaotic sequences.Figure 11
Phase diagrams for the sunspot signal: (a) observation data; (b) the noise reduction result by EMD thresholding; (c) the noise reduction result by wavelet thresholding; (d) the noise reduction result by AMMF; (e) the noise reduction result by ASMMF.
(a)
(b)
(c)
(d)
(e)In order to intuitively verify the advantages and disadvantages of noise reduction methods, MPE is used to describe the complexity and regularity of time series on multiple time scales. Figure12 shows the analysis results by MPE.Figure 12
The MPE of denoising results.By direct comparison, the four methods all perform obvious effect in chaotic signal denoising. The result of noise reduction by these methods is that the MPE value decreases compared with the signal before noise reduction, and the signal after noise reduction is more regular. If the time scale is small, the PE of the unsmoothed adaptive multiscale morphological filtering is lower than that of the EMD threshold method and the wavelet threshold method. As the time scale increases gradually, its performance index gradually exceeds the wavelet threshold method, but the overall denoising performance is poor. The denoising result using ASMMF algorithm has the lowest complexity in different time scales. It indicates that the intrinsic deterministic components of chaotic sequences are enhanced after noise reduction. The simulation results show that compared with WT, EMD-CIIT, and AMMF, ASMMF has a better denoising effect on the observed sunspot sequence, with a smoother phase diagram, and higher self similarity.


## 4.3. Chaotic Vibration Signal Denoising

Detection of faults of rotating machinery by the complex noisy vibration signals is very important and difficult. In recent years, features of chaotic behavior are very good for fault diagnosis and identification. But the features of the LLE, approximate entropy, and correlation dimension are easily affected by noise. Clear attractor orbits and signals that do not damage nonlinear dynamical structures are helpful to obtain chaotic behavior characteristics accurately. So, detecting chaotic vibration signals in noise environment is a significant task.In the experiment, the Bently RK-4 rotor system test bed was used to simulate the fault; it is equipped with a signal preadapter, speed control device, and bearing and oil pump system for oil film instability, which can realize the oil film whirl and oil film oscillation test. Eddy current sensors were installed on both sides of the single-disk rotor to measure radial vibration displacement of rotating shaft. ZonicBook/618E produced by Iotech Company in the United States collects vibration signals of rotors during rotation with sampling frequency of 1280 Hz. When the rotation speed is maintained at 4560 r/min, the rotor system shows obvious resonance oil film oscillation phenomenon, and the time-domain waveform is shown in Figure13.Figure 13
The time-domain waveform of vibration signals.As mentioned inSection 4.2, the LE is an important feature for the study of chaotic behavior. For a chaotic system, the LLE is positive. Figure 14(a) shows the largest Lyapunov exponents, and phase diagram for the vibration signal is positive. So, the collected vibration signal is chaotic. Figure 14(b) shows the phase diagram for the vibration signal. Due to the influence of noise, the chaotic attractors cannot be clearly displayed.Figure 14
The largest Lyapunov exponents and phase diagram for the vibration signal.
(a)
(b)Due to the influence of noise, the chaotic attractors cannot be clearly displayed. It is necessary to denoise the collected vibration signal. In the experiment, the contaminated chaotic vibration signal is denoised by the 4 noise reduction methods mentioned. Figure15 shows the denoising results.Figure 15
Phase diagrams for the vibration signal: (a) the noise reduction result by EMD thresholding; (b) the noise reduction result by wavelet thresholding; (c) the noise reduction result by AMMF; (d) the noise reduction result by ASMMF.
(a)
(b)
(c)
(d)The phase diagram of the ASMMF method after noise reduction is smoother, clearer, and more regular than other methods. As shown in Figure16, as the time scale increases gradually, with the increase of time scale, the PE values of AMMF, WT, and EMD CIIT are not significantly lower than before noise reduction. The PE value of noise reduction result by WT algorithm on some time scales is even greater than that before noise reduction, and the overall denoising performance is poor. The MPE of the ASMMF method’s denoising result is the minimum, which represents that the result of denoising by this algorithm is the most regular. The experimental results show that this method can also achieve good noise reduction effect on the vibration signal, and the denoising performance is better than the proposed three other methods.Figure 16
The MPE of denoising results for the chaotic vibration signal.


## 5. Conclusions

We have proposed an adaptive smoothing morphological filtering method for denoising contaminated chaotic signals. This algorithm solves the problem of top-weakening distortion in morphological filtering and can filter the signal in multilevel and multiscale adaptively. In this paper, AMMF, EMD-CIIT, and WT are compared with ASMMF. The polluted Lorenz signals, the sunspot chaotic signals, and chaotic vibration signals are denoised by 4 methods, respectively. The ASMMF is superior to the other three methods inSNRout, RMSE, and MPE for the noise suppression results of contaminated chaotic signals. Moreover, when the noise of signal is strong, the attractor trajectory of the EMD-CIIT and the WT denoising results have no regular geometric structure, but the result obtained by the proposed method has smoother and more regular attractor orbits. The comparison results show that the topology structure of the chaotic attractor restored by ASMMF is closer to the real signal. So, the proposed method has a certain reference value for actual engineering application. In the future, denoising of different kinds of chaotic signals will still be a hot topic. We will apply the proposed algorithm to the noise reduction of hyperchaotic signals and fractional chaotic signals.



---

*Source: 7242943-2020-02-17.xml*