Design and Verification of a Digital Controller for a 2-Piece Hemispherical Resonator Gyroscope

A Hemispherical Resonator Gyro (HRG) is the Coriolis Vibratory Gyro (CVG) that measures rotation angle or angular velocity using Coriolis force acting the vibrating mass. A HRG can be used as a rate gyro or integrating gyro without structural modification by simply changing the control scheme. In this paper, differential control algorithms are designed for a 2-piece HRG. To design a precision controller, the electromechanical modelling and signal processing must be pre-performed accurately. Therefore, the equations of motion for the HRG resonator with switched harmonic excitations are derived with the Duhamel Integral method. Electromechanical modeling of the resonator, electric module and charge amplifier is performed by considering the mode shape of a thin hemispherical shell. Further, signal processing and control algorithms are designed. The multi-flexing scheme of sensing, driving cycles and x, y-axis switching cycles is appropriate for high precision and low maneuverability systems. The differential control scheme is easily capable of rejecting the common mode errors of x, y-axis signals and changing the rate integrating mode on basis of these studies. In the rate gyro mode the controller is composed of Phase-Locked Loop (PLL), amplitude, quadrature and rate control loop. All controllers are designed on basis of a digital PI controller. The signal processing and control algorithms are verified through Matlab/Simulink simulations. Finally, a FPGA and DSP board with these algorithms is verified through experiments.


Introduction
A Hemispherical Resonance Gyroscope (HRG) is a Coriolis Vibrating Gyroscope (CVG), which measures the angle or angular velocity using the Coriolis force generated by the rotational motion [1]. HRGs are suitable for miniaturization and high precision due to their simple structure made up of five components and easy production process. Particularly, since it is a sensor made of quartz material, which has excellent material properties, using the solid state wave phenomenon, it guarantees high reliability and a long-term lifespan. In addition, it has the advantage that it can be used in the angular velocity mode or the integral angular velocity mode according to the control technique and the electronic circuit used, without requiring any structural changes to the sensor [1].
The development of HRG technology, which started in the 1970s, has been conducted actively in the advanced countries such as the United States (Northrop Grumman), France (SAGEM), Russia (RDC, Medicon), etc. for the purpose of its use in the attitude controllers of satellites and spaceships, long-term navigation systems for strategic missiles, submarines, etc., where long-term reliability is important [2].
The United States (Northrop Grumman) is currently producing miniature HRGs, the size of golf balls, based on the existing HRG 130P from 2012. This is the results of design improvement that converts the 3-piece system, in which the forcer and pickoff are divided, into a 2-piece system, in which the outer forcer is removed, and the process stage is reduced using the multi-flexing method [3]. The multi-flexing method involves obtaining the sensing and driving signal of the x-axis and y-axis using the same element by switching within the specified cycle and is suitable for high precision and low maneuverability systems like satellite systems. In addition, the SIGMA 20 navigation system, according to data made public by France in 2013, is under development based on the 2-piece HRG, which applies the multi-flexing technique to a resonator of 20 mm diameter, showing that it can be utilized in ground navigation systems used in rough environments focusing on the characteristics of the HRG, which can endure thermal and mechanical stress [4].
In the meantime, the traditional control method controls the x-axis to be major axis of the elliptical trajectory expressing the pendulum variables from the axis where the excitation and sensing electrodes are arranged. On the contrary, the differential control method controls the axis which forms a 45a ngle with the x-axis, to be the major axis of the elliptical trajectory, This method has the advantage that it can remove the common mode error for both the x-axis and the y-axis and is easy to switch to the integral angular velocity mode [5].
As such, next generation HRGs will evolve into representative gyros, which can materialize the demands for subminiature structure, high precision and high reliability by the development of 2-piece systems applying multi-flexing methods, differential control algorithms, etc. To develop such subminiature and high precision HRGs, advanced core processes such as the production of low-loss hemispherical resonators and electrode blocks, low-stress heterojunctions, the balancing and tuning, high degrees of vacuum packaging, etc. and electronic module technologies such as low-noise pre-amplifiers, FPGA-based high precision digital signal processing and control circuits, error modelling and compensation techniques, etc. must be developed.
Therefore, in this article, electromechanical modelling will be performed on a 2-piece system equipped with a multi-flexing technique and a signal processing and control algorithm design based on that. In Section 2, the resonator motion equation with continuous harmonic excitation will be deduced and the major electromechanical gains between the resonator and the electrodes calculated through the modeling. At this time, the mode shape of the resonator will be considered. In addition, the equation of the resonator motion for case that the switched harmonic excitation is applied by the multi-flexing will be induced, and its results will be compared with the case of continuous harmonic excitation. In Section 3, the signal processing and control algorithm will be designed based on the electromechanical modelling results. In this section, the signal processing algorithm based on the multi-flexing method and the differential control algorithm in the rate gyro mode (or FTR mode) will be designed. The designed control algorithm will be tested through Matlab/Simulink SW and the design results will be verified finally by comparing the results of an actually made sensor with the simulation results through suitable experiments.

Electromechanical Modeling of a HRG with Switched Harmonic Excitations
As mentioned in the previous section, a HRG is a sensor to measure the input angular velocity using the precession motion of the elastic standing wave caused by a Coriolis force. The circular section of the hemispherical resonator in the non-vibration state repeats the circle, horizontal elliptical shape, circle and vertical elliptical shape in the secondary mode [6]. The location of the maximum amplitude and the non-vibration location are referred to as the antinode and node, respectively, and generally, the equation of HRG motion is induced through the secondary resonance mode having two nodes and antinodes [7]. Basically, such a HRG system can be modeled with the secondary spring damper system and due to causes such as mass unbalance, etc. anisoelasticity errors and damping mismatch occur, which cause the frequency and Q-factor to split, respectively, which are the most important causes of errors of HRG sensors [8].

Full Equations of Motion with Harmonic Excitations
The general full equations of motion of a HRG considering the influences of frequency coupling and damping imperfection are as follows [8]: where Ω is the angular velocity of system about the vertical axis of x-axis and y-aixs, ω is the mean resonant frequency, k is the Brian coefficient (~0.3), , θ ω is the angle of the unbalance between the x-axis and the major axis of resonance mode, θ τ is the angle between the x-axis and the major axis of the linear damper (readers should refer to Appendix A for further details).
To calculate the nominal amplitude and the time constant of each axis, let's assume that there is no frequency coupling by the damping mismatch and the mass unbalance in the above equation and no input angular velocity. Then, the motion characteristics for the each axis in the above equation can be interpreted as the equation of motion of the secondary spring damper system below: ..
where m n is modal mass, τ n = 2Q/ω, Q is the quality factor. In this moment, the excitation force to be applied continuously can be modeled with the harmonic function as follows: where d f dv is the change of control force per unit voltage, v c ptq " v c cosω t is ac control voltage. With the aid of Euler's formula, the solution, x ptq can be written in the form [9]: x ptq " f 0 Q m n ω 2¨s in ωt´1 c 1´´1 2Q¯2 e´ω 2Q t sinω d t‹ ‹ ‚ where f 0 " dC dx V B v c , V B is bias voltage, dC dx is the change of capacitance per unit displacement, Since ω dω with the assumption that the quality factor is high, the equation above can be arranged as follows. The variables needed for deriving the equation are summarized in Table 1 and it can be schematized as shown in Figure 1.  In Equation (5), as

The Nominal Amplitude and Time Constant of HRG with Switched Harmonic Excitations
In the case of a 2-piece HRG system, it switches the sensing and driving cycles and the x-axis and y-axis using the common element, so the excitation force in Equation (3) is not given continuously   In Equation (5), as 1 b 1´p1{2Qq 2 -1 can be assumed, when expressing the amplitude as x ptq " A ptq sin pωt`φq, A ptq " A n´1´e´t {τ n¯, the nominal amplitude A n and the time constant τ n can be calculated as follows: ω " 318.310 s (7)

The Nominal Amplitude and Time Constant of HRG with Switched Harmonic Excitations
In the case of a 2-piece HRG system, it switches the sensing and driving cycles and the x-axis and y-axis using the common element, so the excitation force in Equation (3) is not given continuously but as much as the time determined within the certain cycle, which can be expressed as Equation (8) as follows: where T 0 " 2π ω s, C 1 " 5 is the total operation cycles, including sensing and driving cycles, C 2 " 1 is the x-axis control cycle.
To obtain the system time response characteristics when the non-continuous excitation force is given, the Duhamel integral (or convolution integral) method is used [9]. The Duhamel integral method is the special form of integral to be applied when obtaining the output signal in a linear system if the input signal and the system impulse response are given: x ptq " f ptq˚h ptq " We get the unit impulse response of a viscously damped SDOF system. By convention, the unit impulse response function is frequently called h ptq [9]: Let's calculate x ptq by inserting above Equations (8) and (10) into Equation (9). Since f ptq " 0 when t ă 0, x ptq " 0. In the section of C 1 T 0 pn´1q ă t ď C 2 T 0`C1 T 0 pn´1q, n " 1, 2,¨¨¨, 8, as the excitation and the unit impulse response show the phase difference as Figure 2a, and in the section of C 2 T 0 pn´1q`C 1 T 0 ă t ď 2C 1 T 0 n, n " 1, 2,¨¨¨, 8, they show the phase difference as Figure 2b, the integral can be arranged as follows: Case a. C 1 T 0 pn´1q ă t ď C 2 T 0`C1 T 0 pn´1q , n " 1, 2,¨¨¨, 8 Case b. C 2 T 0 pn´1q`C 1 T 0 ă t ď 2C 1 T 0 n, n " 1, 2,¨¨¨, 8 x ptq " Sensors 2016, 16, 555 5 of 22 but as much as the time determined within the certain cycle, which can be expressed as Equation (8) as follows: where = s, = 5 is the total operation cycles, including sensing and driving cycles, = 1 is the x-axis control cycle. To obtain the system time response characteristics when the non-continuous excitation force is given, the Duhamel integral (or convolution integral) method is used [9]. The Duhamel integral method is the special form of integral to be applied when obtaining the output signal in a linear system if the input signal and the system impulse response are given: We get the unit impulse response of a viscously damped SDOF system. By convention, the unit impulse response function is frequently called h( ) [9]: Let's calculate x( ) by inserting above Equations (8) and (10) into Equation (9). Since f( ) = 0 when < 0, x( ) = 0. In the section of ( − 1) < ≤ + ( − 1), = 1, 2, ⋯ , ∞, as the excitation and the unit impulse response show the phase difference as Figure 2a, and in the section of ( − 1) + < ≤ 2 , = 1, 2, ⋯ , ∞, they show the phase difference as Figure 2b, the integral can be arranged as follows: Case a.
(a) (b) Assuming that x ct ptq " ş t C 1 T 0 pn´1q f pτq h pt´τq dτ, x pt ptq " f pτq h pt´τq dτ, the calculation results are as follows: as in this moment, the Q-factor is high, we can assume that ω dω: it can be arranged as follows: x pt ptq " Refer to Appendix B for further details. If Equations (14) and (15) are inserted to Equations (11) and (12), the time response of SDOF spring-damper system by the switched excitation can be arranged as follows: The time response characteristics of the mean amplitude by the switched excitation are as follows: Refer to Appendix C for further details. Accordingly, when C 1 " 5, C 2 " 1 , the nominal amplitude and the time constant are calculated as follows: Through the above results, it is observed that the amplitude can be obtained by multiplying the rate of the driving cycle against the entire operation cycle compared with the amplitude when the continuous excitation is applied, and the time constants are almost same. This is the important input data for electromechanical gains and design of controller. The amplitude results are schematized as shown in Figure 3.

Verification of the Analytic Results through Simulations
To test the above interpretation results, the Matlab/Simulink SW was prepared as shown in Figure 4 and the simulation was performed. It is observed that the results are as same as Figure 5 and is identical when comparing with Figure 3.

Verification of the Analytic Results through Simulations
To test the above interpretation results, the Matlab/Simulink SW was prepared as shown in Figure 4 and the simulation was performed. It is observed that the results are as same as Figure 5 and is identical when comparing with Figure 3.

Verification of the Analytic Results through Simulations
To test the above interpretation results, the Matlab/Simulink SW was prepared as shown in Figure 4 and the simulation was performed. It is observed that the results are as same as Figure 5 and is identical when comparing with Figure 3.

Electromechanical Modeling between Resonator and Electrodes
Before designing the signal processing and control circuit, the electromechanical modeling between resonator and electrodes must be done to estimate the capacitance change for different electrostatic forces.
The gap between the resonator and the sensing (or driving) electrode block cannot be assumed as parallel plate simply. For the resonator, the 2nd vibration shape corresponding to a standing wave in the thin hemispherical shell should be considered. Therefore, in this article, the capacitance change and the changes in the electrostatic force are calculated considering the mode shape and modal force.
According to Rayleigh, the mode equation for the mode shape analysis among the second vibration mode equations of the thin hemispherical shell of the resonator can be represented [10]: (2 + cos ) tan 2 cos 2( − ) sin ( − ) where A, B is the 1st and 2nd wave amplitude, is the orientation of the wave relative to the resonator. w, α, φ is clear from Figure 6.

Electromechanical Modeling between Resonator and Electrodes
Before designing the signal processing and control circuit, the electromechanical modeling between resonator and electrodes must be done to estimate the capacitance change for different electrostatic forces.
The gap between the resonator and the sensing (or driving) electrode block cannot be assumed as parallel plate simply. For the resonator, the 2nd vibration shape corresponding to a standing wave in the thin hemispherical shell should be considered. Therefore, in this article, the capacitance change and the changes in the electrostatic force are calculated considering the mode shape and modal force.
According to Rayleigh, the mode equation for the mode shape analysis among the second vibration mode equations of the thin hemispherical shell of the resonator can be represented [10]: w " A 2 p2`cosαq tan 2 α 2 cos r2 pϕ´ϕ 0 qs sin rω pt´t 0 qs B 2 p2`cosαq tan 2 α 2 sin r2 pϕ´ϕ 0 qs cos rω pt´t 0 qs where A, B is the 1st and 2nd wave amplitude, ϕ 0 is the orientation of the wave relative to the resonator. w, α, ϕ is clear from Figure 6.

Electromechanical Modeling between Resonator and Electrodes
Before designing the signal processing and control circuit, the electromechanical modeling between resonator and electrodes must be done to estimate the capacitance change for different electrostatic forces.
The gap between the resonator and the sensing (or driving) electrode block cannot be assumed as parallel plate simply. For the resonator, the 2nd vibration shape corresponding to a standing wave in the thin hemispherical shell should be considered. Therefore, in this article, the capacitance change and the changes in the electrostatic force are calculated considering the mode shape and modal force.
According to Rayleigh, the mode equation for the mode shape analysis among the second vibration mode equations of the thin hemispherical shell of the resonator can be represented [10]: where A, B is the 1st and 2nd wave amplitude, is the orientation of the wave relative to the resonator. w, α, φ is clear from Figure 6.  The mode shape calculated by inserting ϕ 0 " 0, A " 1, B " 0 to Equation (21) is as follows: The electrostatic capacitance C 0 and capacitance changes by the displacement dC dx in the sensing electrode are induced as follows using above equation: Table 2 for the rest of the variables. If Equation (22) is inserted to Equation (24), it is as follows: If Equations (26) and (27) is inserted to Equation (25), it is as follows: The calculating the control force f by the control voltage v ac applied to the electrode block is as follows: Equations (26) and (27) are inserted to Equation (29), it can be arranged as follows: The equation related to the amplitude of resonator by the control force is as follows: The equation of the change amplifier output by the change in the amplitude is as follows: By inserting the values in Table 2 into Equations (23), (28), (30)-(32), the electromechanical gains in Table 3 can be obtained:

Design Variables Values
Radius of Resonator R " 15.3 mm Nominal gap between resonator and electrode block d 0 " 120 µm Dielectric permittivity ε " 8.85ˆ10´1 2 F{m Electrodes angles in azimuth ϕ 1 "´18˝, ϕ 2 " 18E lectrodes angles in elevation α 1 " 78.2˝, α 2 " 90B ias voltage V bias " 200 V Quality factor Q " 7ˆ10 6 Modal mass m n " 0.85ˆ10´3 kg Resonant frequency ω " 4.4ˆ10 4 rad{s Feedback capacitance of the charge amp C f " 22 pF If 5-cycles operation multi-flexing method, where C 1 " 5 , C 2 " 1, is applied and four electrode blocks are assigned to the signals of x-axis and y-axis, the control force, the amplitude and the output voltage of charge amplifier by control voltage of 100 mVac are follows: In the measurement results after making actual sensor, the amplitude by the control voltage of 100 mV was 1.02 µm and the output voltage of the charge amplifier was 143 mV showing 5.1% and 6.9% of differences, respectively, which proves that the design results are valid when considering the error in the process and measurement. Based on the electromechanical gains mentioned above, the proper signal processing and control circuit will be designed.

Design of the Signal Processing and Control Algorithm
A HRG can be operated as an angular velocity sensor or an angular sensor with a single sensor without structural change. In this study, FTR mode (or rate gyro mode), which has dynamic range limitations but excellent noise and resolution characteristics, will be handled. FTR mode is driven by the closed loop. In the FTR mode, the driving-axis is excited to remain as the reference amplitude and the sensing-axis generated by input of angular velocity controls the amplitude to be 0. In this moment, the force required to remove the sensing-axis vibration is referred to as the rebalance force and as it is proportional to the angular velocity input, the angular velocity input is estimated by multiplying this force by the conversion scale factor [11].
As shown in Figure 7, if the analog voltage output comes out from the pre-amplifier, it is transmitted to the FPGA in the form of a digital signal through the filter and ADC. It is converted to the in-phase, quadrature signals of x-and y-axis through the demodulation in the FPGA. These signals are calculated from the pendulum variables required for control and digital PI control outputs for amplitude, quadrature, rate and phase control are generated in the DSP. If the control outputs are delivered to the FPGA, the driving signal is generated by the modulation, which becomes an analog control voltage through the DAC controlling the HRG. In addition, in the DSP, it was designed to output the angular velocity by estimating the input angular velocity with rate control output. signals are calculated from the pendulum variables required for control and digital PI control outputs for amplitude, quadrature, rate and phase control are generated in the DSP. If the control outputs are delivered to the FPGA, the driving signal is generated by the modulation, which becomes an analog control voltage through the DAC controlling the HRG. In addition, in the DSP, it was designed to output the angular velocity by estimating the input angular velocity with rate control output.

Design of the Signal Processing Algorithm
In this study, the ADC sampling frequency was designed with 512 times the resonance frequency. If the FPGA clock frequency is assumed as 2 for the convenience and the sensing and control timing diagram is schematized as shown in Figure 8 and Table 4. In the actual algorithm, = ( is an integer). The reference phase ϕ is assumed to be 0, which means that sin , cos are used for the demodulation signal. In the time diagram, [i] indicates i-th operation period. [i-1] indicates previous operation period. If x( ) is expressed as cos + sin and y( ) is expressed as cos + sin , Direct Digital Synthesis (DDS) is used to obtain , , and . As the phase of DDS and the demodulation reference phase are same, the cosine and sine output of the DDS is multiplied by x( ) as is. As y( ) is reversed in the signal intended for demodulation originally and the phase of DDS has difference of π (180°) from the demodulation reference phase, it is multiplied to y( ) as is like x( ).

Design of the Signal Processing Algorithm
In this study, the ADC sampling frequency f s was designed with 512 times the resonance frequency. If the FPGA clock frequency f clk is assumed as 2 f s for the convenience and the sensing and control timing diagram is schematized as shown in Figure 8 and Table 4. In the actual algorithm, f clk " n s f s (n s is an integer). The reference phase φ is assumed to be 0, which means that sinωt, cosω t are used for the demodulation signal. In the time diagram, [i] indicates i-th operation period. [i-1] indicates previous operation period. If x ptq is expressed as c x cosωt`s x sinωt and y ptq is expressed as c y cosωt`s y sinωt, Direct Digital Synthesis (DDS) is used to obtain c x , s x , c y and s y . As the phase of DDS and the demodulation reference phase are same, the cosine and sine output of the DDS is multiplied by x ptq as is. As y ptq is reversed in the signal intended for demodulation originally and the phase of DDS has difference of π (180˝) from the demodulation reference phase, it is multiplied to y ptq as is like x ptq. signals are calculated from the pendulum variables required for control and digital PI control outputs for amplitude, quadrature, rate and phase control are generated in the DSP. If the control outputs are delivered to the FPGA, the driving signal is generated by the modulation, which becomes an analog control voltage through the DAC controlling the HRG. In addition, in the DSP, it was designed to output the angular velocity by estimating the input angular velocity with rate control output.

Design of the Signal Processing Algorithm
In this study, the ADC sampling frequency was designed with 512 times the resonance frequency. If the FPGA clock frequency is assumed as 2 for the convenience and the sensing and control timing diagram is schematized as shown in Figure 8 and Table 4. In the actual algorithm, = ( is an integer). The reference phase ϕ is assumed to be 0, which means that sin , cos are used for the demodulation signal. In the time diagram, [i] indicates i-th operation period. [i-1] indicates previous operation period. If x( ) is expressed as cos + sin and y( ) is expressed as cos + sin , Direct Digital Synthesis (DDS) is used to obtain , , and . As the phase of DDS and the demodulation reference phase are same, the cosine and sine output of the DDS is multiplied by x( ) as is. As y( ) is reversed in the signal intended for demodulation originally and the phase of DDS has difference of π (180°) from the demodulation reference phase, it is multiplied to y( ) as is like x( ).  Figure 8. HRG sensing and control signal time diagram.  Since x ptq and y ptq are measured discontinuously by the multi-flexing, c x , s x , c y and s y cannot be calculated using conventional LPF. Therefore they will be calculated averaging the sensing signal as in the equations below: where N s is the number of FPGA clocks for sampling cycle (= 512). It was designed that the phase delay φ corr of the electronic circuit is compensated in the modulation stage of control command calculated through the PI controller. DDS_Mod cosine and sine outputs correspond to cos pωt`φ corr q and´sin pωt`φ corr q, respectively. The internal frequency of FPGA is generated in the Numerical Controlled Oscillator (NCO). This NCO and the block including the sine and cosine wave Look-Up Table (LUT) are referred to as the DDS. DDS, which generates the signal required for demodulation and modulation, is composed of phase accumulator, phase quantizer and LUT as shown in Figure 9. In Figure 9, the equation related of DDS clock frequency f clk , accumulator bit number N and phase increment value ∆θ and the calculus of the DDS frequency resolution ∆f are as follows: f out " f clk ∆θ 2 N Hz, ∆f " f clk 2 N Hz (38) Since x( ) and y( ) are measured discontinuously by the multi-flexing, , , and cannot be calculated using conventional LPF. Therefore they will be calculated averaging the sensing signal as in the equations below: where is the number of FPGA clocks for sampling cycle (= 512). It was designed that the phase delay of the electronic circuit is compensated in the modulation stage of control command calculated through the PI controller. DDS_Mod cosine and sine outputs correspond to cos( + ) and − sin( + ) , respectively. The internal frequency of FPGA is generated in the Numerical Controlled Oscillator (NCO). This NCO and the block including the sine and cosine wave Look-Up Table (LUT) are referred to as the DDS. DDS, which generates the signal required for demodulation and modulation, is composed of phase accumulator, phase quantizer and LUT as shown in Figure 9. In Figure 9, the equation related of DDS clock frequency , accumulator bit number N and phase increment value Δθ and the calculus of the DDS frequency resolution Δf are as follows:

=
, Δf = (38) The data stored in LUT are sin and cos and n is the integer satisfying ∈ 0, 2 − 1 . For DDS output for demodulation, the phase offset is applied and it was designed to add obtained by below acquisition equation to the address bit number M: Figure 10 shows the entire motion simulation results of FPGA signal processing algorithm explained so far. The data stored in LUT are sin´n 2π 2 M¯a nd cos´n 2π 2 M¯a nd n is the integer satisfying n P " 0, 2 M´1 ‰ . For DDS output for demodulation, the phase offset is applied and it was designed to add N corr obtained by below acquisition equation to the address bit number M: Figure 10 shows the entire motion simulation results of FPGA signal processing algorithm explained so far.

Design of the Control Algorithm
The control algorithm is designed to pursue the ideal gyro motion without flaw. In this design, the definition of IEEE Std 1431 Annex B (D. Lynch) is used [12]. In the case of a shell resonator, which has no damping and the elasticity is axially symmetrical, since E and Q are not changed regardless of pattern angle θ and orbit phase , these two invariants are used as amplitude and quadrature control variable, respectively. In the FTR mode, S and Q are used as rate and phase control variable, respectively. The control variables are shown in Table 5, the control commands are PI controller outputs and the unit is bits. The relation between the demodulated signal ( , , , ) and the control variables in Table 5 is follows [12]:

Design of the Control Algorithm
The control algorithm is designed to pursue the ideal gyro motion without flaw. In this design, the definition of IEEE Std 1431 Annex B (D. Lynch) is used [12]. In the case of a shell resonator, which has no damping and the elasticity is axially symmetrical, since E and Q are not changed regardless of pattern angle θ and orbit phase φ 1 , these two invariants are used as amplitude and quadrature control variable, respectively. In the FTR mode, S and Q are used as rate and phase control variable, respectively. The control variables are shown in Table 5, the control commands are PI controller outputs and the unit is bits.
The relation between the demodulated signal (c x , s x , c y , s y ) and the control variables in Table 5 is follows [12]: 2`c x s y´cy s x˘" 2aq " Q (41) c 2 x`s 2 x´c 2 y´s 2 y "´a 2´q2¯c os2θ " R 2`c x c y`sx s y˘"´a 2´q2¯s in2θ " S 2`c x s x`cy s y˘"´´a 2´q2¯s in2φ 1 " L (44) As above, if the control variables are calculated, the resonant frequency must be sought through PLL and the amplitude, quadrature and rate digital control force must be calculated through PI control. The calculation results by referring to Figure 11 and Table 6 are Equations (47) and (48). When the DAC scale factor is k DA (V/bit) and the voltage-to-force scale factor is k f c (force/volt.), F x " k f c k DA N x , F y " k f c k DA N y . As above, if the control variables are calculated, the resonant frequency must be sought through PLL and the amplitude, quadrature and rate digital control force must be calculated through PI control. The calculation results by referring to Figure 11 and Table 6 are Equations (47) and (48). When the DAC scale factor is (V/bit) and the voltage-to-force scale factor is (force/volt.), = , = . Figure 11. Components of the control forces in the orbit trajectory.
The error for PI control is calculated as shown in Table 7. For integral, the trapezoidal rule is used. It should be noted that for the phase control, the frequency input must be reduced when the phase error is negative. For improving stability and compensating the truncation error, the limiter and summation block are applied to the controller.

Numerical Verification of the Algorithm through Simulations
Based on the design results of Sections 2 and 3, the Matlab/Simulink simulation SW is employed as shown in Figure 12. The target bandwidth by control loop is as shown in Table 8. The process of Figure 11. Components of the control forces in the orbit trajectory. E F a pN a q´sin pωt`φ corr q N a cosθ N a sinθ Q F q`Nq˘c os pωt`φ corr q´N q sinθ N q cosθ S F r pN r q´sin pωt`φ corr q´N r sinθ N r cosθ N x ptq "´pN a cosθ´N r sinθq sin pωt`φ corr q´N q sinθcos pωt`φ corr q (47) N x ptq "´pN a sinθ`N r cosθq sin pωt`φ corr q`N q cosθcos pωt`φ corr q (48) The error for PI control is calculated as shown in Table 7. For integral, the trapezoidal rule is used. It should be noted that for the phase control, the frequency input must be reduced when the phase error is negative. For improving stability and compensating the truncation error, the limiter and summation block are applied to the controller.

Numerical Verification of the Algorithm through Simulations
Based on the design results of Sections 2 and 3 the Matlab/Simulink simulation SW is employed as shown in Figure 12. The target bandwidth by control loop is as shown in Table 8. The process of control variables is approximated by a first-or second-order plus delay model (E & S: first-order, Q & L: second-order). The controller tunings are based on these models by using the recommended SIMC-PID method [13]. Table 9 shows the comparison between the control gains satisfying the bandwidth in the simulation program and the control gains in the actual DSP. control variables is approximated by a first-or second-order plus delay model ( & : first-order, & : second-order). The controller tunings are based on these models by using the recommended SIMC-PID method [13]. Table 9 shows the comparison between the control gains satisfying the bandwidth in the simulation program and the control gains in the actual DSP.  The simulation results are shown in Figures 13 and 14. The amplitude control variable is converged to the target amplitude = 512 bits (= 1 μm ) and the rate control variable is converged The simulation results are shown in Figures 13 and 14. The amplitude control variable is converged to the target amplitude E 0 " 512 bits`" 1 µm 2˘a nd the rate control variable is converged well to the target azimuth angle θ 0 " 45˝pS 0 " 512 bitsq. An estimate of the input rate is obtained by taking the difference of demodulated two forces, while the quadrature is nulled out [13]. well to the target azimuth angle = 45° ( = 512 bits). An estimate of the input rate is obtained by taking the difference of demodulated two forces, while the quadrature is nulled out [13].
(a) (b)  Figure 15 shows the test set to test the design and to conduct the gyro performance test linking with electronic board equipped with sensor, signal processing and control algorithm. First of all, the gyro parameters (Q, ΔQ, f, Δf) were measured after tuning the PI control gain suitable for sensor. When the control is stabilized, the amplitude control is turned off and the Q-factor is calculated by measuring the time constant by the target azimuth angle θ. At this time, ΔQ can be estimated from the difference between the maximum value and the minimum value. The test results are shown in Table 10 and Figure 16a. Then, if only the quadrature and rate control are turned off while maintaining the phase control and the amplitude control, θ is vibrated with certain cycle as shown in Figure 16b and this vibration frequency is the resonance frequency split Δf. As shown in Figure 16b, in case of test sample, the vibration cycle is approximately 55 s and the frequency split Δf is approximately 18 mHz.  Figure 15 shows the test set to test the design and to conduct the gyro performance test linking with electronic board equipped with sensor, signal processing and control algorithm. First of all, the gyro parameters (Q, ΔQ, f, Δf) were measured after tuning the PI control gain suitable for sensor. When the control is stabilized, the amplitude control is turned off and the Q-factor is calculated by measuring the time constant by the target azimuth angle θ. At this time, ΔQ can be estimated from the difference between the maximum value and the minimum value. The test results are shown in Table 10 and Figure 16a. Then, if only the quadrature and rate control are turned off while maintaining the phase control and the amplitude control, θ is vibrated with certain cycle as shown in Figure 16b and this vibration frequency is the resonance frequency split Δf. As shown in Figure 16b, in case of test sample, the vibration cycle is approximately 55 s and the frequency split Δf is approximately 18 mHz.  Figure 15 shows the test set to test the design and to conduct the gyro performance test linking with electronic board equipped with sensor, signal processing and control algorithm. First of all, the gyro parameters (Q, ∆Q, f, ∆f) were measured after tuning the PI control gain suitable for sensor. When the control is stabilized, the amplitude control is turned off and the Q-factor is calculated by measuring the time constant by the target azimuth angle θ. At this time, ∆Q can be estimated from the difference between the maximum value and the minimum value. The test results are shown in Table 10 and Figure 16a. Then, if only the quadrature and rate control are turned off while maintaining the phase control and the amplitude control, θ is vibrated with certain cycle as shown in Figure 16b and this vibration frequency is the resonance frequency split ∆f. As shown in Figure 16b, in case of test sample, the vibration cycle is approximately 55 s and the frequency split ∆f is approximately 18 mHz.   Next, the gyro scale factor and the bias were measured using the rate table and the bias instability and Angle Random Walk (ARW) measuring test were performed within a sound absorbing chamber. Bias instability and ARW are calculated through Allan variance using the data that recorded the gyro output for more than 8 h after stabilizing the gyro in a state isolated from the disturbance. As shown in Figure 17, in case of the sample used in this study, the bias instability is 0.07°/h and ARW is 0.006°/(rt·h).   Next, the gyro scale factor and the bias were measured using the rate table and the bias instability and Angle Random Walk (ARW) measuring test were performed within a sound absorbing chamber. Bias instability and ARW are calculated through Allan variance using the data that recorded the gyro output for more than 8 h after stabilizing the gyro in a state isolated from the disturbance. As shown in Figure 17, in case of the sample used in this study, the bias instability is 0.07°/h and ARW is 0.006°/(rt·h). Next, the gyro scale factor and the bias were measured using the rate table and the bias instability and Angle Random Walk (ARW) measuring test were performed within a sound absorbing chamber. Bias instability and ARW are calculated through Allan variance using the data that recorded the gyro output for more than 8 h after stabilizing the gyro in a state isolated from the disturbance. As shown in Figure 17, in case of the sample used in this study, the bias instability is 0.07˝/h and ARW is 0.006˝/(rt¨h).

Conclusions
HRG is one of the kinds of CVG, which measures the angle or angular velocity using the Coriolis force produced by the rotational motion. The rotational angle or angular velocity can be measured using the principle that the standing wave generated in the quartz resonator in the form of hemisphere shell performs the procession motion.
HRG technology development is underway in many advanced countries and the next generation HRGs, which can materialize the objectives of subminiature size, high precision and high reliability with the 2-piece HRG system is being applied with multi-flexing method and differential control development in the United States, France, etc.
Therefore, in this article, a controller design suitable for a 2-piece HRG system was performed. To design the controller, the electromechanical modeling of the 2-piece HRG system was pre-performed. To interpret the vibration characteristics due to the switched discontinuous excitation force, the Duhamel integral method was applied. In addition the electromechanical gains were calculated considering the mode shape of the thin hemispherical shell. It was proven that the design results are valid by showing an error within 7% in the comparison between design results and measurement results of the amplitude and charge amplifier output.
Based on such modeling, the signal processing based on the multi-flexing method and differential control algorithm were designed. The sensing and driving cycles of x-and y-axis were divided by time with five operation cycles using a common element. The sensing signals generate the control inputs signal through the sampling and demodulation processes. The controller output generates the final control voltages through the modulation processes and the phase delay compensation algorithm was applied in the modulation process. In FTR mode, Control is composed of phase, amplitude, quadrature and rate control and the pendulum variables are used as control variables. The designed algorithm was verified through Matlab/Simulink simulation and in the results, the bandwidth of the amplitude and quadrature control were satisfactory, at 1-5 Hz, and the bandwidth of the rate and phase control were satisfactory at 7.5-12.5 Hz. Finally, the electronic circuit was made by equipping the algorithm in FPGA and DSP and that it satisfied with the target bandwidth was verified through the experiment. In addition, through the sensor linking test, the error identification was performed and it was conformed that the bias instability and ARW are approximately 0.07°/h and 0.006°/(rt·h), respectively by performing the gyro performance test.
Author Contributions: Jungshin Lee performed the electromechanical modelling, signal processing and digital controller simulation. Sung Wook Yun designed the PLL control algorithm and performed the experiments. Jaewook Rhim designed the digital algorithm and experiment. Jungshin Lee wrote the paper.

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

Conclusions
HRG is one of the kinds of CVG, which measures the angle or angular velocity using the Coriolis force produced by the rotational motion. The rotational angle or angular velocity can be measured using the principle that the standing wave generated in the quartz resonator in the form of hemisphere shell performs the procession motion.
HRG technology development is underway in many advanced countries and the next generation HRGs, which can materialize the objectives of subminiature size, high precision and high reliability with the 2-piece HRG system is being applied with multi-flexing method and differential control development in the United States, France, etc.
Therefore, in this article, a controller design suitable for a 2-piece HRG system was performed. To design the controller, the electromechanical modeling of the 2-piece HRG system was pre-performed. To interpret the vibration characteristics due to the switched discontinuous excitation force, the Duhamel integral method was applied. In addition the electromechanical gains were calculated considering the mode shape of the thin hemispherical shell. It was proven that the design results are valid by showing an error within 7% in the comparison between design results and measurement results of the amplitude and charge amplifier output.
Based on such modeling, the signal processing based on the multi-flexing method and differential control algorithm were designed. The sensing and driving cycles of x-and y-axis were divided by time with five operation cycles using a common element. The sensing signals generate the control inputs signal through the sampling and demodulation processes. The controller output generates the final control voltages through the modulation processes and the phase delay compensation algorithm was applied in the modulation process. In FTR mode, Control is composed of phase, amplitude, quadrature and rate control and the pendulum variables are used as control variables. The designed algorithm was verified through Matlab/Simulink simulation and in the results, the bandwidth of the amplitude and quadrature control were satisfactory, at 1-5 Hz, and the bandwidth of the rate and phase control were satisfactory at 7.5-12.5 Hz. Finally, the electronic circuit was made by equipping the algorithm in FPGA and DSP and that it satisfied with the target bandwidth was verified through the experiment. In addition, through the sensor linking test, the error identification was performed and it was conformed that the bias instability and ARW are approximately 0.07˝/h and 0.006˝/(rt¨h), respectively by performing the gyro performance test.
Author Contributions: Jungshin Lee performed the electromechanical modelling, signal processing and digital controller simulation. Sung Wook Yun designed the PLL control algorithm and performed the experiments. Jaewook Rhim designed the digital algorithm and experiment. Jungshin Lee wrote the paper.

Appendix B
Assuming that x ct ptq " t C 1 T 0 pn´1q f pτq h pt´τq dτ, x pt ptq " f pτq h pt´τq dτ, the calculation results are follows. As in this moment, the Q-factor is high, we can assume that ω dω: x ct ptq " x ct ptq " x pt ptq "

Appendix C
The time response of SDOF spring-damper system by the switched excitation can be arranged as follows: With the same method, the mean amplitude by operation cycle is calculated as follows: A pt ď 3C 1 T 0 q " C x ω a