MEMS Gyroscope Bias Drift Self-Calibration Based on Noise-Suppressed Mode Reversal

This paper presents a bias drift self-calibration method for micro-electromechanical systems (MEMS) gyroscopes based on noise-suppressed mode reversal without the modeling of bias drift signal. At first, the bias drift cancellation is accomplished by periodic switching between operation mode of two collinear gyroscopes and subtracting the bias error which is estimated by the rate outputs from a consecutive period interval; then a novel filtering algorithm based on improved complete ensemble empirical mode decomposition (improved complete ensemble empirical mode decomposition with adaptive noise—CEEMDAN) is applied to eliminate the noise in the calibrated signal. A set of intrinsic mode functions (IMFs) is obtained by the decomposition of the calibrated signal using improved CEEMDAN method, and the threshold denoising method is utilized; finally, the de-noised IMFs are reconstructed into the desired signal. To verify the proposed method, the hardware circuit with an embedded field-programmable gate array (FPGA) was implemented and applied in bias drift calibration for the two MEMS gyroscopes manufactured in our laboratory. The experimental results indicate that the proposed method is feasible, and it achieved a better performance than the typical mode reversal. The bias instability of the two gyroscopes decreased from 0.0066°/s and 0.0055°/s to 0.0011°/s; and, benefiting from the threshold denoising based on improved CEEMDAN, the angle random walks decreased from 1.18 ×10−4°/s1/2 and 2.04 ×10−4°/s1/2 to 2.19 ×10−5°/s1/2, respectively.


Introduction
Micro-electromechanical system (MEMS) gyroscopes are developing rapidly due to their extensive applications in the field of consumer electronics and in the military, such as game stations, cameras, vehicle control systems and weapon guidance. Compared with traditional gyroscopes, the MEMS gyroscopes have the advantages of low costs, low power consumption and small size [1].
The major problem for MEMS gyroscopes is to reach the requirement of inertial grade which means a greatly improvement in the performance of bias stability. Therefore, a large number of studies have been conducted, including on the vibrating poly-silicon ring gyroscope, the tuning fork SOI (silicon on insulator) gyroscope with mode-matching and the fully-decoupled SOI/SOG (silicon on glass) gyroscope, which mainly focuses on the MEMS gyroscope-sensing element [2][3][4][5][6][7][8]. Previous published works mainly focus on the superb mechanical characteristics of the MEMS gyroscope sensing element, which has a high requirement for processing technology. In this work, only the bias drift self-calibration method which can improve the accuracy of MEMS gyroscope at a general processing level is concerned.

The Principle and Theoretical Analysis of Mode Reversal
A simplified dynamic model of the MEMS gyroscope having two degrees of freedom can be described as: y + k 2 y = τ y − 2mΩ z .
x . (1) As shown in Equation (1), only the component of the input angular rate along the z-axis causes a coupling between the x-axis and y-axis in an ideal gyroscope. However, manufacture imperfections always exist, which lead to an anisotropic and asymmetric gyroscope structure, resulting in dynamic coupling between the x-axis and y-axis. The asymmetry and anisoelasticity of the structure can be captured as misalignment of principal axis of spring from the physical axis of gyroscope structure. In addition to the structural asymmetry, misalignment of the actuators also contribute to the asymmetric spring terms.
Assume that the principal spring constants are k 1 and k 2 , and the principal spring axes are tilted by angle ε with respect to the drive and sense axis-caused by fabrication imperfections, as shown in Figure 1.
Micromachines 2019, 19, 3 of 17 A simplified dynamic model of the MEMS gyroscope having two degrees of freedom can be described as: As shown in Equation (1), only the component of the input angular rate along the z-axis causes a coupling between the x-axis and y-axis in an ideal gyroscope. However, manufacture imperfections always exist, which lead to an anisotropic and asymmetric gyroscope structure, resulting in dynamic coupling between the x-axis and y-axis. The asymmetry and anisoelasticity of the structure can be captured as misalignment of principal axis of spring from the physical axis of gyroscope structure. In addition to the structural asymmetry, misalignment of the actuators also contribute to the asymmetric spring terms.
Assume that the principal spring constants are and , and the principal spring axes are tilted by angle with respect to the drive and sense axis-caused by fabrication imperfections, as shown in Figure 1. Let {x, y} be the physical drive and sense axis, and {x', y'} be the principal axis of drive and sense resonator. Then, the relationship between spring forces and spring constants in the principal axes is given by: Since the principal spring constants and physical axes are related by: the spring constants in physical drive and sense axis are given by: As shown in Figure 1, the spring forces of principal axes and the actual spring force along the physical drive and sense axis are related by: From Equations (2)-(5), we obtain: Let {x, y} be the physical drive and sense axis, and {x , y } be the principal axis of drive and sense resonator. Then, the relationship between spring forces and spring constants in the principal axes is given by: Since the principal spring constants and physical axes are related by: the spring constants in physical drive and sense axis are given by: Micromachines 2019, 10, 823

of 17
As shown in Figure 1, the spring forces of principal axes and the actual spring force along the physical drive and sense axis are related by: From Equations (2)-(5), we obtain: since the matrix of the actual spring constant can be represented as: where Asymmetric damping can also be represented in terms of principal axes. The coordinate of the principal axis and the magnitude of the damping are determined by an averaging of the damping asymmetries. Similarly to the spring stiffness coefficients case, if damping constants along the principal axes are given by d 1 and d 2 , and the axes are tilted by an angle θ from physical axes of gyroscope, then the damping constant coefficients in physical axes are shown as: Taking into account fabrication imperfections, the dynamic equation, Equation (1), is modified as follows: Equation (14) is the governing equation for the gyroscope in the normal operation mode. Assume that drive displacement x = q 0 cos(ω x t), where the q 0 is the magnitude of the drive displacement and ω x is the natural frequency of drive mode (x-axis). By substituting the stable state condition of feedback loop and the drive displacement into Equation (14), we obtain the sense force τ y as: Through the interface module, demodulation module and low-pass filter, the output signal which contains the angular rate is represented as: Micromachines 2019, 10, 823 (17) and ϕ denotes the total phase delay of the interface module, demodulation module, and low-pass filter. SF and Bias denote the scale factor and bias of the gyroscope in normal operation mode respectively. When the gyroscope in the reverse operation mode, the drive resonator and sense resonator are switched. The location of principal axes of drive and sense mode remains the same, and the coordinates of physical axes rotate 90 • which means that ε = π 2 − ε and θ = π 2 − θ, as shown in Figure 2.
Micromachines 2019, 19, 5 of 17 and denotes the total phase delay of the interface module, demodulation module, and low-pass filter. and denote the scale factor and bias of the gyroscope in normal operation mode respectively.
we obtained the the actual spring constants matrix of mode reversal as follows: where Similarly, the damping constant coefficients in physical axes are shown as: Since Equation (3) and Equation (5) are changed to: we obtained the the actual spring constants matrix of mode reversal as follows: where Similarly, the damping constant coefficients in physical axes are shown as: Therefore, the governing equation for the gyroscope in the mode reversal operation mode can be described as: x + k yx x + k yy y = τ y − 2mΩ .
x . ( Assume that drive displacement y = q 0 sin ω y t , where the q 0 is the magnitude of the drive displacement and ω y is the natural frequency of drive mode (y-axis). By substituting the stable state condition of feedback loop and the drive displacement into Equation (20), we obtain the sense force τ x as: Through the signal processing module includes that interface circuit, demodulation circuit, and low-pass filter, the output signal which contains the angular rate is written as: where and ϕ denotes the phase delay of the signal processing module which has a value identical to normal in operation mode. SF and Bias denote the scale factor and bias of the gyroscope in mode reversal operation mode respectively. From Equations (17) and (30), the estimation of bias drift and measured angular rate can be obtained as: As shown in Equation (33), the bias drift can be eliminated by utilizing the two rate outputs of normal operation mode and reverse operation mode, and double the sensitivity of angle rate.
However, the typical mode reversal scheme degrades the overall noise performance by folding uncorrelated noise from the upper bands down to the baseband and reduced measurement bandwidth, which is limited, just like correlated double sampling method in circuits. Hence, the threshold denoising method was adopted to enhance the noise performance, and an additional gyroscope was applied to increase the bandwidth.

The Threshold Denoising Algorithm Based on Improved CEEMDAN
Empirical mode decomposition (EMD) is an adaptive analysis method for non-stationary signals in nonlinear systems [15]. A series of functions which are called intrinsic mode functions (IMFs) can be obtained by the decomposition of original signal using EMD. However, the original algorithm of EMD exists as an un-desirable phenomenon which is called "mode mixing." In order to solve this issue, the ensemble empirical mode decomposition (EEMD) was presented in [16]. The EEMD utilized the additional white Gaussian noise to suppress the mode mixing by making use of the dyadic filter bank behavior. But EEMD also created new difficulties such that a lot of residual noise was introduced in the reconstructed signal, which is the sum of the modes and the final trend. The Complementary EEMD, which is presented in [17], can reduce the residual noise in signal reconstruction by using complementary pairs of noise expressively. Nevertheless, the complementary pairs of noise may generate a large number of additional modes.
The complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) achieves a negligible residual noise, and the problem of additional modes caused by complementary pairs of noise has also been solved. However, CEEMDAN still has some problems to be solved: there are some residual noise in the modes and spurious modes occur in the early phase of decomposition. In order to avoid the problems, in previous works, the improved CEEMDAN algorithm [18] was adopted. The flowchart of this algorithm is shown as Figure 3. In Figure 3, E k (·) is the operator of the k-th mode of EMD extraction, M(·) is the operator of local mean, w (i) is the Gaussian white noise with zero mean unit variance, and the constant β k is set as β k = ε k std(r k ).
Micromachines 2019, 19, 7 of 17 issue, the ensemble empirical mode decomposition (EEMD) was presented in [16]. The EEMD utilized the additional white Gaussian noise to suppress the mode mixing by making use of the dyadic filter bank behavior. But EEMD also created new difficulties such that a lot of residual noise was introduced in the reconstructed signal, which is the sum of the modes and the final trend. The Complementary EEMD, which is presented in [17], can reduce the residual noise in signal reconstruction by using complementary pairs of noise expressively. Nevertheless, the complementary pairs of noise may generate a large number of additional modes. The complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) achieves a negligible residual noise, and the problem of additional modes caused by complementary pairs of noise has also been solved. However, CEEMDAN still has some problems to be solved: there are some residual noise in the modes and spurious modes occur in the early phase of decomposition. In order to avoid the problems, in previous works, the improved CEEMDAN algorithm [18] was adopted. The flowchart of this algorithm is shown as Figure 3. In Figure 3, (•) is the operator of the k-th mode of EMD extraction, (•) is the operator of local mean, ( ) is the Gaussian white noise with zero mean unit variance, and the constant is set as = ( ).
The set of , , … , is the final result of the decomposition of the original signal. A set of IMFs is obtained by the decomposition of the original signal using improved CEEMDAN method. Both useful components and color noise are contained in the set of IMFs; thus, the threshold denoising method should be adopted [19]. The which is the points in the k-th IMF, can be expressed as ( ) , = 1, … , , and the estimation of empirical variance of is shown in Equation (34): The information of desired and noise signal is all contained in ( ), and the variance of pure noise should be calculated, which is to determine the threshold of denoising. Therefore, the sample entropy (SE) has been adopted [20]. The SE value of signal is positively correlated with the amount The set of IMF 1 , IMF 2 , . . . , IMF k is the final result of the decomposition of the original signal. A set of IMFs is obtained by the decomposition of the original signal using improved CEEMDAN method. Both useful components and color noise are contained in the set of IMFs; thus, the threshold denoising method should be adopted [19]. The im f k which is the points in the k-th IMF, can be expressed as d ( j) k , j = 1, . . . , N , and the estimation of empirical variance of im f k is shown in Equation (34): The information of desired and noise signal is all contained inV(im f k ), and the variance of pure noise should be calculated, which is to determine the threshold of denoising. Therefore, the sample entropy (SE) has been adopted [20]. The SE value of signal is positively correlated with the amount of noise in signal. SE is primarily determined by tow parameters: the embedding dimension m and the tolerance value r. The formula for calculating SE can be shown as follows: After the calculation of sample entropy of each IMF, the first two IMFs which have the highest sample entropy are selected as pure noise signal. The robust estimator shown in Equation (36) can be applied to estimate the standard deviations of noise contained in the two selected IMFs: The noise of MEMS gyroscope can be modeled by fractional Gaussian noise; therefore, the variances of the noise that contained in IMFs can be calculated by Equations (36) and (37). Similarly to the wavelet threshold denoising, the threshold of the denoising algorithm based on improved CEEMDAN can be obtained by Equation (38).
The commonly used threshold denoising algorithm mainly includes two types: the hard-threshold and soft-threshold. The hard-threshold denoising is shown as below: Similarly, the soft-threshold denoising can be expressed as: Since it can make the de-noised waveform smoother, the soft-threshold denoising is adopted in this paper.

The Self-Calibration System Based on Noise-Suppressed Mode Reversal
The single gyroscope bias drift self-calibration model described in Section 2.1 is helpful for principle analysis and performance projections; nevertheless, several difficulties have to be solved to realize a practical real-time self-calibrated system with uninterrupted output. Such a system can be designed using an additional gyroscope for a single axis system. The single axis system in this work consists of two collinear gyroscopes, a state machine sequencer, a bias estimator, and the threshold denoising algorithm based on improved CEEMDAN; it is shown in Figure 4.
The state machine is illustrated in Figure 5. As shown in Figure 5, bias error signal polarities are associated with gyroscopes A and B respectively. According to the analysis in the Section 2.1, the detected input rate signal has a positive relationship with the bias error signal when the gyroscope is in normal operation mode, and a negative relationship when the gyroscope is in reverse operation mode. The state machine is illustrated in Figure 5. As shown in Figure 5, bias error signal polarities are associated with gyroscopes A and B respectively. According to the analysis in the Section 2.1, the detected input rate signal has a positive relationship with the bias error signal when the gyroscope is in normal operation mode, and a negative relationship when the gyroscope is in reverse operation mode.   The state machine is illustrated in Figure 5. As shown in Figure 5, bias error signal polarities are associated with gyroscopes A and B respectively. According to the analysis in the Section 2.1, the detected input rate signal has a positive relationship with the bias error signal when the gyroscope is in normal operation mode, and a negative relationship when the gyroscope is in reverse operation mode. The graph shows four different measurement intervals, Meas1 to Meas4. The operation mode of the gyroscopes A and B switch alternately between normal and reverse mode according to the sequence of state, which enables uninterrupted real-time bias drift calibration. The respective bias error signals of the gyroscopes ( ) and ( ) also switch between the positive and negative polarities, and make bias drifts of both gyroscopes observable by the self-calibration bias estimator. The measurement signal , consists of detecting angular rate signal and bias error signal of Gyro A ( ). Similarly, measurement signal , consists of detecting angular rate signal and bias error signal of Gyro B ( ) . The bias estimator receives the measurement signals , and , and the calibrated outputs , and , which are used to estimate bias error in next measurement interval. The estimated bias error signal ( ) and ( ) are generated by the bias estimator, which are used for error calibration to obtain more accurate measurement. The graph shows four different measurement intervals, Meas1 to Meas4. The operation mode of the gyroscopes A and B switch alternately between normal and reverse mode according to the sequence of state, which enables uninterrupted real-time bias drift calibration. The respective bias error signals of the gyroscopes Bias(A) and Bias(B) also switch between the positive and negative polarities, and make bias drifts of both gyroscopes observable by the self-calibration bias estimator.
The measurement signal V Meas,A consists of detecting angular rate signal V in and bias error signal of Gyro A Bias(A). Similarly, measurement signal V Meas,B consists of detecting angular rate signal V in and bias error signal of Gyro B Bias(B). The bias estimator receives the measurement signals V Meas,A and V Meas,B and the calibrated outputs V out,A and V out,B which are used to estimate bias error in next measurement interval. The estimated bias error signalB ias(A) andB ias(B) are generated by the bias estimator, which are used for error calibration to obtain more accurate measurement.
Both gyroscopes are in a stable operational state during each measurement interval, and the sense axes of the input angle rates are identical. The measurement signals of Gyro A and Gyro B in the i-th measurement interval MeasA(i) and MeasB(i) are shown below, respectively: For the first and second intervals: The above Equations (43)-(46) are transformed into a matrix representation: where H is shown as: The inverse matrix of matrix H is given in Equation (49): The four variables V out (1), V out (2),Bias(A), andBias(B) are, therefore, individually observable, and can be obtained as: (50) The final output signal of the bias estimator V out is calibrated output signal, which can be obtained by Equation (50).
After the calibration by mode reversal, the threshold denoising is used to eliminate the noise in the calibrated signal. The process of threshold denoising algorithm is shown in Figure 4. A set of IMFs is obtained by the decomposition of the calibrated signal using improved CEEMDAN method. Then, we can determine the IMF with dominated noise by calculating the sample entropy of each IMFs, and substitute the variance of noise dominant IMF into the factional Gaussian noise model given in Equation (38) to obtain the threshold of denoising algorithm. Finally, the desired signal can be obtained by reconstruction of the de-noised IMFs.

The Implemention of the Proposed Self-Calibration System
To verify the feasibility of the proposed self-calibration system, the relevant experiment was carried out. The gyroscope we tested which has a symmetrical structure is fabricated by SOI processing, as shown in Figure 6, and the total dimensions of the sensing element are 4 mm × 4 mm [21]. To verify the feasibility of the proposed self-calibration system, the relevant experiment was carried out. The gyroscope we tested which has a symmetrical structure is fabricated by SOI processing, as shown in Figure 6, and the total dimensions of the sensing element are 4 mm × 4 mm [21]. The key parameters of the symmetric gyroscope are shown in Table 1. The two packaged gyroscopes with identical scale factors were arranged on two, identical fourlayer print circuit boards (PCB), and the rest of electronic components were mounted on a 100 mm × 80 mm, six-layer PCB, as shown in Figure 7. The key parameters of the symmetric gyroscope are shown in Table 1. The two packaged gyroscopes with identical scale factors were arranged on two, identical four-layer print circuit boards (PCB), and the rest of electronic components were mounted on a 100 mm × 80 mm, six-layer PCB, as shown in Figure 7.
carried out. The gyroscope we tested which has a symmetrical structure is fabricated by SOI processing, as shown in Figure 6, and the total dimensions of the sensing element are 4 mm × 4 mm [21]. The key parameters of the symmetric gyroscope are shown in Table 1. The two packaged gyroscopes with identical scale factors were arranged on two, identical fourlayer print circuit boards (PCB), and the rest of electronic components were mounted on a 100 mm × 80 mm, six-layer PCB, as shown in Figure 7. The two identical PCBs connected the structures of the weak signals. The third PCB was composed of the 24-bit resolution analog-to-digital converters ADS131 (Texas Instruments Company, Dallas, TX, USA); a FPGA chip XC5VL (Xilinx Company, San Jose, CA, USA) in which the control loops, the demodulation module, state machine sequencer, bias error estimator, and threshold denoising algorithm were implemented; and the 16-bit resolution digital-to-analog converters AD5754 (Analog Device, Norwood, MA, USA).

Experimental Testing
The output signal data was collected with 50 Hz sampling frequency and was recorded in PC with 1 s average value. The outputs of two collinear gyroscopes were recorded for 1000 s at 25 • C controlled temperature, and each measurement interval contains 150 points, as shown in Figure 8. The two identical PCBs connected the structures of the weak signals. The third PCB was composed of the 24-bit resolution analog-to-digital converters ADS131 (Texas Instruments Company, Dallas, TX, USA); a FPGA chip XC5VL (Xilinx Company, San Jose, CA, USA) in which the control loops, the demodulation module, state machine sequencer, bias error estimator, and threshold denoising algorithm were implemented; and the 16-bit resolution digital-to-analog converters AD5754 (Analog Device, Norwood, MA, USA).

Experimental Testing
The output signal data was collected with 50 Hz sampling frequency and was recorded in PC with 1 s average value. The outputs of two collinear gyroscopes were recorded for 1000 s at 25 ℃ controlled temperature, and each measurement interval contains 150 points, as shown in Figure 8.  From Figure 9, it is obvious that there is a large amount of irrelevant noise in the calibrated signal of mode reversal. In order to eliminate the influence of the additional noise in the calibrated signal on the measurement accuracy of MEMS gyroscopes, the threshold denoising method was adopted here. As described in Section 3, the first step of threshold denoising is to decompose the calibrated signal of mode reversal into a set of IMFs using the improved CEEMDAN. The decomposition results From Figure 9, it is obvious that there is a large amount of irrelevant noise in the calibrated signal of mode reversal. In order to eliminate the influence of the additional noise in the calibrated signal on the measurement accuracy of MEMS gyroscopes, the threshold denoising method was adopted here. As described in Section 3, the first step of threshold denoising is to decompose the calibrated signal of mode reversal into a set of IMFs using the improved CEEMDAN. The decomposition results of calibrated signal of mode reversal are given in Figure 10. From Figure 9, it is obvious that there is a large amount of irrelevant noise in the calibrated signal of mode reversal. In order to eliminate the influence of the additional noise in the calibrated signal on the measurement accuracy of MEMS gyroscopes, the threshold denoising method was adopted here. As described in Section 3, the first step of threshold denoising is to decompose the calibrated signal of mode reversal into a set of IMFs using the improved CEEMDAN. The decomposition results of calibrated signal of mode reversal are given in Figure 10.  In order to find the noise dominant IMF from the set of IMFs, the sample entropy of each IMFs is calculated according to Equation (35). The calculated result of the sample entropy of each IMFs is shown in Figure 11. In order to find the noise dominant IMF from the set of IMFs, the sample entropy of each IMFs is calculated according to Equation (35). The calculated result of the sample entropy of each IMFs is shown in Figure 11. From Figure 11, the IMF2 with local maximum value of sample entropy is selected as noise dominant component, and the variance of the IMF2 determines the denoising threshold of the filtering algorithm. Figure 12 gives the filtered results of IMFs after threshold denoising. The de-noised calibrated signal is obtained by the reconstruction of the denoised IMFs. From Figure 11, the IMF2 with local maximum value of sample entropy is selected as noise dominant component, and the variance of the IMF2 determines the denoising threshold of the filtering algorithm. Figure 12 gives the filtered results of IMFs after threshold denoising. The de-noised calibrated signal is obtained by the reconstruction of the denoised IMFs. Figure 11. The sample entropy of the set of intrinsic mode functions (IMFs).
From Figure 11, the IMF2 with local maximum value of sample entropy is selected as noise dominant component, and the variance of the IMF2 determines the denoising threshold of the filtering algorithm. Figure 12 gives the filtered results of IMFs after threshold denoising. The de-noised calibrated signal is obtained by the reconstruction of the denoised IMFs.
As shown in Figure 13, by comparing the calibrated signal before and after denoising, it can be concluded that the threshold denoising algorithm proposed in this paper can effectively filter out the noise in the calibrated signal. As shown in Figure 13, by comparing the calibrated signal before and after denoising, it can be concluded that the threshold denoising algorithm proposed in this paper can effectively filter out the noise in the calibrated signal.  The Allan variance method, which is a widely used quantitative performance evaluation method for inertial sensors [22], was applied for the quantitative comparison of calibration results of different methods. Figure 14 shows the Allan variance analysis of un-calibrated output signals of Gyro A and B, and the calibration results of the typical mode reversal and the proposed improved CEEMDANmode reversal method; and Table 1 shows the quantitative results in detail. The Allan variance method, which is a widely used quantitative performance evaluation method for inertial sensors [22], was applied for the quantitative comparison of calibration results of different methods. Figure 14 shows the Allan variance analysis of un-calibrated output signals of Gyro A and B, and the calibration results of the typical mode reversal and the proposed improved CEEMDAN-mode reversal method; and Table 1 shows the quantitative results in detail. The Allan variance method, which is a widely used quantitative performance evaluation method for inertial sensors [22], was applied for the quantitative comparison of calibration results of different methods. Figure 14 shows the Allan variance analysis of un-calibrated output signals of Gyro A and B, and the calibration results of the typical mode reversal and the proposed improved CEEMDANmode reversal method; and Table 1 shows the quantitative results in detail. From Table 2, we can see that the proposed self-calibration method improves bias instability from 0.0066 °/s and 0.0055 °/s to 0.0011 °/s , and the angle random walk decreases from 1.18 × 10 °/s ⁄ and 2.04 × 10 °/s ⁄ to 2.19 × 10 °/s ⁄ . From Table 2, we can see that the proposed self-calibration method improves bias instability from 0.0066 • /s and 0.0055 • /s to 0.0011 • /s, and the angle random walk decreases from 1.18 × 10 −4 • /s 1/2 and 2.04 × 10 −4 • /s 1/2 to 2.19 × 10 −5 • /s 1/2 . Both Figure 14 and Table 2 show that compared with the typical method, the proposed mode reversal method with improved CEEMDAN can enhance the bias drift stability and the performance of noise suppression significantly.

Conclusions
In summary, a novel bias drift self-calibration method based on noise-suppressed mode reversal for MEMS gyroscopes is proposed. Different from other approaches, the proposed method does not require modeling of bias drift signal. Moreover, benefiting from the improved CEEMDAN threshold denoising, this method eliminates the the noise folding from upper bands down to baseband, which is caused by the mode reversal, and can improve the overall noise performance of the calibration system. The proposed mode reversal self-calibration method does not decrease system bandwidth by adding additional collinear gyroscope. Compared to the typical mode reversal method, the proposed method decreases the bias instability and angle random walk effectively, specifically in Allan variance coefficients, B from 0.0022 • /s to 0.0011 • /s and N from 5.67 ×10 −4 • /s 1/2 to 2.19 ×10 −5 • /s 1/2 . The proposed mode reversal method may enable the MEMS gyroscopes to enter into tactical grade application where low-cost fiberoptic gyroscopes have traditionally dominated.