Turntable IMU Calibration Algorithm Based on the Fourier Transform Technique

The paper suggests a new approach to calibration of a micromechanical inertial measurement unit. The data are collected on a simple rotating turntable with horizontal (or close to) rotation axis. For such a turntable, an electric screwdriver with fairly low rotation rate can be used. The algorithm is based on the Fourier transform applied to the rotation experimental data, implemented as FFT. The frequencies and amplitudes of the spectral peaks are calculated and collected in a small set of data, and calibration is done explicitly with these data. Calibration of an accelerometer triad and choosing the IMU coordinate frame are reduced to approximating the collected data with an ellipsoid in three dimensions. With rotation frequency calculated as the peak frequency of accelerometer readings, calibration of the gyros is a straightforward linear least square problem. The algorithm is purely algebraic, requires no iterations and no initial guess on the parameters, and thus encounters no convergence problems. The algorithm was tested both with simulated and experimental data, with some promising results.


Introduction
Calibration of an inertial measurement unit (IMU) is a necessary step in preparing the IMU for operation, that being either inertial navigation or some supplementary task in biomechanics or robotics. The purpose of calibration is to provide formulas and estimate parameters needed to convert the raw accelerometer and gyro readings into the specific force and the angular rate of the unit. The majority of calibration models found in the literature are linear and include scale factors and biases; similar models are used in this paper. The calibration procedure consists of experiments with the IMU and algorithms to process the collected data.
Calibration of an IMU is fairly straightforward if done on a precise turntable with two or three degrees of freedom, measuring rotation angles with high accuracy. The table is usually well-calibrated, and its azimuth angle relative to the North is known, thus the true angular velocity and specific force are also known. In such a case, the calibration procedure usually consists of several steps-different static positions and rotations with prescribed angular velocities. Therefore, the data processing algorithms consist of several respective steps to calibrate biases and scale factors of accelerometers and angular velocity sensors. Mathematical formulas for these steps are often straightforward and can be reduced to linear least squares [1,2], or Kalman filtering [3,4].
A series of papers considers IMU calibration when a so-called low-grade turntable [5] is used. Such a turntable delivers rotation, but either it does not provide information about rotation rate or angles, or this information is inaccurate. Thus, the only information available is the data taken from the IMU sensors. In such settings, IMU calibration cannot be done step-by-step, and dynamic equations of motion must be taken into account. In Refs. [3,5], an algorithm for calibration on a low-grade turntable based on the extended are calibrated, the angular rate sensor calibration is reduced to a straightforward linear algebraic problem.
We have not found many papers about IMU calibration by applying the Fourier technique to the turntable experiment. One such paper is Ref. [18]. The calibration model in Ref. [18] is the same as ours. However, the big difference is that in Ref. [18], the true angular velocity and the true specific force are assumed known from the turntable data. With this assumption, the equations are completely linear and can be solved without iterations both in the time domain and in the spectral domain. In our case where the true sensors readings are unknown, the initial equations are nonlinear, thus, as far as we know, non-iterative algorithms in the time domain do not exist, which leaves the spectral approach we propose the only non-iterative one. Some preliminary results were published by us in Ref. [19].
The paper is organized as follows. First, we state the problem. Next, the necessary mathematical formulas are derived and analyzed. Further on, we present results of simulation and experiments with a micromechanical IMU that support our reasoning. The bold letters such as f denote column vectors in the three-dimensional space, the bold capital letters such as F denote column vectors in the three-dimensional space of Fourier transforms. When the coordinate frame the column vectors are written in is of importance, we write them as f .z , where z 1 z 2 z 3 are the coordinate frame axes.

Calibration Model
The IMU measures the angular rate ω of the sensor relative to the inertial frame and the specific force f acting on the proof mass M of the accelerometer. In some cases, especially for IMU operating in the environment with high angular rates, it is necessary to assume several proof masses. However, here, we assume just one proof mass. Both the angular rate and the specific force are given by the column vectors of their coordinates ω z and f z in the so-called instrument frame Mz 1 z 2 z 3 . The origin of this frame is at the proof mass, its axes z 1 , z 2 , z 3 are rigidly connected with the instrument body.
Under the linear calibration model, the raw accelerometer data a(t) and the raw gyro data w(t) are connected with the true specific force f .z and angular velocity ω .z written in the instrument frame Mz 1 z 2 z 3 by the calibration formulas [3,7]: where δ f (t), δω(t) are the unmodeled measurement errors. The goal of calibration is to estimate the 3 × 3 matrices S f , S ω and the 3 × 1 column vectors b f , b ω . Many calibration procedures can be found in the literature; some are mentioned in the introduction section.
In this paper, we consider the following calibration procedure settings. The IMU is mounted on the shaft of a turntable or a screwdriver in several positions and, in each position, it is rotated for several tens of revolutions using the screwdriver motor. These rotations are called "experiments". The screwdriver body does not move relative to the Earth during all the experiments.
We use the following coordinate frames (Figure 1). The navigation frame On 1 n 2 n 3 is connected with the Earth, its axis n 3 points up, the axis n 1 is the projection of the shaft rotation axis on the horizon. The axes origin O lies on the shaft rotation axis. The azimuth angle between n 1 and the East direction is denoted as ψ. The frame Os 1 s 2 s 3 is rigidly tied to the body of the screwdriver, its s 1 axis is directed along the shaft, the angle 0 < α < π/2 of this axis with the horizon is constant. The axis s 3 lies in the vertical plane, the axis s 2 lies in the horizontal plane. The frame Oe 1 e 2 e 3 is rigidly tied to the shaft: the axis Oe 1 points along the shaft, the axes Oe 2 , Oe 3 rotate in the plane perpendicular to the shaft during experiments. The IMU instrument frame Mz 1 z 2 z 3 is rigidly tied to the shaft in each experiment but changes its orientation relative to the shaft from experiment to experiment. The proof mass M of the accelerometer can lie out of the shaft rotation axis: the centrifugal forces induced by the shaft rotation are accounted for by the algorithm. Figure 1. Coordinate frames. The axes e 1 , e 2 , e 3 form the shaft frame; the axes z 1 , z 2 , z 3 form the IMU instrument frame. The axis n 1 is horizontal and lies in the vertical plane of the screwdriver shaft and diverges by the angle ψ from the East direction, the axis n 3 points up. The angle α is the angle of the shaft to the horizon.

Experiments and Data Models
Let P be the number of experiments (different rotations) numbered by p = 1, . . . , P: in each experiment, the shaft axes e 2 , e 3 rotate around the axis e 1 , the angle of the axis e 2 relative to the stand axis s 2 is denoted by θ p (t): Here, T p is the duration of the experiment, ω p is the angular velocity of the shaft, γ p is the phase shift, δθ p (t) is the disturbance due to imperfection of the turntable. For brevity, we set δθ p (t) to zero in the theoretical section of this paper (influence of this term can be compensated for by analyzing the higher frequency spectrum of the sensors). Let us complement experiments-rotations with static experiments were the IMU stands still at different positions, denoting these experiments with numbers q = 1, . . . , Q: Let g and u be the gravity and the Earth angular velocity vectors, g and u be their absolute values. In the navigation frame, these vectors are written as g = −gn 3 , u = n 1 sin ψ cos φ + n 2 cos ψ cos φ + n 3 sin φ, where φ is the latitude, ψ is the azimuth angle of the screwdriver shaft. In the screwdriver frame Os 1 s 2 s 3 , these vectors can be written as Here, we use the notation g 1 = g sin α, g 2 = g cos α, u 1 = cos α sin ψ cos φ + sin α sin φ, u 2 = cos ψ cos φ, u 3 = cos α sin φ − sin α sin ψ cos θ.

of 16
Note that ψ, α are unknown before calibration, thus the constants in (7) are also unknown. In the shaft frame, the gravity g and the Earth angular velocity u during the p-th experiment can be written as −g p (t) = g 1 e 1 + g 2 [e 2p sin θ p (t) + e 3p cos θ p (t)], Now, we remember that the sensor angular velocity and the specific force acting on the proof mass in the p-th experiment can be written as where ρ p = const is the vector of displacement of the proof mass relative to the shaft. In the instrument frame in the p-th experiment, the above vectors can be written as Here, ρ 1p , ρ 2p , ρ 3p are the components of ρ p along the e 1p , e 2p , e 3p axes. Writing down the centrifugal force, we have neglected its component due to the Earth rotation. Turning e 3p , e 2p , e 3q , e 2q appropriately in the plane perpendicular to the shaft, we can assume without loss of generality that γ p = 0, γ q = 0 in (3), (4). Thus, we can write For the static experiments (4), where the screwdriver motor is switched off, the IMU measurements can be written as Equations (10)- (13) were written in the invariant vector form. Further on, we assume that they are written as column vectors in the instrument frame: f = f .z , etc. However, since the instrument frame is the only frame used below, we omit the superscript .z. Note that the column vectors e 1p , e 2p , e 3p and e 1p , e 2p , e 3p written in the instrument frame are constant during each experiment but change from experiment to experiment.

Calibration Algorithm
This section covers mathematical foundations of the proposed algorithm; some technical details of implementation are discussed in the results section.

Fourier Transform for Accelerometer Calibration Formulas
In this section, we neglect imperfections of the motor, setting δθ p (t) = 0, thus setting θ p (t) = ω p t. Note that these imperfections do not influence the gyro calibration due to their high-frequency spectrum, but can cause a bias in accelerometer calibration due to the centrifugal force bias. However, these biases can be calculated by analyzing the oscillations in the angular velocity measurements. We do not pursue this matter further here. Let us apply the Fourier transform [20] in (1), (10) and denote the Fourier images as f p (t) → F p (ω), etc.: The Fourier transform in (1) for both rotations and static experiments takes the form where δ(ω) is the delta-function. The Fourier transform in (10) takes the form We have used the following well-known Fourier transform formulas [20]: Comparing (17), (18) and (15), (16), we see that the Fourier transform A p (ω) has delta-function-type peaks at ω = 0, ω = ±ω p , while the Fourier transform A q (ω) has delta-function-type peaks at ω = 0. Let us denote the amplitudes of these peaks as The Fourier transforms A p (ω), A p (ω) take the form Note that ω p here can be calculated as the position of a peak of A p (ω) in the interval ω > 0, if the rotation was clockwise, or in the interval ω < 0, if the rotation was counterclockwise. Equations (15) and (16) can now be rewritten as (22) Collecting coefficients at delta-functions in the equation above and in (17), (18), we obtain g 1 e 1q + g 2 e 3q = S fĀq + b f , (23) Equation (23) are the base for calibrating the accelerometers. They can be resolved in the unknowns S f , b f , g 1 , g 2 , e 1p , e 2p , e 1p , e 2p in several ways. We take the simplest way here, not claiming it to be the most accurate one.
Let us exclude the unknown constants g 1 , g 2 , e 1p , e 2p , e 3p , e 1q , e 2q from (23) by algebraic manipulations: (24) Equation (24) can be resolved to find S f , b f . Note that the choice of the sensor frame is arbitrary provided it is fixed in the instrument body; thus, S f includes only six independent variables [3]. For a given number P, Q of experiments, we have Q + 2P equations. It looks like that to calculate 6 + 3 = 9 calibration parameters, it suffices to set Q = 3, P = 3. However, a more careful look shows that we must take Q ≥ 4, see below.

Fitting Accelerometer Data to an Ellipsoid
Equation (24) is quadratic. To reduce these to linear equations, we use the well-known trick [17]: introduce the symmetric positive matrix M, the vector m, and the constant m 0 as Equation (24) can then be rewritten as Equation (26) defines an ellipsoid in the three-dimensional space. To obtain its parameters, let us divide the first equation in (26) with g 2 − m 2 0 , and use the notation We get Note that if Q = 3, and the column vectorsĀ 1 ,Ā 2 ,Ā 3 are linearly independent, then the system (28) allows the false trivial solution which will follow with the estimate of the scaling matrix as S f = 0. Thus, for reliable calibration, we must set Q ≥ 4.
Expanding coefficients forM,m to a vector as and expanding the unknowns to the vector we obtain the system of linear equations in x (32) The overdetermined system (32) can be solved in the mean square sense to obtain x, thenM,m from (31), and further on Finally, the calibration parameters for the accelerometers S f , b f are obtained by factorizing M in upper triangular, lower triangular or symmetric form, depending on the situation: Knowing the calibration parameters, we can find g 1 , g 2 , α from (23) as: Next, we calculate the orts e 1p , e 2p , e 3p as column vectors in the instrument frame as: Now, we have to discuss one question. If the directions of rotation are known beforehand, then the signs of ω p are also known, and the formulas (36) are quite correct. Now suppose that we do not know these signs, and have taken the wrong sign of some ω p . Then, the sign of I p will be wrong; consequently, the signs of e 1p , e 2p will be also wrong. To find the correct sign, let us multiply the second equation in (23) with e T 1p taken from (36): we get g 1 = e T 1 (S fĀp + b f ). We know that the value g 1 = g sin α is positive; thus, the scalar product on the right must be also positive. If it is negative, we must change the sign of ω p and change the signs of e 1p , e 2p before proceeding to calibrate the gyros.
Note that the vectors e 1q , e 2q , e 3q cannot be explicitly calculated from the above equations. However, if the experiments were done in such a way that each q-th static experiment precedes to or follows after some p-th rotation experiment without dismounting the IMU from the screwdriver shaft, then we can set e 1q = e 1p . However, we cannot do the same with e 2q , e 3q , because the rotation angle of these unit vectors around the shaft can be different from the rotation angles of e 2p , e 3p .

Calibration of Angular Rate Sensors
When calibrating accelerometers, we have obtained the axes e 1p , e 2p , e 3p of the shaft frame in the sensor frame. We have also obtained the rotation angular rates ω p as frequencies of the peaks of the Fourier transform of the accelerometer data.
(37) Here, Ω q (ω) = (u 1 e 1q + u 2 e 2q + u 3 e 3q )δ(ω), As earlier mentioned, introducing the peak values of the angular velocity spectrum we can write When we collect factors before the delta-functions in (37), (38), we obtain the system of 9P + 3Q equations We have 21 unknowns here: besides the twelve calibration parameters S ω , b ω , we do not know nine more: u 1 , u 2 , u 3 , e q2 , e q3 . Further on, we can select one of two ways.
If we are calibrating a low-accuracy IMU, we can neglect the Earth angular velocity, and get the system of 3P + 3Q linear equations in the twelve unknowns: If P ≥ 3, Q ≥ 1, these equations allow us to calculate the elements of S w , b ω . If we are calibrating a high-accuracy IMU, then we must take into account the Earth angular velocity. To get rid of the unknown vectors e 2q , e 3q , we multiply (41) with e T 1q to get the system of equations This is a system of Q + 9P linear equations, with additional unknowns u 1 , u 2 , u 3 , 18 unknowns in total. If P = 3, Q = 4, which is the minimal number of experiments for calibrating the accelerometers, we get 31 equations-more than enough. Note that we need to know only the absolute value of the Earth angular velocity here-no need to know the azimuth angle or even the latitude.

Results
The algorithm was tested both with simulated data and with the commercial micro-IMU named x-IMU (Figure 1) from the x-io Technologies Limited company [21]. The goal of simulation was to estimate the potential accuracy of calibration. Since x-IMU is a lowaccuracy IMU, the Earth angular velocity was not taken into account in data processing.

Calculating the Fourier Transform with FFT
In the above mathematical formulas, we have assumed infinite duration of experiments; thus, the Fourier transforms of the sine wave and of the constant signal were delta-functions. In the real world, experimental data are recorded at some finite time intervals of duration T p , p = 1, . . . , P for rotations, or T q , q = 1, . . . , Q for static experiments. Thus, the Fourier transform is calculated for a finite support function weighted with some window function h(t) It is well known that this provides spectral leakage [20]: the Fourier transform for a windowed exponential e iω p t is H(ω − ω p ), where H(ω) is the Fourier transform for h(t). The Fourier transform for windowed f p (t) is The choice of the window h(t) depends on the purpose. For maximal resolution of the spectrum, the rectangular window is the optimum. In this case, H(ω) = sinc(ωT/2) which is a fast oscillating function with minimal width of the central peak. If it is important to avoid false peaks of the spectrum, windows with nearly monotonically decreasing H(iω) are preferred. When the rotation rate is not very stable (as was the case in our experiments, see below), it is better to widen the spectral window so that it would cover the varied main frequencies.
Another consideration arises when we remember that the spectrum is calculated by FFT at discrete set of frequencies with the frequency step ∆ω = 2π/T p [20]. In addition, we need to estimate not just the peak frequency and amplitude, but also the real and imaginary parts of the spectrum which can change fast near the peak frequency. To calculate these values with discrete frequencies, it is preferable to use a window function with a fairly wide central lobe.
As a good compromise, accounting for most of the factors mentioned above, we use here the Hann window function To increase the frequency sampling rate of FFT, we use zero padding [20], adding zeros to the signal on the left and right of the time interval. When N zero signal intervals each of duration T p are padded to the signal both on the left and on the right of the recorded data, the frequency sampling rate increases 2N + 1 times, to ∆ω = 2π/T p /(2N + 1). To preserve the spectral amplitudes with zero padding, appropriate scaling of the spectrum is done. In the experiments described in the next subsections, we set N = 10.

Calibration with Simulated Data
When generating the simulated data, their sampling rate was set to 100 Hz, the shaft rotation main frequency ω p was chosen as 2.1 rotations per second, the number of rotation experiments was P = 3, while the number of static experiments was set to Q = 4. In each of the three rotation experiments, one of the axes of the instrument body pointed along the rotation shaft. The generated noise was a Gaussian white noise; RMS of the accelerometer noise was set to σ f = 0.01m/s 2 ; RMS of the angular rate sensor noise was set to σ ω = 0.1 rad/s. Duration T p = T q of the experiments was varied. The results are shown in Table 1. The errors in S f , S ω are relative, the errors in b f , b ω are in m/s/s and rad/s, respectively. When rotations are done for 10 s at each position (about 20 full revolutions of the body), the calibration accuracy is mediocre. When rotations are done for 80 s (about 160 full revolutions of the body), calibration results become very good. Note that the accuracy of angular rate sensor calibration is considerably higher than that of accelerometer calibration. Probably this is due to using the ellipsoid fitting method to calibrate the accelerometers, which is known not to be the most accurate. Table 1. Calibration errors depending on duration T p = T q of each experiment. The errors in S f , S ω are relative, the errors in b f , b ω are absolute, in m/s/s and rad/s, respectively.

Calibration with Experimental Data
We show here the results of just one test with x-IMU [21], with the data recording rate of 256 Hz. Data recording was done via bluetooth connection. We used a consumer screwdriver as a turntable (Figure 2). This screwdriver is two-speed, the revolution frequency in the slow mode is about 4 rotations per second. The IMU was fastened to the shaft in three different positions: p = 1, 2, 3 = P. The static experiments were just the positions of the IMU on the shaft before and after rotations: q = 1, 2, 3, 4 = Q. No data were collected when the IMU was not fastened to the shaft. All rotations were done clockwise. Figure 2. The IMU mounted on the shaft of a screwdriver. The IMU is inside a small wooden box which is fixed to the shaft in several positions by the adhesive tape. The trigger of the screwdriver is held by a rope in the partly pushed position to limit the angular rate of rotations (otherwise, the gyro readings are saturated). Unfortunately, for the partly pushed trigger, the low angular rates of the screwdriver we use are fairly unstable. Books on inertial navigation are of great help in the experiments.
We used a factory pre-calibrated IMU. Thus, the raw data we obtained were actually close to the true data; the accelerometer readings were in fractions of g, while the angular rate sensor data were in degrees per second. Duration of experiments, including the time taken by mounting and dismounting the IMU on the shaft of the screwdriver, took approximately 15 min, as can be seen from Figure 3. Figure 3 shows the graph of angular sensor raw reading records during the whole sequence of the experiments. The trigger of the screwdriver was pressed only partially to limit the angular rate of rotations to be within the working range of the IMU angular velocity sensor which is about three rotations per second (see the rope holding the trigger). With the half-pushed trigger, the angular velocity can be seen to be unstable; this fact is seriously decreasing the accuracy of our spectral method, as the main spectral peak at revolution frequency tends to spread wider. Figure 4 shows a fragment of raw accelerometer readings. We see high-frequency oscillations in the data which are probably due to the imperfections of the gearbox. Figure 5 shows a fragment of the FFT of the accelerometer data near the peak at the revolution angular rate of approximately 700 deg/s. We see that the spectrum is spread near the central frequency, partly due to the instability of rotations, partly due to the spectral leakage discussed above [20]. To make the leaked spectrum more smooth, time weighting of the data with the Hann window was done. To increase the sampling frequency in the frequency domain to obtain better estimates of positions of the peaks, the data in time were zero-padded, increasing the length of the data vector 21 times. Figure 6 shows a fragment of FFT of the angular rate sensor data near the zero frequency.     Figure 5. The raw accelerometer data spectrum F p (ω) for p = 1 at the zero frequency. The blue line is the absolute value of A p (ω), the red and green lines denote the real and imaginary parts of A p (ω). The marker shows the estimated amplitude peaks. The frequency along the x-axis is in rad/s.  Figure 6. The raw gyro data spectrum W p (ω) for p = 1 near the zero frequency. The spectrum is normalized in such a way that the peak amplitudeW p is the rotation frequency in deg/s. The blue line is the absolute value of W p (ω), the red and green lines denote the real and imaginary part of W p (ω). The marker shows the estimated amplitude peaks. The frequency along the x-axis is in rad/s. The factory set scaling matrices S f , S ω and these produced by our calibration algorithm for the x-IMU instrument are shown below: The differences in alignment (non-diagonal elements in the matrices) are not more than 0.06 degrees; the differences in scaling factors (diagonal elements of the matrices) are not more than 0.2%. We consider this as a sound verification of our algorithm, since factory calibration of x-IMU has proved to be very accurate in our experiments in pedestrian navigation [10].
We do not discuss here the estimates of the sensor biases since for the micromechanical IMU, the biases are usually fairly unstable and change from switch-on to switch-on.

Discussion
The paper presented an approach for calibration of a micro-mechanical IMU using a simple turntable, such as a domestic screwdriver, with a calibration algorithm based on the Fourier transform of the data. Verification of our calibration algorithm was done in two ways. First, we used the simulation data and showed that when the duration of experimental rotations was taken more than 40 s (80 revolutions at two rotations per second), the calibration accuracy became sufficient for any purpose. When duration of rotations was below 10 s (20 revolutions at two rotations per second), calibration results were inaccurate. The latter can be explained by the effect of "spectrum leakage" when using FFT on a short time interval.
Second, we performed experiments with the micromechanical inertial measurement unit called x-IMU [21] and compared the results of calibration using our algorithm to the x-IMU factory calibration data, which had proved to be very accurate in our experiments in pedestrian navigation [10]. The differences in the scale factors were around 0.2%, the errors in estimating directions of sensitivity axes were around 0.06 degrees. However, the calibration accuracy was not as high as was expected from simulations. We see the reason for this in that, in our experiments, the rotation angular velocity drifted considerably due to instability of screwdriver RPM. We expect the results to be much better when we use a better turntable than our screwdriver. Another way to increase accuracy is to use reverse rotations, which we did not do.
Some rough estimates not included here show that the accuracy of our algorithm is considerably lower than that of the well-known Kalman filter approach to calibration on a rotating turntable [3], which can be seen to be asymptotically optimal under the assumption of the sensor errors being Gaussian white noise. Nevertheless, a big advantage of our algorithm, which comes from purely algebraic formulas used, is its guaranteed convergence without any need for an initial guess. Without such a guess, the Kalman filter algorithm often diverges. The results of calibration with our algorithm can be used as an initial guess for the Kalman filter algorithm. Joining together the Fourier approach and the Kalman filter approach, with an initial guess provided by the former, is an aim of our future work.
Our algorithm has two versions-for high-grade and low-grade IMU. For the case of low-grade IMU, which is discussed in this paper in more detail, we do not account for the Earth angular velocity, and the algorithm becomes much simpler. Note that when rotation rates become about one turn per second, the Earth angular velocity influence on estimating the gyro scaling matrices is negligible since the rotation rates are several orders of magnitude higher than the Earth rotation rate of 13 degrees per hour. We think that the high-grade IMU version of our algorithm should be used only when a turntable with fairly constant rotation rate is used.
One of the popular micromechanical IMU calibration methods is the so-called calibration without external equipment [13,14], etc., which we have used a lot in our experiments in pedestrian navigation [9]. We have found (the results are not included here) that the new algorithm performs considerably better when it comes to the angular rate calibration of the sensors because their scale factors are estimated better with relatively fast and lengthy rotations, which are impossible to perform by hand. Moreover, contrary to our algorithm, the rotation-by-hand angular rate calibration requires an initial guess on the parameters. Accuracy of accelerometer calibration for the two approaches is comparable. We note also a general advantage of the spectral approach which is that the spectrum provides additional information about the IMU, helping to find out certain effects such as proof mass separation the initial mathematical calibration model has not accounted for.
As already mentioned in the Introduction, we have found just one paper [18] about IMU calibration by applying the Fourier technique which is in relation to ours. The big difference is that in Ref. [18], the true angular velocity and the true specific force are assumed known from the turntable data; thus, the equations are completely linear and can be solved without iterations both in the time domain and in the spectral domain.
In our case where the true sensor readings are unknown, the initial equations are nonlinear, which makes the spectral approach we propose look special from the point of view of guaranteed convergence.