# A Method Combining Fractal Analysis and Single Channel ICA for Vibration Noise Reduction

Authors: Quanbo Lu; Mei Li
Journal: Shock and Vibration
 (2021)

---

## Abstract

Aiming at the problem that real engineering vibration signals are interfered by strong noise, this paper proposes a method combining single channel-independent component analysis (SCICA) and fractal analysis (FD) to reduce the effect of noise on the time-frequency analysis of vibration signals. First, phase space reconstruction is performed on the vibration signal to make the proper input for ICA algorithm. The original is then decomposed into several component signals. The fractal dimension of each component signals is calculated to determine whether the signal should be considered noise. Noisy component signals are then processed by wavelet denoising. Finally, the output signal after noise reduction is reconstructed using the filtered “right” component signals. This paper uses the method to analyze real noisy vibration signal. Experimental results show the effectiveness of the proposed method.


---

## Body

## 1. Introduction

The measured vibration signals contain many high-frequency and low-frequency noises because of the test environment. Thus, the spectral components of the signal are complex. In addition, extracting the envelope of the signal is difficult for the reason of noise interference. Noise causes the inaccurate decomposition of the signal envelope, particularly, when using time-frequency analysis technology. An effective noise reduction method for time-frequency analysis would improve time-frequency resolution and algorithm performance, which are required both for practical applications and theoretical analysis. Wavelet denoising is a vibration signal, Qin et al. used the wavelet denoising method to improving the detection performance of distributed vibration sensor [1], and Xie et al. proposed an improved wavelet denoise method for vibration signal processing [2]. Other methods proposed an improved EMD method for time-frequency feature extraction of telemetry vibration [3]. Li et al. proposed a novel noise reduction technique for underwater acoustic signals based on complete ensemble empirical mode decomposition with adaptive noise, minimum mean square variance criterion, and least mean square adaptive filter [4]. For warship radio noise, Li et al. proposed a novel linear spectrum frequency feature extraction technique for warship radio noise based on complete ensemble empirical mode decomposition with adaptive noise, duffing chaotic oscillator, and weighted-permutation entropy [5]. Wang et al. proposed a new fault diagnosis method based on adaptive spectrum mode extraction [6]. For the data-driven fault diagnosis, Wang et al. proposed a novel technique based on AlexNet Convolutional neural network [7].ICA is a common blind source separation method but generally requires at least two input signals [8]. In 2001, James and Lowe proposed an algorithm based on a dynamic embedded framework to apply an ICA algorithm to extract multidimensional brain activity signals from single brain electromagnetic signals [9]. Davies and James suggested that SCICA can sometimes be used to separate important components from a time series [10]. And, in [11], the author reviewed the concept and current practice of BSS based on second-order statistics and on higher-order statistics, which is known as ICA [12], and applied ensemble empirical mode decomposition to make use of ICA on single-channel neuromuscular signals [13]. Guerrero has proposed an adaptive noise reduction method that combines ICA and recursive least square (RLS) for processing EEG signals [14]. However, this method filters all source signals after ICA decomposition. In other words, part of the useful signal will be likely filtered out as noise, which is unacceptable for vibration signal processing. If useful information is lost because of filtering, the analysis and identification of vibration signals will be ambiguous.To solve this problem, this paper proposes a vibration signal noise reduction method based on SCICA and FD. First, the Newton iteration ICA method is performed using real vibration signal based on kurtosis maximization of the vibration signal for blind source separation. Second, the component signal required for filtering is selected using fractal dimensions, different denoising methods are performed on the selected component signals, and finally, the signals are reconstructed. The results show that the proposed algorithm not only greatly suppresses noise but also effectively avoids the loss of useful information.


## 2. Methods

Abnormal signal filtering is especially problematic because noise greatly reduces the vibration signals. If the method filters the entire source signal, some information will be lost. To address this problem, this paper uses the signals fractal dimension to constrain the ICA filtering process, avoiding the loss of useful information. As shown in Figure1, the algorithm is described below:(1)
The collected vibration signal is preprocessed using zero drift correction and term trend elimination. Phase space reconstruction of the signal is then performed.(2)
Blind source separation is performed on the signal using kurtosis maximization via the Newton iterative ICA algorithm [15]. The source signal component S=S^1,S^2,…,S^M and separating matrix W are then calculated.(3)
The grid dimensionDi, i=1,2,…,M, and M of each source signal are calculated. The source signal is screen filtered using the grid dimension. According to the preliminary analysis of many experiments, the grid dimension of source signals containing an abnormal signal is greater than 1.5. As a conservative screening condition, this paper sets the maximum step point of the grid dimension curve close to 1.5. Specific operational procedures: sequence grid dimension Di is in ascending order. The most significant step point of the arrangement curve is the screening conditions. The source signal component S^ can be divided into one part S^f to be filtered and the other part S^b to be retained.(4)
Next, wavelet denoising is performed on the filtering source signalS^ff.(5)
The filtered source signalS^f and S^b are combined to produce a recombinant source signal S^′. The same relative position is maintained during the combination. The separating matrix W is used to reconstruct the phase space signal X=WS′. The noise reduction vibration signal is obtained from X. A description of the phase space reduction method may be found elsewhere.Figure 1
Total procedure.


### 2.1. Phase Space Reconstruction

The phase space reconstruction method is required when using ICA blind source separation for processing the single-channel vector data. In this paper, the chaotic phase space reconstruction method is performed on the signal to extend to multidimensional phase space. Thus, each aliasing source signal is reopened in a multidimensional space, and the ICA method is used to separate the source signals. In the 1980s, Packard et al. proposed the use of chaotic time series to reconstruct the phase space to reveal the dynamics of the system [16]. The basic idea is that the evolution of any component in the system is related to the change of all interacting components. By the Takens embedding theory for one dimension signal, a phase space equivalent to the original dynamic system can be constructed by determining its embedding dimension and adding time delay. The phase space contains the required information and characteristics of the original signal to obtain enough hidden information in the observation signal. As the dynamic embedding principle, the phase space reconstruction method construct delay vectors by delaying one-dimensional observation signal, and combining these vectors to build a multidimensional phase space can reveal the dynamic system in the chaotic attractor and chaotic time series. In order to obtain enough information content, the reconstructed phase space embedding dimension m should be at least greater than system freedom D, i.e., m≥2D+1.First, determine the embedding dimensionm and the time delay τ; for a complex signal xt, the time delay τ is added to obtain the delay vector at a time t:(1)xt=xt,xt+τ,…,xt+m−1τT,where t=1,2,…,N−m−1τ and N is the length of the observed signal. These delay vectors are combined to form an embedding matrix, which is a multidimensional delay signal:(2)x=x1x2⋯xN−m−1τx1+τx2+τ⋯xN−m−1τ+τ⋮⋮⋱⋮x1+m−1τx2+m−1τ⋯xN.Without prior information about vibration signals, selecting an embedding dimension is difficult. Instead, a delay time is selected before determining the embedding dimension according to the delay time. Many algorithms may be used to calculate the delay time, including the C-C algorithm [17], G-P algorithm, and mutual information method [18]. The algorithm processes are detailed elsewhere. This paper compared the performance and computational efficiency and used the mutual information method to obtain the delay time. Finally, using the method proposed by Takens [19], the Cao algorithm [20] was implemented to calculate the embedding dimension.


### 2.2. Newton Iteration ICA Algorithm Based on Kurtosis Maximization

The phase space reconstruction expended the single-channel signal into multiple channels. Thus, ICA can be applied. The ICA algorithm is a blind source separation method that uses the assumption of statistical independence of alias components. A separating matrixW can be constructed, and blind source separation can be realized as follows:(3)X=AS,S^=WX,where X=x1,x2,…,xNT is the aliasing signal vector, A is the mixing matrix, S=S1,S2,…,SMT is the source signal vector, and S^=S^1,S^2,…,S^MT is the estimated value of S. The traditional method, which is based on ICA kurtosis and uses the fixed point iterative method to calculate the separating matrix, has many shortcomings, such as low computational efficiency, the possibility that the separating matrix fails to converge, and poor separation effects. Li et al. proposed an alternative method based on the Newton iteration to calculate the separating matrix [21]. However, the results of the algorithm are unstable when the mean of the source signal is zero. To address these problems, this paper uses Newton iteration ICA by maximizing the kurtosis. The constraint matrix separation operation using Newton iterative algorithm is based on the kurtosis formula. The source signals may be obtained by reducing the Gaussian of the separation results. Simulation results show that this method improves the computational efficiency and separation effects. First, the kurtosis of the signal is defined as follows [22]:(4)kurts^=Es^4−3Es^22.Lets^=wTx, where x and w is a component of the separating matrix W. The gradient kurts^ of w is(5)∇wkurtwTx=4EwTx3−3EwTx2.To maximize the kurtosis, the extreme solution is calculated under the constraintsEwTx2=w=1 using the Lagrange algorithm as shown below:(6)l=kurtwTx=4EwTx3+1.5EwTx2−1.The iterative formula forw is calculated using the Newton iterative method:(7)wk+1=wk−ExwkTx33EwkTx2.After each iteration, the decorrelation extraction and normalization are performed on the linear combination of independent components to prevent these components from converging to the same extreme point, which leads to the instability of the algorithm because of the zero-mean of the source signal. The formulas for this method are shown as follows:(8)wp+1k+1=wp+1k−Exwp+1kTx33Ewp+1kTx2,wp+1k+1=wp+1k+1−∑j=1pwp+1k+1Twjwj,wp+1k+1=wp+1k+1wp+1k+1Twp+1k+1.The separating matrixW is calculated by estimating all the components w with the method described above.


### 2.3. Fractal Dimension

Fractal theory, an active area of signal processing research, is a mathematical tool that quantitatively describes various complex signals in nature. The fractal dimension of the vibration signal may be calculated to measure the complexity and irregularity of the signal. In fractal theory, larger fractal dimensions indicate that the signal is complex. For vibration signal processing, noise often causes the vibration signal components to be complicated, nonlinear, and nonstationary. Thus, the fractal dimension of the vibration signal is large. Many experimental results show the fractal dimension (grid dimension) of abnormal signals is greater than 1.5. Hence, the source signals with a larger fractal dimension may be filtered out avoiding meaningful information loss. The noise may be suppressed using a wavelet filter on the source signal, which has smaller fractal dimension. The fractal dimension contains the box dimension, correlation dimension, Hausdorff dimension, and grid dimension. For one-dimensional signals, the box dimension is equivalent to the grid dimension, which simplifies the computer processing. For an signal channel signalx|xi,i=1,2,…,N, the grid dimension is defined as(9)D=−lnMjlnΔt,where Δt is the sampling interval, Mj=∑xi+1−xi/Δt, and x represents the closest integer that is greater than x. First, the space of discrete signals is divided into a multidimensional grid, and the grid dimension is calculated by counting the grid points. The Euclidean space of discrete signals RN is divided as small as possible into grids with the area of Δ. Points’ count NΔK of grid space of the signal x|xi,i=1,2,…,N is calculated. Different areas NΔK are used to calculate the NΔK point, yielding a point set xk,yk:(10)xk=−lnΔK,yk=−lnNΔK.The slope of the line formed by the set of pointsxk,yk is the grid dimension of the signal. The line is not straight, so the least-squares method is often used to fit a line for practical applications. Finally, the grid dimension D may be obtained.


## 2.1. Phase Space Reconstruction

The phase space reconstruction method is required when using ICA blind source separation for processing the single-channel vector data. In this paper, the chaotic phase space reconstruction method is performed on the signal to extend to multidimensional phase space. Thus, each aliasing source signal is reopened in a multidimensional space, and the ICA method is used to separate the source signals. In the 1980s, Packard et al. proposed the use of chaotic time series to reconstruct the phase space to reveal the dynamics of the system [16]. The basic idea is that the evolution of any component in the system is related to the change of all interacting components. By the Takens embedding theory for one dimension signal, a phase space equivalent to the original dynamic system can be constructed by determining its embedding dimension and adding time delay. The phase space contains the required information and characteristics of the original signal to obtain enough hidden information in the observation signal. As the dynamic embedding principle, the phase space reconstruction method construct delay vectors by delaying one-dimensional observation signal, and combining these vectors to build a multidimensional phase space can reveal the dynamic system in the chaotic attractor and chaotic time series. In order to obtain enough information content, the reconstructed phase space embedding dimension m should be at least greater than system freedom D, i.e., m≥2D+1.First, determine the embedding dimensionm and the time delay τ; for a complex signal xt, the time delay τ is added to obtain the delay vector at a time t:(1)xt=xt,xt+τ,…,xt+m−1τT,where t=1,2,…,N−m−1τ and N is the length of the observed signal. These delay vectors are combined to form an embedding matrix, which is a multidimensional delay signal:(2)x=x1x2⋯xN−m−1τx1+τx2+τ⋯xN−m−1τ+τ⋮⋮⋱⋮x1+m−1τx2+m−1τ⋯xN.Without prior information about vibration signals, selecting an embedding dimension is difficult. Instead, a delay time is selected before determining the embedding dimension according to the delay time. Many algorithms may be used to calculate the delay time, including the C-C algorithm [17], G-P algorithm, and mutual information method [18]. The algorithm processes are detailed elsewhere. This paper compared the performance and computational efficiency and used the mutual information method to obtain the delay time. Finally, using the method proposed by Takens [19], the Cao algorithm [20] was implemented to calculate the embedding dimension.


## 2.2. Newton Iteration ICA Algorithm Based on Kurtosis Maximization

The phase space reconstruction expended the single-channel signal into multiple channels. Thus, ICA can be applied. The ICA algorithm is a blind source separation method that uses the assumption of statistical independence of alias components. A separating matrixW can be constructed, and blind source separation can be realized as follows:(3)X=AS,S^=WX,where X=x1,x2,…,xNT is the aliasing signal vector, A is the mixing matrix, S=S1,S2,…,SMT is the source signal vector, and S^=S^1,S^2,…,S^MT is the estimated value of S. The traditional method, which is based on ICA kurtosis and uses the fixed point iterative method to calculate the separating matrix, has many shortcomings, such as low computational efficiency, the possibility that the separating matrix fails to converge, and poor separation effects. Li et al. proposed an alternative method based on the Newton iteration to calculate the separating matrix [21]. However, the results of the algorithm are unstable when the mean of the source signal is zero. To address these problems, this paper uses Newton iteration ICA by maximizing the kurtosis. The constraint matrix separation operation using Newton iterative algorithm is based on the kurtosis formula. The source signals may be obtained by reducing the Gaussian of the separation results. Simulation results show that this method improves the computational efficiency and separation effects. First, the kurtosis of the signal is defined as follows [22]:(4)kurts^=Es^4−3Es^22.Lets^=wTx, where x and w is a component of the separating matrix W. The gradient kurts^ of w is(5)∇wkurtwTx=4EwTx3−3EwTx2.To maximize the kurtosis, the extreme solution is calculated under the constraintsEwTx2=w=1 using the Lagrange algorithm as shown below:(6)l=kurtwTx=4EwTx3+1.5EwTx2−1.The iterative formula forw is calculated using the Newton iterative method:(7)wk+1=wk−ExwkTx33EwkTx2.After each iteration, the decorrelation extraction and normalization are performed on the linear combination of independent components to prevent these components from converging to the same extreme point, which leads to the instability of the algorithm because of the zero-mean of the source signal. The formulas for this method are shown as follows:(8)wp+1k+1=wp+1k−Exwp+1kTx33Ewp+1kTx2,wp+1k+1=wp+1k+1−∑j=1pwp+1k+1Twjwj,wp+1k+1=wp+1k+1wp+1k+1Twp+1k+1.The separating matrixW is calculated by estimating all the components w with the method described above.


## 2.3. Fractal Dimension

Fractal theory, an active area of signal processing research, is a mathematical tool that quantitatively describes various complex signals in nature. The fractal dimension of the vibration signal may be calculated to measure the complexity and irregularity of the signal. In fractal theory, larger fractal dimensions indicate that the signal is complex. For vibration signal processing, noise often causes the vibration signal components to be complicated, nonlinear, and nonstationary. Thus, the fractal dimension of the vibration signal is large. Many experimental results show the fractal dimension (grid dimension) of abnormal signals is greater than 1.5. Hence, the source signals with a larger fractal dimension may be filtered out avoiding meaningful information loss. The noise may be suppressed using a wavelet filter on the source signal, which has smaller fractal dimension. The fractal dimension contains the box dimension, correlation dimension, Hausdorff dimension, and grid dimension. For one-dimensional signals, the box dimension is equivalent to the grid dimension, which simplifies the computer processing. For an signal channel signalx|xi,i=1,2,…,N, the grid dimension is defined as(9)D=−lnMjlnΔt,where Δt is the sampling interval, Mj=∑xi+1−xi/Δt, and x represents the closest integer that is greater than x. First, the space of discrete signals is divided into a multidimensional grid, and the grid dimension is calculated by counting the grid points. The Euclidean space of discrete signals RN is divided as small as possible into grids with the area of Δ. Points’ count NΔK of grid space of the signal x|xi,i=1,2,…,N is calculated. Different areas NΔK are used to calculate the NΔK point, yielding a point set xk,yk:(10)xk=−lnΔK,yk=−lnNΔK.The slope of the line formed by the set of pointsxk,yk is the grid dimension of the signal. The line is not straight, so the least-squares method is often used to fit a line for practical applications. Finally, the grid dimension D may be obtained.


## 3. Experiments and Results

This section contains 3 groups of experiments. The first uses general sin wave to show the effectiveness of the SCICA with phase space reconstruction. Two cases are given respective correspondence to the condition with or without noise. The section and third experiments apply the SCICA with FD (SCICA-FD) in noise reduction. Experiment 2 combines SCICA-FD with wavelet denoising, the basic knowledge of wavelet is given, and the experiment result is discussed. Similar to experiment 2, experiment 3 applies smoothing filtering to SCICA-FD to show the proposed method can be applied to different denoising methods.


### 3.1. Validation of SCICA Based on Phase Space Reconstruction

Case 1.
Using general sine waves and performing phase space reconstruction-based SCICA to verify the effect, in the absence of noise, signal 1 and signal 2 and their mixed-signal are shown in Figure2(a). SCICA firstly separates the mixed-signal with time delay calculated as 3 and the embedding dimension 7, which means 7 ICA components are obtained. Because too many signals are obtained by SCICA, principal component analysis (PCA) is applied to reduce the dimension to avoid choosing the signal from the 7 ICA output signals. It can be seen in Figure 2(b) that the three principal components with the largest contribution is added to more than 95%, so PCA is used to reduce the 7 dimension signal matrix to three dimensions. From the result shown in Figure 2(c), it can be seen that the phase space reconstruction-based SCICA can effectively separate the two sine signals from the combined signal.Figure 2
Noise-free condition. (a) Original and combined signal, (b) PCA variance, and (c) separation results.
(a)(b)(c)Case 2.
Similar to Case1, signal 1, signal 2, a random noise, and their mix signal are shown in Figure 3(a). The mixed signal is separated by SCICA with time delay which equals to 2 and embedding dimension which equals to 9 which means there is 9 ICA components. Again, PCA is used to reduce the dimension of the ICA output signal matrix. From Figure 3(b), the matrix should be reduced to 4 dimensions. The result is shown in Figure 3(c). It can be seen that compare to that of Case 1 (without noise), the output waveform has some fluctuation, but the signal obtained is still close enough to the original one.Figure 3
Noise added condition. (a) Original and combined signal, (b) PCA variance (with noise), and (c) separation results.
(a)(b)(c)


### 3.2. Application to Wavelet Denoising

From a signal perspective, wavelet denoising is a signal filtering method. Although wavelet denoising can be viewed as a low pass filter, the signal features can be preserved successfully after denoising, so it is better than the traditional low-pass filter.


#### 3.2.1. Necessary Knowledge of Wavelet Denoising

The model of a nonstationary, noisy, and one-dimensional signal is represented as follows:(11)si=fi+σei,i=0,1,…,n−1,where fi is the true signal, ei is the noise, the noise property is unknown, and si is the noisy signal series. To improve the signal-to-noise ratio (SNR) of the signal, it is necessary to reduce the noise. Taking the simplest noise model as an example, ei is considered as Gauss white noise N0,1. In the actual acquisition of vibration signal, the useful signal is usually low or some relatively stable signal and noise signal is normally high-frequency signal, and the noise reduction process can be processed according to the following methods. First, decompose the noise signal by wavelet transform, such as the three layers’ decomposition, and the decomposition process is shown in Figure 4; the decomposed signals include two parts: approximate signal (CA3) and the detail signal (CD3, CD2, and CD1).Figure 4
Sketch map of wave decomposition.The approximate signal represents the decomposed low-frequency signal. Detail signals stand for high-frequency part. Because the noise belongs to high frequency signals, it is usually included in each detail signal level. In the process, the approximate signal of wavelet transform is reserved, and the detail signal is processed in the form of a threshold value. After the above process, the signal is reconstructed by wavelet transform, thus obtains the denoised signal in which the useful part is well preserved.The noise reduction process of the vibration signals can be divided into three steps:Step 1 (vibration signal wavelet decomposition): select a wavelet, determine the wavelet decomposition layersN, and then take N layer wavelet decomposition to signal s. Also, the appropriate choice of wavelet base has a significant impact on the noise reduction effect. A suitable wavelet basis can make the signal energy be distributed on a small number of substrates.Step 2: after the wavelet decomposition, the high frequency coefficient’s threshold value is quantized. The high-frequency coefficient of each layer of the first layer to the N layer is selected, and a threshold value is selected to process the soft threshold quantization.Step 3 (one-dimensional wavelet reconstruction): according to the low-frequency coefficients of the fourth layer of wavelet decomposition and the high-frequency coefficients of the first layer to theN layer, the wavelet reconstruction of the vibration signal is carried out. In the process of noise reduction, the threshold value of the high-frequency coefficient is selected, and the threshold and the soft threshold are given in two forms of Donoho [23]. Usually, the hard threshold denoising will cause the loss of available information, so the engineering signal analysis is mostly based on the soft threshold denoising technique. The method described in this paper is also a soft thresholding.


#### 3.2.2. SCICA-FD with Wavelet Denoising

Using Case Western Reserve University’s bearing data, Figure5(a) shows the vibration signal in time domain and its spectrum is shown in Figure 5(b). As shown in Figure 5(b), the vibration signal time domain contains large background noise, and the fault information is submerged in the noise.Figure 5
Vibration signal and its spectrum. (a) Time-domain waveform of vibration signal and (b) spectrum of the vibration signal.
(a)(b)Using the proposed method, the vibration data is separated by SCICA. The fractal dimension of each separated signal is shown in Figure6. As mentioned before, this paper chooses signals with fractal dimension close to 1.5, which is number 4 to 13 in Figure 6.Figure 6
Fractal dimension of each ICA output signal.Then, the same wavelet threshold is applied to the selected signal. The spectrum of the restored signal is shown in Figure7(b). Performing only the wavelet denoising, the result is shown in Figure 7(a). Comparing to Figure 7(a), Figure 7(a) can tell that the proposed method reduces the noise and reserves the high-frequency part. Some of the high-frequency parts are a little weak.Figure 7
The spectrum of vibration signal (a) after wavelet denoising and (b) after SCICA-wavelet denoising.
(a)(b)


### 3.3. Application to Smooth Filtering

Wavelet denoising is applied in experiment 2; as mentioned, wavelet denoising can be considered a lowpass filter but retains more useful frequencies. The smoothing filter used here is called Savitzky–Golay filter which is proposed in 1964 [24]. Savitzky–Golay filter is a low-pass filter based on polynomial fitting, compared to wavelet denoising, which is a simpler and more straightforward way. This paper shows in this part that our proposed SCICA-FD can be applied with even simpler methods.


#### 3.3.1. Savitzky–Golay Filter

Let a set of data inxn be xi, i=M,…,0,…,M, and a polynomial of order p is constructed:(12)fi=a0+a1i+⋯+apip=∑k=0pakik,p≤2M.The polynomial is used to fit the data. Due to the existence of fitting errors, the sum of squares of the total error is(13)E=∑−MMfi−xi2=∑−MM∑k=0pakik−xi2.To minimizeE, the derivative of E to each coefficient is set to 0, ∂E/∂ar=0,r=0,1,2,…,p, that is,(14)∂E∂ar=∂∑−MM∑k=0pakik−xi2∂ar=0,and thus,(15)∑k=0pak∑i=−MMik+r=∑i=−MMxiir.Let(16)Fr=∑i=−mMxiir,Sk+r=∑i=−MMik+r,and then, equation (15) can be written as(17)Fr=∑k=0PakSk+r.Given the number of pointM needed to fit, polynomial order P, and the data xi, i=M,…,0,…,M, to be fitted, the coefficients a0,a1,…,ap can be calculated by equation (17) and thus get the polynomial fi.In the actual calculation, the best fitting polynomial to center pointx0 can be obtained as long as the coefficient a0 is obtained, and assume M=2andP=2, and applying equation (17), we obtain(18)a0=−2x−2+12x−1+17x0+12x1−3x235,and consider a0 as a filter factor or template:(19)hn=−3,12,17,12,−335.Moving the template on dataxn, central point f0 can be calculated according to equation (19) and thus achieve fitting of xn. hn is a symmetric FIR filter, which is n=O symmetry, so it has the characteristics of zero phase, and its coefficients sum equals 1, and its frequency characteristic is as follows:(20)Hejw=17+24cosw−6cos2w35.The smooth filter used in this paper set isM=5 and p=3.


#### 3.3.2. SCICA-FD with Smooth Filtering

Using the same data, recall that the time domain signal is plotted in Figure5(a) and its spectrum is plotted in Figure 5(b). Apart from the denoising method, nothing changes, and compare d to experiment 2. That is, the SCICA-FD output is the same, and the fractal dimension of each separated signal is the same as that in Figure 6.The spectrum of the restored signal is shown in Figure8(b). Applying only smooth filtering on the vibration signal, the filtering result is shown in Figure 8(a). After smooth filtering, the spectrum is similar to that of wavelet denoising as both algorithms have the characteristic of lowpass filter. First, carry out the proposed SCICA-FD than filtering only part of the component signals, and the result is shown in Figure 8(b). Again, both the 1000 Hz and 2000 Hz high frequency are retained and frequency elsewhere is depressed. From Figure 8(b), it is clear that the spectrum of vibration signal using SCICA-FD with smooth filtering has obvious fault shock characteristics at 1000 Hz and 2000 Hz.Figure 8
The spectrum of vibration signal (a) using smooth filtering and (b) using SCICA-FD with smooth filtering.
(a)(b)


## 3.1. Validation of SCICA Based on Phase Space Reconstruction

Case 1.
Using general sine waves and performing phase space reconstruction-based SCICA to verify the effect, in the absence of noise, signal 1 and signal 2 and their mixed-signal are shown in Figure2(a). SCICA firstly separates the mixed-signal with time delay calculated as 3 and the embedding dimension 7, which means 7 ICA components are obtained. Because too many signals are obtained by SCICA, principal component analysis (PCA) is applied to reduce the dimension to avoid choosing the signal from the 7 ICA output signals. It can be seen in Figure 2(b) that the three principal components with the largest contribution is added to more than 95%, so PCA is used to reduce the 7 dimension signal matrix to three dimensions. From the result shown in Figure 2(c), it can be seen that the phase space reconstruction-based SCICA can effectively separate the two sine signals from the combined signal.Figure 2
Noise-free condition. (a) Original and combined signal, (b) PCA variance, and (c) separation results.
(a)(b)(c)Case 2.
Similar to Case1, signal 1, signal 2, a random noise, and their mix signal are shown in Figure 3(a). The mixed signal is separated by SCICA with time delay which equals to 2 and embedding dimension which equals to 9 which means there is 9 ICA components. Again, PCA is used to reduce the dimension of the ICA output signal matrix. From Figure 3(b), the matrix should be reduced to 4 dimensions. The result is shown in Figure 3(c). It can be seen that compare to that of Case 1 (without noise), the output waveform has some fluctuation, but the signal obtained is still close enough to the original one.Figure 3
Noise added condition. (a) Original and combined signal, (b) PCA variance (with noise), and (c) separation results.
(a)(b)(c)


## 3.2. Application to Wavelet Denoising

From a signal perspective, wavelet denoising is a signal filtering method. Although wavelet denoising can be viewed as a low pass filter, the signal features can be preserved successfully after denoising, so it is better than the traditional low-pass filter.


### 3.2.1. Necessary Knowledge of Wavelet Denoising

The model of a nonstationary, noisy, and one-dimensional signal is represented as follows:(11)si=fi+σei,i=0,1,…,n−1,where fi is the true signal, ei is the noise, the noise property is unknown, and si is the noisy signal series. To improve the signal-to-noise ratio (SNR) of the signal, it is necessary to reduce the noise. Taking the simplest noise model as an example, ei is considered as Gauss white noise N0,1. In the actual acquisition of vibration signal, the useful signal is usually low or some relatively stable signal and noise signal is normally high-frequency signal, and the noise reduction process can be processed according to the following methods. First, decompose the noise signal by wavelet transform, such as the three layers’ decomposition, and the decomposition process is shown in Figure 4; the decomposed signals include two parts: approximate signal (CA3) and the detail signal (CD3, CD2, and CD1).Figure 4
Sketch map of wave decomposition.The approximate signal represents the decomposed low-frequency signal. Detail signals stand for high-frequency part. Because the noise belongs to high frequency signals, it is usually included in each detail signal level. In the process, the approximate signal of wavelet transform is reserved, and the detail signal is processed in the form of a threshold value. After the above process, the signal is reconstructed by wavelet transform, thus obtains the denoised signal in which the useful part is well preserved.The noise reduction process of the vibration signals can be divided into three steps:Step 1 (vibration signal wavelet decomposition): select a wavelet, determine the wavelet decomposition layersN, and then take N layer wavelet decomposition to signal s. Also, the appropriate choice of wavelet base has a significant impact on the noise reduction effect. A suitable wavelet basis can make the signal energy be distributed on a small number of substrates.Step 2: after the wavelet decomposition, the high frequency coefficient’s threshold value is quantized. The high-frequency coefficient of each layer of the first layer to the N layer is selected, and a threshold value is selected to process the soft threshold quantization.Step 3 (one-dimensional wavelet reconstruction): according to the low-frequency coefficients of the fourth layer of wavelet decomposition and the high-frequency coefficients of the first layer to theN layer, the wavelet reconstruction of the vibration signal is carried out. In the process of noise reduction, the threshold value of the high-frequency coefficient is selected, and the threshold and the soft threshold are given in two forms of Donoho [23]. Usually, the hard threshold denoising will cause the loss of available information, so the engineering signal analysis is mostly based on the soft threshold denoising technique. The method described in this paper is also a soft thresholding.


### 3.2.2. SCICA-FD with Wavelet Denoising

Using Case Western Reserve University’s bearing data, Figure5(a) shows the vibration signal in time domain and its spectrum is shown in Figure 5(b). As shown in Figure 5(b), the vibration signal time domain contains large background noise, and the fault information is submerged in the noise.Figure 5
Vibration signal and its spectrum. (a) Time-domain waveform of vibration signal and (b) spectrum of the vibration signal.
(a)(b)Using the proposed method, the vibration data is separated by SCICA. The fractal dimension of each separated signal is shown in Figure6. As mentioned before, this paper chooses signals with fractal dimension close to 1.5, which is number 4 to 13 in Figure 6.Figure 6
Fractal dimension of each ICA output signal.Then, the same wavelet threshold is applied to the selected signal. The spectrum of the restored signal is shown in Figure7(b). Performing only the wavelet denoising, the result is shown in Figure 7(a). Comparing to Figure 7(a), Figure 7(a) can tell that the proposed method reduces the noise and reserves the high-frequency part. Some of the high-frequency parts are a little weak.Figure 7
The spectrum of vibration signal (a) after wavelet denoising and (b) after SCICA-wavelet denoising.
(a)(b)


## 3.2.1. Necessary Knowledge of Wavelet Denoising

The model of a nonstationary, noisy, and one-dimensional signal is represented as follows:(11)si=fi+σei,i=0,1,…,n−1,where fi is the true signal, ei is the noise, the noise property is unknown, and si is the noisy signal series. To improve the signal-to-noise ratio (SNR) of the signal, it is necessary to reduce the noise. Taking the simplest noise model as an example, ei is considered as Gauss white noise N0,1. In the actual acquisition of vibration signal, the useful signal is usually low or some relatively stable signal and noise signal is normally high-frequency signal, and the noise reduction process can be processed according to the following methods. First, decompose the noise signal by wavelet transform, such as the three layers’ decomposition, and the decomposition process is shown in Figure 4; the decomposed signals include two parts: approximate signal (CA3) and the detail signal (CD3, CD2, and CD1).Figure 4
Sketch map of wave decomposition.The approximate signal represents the decomposed low-frequency signal. Detail signals stand for high-frequency part. Because the noise belongs to high frequency signals, it is usually included in each detail signal level. In the process, the approximate signal of wavelet transform is reserved, and the detail signal is processed in the form of a threshold value. After the above process, the signal is reconstructed by wavelet transform, thus obtains the denoised signal in which the useful part is well preserved.The noise reduction process of the vibration signals can be divided into three steps:Step 1 (vibration signal wavelet decomposition): select a wavelet, determine the wavelet decomposition layersN, and then take N layer wavelet decomposition to signal s. Also, the appropriate choice of wavelet base has a significant impact on the noise reduction effect. A suitable wavelet basis can make the signal energy be distributed on a small number of substrates.Step 2: after the wavelet decomposition, the high frequency coefficient’s threshold value is quantized. The high-frequency coefficient of each layer of the first layer to the N layer is selected, and a threshold value is selected to process the soft threshold quantization.Step 3 (one-dimensional wavelet reconstruction): according to the low-frequency coefficients of the fourth layer of wavelet decomposition and the high-frequency coefficients of the first layer to theN layer, the wavelet reconstruction of the vibration signal is carried out. In the process of noise reduction, the threshold value of the high-frequency coefficient is selected, and the threshold and the soft threshold are given in two forms of Donoho [23]. Usually, the hard threshold denoising will cause the loss of available information, so the engineering signal analysis is mostly based on the soft threshold denoising technique. The method described in this paper is also a soft thresholding.


## 3.2.2. SCICA-FD with Wavelet Denoising

Using Case Western Reserve University’s bearing data, Figure5(a) shows the vibration signal in time domain and its spectrum is shown in Figure 5(b). As shown in Figure 5(b), the vibration signal time domain contains large background noise, and the fault information is submerged in the noise.Figure 5
Vibration signal and its spectrum. (a) Time-domain waveform of vibration signal and (b) spectrum of the vibration signal.
(a)(b)Using the proposed method, the vibration data is separated by SCICA. The fractal dimension of each separated signal is shown in Figure6. As mentioned before, this paper chooses signals with fractal dimension close to 1.5, which is number 4 to 13 in Figure 6.Figure 6
Fractal dimension of each ICA output signal.Then, the same wavelet threshold is applied to the selected signal. The spectrum of the restored signal is shown in Figure7(b). Performing only the wavelet denoising, the result is shown in Figure 7(a). Comparing to Figure 7(a), Figure 7(a) can tell that the proposed method reduces the noise and reserves the high-frequency part. Some of the high-frequency parts are a little weak.Figure 7
The spectrum of vibration signal (a) after wavelet denoising and (b) after SCICA-wavelet denoising.
(a)(b)


## 3.3. Application to Smooth Filtering

Wavelet denoising is applied in experiment 2; as mentioned, wavelet denoising can be considered a lowpass filter but retains more useful frequencies. The smoothing filter used here is called Savitzky–Golay filter which is proposed in 1964 [24]. Savitzky–Golay filter is a low-pass filter based on polynomial fitting, compared to wavelet denoising, which is a simpler and more straightforward way. This paper shows in this part that our proposed SCICA-FD can be applied with even simpler methods.


### 3.3.1. Savitzky–Golay Filter

Let a set of data inxn be xi, i=M,…,0,…,M, and a polynomial of order p is constructed:(12)fi=a0+a1i+⋯+apip=∑k=0pakik,p≤2M.The polynomial is used to fit the data. Due to the existence of fitting errors, the sum of squares of the total error is(13)E=∑−MMfi−xi2=∑−MM∑k=0pakik−xi2.To minimizeE, the derivative of E to each coefficient is set to 0, ∂E/∂ar=0,r=0,1,2,…,p, that is,(14)∂E∂ar=∂∑−MM∑k=0pakik−xi2∂ar=0,and thus,(15)∑k=0pak∑i=−MMik+r=∑i=−MMxiir.Let(16)Fr=∑i=−mMxiir,Sk+r=∑i=−MMik+r,and then, equation (15) can be written as(17)Fr=∑k=0PakSk+r.Given the number of pointM needed to fit, polynomial order P, and the data xi, i=M,…,0,…,M, to be fitted, the coefficients a0,a1,…,ap can be calculated by equation (17) and thus get the polynomial fi.In the actual calculation, the best fitting polynomial to center pointx0 can be obtained as long as the coefficient a0 is obtained, and assume M=2andP=2, and applying equation (17), we obtain(18)a0=−2x−2+12x−1+17x0+12x1−3x235,and consider a0 as a filter factor or template:(19)hn=−3,12,17,12,−335.Moving the template on dataxn, central point f0 can be calculated according to equation (19) and thus achieve fitting of xn. hn is a symmetric FIR filter, which is n=O symmetry, so it has the characteristics of zero phase, and its coefficients sum equals 1, and its frequency characteristic is as follows:(20)Hejw=17+24cosw−6cos2w35.The smooth filter used in this paper set isM=5 and p=3.


### 3.3.2. SCICA-FD with Smooth Filtering

Using the same data, recall that the time domain signal is plotted in Figure5(a) and its spectrum is plotted in Figure 5(b). Apart from the denoising method, nothing changes, and compare d to experiment 2. That is, the SCICA-FD output is the same, and the fractal dimension of each separated signal is the same as that in Figure 6.The spectrum of the restored signal is shown in Figure8(b). Applying only smooth filtering on the vibration signal, the filtering result is shown in Figure 8(a). After smooth filtering, the spectrum is similar to that of wavelet denoising as both algorithms have the characteristic of lowpass filter. First, carry out the proposed SCICA-FD than filtering only part of the component signals, and the result is shown in Figure 8(b). Again, both the 1000 Hz and 2000 Hz high frequency are retained and frequency elsewhere is depressed. From Figure 8(b), it is clear that the spectrum of vibration signal using SCICA-FD with smooth filtering has obvious fault shock characteristics at 1000 Hz and 2000 Hz.Figure 8
The spectrum of vibration signal (a) using smooth filtering and (b) using SCICA-FD with smooth filtering.
(a)(b)


## 3.3.1. Savitzky–Golay Filter

Let a set of data inxn be xi, i=M,…,0,…,M, and a polynomial of order p is constructed:(12)fi=a0+a1i+⋯+apip=∑k=0pakik,p≤2M.The polynomial is used to fit the data. Due to the existence of fitting errors, the sum of squares of the total error is(13)E=∑−MMfi−xi2=∑−MM∑k=0pakik−xi2.To minimizeE, the derivative of E to each coefficient is set to 0, ∂E/∂ar=0,r=0,1,2,…,p, that is,(14)∂E∂ar=∂∑−MM∑k=0pakik−xi2∂ar=0,and thus,(15)∑k=0pak∑i=−MMik+r=∑i=−MMxiir.Let(16)Fr=∑i=−mMxiir,Sk+r=∑i=−MMik+r,and then, equation (15) can be written as(17)Fr=∑k=0PakSk+r.Given the number of pointM needed to fit, polynomial order P, and the data xi, i=M,…,0,…,M, to be fitted, the coefficients a0,a1,…,ap can be calculated by equation (17) and thus get the polynomial fi.In the actual calculation, the best fitting polynomial to center pointx0 can be obtained as long as the coefficient a0 is obtained, and assume M=2andP=2, and applying equation (17), we obtain(18)a0=−2x−2+12x−1+17x0+12x1−3x235,and consider a0 as a filter factor or template:(19)hn=−3,12,17,12,−335.Moving the template on dataxn, central point f0 can be calculated according to equation (19) and thus achieve fitting of xn. hn is a symmetric FIR filter, which is n=O symmetry, so it has the characteristics of zero phase, and its coefficients sum equals 1, and its frequency characteristic is as follows:(20)Hejw=17+24cosw−6cos2w35.The smooth filter used in this paper set isM=5 and p=3.


## 3.3.2. SCICA-FD with Smooth Filtering

Using the same data, recall that the time domain signal is plotted in Figure5(a) and its spectrum is plotted in Figure 5(b). Apart from the denoising method, nothing changes, and compare d to experiment 2. That is, the SCICA-FD output is the same, and the fractal dimension of each separated signal is the same as that in Figure 6.The spectrum of the restored signal is shown in Figure8(b). Applying only smooth filtering on the vibration signal, the filtering result is shown in Figure 8(a). After smooth filtering, the spectrum is similar to that of wavelet denoising as both algorithms have the characteristic of lowpass filter. First, carry out the proposed SCICA-FD than filtering only part of the component signals, and the result is shown in Figure 8(b). Again, both the 1000 Hz and 2000 Hz high frequency are retained and frequency elsewhere is depressed. From Figure 8(b), it is clear that the spectrum of vibration signal using SCICA-FD with smooth filtering has obvious fault shock characteristics at 1000 Hz and 2000 Hz.Figure 8
The spectrum of vibration signal (a) using smooth filtering and (b) using SCICA-FD with smooth filtering.
(a)(b)


## 4. Conclusion

The experimental results demonstrate that vibration signal components are complex. The signals exhibit both nonlinear and nonstationary characteristics. The proposed vibration noise reduction method combines the advantages of fractal analysis, blind source separation, and different denoising algorithms. Using the fractal dimension constraint with the SCICA, the denoising process avoids filtering useful information. Thus, the proposed method can reduce noise while ensuring that useful information is retained. The signal after noise reduction is smoother, and the burr is reduced. Although the amplitude of the noise was significantly attenuated because of the ICA algorithm, this change does not affect the analysis of the frequency domain components because changes in amplitude do not affect the time-frequency distribution. The proposed method effectively reduces noise and may have many theoretical and practical applications. Future work may show that the proposed method may be extended to other areas of research in which noise reduction is required.



---

*Source: 5583587-2021-05-17.xml*