A Gyroscope Signal Denoising Method Based on Empirical Mode Decomposition and Signal Reconstruction

To suppress the random drift error of a gyroscope signal, this paper proposes a novel denoising method, which is based on processing the intrinsic mode functions (IMFs) obtained by empirical mode decomposition (EMD). Considering that a gyroscope signal contains colored noise in addition to Gaussian white noise, fractal Gaussian noise (FGN) was introduced to quantify the noise in the gyroscope data. The proposed denoising method combines the FGN energy model and the modified method of Hausdorff distance (HD) to adaptively divide the IMFs into three categories (pure noise, pure information, and mixed components of noise and information). Then, the information IMFs and the mixed components after thresholding were selected to give the optimal signal reconstruction. Static and dynamic signal tests of the fiber optic gyroscope (FOG) were carried out to illustrate the performance of the proposed method, and compared with other traditional EMD denoising methods, such as the Euclidean norm measure method (EMD-l2-norm) and the sliding average filtering method (EMD-SA). The results of the analysis of both the static and dynamic signal tests indicate the effectiveness of the proposed method.


Introduction
Optical gyroscopes have been widely used in strap-down inertial navigation systems (SINS) [1,2], in which fiber optic gyroscopes (FOGs) have become ideal sensors because of their low cost, light weight, high reliability, easy integration, and high accuracy [3]. FOGs are based on the Sagnac effect, which is mainly composed of a light source, a phase modulator, a fiber coil, a detector, and a coupler, each of which can generate noise [4]. Similar to the core sensor in SINS, the random errors of FOG increase cumulatively with time [5], which further affects the accuracy of the position, velocity, and attitude (PVA). Therefore, suppressing random drift errors is a key issue in gyroscope signal processing.
Two methods are commonly adopted to suppress or eliminate random errors: establishing a random errors model and directly denoising the gyroscope data [6,7]. The former method needs to ensure the accuracy of the random model, and is then combined with a variety of filtering algorithms to supplement the model. However, random errors are susceptible to various uncertainty factors, which have no obvious regularity. Thus, it is difficult to estimate random errors effectively [8]. Since the gyroscope data contains both the actual motion information and the sensor noise, directly denoising gyroscope data can effectively reduce random error. Therefore, the direct denoising method is widely used to reduce random errors of gyroscope sensors [3,6,9,10].
IMFs and a residual [19]. The IMFs and the residual are obtained through a screening process, the detailed steps for which are shown in Figure 1.
There are two judgment processes in Figure 1, which are needed to define the IMFs and the residual. The IMFs have no restrictions on linearity and stationarity, however, two conditions must be satisfied: (1) in the whole signal, the number of extremes must be equal to the number of zero crossings, or at most may differ by one; (2) at any point, the mean value of the envelopes is zero, where the upper envelope is determined by the local maxima and the lower envelope is determined by the local minima. For the residual, this is the stoppage criterion. When the function becomes monotonic or has only one extremum, the screening process finally stops. After decomposition, the signal can be expressed as: where x(t) is the raw signal; L is the number of IMFs; ℎ ( ) ( ) is the i-th order IMF; and ( ) is the residual. Through decomposition, all intrinsic mode functions (IMFs) become zero-mean waveforms, and the lower-order IMFs hold more extreme points than that of higher-order IMFs. In other words, the lower-order IMFs capture fast oscillations, while higher-order IMFs typically represent slow oscillations. The noises are mainly contained in the lower order IMFs, while the useful information is contained in the higher order IMFs.

Quantifying the Noise in the Gyroscope Sensor through Fractal Gaussian Noise
Fractal Gaussian noise is essentially a discrete time process, described as an incremental process of fractal Brownian motion (FBM) [36]. FBM is the only self-similar Gaussian process with fixed increments; therefore, its statistical properties are completely dependent on the second moment of There are two judgment processes in Figure 1, which are needed to define the IMFs and the residual. The IMFs have no restrictions on linearity and stationarity, however, two conditions must be satisfied: (1) in the whole signal, the number of extremes must be equal to the number of zero crossings, or at most may differ by one; (2) at any point, the mean value of the envelopes is zero, where the upper envelope is determined by the local maxima and the lower envelope is determined by the local minima. For the residual, this is the stoppage criterion. When the function becomes monotonic or has only one extremum, the screening process finally stops. After decomposition, the signal can be expressed as: where x(t) is the raw signal; L is the number of IMFs; h (i) (t) is the i-th order IMF; and r L (t) is the residual. Through decomposition, all intrinsic mode functions (IMFs) become zero-mean waveforms, and the lower-order IMFs hold more extreme points than that of higher-order IMFs. In other words, the lower-order IMFs capture fast oscillations, while higher-order IMFs typically represent slow oscillations. The noises are mainly contained in the lower order IMFs, while the useful information is contained in the higher order IMFs.

Quantifying the Noise in the Gyroscope Sensor through Fractal Gaussian Noise
Fractal Gaussian noise is essentially a discrete time process, described as an incremental process of fractal Brownian motion (FBM) [36]. FBM is the only self-similar Gaussian process with fixed increments; therefore, its statistical properties are completely dependent on the second moment of parameter H, which is termed the Hurst parameter, satisfying 0 < H < 1 [37]. Typically, fractal Gaussian noise is defined as a zero-mean Gaussian stationary process; its autocorrelation function is: where r H (τ) is an autocorrelation function; τ is the lag length; H is the Hurst parameter; and σ 2 is the variance, which is equal to r H (0). When τ → ∞ , we can get: For H > 0.5, r H (τ) decays slowly (long-range dependence) and its value is positive. For H < 0.5, r H (τ) decays rapidly (short-range dependence) and its value is negative. Importantly, for H = 0.5, the FGN is the white Gaussian noise: Thus, fractal Gaussian noise is a generalization of white Gaussian noise. White Gaussian noise describes uncorrelated error, and fractal Gaussian noise can describe correlated noise according to different H values. This fact extends the applicable range of the error model.
Since the gyroscope sensor contains both uncorrelated noise and colored noise, it is very suitable to use fractal Gaussian noise to construct the error model instead of white Gaussian noise. As the core parameter of FGN, the Hurst parameter should be estimated before the process to construct a noise model containing colored noise. A simple and practical method to estimate the value of the Hurst parameter [27] is power spectral density (PSD). First, take the discrete Fourier transform of Equation (2) to obtain the PSD of FGN: where S H ( f ) is the PSD of FGN; f is the frequency satisfying f ≤ 0.5; and C is the constant. If H 0.5, f → 0 , and the PSD of FGN is described as: In the Nyquist frequency band, Equation (6) can be established. This is translated into the double logarithm as: log A least squares adjustment of Equation (7) is then used to calculate the value of H. Let k be the slope of the fitted line. Then: The Hurst parameter can expose all the information on a gyroscope's sensor errors. Therefore, after the value of H is estimated, we can quantify the noise in the gyroscope sensor through the FGN.

Proposed Denoising Method
The decomposed IMFs are arranged from high to low according to frequency and amplitude. We divided all IMFs into three categories: noise IMFs, mixed IMFs, and information IMFs. The reconstruction scheme of the optimal signals is shown in Figure 2. The components defined as noise IMFs are discarded, and the mixed IMFs are processed by the hard interval threshold, while the information IMFs and the residual are completely reserved. Based on the above analysis, the denoised signal can be expressed as where x(t) is the denoised signal; M 1 is the boundary of the noise IMFs and mixed IMFs, and M 2 represents the boundary of mixed IMFs and information IMFs; h (i) (t) is the process of thresholding; r L (t) is residual. The mixed IMFs and information IMFs can be expressed as M 2 where ( ) is the denoised signal; is the boundary of the noise IMFs and mixed IMFs, and represents the boundary of mixed IMFs and information IMFs; ℎ ( ) ( ) is the process of thresholding; ( ) is residual. The mixed IMFs and information IMFs can be expressed as ∑ ℎ ( ) ( ) and ∑ ℎ ( ) ( ) , respectively. In summary, the optimal signals for reconstruction are the mixed IMFs after thresholding, the information IMFs, and the residual. The steps of the proposed denoising method are as follows: 1. Estimate the value of the Hurst parameter (H) and quantify the noise through FGN; 2. Implement the EMD algorithm to obtain the IMFs and residual; 3. Determine the pure noise IMFs by comparing the significant difference between the actual and theoretical energy values; 4. Determine the mixed IMFs using an improved Hausdorff distance (HD) method; 5. Threshold the mixed IMFs via a hard interval threshold; 6. Reconstruct the optimal signal. Therefore, the key steps of this method are to classify the IMFs (determine the values of and ) and filter the mixed IMF.

Determining the Value of and
To determine the value of , the energy values of the noise were calculated. Flandrin et al. estimated the empirically observed energy of IMFs in 2004 [27], and proposed that the energy values of pure noise can be estimated as where [ ] is the noise energy estimate of the k-th order IMF; N is the length of the data; H is the In summary, the optimal signals for reconstruction are the mixed IMFs after thresholding, the information IMFs, and the residual. The steps of the proposed denoising method are as follows:

1.
Estimate the value of the Hurst parameter (H) and quantify the noise through FGN; 2.
Implement the EMD algorithm to obtain the IMFs and residual; 3.
Determine the pure noise IMFs by comparing the significant difference between the actual and theoretical energy values; 4.
Determine the mixed IMFs using an improved Hausdorff distance (HD) method; 5.
Threshold the mixed IMFs via a hard interval threshold; 6.
Reconstruct the optimal signal. Therefore, the key steps of this method are to classify the IMFs (determine the values of M 1 and M 2 ) and filter the mixed IMF.  [27], and proposed that the energy values of pure noise can be estimated asŴ whereŴ H [k] is the noise energy estimate of the k-th order IMF; N is the length of the data; H is the Hurst exponent; d 1,H [n] is the first-order IMF sequence; and β H is a constant. For IMF power spectra, the behavior of the first IMF, which has high-pass characteristics, is completely different from that of other modes. For modes k = 2 to 6, the power spectra appear to be nearly the same; the spectrum of IMF 7 gradually evolves from a band-pass to low-pass as the H increases [38]. Based on the characteristics of the spectra, we directly set the first IMF to the pure noise mode, while other noise modes are generated from IMF 2 to IMF 6 . By comparing the actual energy value and the theoretical energy value of the IMFs, the modes that can match the theoretical value without significantly difference are defined as noise IMFs. Here, we set the maximum difference value of the energy between the theoretical noise and the actual IMF as M 1 .
Next, the Hausdorff distance (HD) between the probability density function (pdf) of the input signal and that of each IMF is introduced to determine the value of M 2 . HD is a nonlinear operator, which can be used to measure the degree of similarity between two geometric shapes or two sets [39]. The pdf can denote the distribution shape of the signal [40]. By combining the two algorithms, the difference and similarity of the two signals can calculated. The HD is defined as follows.
Given two sets A = {a 1 , · · · , a M } and B = {b 1 , · · · , b N }, the distance is defined as: where H(A, B) is the HD of A and B. Here, h(A, B) is the one-way distance from set A to set B.
Similarly, h(B, A) is the one-way distance from set B to set A and a − b is the Euclidean distance between the two points. The probability density functions of the IMFs can be estimated using the kernel density estimator. The similarity measurement between the input signal x(t) and each IMF h i (t) expressed by the PDF can be defined as follows: where L(i) is the similarity measurement and HD stands for the calculated distance between the two PDFs estimated by the Hausdorff Distance. M 2 is the next point of the Hausdorff distance maximum; that is, the first decreasing point, which can be identified by In order to prevent useful information in the high order IMFs from being discarded, we set M 2 to be generated in the first half of the IMFs. If there is no decreasing point in the first half of the IMFs, then M 2 = L/2. This improved method can be given by: The improved HD method gives a good choice of mixed mode to prevent removing useful information. In this case, the second half of the IMF does not participate in the thresholding.

Thresholding for Mixed IMFs
After determining the M 1 and M 2 boundaries of the IMFs, the mixed IMFs can be determined as IMF M 1 to IMF M 2 , which should be denoised. Inspired by the wavelet threshold, EMD direct thresholding (EMD-DT) is generated [34]. EMD-DT can be modified by hard and soft thresholding [41], which can be expressed as Equations (20) and (21), respectively: where T i is the threshold of the i-th IMF. However, the direct hard threshold causes discontinuities, and oscillations are likely to occur. Although the soft threshold is continuous, there is permanent bias from the actual signal after processing. This threshold is a good choice to circumvent permanent bias and overcome discontinuities. An interval threshold denoising method based on the mode cell filter is introduced. The signal between two adjacent zero-crossings in the IMFs is defined as a mode unit, which can be described as Thus, the hard threshold is translated into: where Z i j is all the samples from Z i j to Z i j+1 ; and h (i) r  (10) and (11). Other variations in IMFs can be expressed as: The interval threshold is: The framework of our denoising method is shown in Figure 3, where all steps in the framework can be calculated by the above equations.
The framework of our denoising method is shown in Figure 3, where all steps in the framework can be calculated by the above equations.

Experimental Design
The experiments were performed on both static and dynamic test signal, which were collected in a laboratory with a room temperature of 25 °C. The fiber optic gyroscopes (FOG), the construction of the turntable test, and the data transmission equipment are shown in Figure 4. The turntable was made by China National Shipbuilding Corporation; its rotational angular velocity is 0.0001°. The FOG was made by Beijing Aerospace Times Optical-Electronics Technology Co., and is a test product with a zero offset of 0.3 °/h and a sample period of 0.39 ms. The static data collected from the gyroscope remained stationary for more than one hour, and each rotational dynamic datum was collected from the turntable after more than ten minutes. The dynamic test was set to five different uniform rotation

Experimental Design
The experiments were performed on both static and dynamic test signal, which were collected in a laboratory with a room temperature of 25 • C. The fiber optic gyroscopes (FOG), the construction of the turntable test, and the data transmission equipment are shown in Figure 4. The turntable was made by China National Shipbuilding Corporation; its rotational angular velocity is 0.0001 • . The FOG was made by Beijing Aerospace Times Optical-Electronics Technology Co., and is a test product with a zero offset of 0.3 • /h and a sample period of 0.39 ms. The static data collected from the gyroscope remained stationary for more than one hour, and each rotational dynamic datum was collected from the turntable after more than ten minutes. The dynamic test was set to five different uniform rotation speeds: 10 • /s, 20 • /s, 30 • /s, 40 • /s, and 50 • /s, respectively. The FOG and the turntable were restarted for each different test.
Moreover, we compare the denoising effects between the traditional and proposed methods.   To quantitatively illustrate the de-nosing effects of the different methods, we selected the root mean square error (RMSE) as the indicator. The RMSE is calculated as Equation (25), which can reflect the absolute accuracy of the denoising signal, where ( ) and ( ) are the original signal and denoising signal, respectively; N is the total number of original signals. The smaller the RMSE, the better the performance of the denoising method. The amplitude spectrum and Allan variance analysis are used to describe the denoising effect, and the values of bias drift, angle random walking (ARW), quantization noise, and RMSE are also calculated for the four schemes.

Static Test Signal
For the real static data, we used a fiber optic gyroscope (FOG) sensor with a zero offset of 0.3 °/h. The sample period was 0.39 ms, and the total number of data samples was N = 9262,589. To avoid the influence of human operations of the instrument during the startup and shutdown phases, 100,000 data points in the middle part of the samples were chosen. Using these 100,000 data points and Equation (8), the Hurst exponent was estimated as a value of 0.5167. The data were decomposed into 17 IMFs and one residual by EMD processing, and are shown in Figure 5. After determining the pure noise IMFs, the mixed IMFs, and the information IMFs, by thresholding the mixed IMFs and reconstructing the optimal signal, we obtained the denoising results. Figure 6 shows the static test signal and the denoised signals using the three denoising methods. The proposed method (Scheme 4) offers the best denoising effect, with which the noise is significantly reduced. To quantitatively illustrate the de-nosing effects of the different methods, we selected the root mean square error (RMSE) as the indicator. The RMSE is calculated as Equation (25), which can reflect the absolute accuracy of the denoising signal, where y(t) and y(t) are the original signal and denoising signal, respectively; N is the total number of original signals. The smaller the RMSE, the better the performance of the denoising method. The amplitude spectrum and Allan variance analysis are used to describe the denoising effect, and the values of bias drift, angle random walking (ARW), quantization noise, and RMSE are also calculated for the four schemes.

Static Test Signal
For the real static data, we used a fiber optic gyroscope (FOG) sensor with a zero offset of 0.3 • /h. The sample period was 0.39 ms, and the total number of data samples was N = 9262,589. To avoid the influence of human operations of the instrument during the startup and shutdown phases, 100,000 data points in the middle part of the samples were chosen. Using these 100,000 data points and Equation (8), the Hurst exponent was estimated as a value of 0.5167. The data were decomposed into 17 IMFs and one residual by EMD processing, and are shown in Figure 5. After determining the pure noise IMFs, the mixed IMFs, and the information IMFs, by thresholding the mixed IMFs and reconstructing the optimal signal, we obtained the denoising results. Figure 6 shows the static test signal and the denoised signals using the three denoising methods. The proposed method (Scheme 4) offers the best denoising effect, with which the noise is significantly reduced.
The amplitude spectra of the signals are obtained by fast Fourier transform (FFT), and are shown in Figure 7. The real test signal of FOG contains significant noise with a high frequency (Scheme 1). The range of amplitude is significantly reduced for high frequencies via the three denoising methods, which means that the higher frequency IMFs of the real test signal have been removed. For the EMD-SA method (Scheme 2), each IMF is filtered to remove a portion of the noise, but the high-frequency noise is retained. The EMD-l 2 -norm (Scheme 3) and the proposed methods (Scheme 4) can deal with high-frequency noise thoroughly, leaving only the low-frequency parts. The proposed denoising method retains the ability to handle low frequency signals.     The amplitude spectra of the signals are obtained by fast Fourier transform (FFT), and are shown in Figure 7. The real test signal of FOG contains significant noise with a high frequency (Scheme 1). The range of amplitude is significantly reduced for high frequencies via the three denoising methods, which means that the higher frequency IMFs of the real test signal have been removed. For the EMD-SA     The amplitude spectra of the signals are obtained by fast Fourier transform (FFT), and are shown in Figure 7. The real test signal of FOG contains significant noise with a high frequency (Scheme 1). The range of amplitude is significantly reduced for high frequencies via the three denoising methods, which means that the higher frequency IMFs of the real test signal have been removed. For the EMD-SA Figure 6. Gyroscope signal and its denoised results for three denoising methods. The red line is the gyroscope signal; the blue, orange, and green lines are the denoising results of Euclidean norm measure method (EMD-l 2 -norm), the sliding average filtering method (EMD-SA), and the proposed method, respectively. method (Scheme 2), each IMF is filtered to remove a portion of the noise, but the high-frequency noise is retained. The EMD--norm (Scheme 3) and the proposed methods (Scheme 4) can deal with highfrequency noise thoroughly, leaving only the low-frequency parts. The proposed denoising method retains the ability to handle low frequency signals. Allan variance is recommended by the Institute of Electrical and Electronics Engineers (IEEE) to analyze random gyroscope error. Thus, we adopted Allan variance to further compare the effects of the three denoising methods. The results are shown in Figure 8. The three de-nosing methods can reduce the noise, and the proposed method offers the best performance. Allan variance is recommended by the Institute of Electrical and Electronics Engineers (IEEE) to analyze random gyroscope error. Thus, we adopted Allan variance to further compare the effects of the three denoising methods. The results are shown in Figure 8. The three de-nosing methods can reduce the noise, and the proposed method offers the best performance. To further accurately describe the denoising effects, we collected a six-hour output signal and calculated the denoising results for different schemes. In Table 1, we can see that the values of bias drift, angle random walking (ARW), quantization noise, and RMSE for the proposed method in this paper are 0.0038, 0.0191, 0.1622, and 0.0045, respectively, which are lower than the values for other denoising methods. Compared with Scheme 1, the proposed method (Scheme 4) improved the denoising effects by 87.7%, 63.4%, 68.4%, and 85.6%, respectively. To further accurately describe the denoising effects, we collected a six-hour output signal and calculated the denoising results for different schemes. In Table 1, we can see that the values of bias drift, angle random walking (ARW), quantization noise, and RMSE for the proposed method in this paper are 0.0038, 0.0191, 0.1622, and 0.0045, respectively, which are lower than the values for other denoising methods. Compared with Scheme 1, the proposed method (Scheme 4) improved the denoising effects by 87.7%, 63.4%, 68.4%, and 85.6%, respectively.

Dynamic Test Signal
The FOG in the dynamic test used for data acquisition is the same sensor as the static test; from these tests, five sets of dynamic rotation data were collected. In each dynamic test (the total number of data samples was N = 1800,000); we chose 100,000 stable rotation data points in the middle part of the samples for analysis. After the denoising process, the five sets of signals were denoised using the three denoising methods, which are shown in Figure 9. From this figure, we can see that the amplitude of the noise increases with the increase in the rotational speed. The three denoising methods (Schemes 2-4) can reduce the noise in all dynamic data, and the proposed method (Scheme 4) offers the better performance, as the noise can be significantly reduced.  Figure 9. The gyroscope dynamic signal and its denoised results for three denoising methods in the five dynamic tests. The red line is the gyroscope's dynamic signal; the blue, orange, and green lines are the denoising results of the EMD-l2-norm, EMD-SA, and the proposed method, respectively. (a-e) Dynamic gyroscope signals and the denoised results for the three denoising methods, with speeds of 10 °/s, 20 °/s, 30 °/s, 40 °/s, and 50 °/s, respectively. Table 2 shows the results of the dynamic signal with different rotational speeds. In Table 2, we can see that all Hurst values for the dynamic data are less than 0.5. Obviously, the random errors of this gyroscope present a negative correlation in this test. The RMSE and standard deviation (STD) are improved by the three denoising methods (Schemes 2-4), and the proposed method (Scheme 4) has the minimum value. The proposed method improves the RMSE by 4.5 times. Table 2. Comparison of root mean square error (RMSE) and standard deviation (STD) for the  Table 2 shows the results of the dynamic signal with different rotational speeds. In Table 2, we can see that all Hurst values for the dynamic data are less than 0.5. Obviously, the random errors of this gyroscope present a negative correlation in this test. The RMSE and standard deviation (STD) are improved by the three denoising methods (Schemes 2-4), and the proposed method (Scheme 4) has the minimum value. The proposed method improves the RMSE by 4.5 times.

Conclusions
Traditional denoising methods based on empirical mode decomposition (EMD) are mainly classified into two categories: the partial reconstruction of relevant modes and the whole reconstruction of all filtered modes [26,27]. However, when the signal-to-noise ratio (SNR) of the signal is high, the useful signal is also decomposed into lower-order intrinsic mode functions (IMFs), in which case the useful information may be mistaken for the discarded irrelevant information. To avoid losing useful information, the irrelevant IMFs should be defined as accurately as possible. To this end, we proposed a novel method to denoise the noisy signal; our focus was on discriminating the IMF type (pure noise IMF, mixed IMF, or information IMF). Then, the IMFs in different categories were processed separately.
The proposed method is fully data-driven and adaptive. Firstly, the original signal was decomposed by EMD to obtain IMFs according to the characteristics of the signal itself. Next, in the process of IMF classification and for the interval threshold, no prior parameter settings were required to obtain the optimal signal reconstruction. Then, the values of the threshold and the mode demarcation points required for denoising were also calculated entirely from the IMFs.
This improved method can reduce high frequency noise from a noisy signal, which offers enhanced performance compared to other denoising methods in our experiment. For the real signal test in this paper, the denoising results of the static and dynamic data indicate the effectiveness of the new method, which improves the values of bias drift, angle random walking (ARW), quantization noise, and root mean square error (RMSE) by 87.7%, 63.4%, 68.4%, and 85.6%, respectively. However, the value of ARW ranges from 0.0526 to 0.0192 • / √ h and the value of the bias drift ranges from 0.0308 to 0.0038, which is not good compared to state-of-the-art models [1,4]. The FOG we used is not highly accurate. In addition, improved EMD algorithms [42,43] and different thresholding methods for mixed IMFs will be adopted in our next study, to further improve denoising performance.