Signal Denoising Based on Duffing Oscillators System

Signal denoising is an important aspect of signal processing. It is a meaningful direction to solve this problem by using chaotic oscillator. As a kind of chaotic oscillator, Duffing oscillator is often used in the field of periodic signal or pulse signal detection, but it has not been used in signal denoising. In this paper, a new Duffing oscillators system for signal denoising is proposed. The system utilizes the coupling between linear and nonlinear restoring forces of the oscillators to achieve complete synchronization of the system. Because the external signal will cause a large deviation of the oscillator’s trajectory, while the Gauss white noise will only produce a small deviation. Based on this phenomenon, signal denoising can be realized. The simulation results show that, compared with the existing technology, our method is effective for the denoising of complex irregular signal, and can recover the mutation points of signal well. This method has less computational complexity and is easy to implement parallel computation. It can be used in real-time signal processing.


I. INTRODUCTION
Signal denoising is an important aspect of signal processing. For different types of signals and noises, finding the best denoising method has always been one of the main problems in signal processing. According to the characteristics of different signals and noises, many basic denoising methods have been proposed, such as Fourier transform, wavelet transform [1], [2], singular value decomposition (SVD) [3], [4], empirical mode decomposition (EMD) [5], [6], blind source separation [7], sparse decomposition [8], stochastic resonance [9], and etc. These basic denoising methods have their own limitations, and cannot get ideal signal denoising result generally. So a research hotspot is to combine them to get better performance nowadays [10]- [13]. On the other hand, introducing deep learning into signal processing [14]- [16] is another hot topic at present. However, deep learning must have some prior knowledge, such as speech denoising and enhancement. Therefore, without prior knowledge, deep learning cannot be used.
The associate editor coordinating the review of this manuscript and approving it for publication was Jun Shi .
So far, there is no new basic signal denoising method proposed. In this paper, we will try to propose a new basic signal denoising method using chaotic oscillator.
Chaotic oscillator has been applied in sinusoidal signal detection [17], [18], aperiodic pulse signal detection and other fields, and has its unique advantages compared with other signal detection methods. Duffing oscillator is most widely used in signal processing because it has a good ability to suppress white noise and colored noise [19]. Yang et al. [20] pointed out that multi Duffing oscillators system can be applied to detect aperiodic signal, and proposed a double Duffing oscillators system based on the coupling of linear restoring forces. After that a bidirectional ring coupled Duffing oscillators system [21] is proposed to improve the detection performance by the phenomenon of transient synchronization mutation between the oscillators. Recently, a Duffing oscillators system [22] based on the coupling of linear restoring forces and damping forces is proposed to improve the performance of signal detection.
All of the above literatures applied Duffing oscillator to signal detection, and the purpose of signal detection is not to extract waveform from signal.
Recently, a new filtering algorithm based on cascade of driven chaotic oscillators is proposed in [23]. That paper pointed out that chaotic system can work as filter, and proposed a non-linear filter construction as a cascade of driven Rössler oscillators. In that paper, the chaotic system is used as a filter to process the digital chaotic signal, which can improve the SNR. The purpose of that paper is more in the field of digital signal detection than in the field of signal denoising. The application of Duffing oscillator in general signal denoising has not been discussed so far.
In order to extract complex and irregular signal from noise, a denoising method based on a coupled Duffing oscillators system is proposed, which is characterized by the coupling of linear and nonlinear restoring forces. In order to preserve the waveform features as much as possible, the selection of system parameters is studied. Simulation results show that, compared with existing technology, Duffing oscillator denosing can better restore the mutation parts and overall trend of complex irregular signal at lower SNR environment, but there are some small high-frequency parasitic oscillations in the denoised signal. Therefore, after the Duffing oscillator denoising, wavelet denoising is used to suppress the small parasitic oscillation of the signal waveform.
In addition, the calculation cost of this method is low, and it is an easy signal denoising method.

II. CONSTRUCTION OF DUFFING OSCILLATORS SYSTEM
In the application of Duffing oscillator in signal processing, a single oscillator can only use its phase transition to detect the sinusoidal signal, while the coupled oscillators can use synchronous burst between oscillators to detect the pulse signal [21].
In order to perform signal denoising, a double Holmes-Duffing oscillators system based on the coupling of linear and nonlinear restoring forces is proposed. It can be expressed as where ξ is the damping ratio, f denotes the amplitude of periodic driving force, α is the coupling coefficient of the linear restoring forces, β is the coupling coefficient of the nonlinear restoring forces and s(t) represents signal to be denoised. The system output can be obtained through (x 1 − x 2 ).
As shown in Eq. (1), the angular frequency of the periodic driving force is 1, so it can only be adopted to handle digital signal with low sampling frequency. In order to handle signal with any sampling frequency, the Eq. (1) can be transformed to a generalized time scale form as where t is replaced by ωτ firstly, then τ is substituted by t. Time scale transformation does not change the characteristics and performance of the system. It should be noted that only Eq. (1) will be discussed in this paper. Runge-Kutta (RK) method is a commonly used numerical solution for nonlinear system. In this paper, the fixed-step fourth-order Runge-Kutta (RK4) method is used for basic numerical simulation.
In this paper, without special statement, the system solving parameters are set as typical values: ξ is 8, α is 3, β is 0.5, f is 0.2, and the total simulation time and recursive step size of RK4 are 2000s and 0.2s respectively.

A. THE PRINCIPLE OF SIGNAL DENOISING
The proposed system is a fully synchronous system with strong coupling strength, and the two oscillators are always in the same potential well after synchronization. When external signal with small amplitude is added into a certain oscillator, the oscillators will be out of synchronization accordingly, but will not jump to different potential wells. The phenomenon can be described as ''out of synchronization in one potential well''. While some weakly coupled systems will jump to different potential wells. Because the response of the excited oscillator and the coupled oscillator to external signal is different, the external signal can be recovered by comparing the different oscillation orbits of the two oscillators.
When noise with low intensity is added into the system, the energy at the resonance frequency of the system is weak due to the wide-band characteristics of noise. Noise will not disturb the phase orbit of the oscillators greatly, but will leave some small rough marks on the phase orbits. Therefore, the system has anti-noise performance.
Above is the principle of signal denoising. The following numerical simulation illustrates this conception.
A square wave is added to the first oscillator of the system, which lasts from 1000s to 1020s and has an amplitude of 0.5. The system solving parameters are selected as typical values. Fig. 1(a) shows the waveforms of the state variables x 1 and x 2 . It can be found that the oscillation trajectories of the two oscillators are completely synchronized without external interference, and the oscillation period is 2π which is equal to the periodic driving force of the system. At the position of the external square wave, the trajectory of the first oscillator has a large deviation, and the trajectory of the second oscillator VOLUME 8, 2020 also migrates almost at the same time. The deviation of the two oscillators is in the same direction. The waveform of the external signal can be recovered by the difference of the state variable, e.g. x 1 − x 2 , as shown in Fig. 1(b). In time domain, the system can respond to external signal with 2π width.
If β is set to 1 and other conditions remain unchanged, the calculation results are shown in Fig. 2.
It can be seen that the orbit deviation of the first oscillator in Fig. 2(a) is larger than that in Fig. 1(a), while the orbit deviation direction of the second oscillator has changed. The main restraint of the above two cases is that the total potential energy of the system remains unchanged. Compared with Fig. 1(b) and Fig. 2(b), the larger the coupling strength, the larger the output amplitude, that is, the higher the output signal-to-noise ratio (SNR) of the system. But the output waveform in Fig. 2(b) deviates greatly from the original signal. In the aspect of signal detection, it is more inclined to improve the output SNR of the system. However, this paper pays more attention to the waveform distortion after signal denoising, and subsequent discussions focus on waveform distortion.

B. THE INFLUENCE OF PARAMETERS ON SYSTEM STATE
The initial values of oscillator variables, damping coefficient ξ , coupling coefficient α and β, amplitude of periodic driving force f and recursive step size of RK4 may affect the state of the system. The system parameters must be carefully selected to enable the system to work in the required state.
As an example, the influence of different β values on the state of the system is analyzed, where other system parameters are set to typical values.
When β is less than 1.336, the system is in synchronous state, and the two oscillators oscillate periodically around focal point (+1, 0). While β is in the range of 1.337∼1.9732, the two oscillators are out of synchronization and the orbits of the oscillators are bifurcated. Each oscillator has two focal points, one of which is (+1, 0). When β is greater than 1.9733, the solution will diverge. Fig. 3(a) is the phase diagram of the two oscillators when β is 0.5, from which we can see that the two oscillators enter synchronization from different initial values. Fig. 3(b) is the phase diagram of the two oscillators when β equals 1.5, and it can be found that the two oscillators have two focal points and are asynchronous.

C. THE INFLUENCE OF PARAMETERS ON SYSTEM PERFORMANCE
The following describes the principles for system parameters selection.

1) INITIAL VALUES OF SYSTEM STATE VARIABLES
Because of high coupling degree of the system, small initial values of system will not affect the response of the system to external signal. Therefore, the initial values of system can be chosen as any smaller values. In order to synchronize at the beginning of the system, it is recommended to set the initial values of the system variable to zeroes.

2) SOLVING PARAMETERS OF RK4
The total computation time is suggested to be about 2000s. If it is too long, the solution will diverge.
The recursive step size is best chosen to be about 0.2s. If it is too large or too small, the solution will diverge. The system solution is stable when that value is approximate to the period of the periodic driving force of the oscillators. If the solution step is small, the number of sampling points of the external input signal will be insufficient. That will lead to a large deviation in the three-stage slope estimation of RK4 calculation, and the error will be very large after iteration accumulation, which will lead to divergence of system solution. On the other hand, if the solution step is large, the periodic motion of the system will not be sampled enough. The estimation deviation of the three-stage slope in RK4 recursive calculation will also accumulate to a very large extent, and it will lead to divergence of system solution.

3) ξ AND f
When ξ is small, the oscillation amplitude of the system is large and the synchronization process is long because of the small damping force. Under the external excitation, the output of the system has more parasitic oscillations, which is undesirable for the recovery of original signal waveform. When ξ is larger, the system will be out of synchronization. ξ is recommended as 8, which corresponds to strong damping force.
When the damping force is large, the value of f has little effect on the performance of the system. The output waveform is smoother when the value of f is smaller. The value of f is suggested to be 0.2. When ξ and f are set to 8 and 0.2 respectively, the oscillators are in the limit circle state, and the system has the best performance in that state. In fact, the coupling system can achieve synchronization and signal denoising in limit cycle state, large period state and even chaos state. However, the denoising effect is not very good in large periodic state and chaotic state.

4) α AND β
When the value of α is larger or smaller, the numerical solution will diverge. When the value of β is smaller, the oscillators are asynchronous. When the value of β is larger, the numerical solution will diverge. Typical values of α and β are suggested to be 3 and 0.5 respectively.
For the selection of the above parameters, we need an objective evaluation method. Hence Pearson correlation coefficient is chosen to define the similarity between the system output signal and original signal, which can be expressed as where x i and y i are discrete input and output signal sequences respectively,x andȳ are the mean value of x i and y i respectively, r is the correlation coefficient of x i and y i with the range of −1∼1. The larger r, the higher the positive correlation. This expression can be found in the formula (1) of reference [24]. Taking r versus ξ as an example, the selection method of system parameters is illustrated. In the simulation, other parameters take typical values except ξ . A square wave of 40s duration and 0.5 amplitude with a SNR of −10dB is input into the system. ξ increases from 0.3 to 10.3, and the correlation coefficient of the output and input signal is shown in Fig. 4. We can find that when ξ is in the range of about 7∼10, r is at the extreme value. Similarly, other parameters can be optimized by the same method.

D. FREQUENCY-DOMAIN CHARACTERISTICS OF THE SYSTEM
A series of single-frequency signals with the same amplitude and different frequencies are input into the system, and the amplitudes of output signals are recorded. Then the frequency response characteristics of the system can be obtained by fitting those recorded amplitudes. Through such experiments, the normalized amplitude-frequency response function of the system can be expressed as a Gauss function: According to the equation, digital filters can also be designed. Fig. 5 is the processing result of a white Gaussian noise with 10000 points by the oscillators system. Fig. 5(a) is the time domain waveform of the noise, and Fig. 5(b) is the output of the system. Fig. 5(c) is the normalized spectrum of the output signal. It can be seen from Fig. 5(c) that the response of the system to white Gauss noise is similar to that of a Gauss low pass filter. Fig. 5(d) is the normalized amplitude-frequency response curve of the system represented by Eq. (4). Comparing Fig. 5(c) and Fig. 5(d), it can be found that the frequency response of the system to Gauss white noise accords with the description of Eq. (4).

III. PERFORMANCE ANALYSIS OF SIGNAL DENOISING A. COMPARISON WITH WAVELET DENOISING
Due to its good time-frequency characteristics, wavelet has been widely used in the field of signal denoising. However, the performance of simple wavelet denoising method is not good, and it must be combined with other methods to get satisfactory results.
In MATLAB's official wavelet toolbox, there is an example named ''multivariate wavelet denoising'', and that VOLUME 8, 2020  [25]. That method is one of the best denoising methods for complex and irregular signals.
In this part, the proposed oscillators system is compared with that wavelet method. The third signal in that MAT-LAB's example is taken out for simulation, which is shown in Fig. 6(a). Fig. 6(b) is the signal contaminated with white Gaussian noise at a SNR of 12dB, and Fig. 6(c) is the denoising result of the wavelet method. The denoising result of the oscillators system is shown in Fig. 6(d). The following simulation results of the proposed method are all linearly magnified by five times.
Comparing Fig. 6(c) and Fig. 6(d), it can be found that the result of the proposed method is similar to that of the wavelet method, and the result of wavelet is smoother because it removes the high frequency part of signal. As shown in the amplification part in Fig. 6(d), the denoised waveform of the proposed system has small ripples, because it does not completely eliminate the high frequency part of noise.
It should be noted that the SNR of the paper is relatively high because existing technologies cannot extract complex irregular waveform from strong noise at present.
The fourth signal in that MATLAB's example is taken out for further simulation, which is shown in Fig. 7(a). Fig. 7(b) is the signal contaminated with white Gaussian noise at a SNR of 3dB, and Fig. 7(c) is the denoising result of the wavelet method. The denoising result of the proposed system is shown in Fig. 7(d).
The wavelet method can correctly denoise the signal with only about 50% probability under that SNR, and Fig. 7(c) shows a failure case. As can be seen from Fig. 7(d), the proposed method can restore more details of the original signal at lower SNR.  Next, a more complex irregular signal with 1200 sampling points is selected, which is shown in Fig. 8(a). Fig. 8(b) is the signal contaminated with white Gaussian noise at a SNR of 5dB, and Fig. 8(c) is the denoising result of the wavelet method. The denoising result of the proposed system is shown in Fig. 8(d).
Comparing Fig. 8(c) and Fig. 8(d), it can be found that the proposed method can better preserve the mutation point and trend of the original signal, while the wavelet method may leave out some signal details because of its spectrum truncation operation.
Furthermore, another complex and irregular signal with 3000 sampling points is used for further experiments, which is shown in Fig. 9(a). Fig. 9(b) is the signal contaminated with white Gaussian noise at a SNR of 10dB, and Fig. 9(c) is the denoising result of the wavelet method. The denoising result of the proposed system is shown in Fig. 9(d).
Comparing Fig. 9(c) and Fig. 9(d), it can be found that the proposed method is a little more accurate to the detail restoration of the test signal.
A signal like noise with 2000 sampling points is also processed, which is shown in Fig. 10(a). Fig. 10(b) is the signal contaminated with white Gaussian noise at a SNR of 6dB, and Fig. 10(c) is the denoising result of the wavelet  method. The denoising result of the proposed system is shown in Fig. 10(d).
Comparing Fig. 10(c) and Fig. 10(d), it can be found that the denoising result of the proposed method is similar to that of the wavelet method. And at the circle position of the figure, the proposed method is a little more accurate than the wavelet method.
In addition, because the chaotic signal is a kind of wideband signal, and its spectrum is close to noise spectrum. It is also a problem to be studied whether the chaotic signal can be effectively denoised after being disturbed by noise. We simulate the chaotic signal of Lorenz attractor with 2000 sampling points, as shown in Fig. 11(a). Fig. 11(b) is the signal contaminated with white Gaussian noise at a SNR of 0dB, and Fig. 11(c) is the denoising result of the wavelet method. The denoising result of the proposed system is shown in Fig. 11(d). It can be seen from the simulation results that the proposed method is also effective for denoising chaotic signal, but the effect of wavelet method is not good.
Finally, a stable deterministic chirp signal is used for analysis, which is expressed as where t is time.  When the signal lasts from 0s to 1s and the sampling frequency is 2000 Hz, the waveform of the signal is shown in Fig. 12(a). Fig. 12(b) is the signal contaminated with white Gaussian noise at a SNR of −6dB, and Fig. 12(c) is the denoising result of the wavelet method. The denoising result of the proposed system is shown in Fig. 12(d).
Through Fig. 12(c), it's easy to find out that the denoising result of the wavelet method cannot restore every signal period very well. Through Fig. 12(d), we can find that the denoising result of the proposed method is better than that of the wavelet method.
Through the above simulation examples, we can see that the proposed method can get better denoising results for both unstable and stable signals.
Generally speaking, the denoising result of the proposed method can recover the overall trend and mutations of irregular signal accurately with different SNR. On the other hand, because the proposed method retains more high-frequency components of the signal, the denoising of the signal with more low-frequency components is not as ideal as the wavelet method. Therefore, this method is suitable for complex irregular waveform denoising, but not for slowly varying waveform denoising. VOLUME 8, 2020 FIGURE 13. The EMD denoising result of the signal shown in Fig. 8(b).

FIGURE 14.
The EMD denoising result of the signal shown in Fig. 7(b).
Moreover, the proposed method does not need to adjust the system parameters artificially, and it is also a simple and operational signal denoising method.

B. COMPARISON WITH EMD DENOISING
EMD is an adaptive signal time-frequency processing method, which is especially suitable for the analysis of nonlinear and non-stationary signals. EMD decomposes the signal according to the time scale characteristics of the data itself, without setting any basis function in advance. This is essentially different from the wavelet decomposition based on the prior wavelet basis function.
Kopsinis and McLaughlin [26] proposed a wavelet threshold denoising method suitable for EMD mode. That paper improved energy estimation method of IMFs and wavelet threshold condition. That method has better performance than the traditional EMD denoising method.
Kopsinis' method is used to denoise the signal in Fig. 8(b), and the result is as shown in Fig. 13. From Fig. 13, it can be seen that the denoising result of EMD method is similar to that of wavelet method, and the mutation details of signals cannot be well preserved.
In order to investigate the anti-noise performance of EMD method, the signal in Fig. 7 (b) is denoised by kopsinis' method, and the result is shown in Fig. 14. From Fig. 14, it can be seen that EMD denoising still cannot attain good anti-noise performance.
That method uses the first intrinsic mode function (IMF1) of EMD decomposition to the energy of noise, and the estimated value is not completely accurate because partial energy of noise will leak into other modes. When the noise is strong, the estimation is more inaccurate. Therefore, the anti-noise performance of that method is general.

C. IMPLEMENTATION OF THE PROPOSED METHOD
The above simulation is based on RK4 method, which contains 130 × (N-1) real additions and 170 × (N-1) real multiplications, where N is the length of input signal. The typical It can be seen from the comparison that the wavelet method is faster than the proposed method, but the proposed method is also sufficient for real-time calculation. EMD method is the slowest, which is mainly due to the slow parameter search process of that algorithm.
Although RK4 is a common method to solve ordinary differential equations, its computational efficiency is not high. Semi-implicit method can be used to improve the speed of the algorithm. The semi-implicit numerical integration method is an effective tradeoff between the weakly stable explicit method and the expensive implicit ODE solver. It allows the construction of computationally cheap symmetric methods, and it is easy to construct parallel algorithms [27]. Recently it has been discussed for the solution of various chaotic systems [28]- [30].
We choose the semi-implicit D-method to solve differential Eq. (1), because it has superior performance in serial and parallel computer architectures [27].
First, we change Eq. (1) to the first order form: The semi-implicit extrapolation solver (CD method [27]) for Eq. (6) can be written as: + β x n+0.5 where h is the local step size, and it is one-half of the global step size H. Note that H is the same as the step size of RK4 discussed in section II C. The superscript n represents the value of the variable at the current time. The superscript n+0.5 indicates the value after half step. The superscript n+1 indicates the next value after a global step.
Eq. (7) is of 2 nd order accuracy, which is lower than RK4, but it does not affect the denoising result of the system. For example, the signal in Fig. 8(b) is calculated by RK4 and semi-implicit method, and the calculation results are shown in Fig. 15. It can be found that the output signals of the two methods are very similar, and the two methods are also very similar for other signals. So, the semi-implicit method can be used to solve the problem of coupled Duffing oscillator denoising.
In order to clarify the accuracy of RK4 and semi-implicit method, the synchronization error signal x 1 − x 2 of Eq. (1) without external signal is simulated. The synchronization error signal obtained by RK4 method is shown in Fig. 16(a), while the result obtained by semi-implicit method is shown in Fig. 16(b). It can be seen from the figure that when RK4 method is used to solve the problem, the error signal after synchronization is zero, which shows that the method has enough accuracy. When the semi-implicit method is used, the synchronization error signal has small oscillations with amplitude of about 1.25 × 10 −4 , which is also relatively low.
The CD method reduces the computation by about half compared with the explicit RK4 method. Speed is the advantage of semi-implicit method. In addition, the input signal  can be divided into several segments for parallel processing, which improves the efficiency of operation.
The proposed method has the advantages of low computational complexity and is suitable for parallel computing with semi-implicit numerical integration method, so it can be used in real-time signal processing. Besides, the system can be realized simply either by independent circuit structure [28] as shown in [31] or by recursive algorithm with universal digital integrated circuit.

IV. FURTHER DENOISING
The simulation results in the above part show that there are some small parasitic oscillations in the output waveform when the Duffing oscillator system is used only. As mentioned before, the reason is that the high frequency part of the signal is not completely removed by Duffing oscillator system.
In order to further remove the high-frequency part of the signal and ensure the operation speed, a wavelet denoising after Duffing oscillator denoising is proposed, as shown in Fig. 17.
The wavelet denoising method has no parameter search process and runs very fast. Adding this step after Duffing denoising does not affect the speed of the whole algorithm.
The signals shown in Fig. 6 (d), Fig. 7 (d) and Fig. 8 (d) are respectively input to the above ''wden'' function for denoising, and the corresponding results are obtained as shown in Fig. 18.
As can be seen from Fig. 18, the waveforms after wavelet filtering become smooth, and are more similar to the original signal waveforms.
It can be seen that the Duffing oscillator plus wavelet denoising method can recover the details of the signal more  effectively and has stronger anti-noise ability than other methods.
In order to illustrate the similarity between the denoised waveform and the original waveform of the above simulation examples, the correlation distance is used to measure it. Specifically, the Matlab's own function ''pdist (x, 'correlation')'' is used to calculate the correlation distance, and the calculation results are shown in Tab. 2.
It can be seen from Tab. 2 that when Duffing system is used alone, the denoising output waveform is not smooth and the correlation distance from the original signal is larger. The correlation distance between the output of multivariate wavelet denoising and original signal is smaller, and the result of multivariate wavelet denoising is closer to the original signal. It can also be found that the correlation distance of Duffing plus wavelet denoising is the smallest; that is to say, its denoising effect is the best.

V. CONCLUSION
In this paper, a new basic signal denoising method is constructed with a linear and nonlinear restoring forces coupled Duffing oscillators system, and the selection of system parameters is analyzed. In addition, wavelet denoising after Duffing oscillator denoising can get good signal denoising results. Compared with other methods, the simulation results show that the method has good denoising performance for complex and irregular signals in different SNR and fast computing speed. In addition, when the semi-implicit method is used, the structure of numerical solution of the proposed method is simple and easy to construct parallel computing, and the proposed method is suitable for real-time signal processing.