Characterization of Modal Frequencies and Orientation of Axisymmetric Resonators in Coriolis Vibratory Gyroscopes

This paper presents the characterization of the modal frequencies and the modal orientation of the axisymmetric resonators in Coriolis vibratory gyroscopes based on the approaches of the frequency sweep and the ring down. The modal frequencies and the orientation of the stiffness axis are the key parameters for the mechanical correction of the stiffness imperfections. The frequency sweep method utilizes the zero and the poles in the magnitude-frequency responses of the two-dimensional transfer function to extract the modal orientation information within the frequency domain. The ring down method makes use of the peak and the valley values of the beat signals at the readout electrodes to obtain the modal orientation and the coefficient of the nonlinear stiffness directly within the time domain. The proposed approaches were verified via a silicon ring resonator designed for gyroscopic sensing and the modal information from the experiments exhibited a good agreement between the methods of the frequency sweep and the ring down.


Introduction
Microelectromechanical system (MEMS) Coriolis vibratory gyroscopes (CVGs) have been revealed enormous potential for high-end applications, even for the applications of near-navigation grade [1][2][3][4][5]. Among the MEMS gyroscopes with the top performances, the MEMS CVGs with axisymmetric resonators, for instance, ring [6] or disk [7][8][9], hemispherical [10][11][12], and quad-mass [13][14][15] resonators, are particularly appealing by virtue of their angular rate-integrating ability, which offers benefits of high dynamics, wide ranges, and shock and vibration resistance. Besides, the symmetry in the structure of the resonators provides possibilities of various in-run calibration and compensation algorithms in the control system [16][17][18][19], further improving the bias and the scale factor performances for the gyroscopes.
The most significant property of the axisymmetric resonators is the symmetry between the two working modes of the gyroscopic operation. The efficiency of the energy transfer between the working modes primarily depends on the modal symmetry. However, the errors and imperfections in the fabrication and assembly processes are inevitable and will introduce the modal frequency split. It is essential to correct these errors after the fabrication for improving the gyroscope performances. Generally, the post-fabricated correction can be performed electrostatically or mechanically. Commonly, the electrostatic correct is actively and temporarily applied during the gyroscopic operation. Nonetheless, the mechanical correction is performed before the resonator is packaged and is permanent. Researchers have developed numerous approaches, such as mass deposition [20][21][22], laser ablation [23][24][25], ion beam trimming [14] and chemical etching [26], to mechanically correct the modal errors to diminish the frequency split. To accurately perform the mechanical correction, it is vital to find the precise position, namely, the modal orientation, where the deposition, the ablation, the trimming, or the etching should be carried out.
However, as far as we know, very few literatures particularly address the characterization of the modal orientation of the axisymmetric resonators. Oriented to the Jet Propulsion Laboratory MEMS vibratory rate gyroscope, the researchers reported an experimental system identification method to estimate the orientation of the stiffness matrix principal axes [27]. The proposed procedure includes first applying modern system identification algorithms to determine the multivariable input-output model and then the orientation of the principal axes was determined based on the identified model. Different from the methods in [27], this paper provides detailed approaches to characterize the modal frequencies and the modal orientation based upon the frequency sweep and upon the ring down method, which can be respectively performed in the frequency domain and the time domain, without any sophisticated modern system identification algorithms. The proposed approaches provide the necessary information utilized in the mechanical trimming processing.
The coming sections of this paper are organized as follows. Section 2 first reviews the dynamics of the axisymmetric resonators with respective to the displacements along the orthogonal working modes and those along the directions of capacitive forcing and pick-off, respectively. The second part of Section 2 introduces the two-dimensional transfer functions, in terms of the forcing frequency and the forcing direction, of the axisymmetric resonators. Section 3 reveals the relationships between the modal frequencies, the modal orientation, and the poles and the zero in the two-dimensional transfer function described in Section 2. Through the obtained relationships, the method based on the continuous linear frequency sweep to extract the modal characters is analyzed in detail. Section 4 gives the characterization approach based on the ring down method, which may be more effective as for resonators with high Q values. In addition, the ring down method can extract not only the modal frequencies and the orientation but also the nonlinear coefficients of the stiffness. Section 5 demonstrates the experimental verification of the methods discussed in Section 3 and Section 4 by means of a ring resonator. Finally, Section 6 briefly concludes this paper.

Equations of Motion of Axisymmetric Resonators' Working Modes
The mode shape with the wave number of n = 2, shown in Figure 1 where a micro ring resonator and a micro hemispherical resonator are taken as examples, is commonly implemented as working modes in MEMS axisymmetric gyroscopes. From Figure 1, we can observe that there are two orthogonal modes for the mode shape of n = 2. The angle between the principal axes of the working modes is 45 • in geometry, along the directions of p and q.
The fundamental equations of motion of the working modes can be described as two coupled oscillators where p(t) and q(t) are displacements of the oscillator along p and q, k is an angular gain factor determined by the geometry of the oscillator, Ω is the applied angular rate, ω 1 and ω 2 are the modal frequencies of the working modes. The external forces and the damping forces are not presented in Equation (1). For convenience, we assign ω 1 ≤ ω 2 in this paper. According to different control approaches for the working modes, CVGs with axisymmetric resonators can operate under the force-to-rebalance (FTR) scheme or the whole angle (WA) scheme. Under the FTR scheme, one of the working modes (called the primary mode) will be driven into resonance with a constant vibration amplitude. Once the gyroscope is subject to any rotation, the energy will transfer from the primary mode to the other (the secondary mode), making the latter vibrate as well. The control loop picks off the vibration induced by the Coriolis effect and forces the secondary mode stationary again. The force utilized to rebalance the secondary mode is taken as the measurement of the angular rate information.
Under the WA scheme, both of the working modes participate in in-phase oscillations, forming a standing wave whose azimuth will proportionally rotate with the input rotation. Therefore, the azimuth of the standing wave can be a measure of the angle rotated.
For both the FTR scheme and the WA scheme, the resonator is expected perfectly axisymmetry, presenting the mode-matching condition, namely, ω 1 = ω 2 . However, as a result of structure imperfections introduced by fabrication errors and material defects, the parameters between the two working modes slightly differ from each. The imperfections will cause the azimuth of the working modes' principal axes towards a specific direction, θ ω , which is regularly misaligned with the center of the discrete electrodes, as illustrated in Figure 2a. In this paper, θ ω represents the orientation of the principal axis with the lower frequency, namely ω 1 , and is defined within the range from −π/4 to π/4. Although the motions of the resonator can be simply decomposed into the vibrations along the principal axes, only the displacements along the pick-off electrodes can be readout effectively. Figure 2b shows the relationship between the displacements along the principal axes (p and q) and those along the electrode axes, 0 • (x) and 45 • (y). It should be noted that the azimuth of the principal axis in the Cartesian coordinate is twice of that in the geometry, due to the mapping of the coordinate y from 45 • to 90 • . Therefore, we can readily have where x(t) and y(t) are the displacements along 0 • and 45 • , and R ω is a rotation matrix defined as By combining Equation (1) and Equation (2), from the perspective of the capacitive readout and drive, we can write the equations of motion as where Ω is omitted for the purpose of parameter characterizations in this paper. If we further include the damping forces and the external electrostatic forces, the dynamics of the resonator can be described as where τ 1 and τ 2 are time constants along the damping axes, K T is the stiffness matrix introduced by the electrostatic forces, m is the effective mass of the oscillator, f x (t) and f y (t) are electrostatic forces along the x and the y electrodes, respectively. In Equation (5), R τ is a rotation matrix with respect to the azimuth of the damping axes and is in the form of where θ τ is the azimuth of the damping axes. The electrostatic stiffness matrix can be represented as [28] where η k , in N/m/V 2 , is an efficiency coefficient from voltage to stiffness, V x, rms and V y, rms are the voltages in RMS applied on the x and the y electrodes. The electrostatic forces exerted on the electrodes are in the form of where η f , in N/V 2 , is an efficiency coefficient from voltage to force, V x (t) and V y (t) are the voltages on the electrodes. Equation (5) gives the basic dynamics of the working modes of the axisymmetric resonators in MEMS CVGs, from which the characterizations of the modal parameters can start.

Two-Dimensional Transfer Function
Equation (5) describes the axisymmetric resonator's equations of motion along the electrodes of capacitive forcing and pick-off. However, in some situations, especially when the gyroscope operates under the WA scheme, the dynamics along particular forcing directions are of more interest. Under the WA scheme, the control system of the CVG picks off the azimuth of the standing wave and exerts the electrostatic force just towards the direction of the motion to maintain a stable vibration. Meanwhile, the control system also suppresses any orthogonal vibration which undermines the standing wave. To better understand the motions of the axisymmetric resonator vibrating in different directions, we can examine the dynamics in Equation (5) through different forcing angles.
By utilizing the orthogonality of f x and f y , we can reconstruct the electrostatic forces as where f θ and f θ+π/4 are the electrostatic forces towards the forcing angles of θ and θ + π/4, and Similarly, the displacements along the directions of θ and θ + π/4 relate to x and y through By substituting Equation (9) and Equation (11) into Equation (5), we will arrive at . (12) Applying the Laplace transform to the both sides of Equation (12) yields where R θ (s), R θ+π/4 (s), F θ (s), and F θ+π/4 (s) respectively are the Laplace transforms of r θ (t), r θ+π/4 (t), f θ (t), and f θ+π/4 (t), and By setting F θ+π/4 (s) = 0 in Equation (13), we can get where T 11 (s, θ) is the transfer function from the electrostatic force f θ (t) to the displacement r θ (t), and T 21 (s, θ) actually reflects the cross-coupling from f θ (t) to the displacement along θ + π/4. T 11 (s, θ) and T 21 (s, θ) can be obtained by matrix manipulations in Equation (14) T where a(s, θ) = s 2 + 2 and To gain some insights, if we ignore the damping mismatch and assume that V x,rms = V y,rms = V rms , then Equations (16) and (17) will be reduced to where τ 0 is a time constant for the ideal damped decay. Figure 3 visualizes the magnitude-frequency responses of T 11 and T 21 as functions of both the forcing frequency, f , and the forcing angle, θ. The parameters used in Figure 3 are listed in Table 1. Figure 3 clearly exhibits the changes of the magnitude responses over different forcing angles. In T 11 , as shown in Figure 3a,b, the response will experience two peaks and one valley, except in the directions of θ + nπ/4 (0 ≤ n ≤ 7, n ∈ Z) where only one peak will occur at the frequency of ω 1 or ω 2 . The responses of cross-couplings, in Figure 3c,d, display two resonances in most directions expect in θ + nπ/4 where no couplings occur at all, which can be verified by examining the numerator of Equation (22). Also, the most pronounced difference between T 11 and T 21 is the valley between the two peaks.
The peaks and the valley in the transfer function of T 11 can be utilized to reveal the modal frequencies and the orientation of the stiffness axes, which will be discussed in Section 3.    Figure 3.

Parameters Values Units
f 1 9995 Hz f 2 10,005 Hz

Poles and Zeros in Magnitude-Frequency Response of the Transfer Function
It is clear to observe, in Figure 3, that the poles and the zero in the magnitude-frequency responses of the two-dimensional transfer function contain sufficient information about the modal frequencies and the orientation of the stiffness axes.
If the resonator is driven towards the electrodes at 0 • , the transfer function of T 11 (s, 0) can be obtained as In Equation (23), the errors in damping forces are still ignored because the damping mismatch will not significantly affect the modal frequencies and their directions. Figure 4 manifests the transfer function of T 11 (s, 0) with different damping mismatches. The Qfactor of the working mode with the frequency of ω 2 is fixed as 100,000. The corresponding poles and zeros of the magnitude-frequency responses almost stay at the same frequencies while the Q-factor of the frequency of ω 1 changes from 50,000 to 200,000.  According to the fact demonstrated in Figure 4, regardless of the impacts of the damping errors on the modal frequencies, the locations of the poles and the zero in the transfer function can be found as from which we can attain the relation between the poles, the zero, and the modal orientation as z 2 = p 2 1 sin 2 2θ ω + p 2 2 cos 2 2θ ω However, by definition, the modal orientation θ ω is between −π/4 to π/4. Hence, Equation (25) only gives the absolute value of θ ω To determine the direction of θ ω , once the absolute value of θ ω is obtained, we can further drive the resonator towards θ ω and −θ ω . The resultant magnitude-frequency response can be found as and There will be only one resonance in the response along the orientation of θ ω and will be still two resonances along −θ ω , which is depicted in Figure 5.

Continuous Linear Frequency Sweep
Voltage with continuous linear frequency sweep, also named as the linear chirp signal, defined in Equation (29), can be performed to effectively extract the poles and the zero in the transfer function.
where V x (t) is the voltage applied on the electrodes along 0 • direction, first introduced in Equation (8), V d is the direct current voltage, and where A is the amplitude of the chirp signal, T is the duration time of the sweep, ϑ(t) is the instantaneous phase, Π(x) is the standard rectangle function and is defined as The instantaneous angular frequency can be calculated as where ω L is the start angular frequency of the linear sweep, and ω H is the end angular frequency. Figure 6 illustrates the chirp signal, with the parameters listed in Table 2, in the time-frequency domain. The instantaneous frequency of the signal linearly raises from 8 kHz to 12 kHz within 100 s. Table 2. Parameters of the chirp signal shown in Figure 6.   Table 2.

Parameters
By omitting the low amplitude oscillations and assuming that t T, we can approximate the autocorrelation function of Equation (30) as [29] where W is the signal bandwidth in Hz, ω c is the central angular frequency. We have the following definitions: The Fourier transform of Equation (33) gives the energy spectral density of the chirp signal V a (t) where Π(x) is defined in Equation (31). Figure 7 demonstrates the approximated energy spectral density of the chirp signal from Equation (34) and the exact curve from the numerical calculation. Except the oscillatory ripples near the edge frequencies, Equation (34) can retain the essential properties of the chirp signal's energy spectral density. To ensure the poles and the zero in the magnitude-frequency response can be covered, we set the relationship between the frequencies as The amplitude spectral density of the displacement along 0 • direction can be acquired by Equation (36) indicates that, if ignoring the ripples at the edge frequencies, the amplitude spectral density of the chirp-driven resonator's outcome just gives the curve of the frequency response. It should be noted that the strength of the resultant amplitude spectral density can be adjusted by the chirp input under the law of A √ T/W. On the other hand, the resultant frequency resolution in Φ x (j2π f ) is inversely proportional to the sweep time T. Figure 8 presents the simulated response to the chirp excitation in the time domain and in the frequency domain, respectively. The numerically simulated amplitude spectral density and that calculated from Equation (36) are also compared in Figure 8b, showing the full consistence. Once the peaks and the valley in the transfer function of T 11 (s, 0) are extracted through the chirp excitation along 0 • , the absolute value of the modal orientation can be derived from Equation (26). The sign of the orientation can be conformed by the chirp excitation towards θ ω and −θ ω , from which the magnitude responses of T 11 (s, θ ω ) and T 11 (s, −θ ω ) can be computed.
However, for resonators vibrating with considerable amplitudes, the stiffness force will exhibit the nonlinear property. The pronounced consequence of the nonlinearity on the frequency sweep is the hysteresis in the magnitude-frequency response, presenting sharp jumps at certain frequencies and, thus, making the positions of the poles indistinct in the frequency domain.

Ring-Down and Nonlinear Stiffness Coefficient
The modal frequencies, the orientation, and even the nonlinear stiffness coefficient, are hidden in the free oscillation signals of the resonator. First we introduce the undamped free-oscillation of the resonator along the p and q modes where A 1 and A 2 are the amplitudes of p and q modes' oscillations, A 1,2 ≥ 0, and φ is the initial phase difference between them. By combining Equation (37) and Equation (2), we get the undamped free-vibrations along the forcer electrodes and the readout electrodes as There distinctly exhibits beats in X(t) and Y(t) since ω 1 and ω 2 are pretty close, as demonstrated in Figure 9. To study the envelope of the beat signal, we can utilize the concept of the analytic signal which is defined as the Hilbert transform of X(t). By the definition in Equation (40), the Hilbert transform of X(t) can be calculated as Therefore, the envelope of X(t) can be attained by (42) Figure 9 also shows the Hilbert transform of X(t) and the magnitude of the analytic signal X A (t), which exactly is the envelope of X(t). Equation (42) implies that the fundamental frequency of the signal envelope reflects the modal frequency split, ω 2 − ω 1 . Similarly, the envelope of Y(t) can be expressed by By observing the beats in both X(t) and Y(t), demonstrated in Figure 10, we can find that, if θ ω > 0, when t 1 satisfies sin(ω 1 t 1 ) = 1 and sin(ω 2 t 1 + φ) = −1, E X (t) will reach the peak, that is A 1 cos 2θ ω + A 2 sin 2θ ω = X(t 1 ).
Meanwhile, at t = t 1 the value of Y(t) will be where |Y(t 1 )| is indeed the minimum of E Y (t) but the sign of Y(t 1 ) should be determined by checking the instantaneous phase of Y(t). On the other hand, when t = t 2 satisfies sin(ω 1 t 2 ) = 1 and sin(ω 2 t 2 + φ) = 1, E Y (t) will reach the peak and X(t) is The constraints described by Equations (44)-(47) are sufficient for resonators with extremely high Q-factor, whose time constant of the ring-down is large enough to perform the above analysis. However, if the ring-down of the resonators decays notably, the peak and the valley values of the beats are not easy to maintain, as shown in Figure 11.
The envelope of the ring-down beats, calculated from the magnitude of the analytic signal, are also depicted in Figure 11, where the exponential fitting curve of the envelope is also presented. It should be noted that the start point and the end point of the envelope exist serious errors due to the end effect of the Hilbert transform. Therefore, the fitting processing of the ring-down should exclude the start and the end of the envelope.
Based on the fitting curve we can obtain the time constant of the ring-down for both x(t) and y(t) as τ x and τ y , respectively. Then, the undamped signals of X(t), Y(t), and their corresponding envelopes can be obtained by According to the relations in Equation (48), the peak and the valley values of the beats can be acquired as the undamped cases. Once the essential values are determined effectively, we rewrite Equations (44)-(47) as where ε 1 to ε 4 can be interpreted as the residual errors, and the case for θ ω < 0 is also considered.
By minimizing the sum of the residual errors we can find the estimation of the modal orientation aŝ where Similar to the method of the frequency sweep, the value of θ ω still cannot be figured out by one single constraint. However, by virtue of the obtained θ 1 , we can rotate the ring-down vibrations along xy to those along pq through Equation (2), as illustrated in Figure 12.
Considering the rotated ring-down oscillation is a single tone, almost is monocomponent, along the principal axis, we can utilize the Hilbert transform again to extract the precise frequencies. We construct the analytic signal of p(t) as well as Equation (39) where E p (t) is the instantaneous amplitude of p(t), and ψ 1 (t) is the corresponding instantaneous phase, which relates to the instantaneous frequency, ω 1 (t), by Meanwhile, the retrieved instantaneous frequency against the instantaneous amplitude can reveal the nonlinearity in the stiffness. Figure 13 illustrates the evolution of the amplitude and the instantaneous frequency with respect to the time, showing the frequency-amplitude dependency in the ring-down oscillation. A resonator with the nonlinear stiffness commonly can be described as a duffing oscillatorp where γ n is the coefficient of the n-th order nonlinearity. Generally, from observation the most dominate nonlinearity is the cubic term in Equation (54), which will lead to the magnification of the force amplitude to the displacement amplitude roughly as Equation (55) shows the strong dependency between the magnification M and the vibration amplitude E p . The resonant frequency can be approximated as Generally, the relation between the amplitude E p , in Volts, and the resonant frequency f r , in Hz, can be obtained through experiments, then we can fit the relation as where a and b are the corresponding coefficients of the fitting. Therefore, the nonlinear coefficient of the cubic stiffness term can be estimated aŝ where K pre is the gain of the frond-end circuit which converts the displacement of the resonator to the measured voltage.

Experiments
The methods discussed in Sections 3 and 4 were verified by a vacuum-packaged ring resonator shown in Figure 14. The resonator was fabricated by a standard silicon-on-glass processing, which was previously reported in [8]. The chirp forcing voltage for the frequency sweep was generated by a digital-to-analog convertor (DAC) of AD5791 from driven by a Cyclone ® IV FPGA of EP4CE from Altera ® . The readout interface circuit was based on a ring-diode capacitance-to-voltage convertor discussed in [30]. The readout signals were acquired by a dynamic signal analyzer, PXI-4461, from National Instruments.
The DAC output a chirp signal with the start frequency of 7280 Hz and the end frequency of 7310 Hz at the duration time of 30 s. The analyzer sampled the signals from the chirp-driven resonator at a rate of 204.8 kHz. The resultant amplitude spectral density is presented in Figure 15a, giving that p 1 = 7288.9 Hz, z = 7290.7 Hz, and p 2 = 7299 Hz. From Equation (26), we can conclude that |θ ω | ≈ 32.52 • . To further determine the sign of θ ω , we additionally excited the resonator along 32 • and −32 • , respectively. The obtained spectrum was demonstrated in Figure 15b, which confirmed that θ ω ≈ −32.52 • .
The resonator was also characterized by the ring-down method, showing in Figure 16. The resonator was first excited with a phase-locked loop and then released by the loop, yielding a free-oscillation of the working modes. The readout signals from the xy directions, x(t) and y(t), were simultaneously sampled and, then, filtered by a narrow bandpass filter to remove the disturbance and noise.
The envelopes of the filtered x(t) and y(t) were obtained via the magnitudes the analytic signals. The time constants of the ring-down decays were given by the exponentially fitted curves of the envelopes: τ x = 7.652 s and τ y = 7.662 s. By applying Equation (48), we flattened the ring-down signals and attained the peak and the valley values of X(t) and Y(t) as X(t 1 ) = 0.86812, Y(t 1 ) = 0.27495, Y(t 2 ) = 0.46726, and X(t 2 ) = 0.77464, from which θ ω would be estimated as 12.17 • or −32.83 • .
Once we got the candidates of the modal orientation, we tentatively rotated the oscillations of x(t) and y(t) by −32.83 • , resulting in p(t) and q(t) presented in Figure 17.
To utilize the method described by Equations (52) and (53), the variational mode decomposition (VMD) was respectively performed to p(t) and q(t) to verify the property of monocomponent. It can be readily found that, in both p(t) and q(t), there was only one VMD component whose amplitude significantly greater than the others, suggesting the Hilbert transform can be effectively applied.
Finally, the instantaneous frequencies of the ring-down signals were extracted by the numerical differences of the instantaneous phase of the analytic signals constructed by the Hilbert transform, as shown in Figure 18. The attained frequency of p(t) was lower than that of q(t), confirming the modal orientation was −32.83 • . Otherwise, if we rotated x(t) and y(t) by 12.17 • in the first place, the attained frequency order would be opposite, contradicting with the definition.  From Figure 18, the modal frequencies of the working modes were measured as 7288.8 Hz and 7298.9 Hz. In addition, Figure 18 remarkably demonstrated the quadratic relationship between the vibration magnitude and the resonant frequency. From the fitting curves, together with Equation (58), we can calculate the nonlinear coefficients as γ 31 /K 2 pre = 2 × 10 −3 V −2 for p mode and γ 32 /K 2 pre = 2.67 × 10 −5 V −2 for q mode. The measured modal frequencies and the orientation by the frequency sweep and the ring-down analysis are summarized in Table 3, where the measured parameters demonstrated good consistence between the two approaches.

Conclusions
This paper demonstrated the frequency sweep and the ring down approaches to characterize the modal frequencies and the modal orientation of the axisymmetric resonators utilized in the Coriolis vibratory gyroscopes. In the frequency sweep method, the resonator was excited by chirp voltages and the responses were transformed into the frequency domain, within which the zero and the poles imply the modal information of interest. As for resonators with extraordinarily high Q values, it was more convenient to adopt the ring down method, which extracted the modal characterizations, including the modal frequencies, the modal orientation, and the nonlinear stiffness coefficient, by the peaks and the valleys of the beat signals. The results from the experiments were coincident with each other in the two approaches.

Conflicts of Interest:
The authors declare no conflict of interest.