{'fraction_non_alphanumeric': 0.026908946479992056, 'fraction_numerical': 0.011170688114387846, 'mean_word_length': 5.194034440344404, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Attitude Reconstruction from Inertial Measurement: Mitigating Runge Effect for Dynamic Applications

Time-equispaced inertial measurements are practically used as inputs for motion determination. Polynomial interpolation is a common technique of recovering the gyroscope signal but is subject to a fundamentally numerical stability problem due to the Runge effect on equispaced samples. This paper reviews the theoretical results of Runge phenomenon in related areas and proposes a straightforward borrowing-and-cutting (BAC) strategy to depress it. It employs the neighboring samples for higher-order polynomial interpolation but only uses the middle polynomial segment in the actual time interval. The BAC strategy has been incorporated into attitude computation by functional iteration, leading to accuracy benefit of several orders of magnitude under the classical coning motion. It would potentially bring significant benefits to the inertial navigation computation under sustained dynamic motions.


I. INTRODUCTION
Inertial sensors, i.e., gyroscopes and accelerometers, produce discrete measurements mostly at equally-spaced time instants, which are widely used to design algorithms for motion determination, e.g., in various applications of attitude estimation or inertial navigation [1,2]. Usually, the inertial sensor signals are reconstructed by means of popular polynomial interpolation methods [3][4][5][6][7][8][9][10][11], as polynomials are supposed to approximate a smooth function as closely as possible if the number of samples are sufficiently large. The smooth interpolant of the underlying smooth functions at equispaced samples is a ubiquitous problem [12]. However, the close approximation of polynomials to a smooth function is preconditioned on the distribution of samples in addition to the function property [13]. It is well known that polynomial interpolants in Chebyshev samples, which are clustered towards the interval boundary, are well-behaved and have been widely used. For interpolants of a function at equispaced samples to converge along with the growing number of samples, the function must be analytic, not just on the considered interval but throughout the so-called Runge region [13]. Furthermore, computer rounding errors might be considerably amplified or even cause divergence in practice, depending on the distribution of the samples. Unfortunately, the equispaced samples are among those that are exponentially ill-conditioned in term of numerical stability [13]. The Runge phenomenon, in which wild oscillations are observed near the interval boundary, is not unusual for an equispaced high-order interpolant.
The work was supported in part by National Key R&D Program of China (2018YFB1305103). An abstract was submitted to DGON ISS 2021, Germany.
A number of methods have been proposed to approximate the underlying functions from equispaced samples, but unfortunately it was proven that even for analytic functions no stable method can be simultaneously exponentially convergent and well-conditioned [14]. The linear rational interpolant proposed in [15] has been well known as a significant development in the field of numerical analysis [12], and was further extended in [16] making use of additional samples beyond the interval boundary for much improved error mitigation in the original interval. Notably, it possesses error characteristics comparable to the excellent polynomial interpolant in Chebyshev samples and is very close to be optimal for equispaced samples [12]. Recent advance in the inertial navigation computation has been achieved by the functional iterative integration approach to solve the kinematics differential equations up to the machine precision [8][9][10][11], which would potentially bring significant benefits to the inertial navigation computation under sustained dynamic motions [17,18]. As it takes as inputs the angular velocity/specific force polynomials fitted from timeequispaced gyroscope/accelerometer measurements, in principle any approximation error in the fitted polynomials would inevitably decay the computation accuracy. For instance, the Runge phenomenon is recently found having an obvious adverse effect on the accuracy and robustness of attitude computation for dynamic applications [19]. For a fixed sampling rate restricted by any practical inertial navigation system, more gyroscope samples might be desired for more accurate polynomial approximation of angular velocity, but it is fundamentally prohibited by the Runge phenomenon, as reported in [19]. The technical contribution of the current paper is proposing a straightforward strategy of suppressing the Runge effect in equispaced interpolation and incorporating it into attitude computation for enhanced accuracy in highly dynamic applications. It helps break the barrier of high-accuracy senor measurement interpolation on equispaced samples for rigid motion reconstruction. The content is organized as follows. Section II briefly reviews the polynomial and rational interpolation advances in the applied mathematics community, and specifically underlines their numerical properties in terms of Lebesgue constant by a simple example. Section III investigates polynomial and rational interpolations in angular velocity reconstruction and attitude computation. Conclusions are drawn in the last section.

II. POLYNOMIAL/RATIONAL INTERPOLATION AND RUNGE EFFECT
For the ordered distinct nodes and   the corresponding signal values 0 1 , , , N f f f  , this section reviews the technique of polynomial and rational interpolations for obtaining functional reconstruction of the original signal. Specifically, the underlying reason of the Runge effect is highlighted.

A. Polynomial Interpolation
It can be readily checked that . For t other than the above distinct nodes, the polynomial (1) can be reformulated as from the uniqueness of the polynomial interpolation. Dividing (2) by (3) gives the well-known barycentric form of the polynomial interpolant [12]   with the special case (4) is provably numerically stable [12,13]. Seemingly like a rational function, (4) is indeed a polynomial with the specially defined weights given as in (2). For other weights, however, it will in general be a rational function instead of a polynomial, which is to be discussed in the sequel. The above barycentric form has the advantage of cancelling common factors and resulting in simple weights [13]. Regarding the equispaced nodes, the weights are   The error of an interpolant is governed by not only the theoretical interpolation error in exact arithmetic but the amplification of rounding/measurement error [21]. The latter is characterized by the Lebesgue constant in which the first term is related to the theoretical interpolation error in exact arithmetic and the second term is related to the interpolation error in presence of rounding/measurement error. It has been proved that the Lebesgue constant for equispaced nodes and Chebyshev nodes [13]. No interpolation method for equispaced nodes exists that converges geometrically and has a Lebesgue constant decreasing geometrically [14,21]. As a matter of fact, the polynomial interpolation in equispaced points is exponentially ill-conditioned, which leads to the wellknown Runge phenomenon near the interval boundary. See a demo example (polynomial interpolation in red lines in Fig. 2) illustrated in Sec. II.C below.

B. Rational Interpolation
The rational interpolation methods can be classified in two classes: nonlinear ones and linear ones [12]. In the latter case, the interpolant is linear in f , resembling the form in (4) of which the denominator depends on the nodes but not on f .
Notably, the famous family of Floater and Hormann (FH) rational interpolant [15], for a design parameter d , is given by As long as f is smooth enough, the FH interpolant has the convergence rate independent of how the nodes are distributed, where h is the maximum spacing between consecutive nodes. Larger value of d is supposed to lead to faster convergence, but rounding or measurement errors also come into effect. For equispaced nodes, it was shown in [22] that the Lebesgue constant for the FH interpolant (8) [16]. It was conjectured that the extended FH interpolant is likely very close to optimality for equispaced nodes [12]. The strategy of using the neighbouring samples to assist signal reconstruction is not uncommon in the navigation community. For example, actual gyroscope samples in the previous computation interval, in contrast to the Taylor-expansion extrapolated points in EFH [16], have been frequently used in designing attitude algorithms [23][24][25]. For the attitude computation of main interest in this paper, the required additional signal values on both ends of the current interval would be actual gyroscope samples at the cost of moderate computational latency. To use the neighbouring samples is specifically named as the borrowing-and-cutting (BAC) strategy for interpolation. This name is used in this paper to intuitively underline the main idea of the approach.   Table I lists the Lebesgue constant bounds for polynomial/rational interpolation on equispaced and Chebyshev nodes. The Lebesgue constant bounds at equispaced samples are plotted in Fig. 1. The bound trend tells that the Lebesgue constant of the polynomial interpolation is by far larger than that of the EFH rational interpolation when the number of samples surpasses ten. Note that the lower bound for the polynomial seems quite overoptimistic, as will be seen from Fig. 2 below. Generally speaking, it is safe for us to use the polynomial interpolation for 10 N  . This observation is consistent with the selection of eight samples in attitude/navigation computation in [8][9][10][11]. Next, we examine the polynomial and EFH rational interpolation reconstruction quality for a trigonometric signal . It appears that the polynomial interpolation also has the same feature as reported in [16], namely, its errors or Lebesgue constants are big near both ends but much smaller in the middle part. Therefore, we apply the EFH's BAC strategy to the polynomial interpolation and plot the result in Fig. 2    the EFH rational and polynomial-BAC interpolations become significantly better than the normal polynomial for 20 N d   . Compared with the EFH rational interpolation, the polynomial-BAC interpolation has comparable Lebesgue function values across the time interval and is slightly worse in accuracy by no more than one order of magnitude. Figure 3 further assesses the interpolation performance with samples error-contaminated by a Gaussian noise of standard deviation 0.01. It can be seen from Figs. 2-3 that the Lebesgue constants well characterize the error behaviors of the polynomial and EFH rational interpolations, and the accuracy advantage of EFH rational interpolation for 20 N d   is still prominent in this error-contaminated scenario, so is the polynomial-BAC interpolation. In other words, the signal reconstruction quality at the presence of rounding/measurement errors is mostly governed by the Lebesgue constant as well. The advantage of EFH rational and polynomial-BAC interpolations would not be overwhelmed by the rounding/measurement errors.  Fig. 2, while a great discrepancy is observed for that of 8 N d   . The normal polynomial seems unable to cope with the signal of higher frequency by using only 1 9 N   samples, i.e., the polynomial order is 8, while the polynomial-BAC interpolation reaches a high order of 24 by using 2 1 25 N d    samples. As the EFH rational interpolant is a blend of polynomials of order 8, the insufficient representation of the normal polynomial is also shared by the EFH rational interpolant. Note the overlapped spikes between polynomial and EFH rational interpolations.

III. APPLICATION TO ANGULAR VELOCITY RECONSTRUCTION
AND ATTITUDE COMPUTATION This section applies the EFH rational and polynomial-BAC interpolations to the functional iterative attitude computation [8][9][10]19]. Gyroscopes have two common output forms: angular velocity measurement and angular increment (integrated angular velocity) measurement. This paper takes the latter case as a demonstration.

A. Angular Velocity Reconstruction
Assume discrete angular increment measurements As the reconstructed integrated angular velocity are normal polynomial or rational functions, we need to map the current time interval and then approximate the reconstructed function via Chebyshev polynomials, as required by attitude computation with the functional iteration approach [8][9][10][11]. The Chebyshev polynomial approximations of the reconstructed integrated angular velocity can be readily given as where   i F x is the i th -degree Chebyshev polynomial of the first kind, and n  is the maximum order of Chebyshev polynomials [20]  Polynomial Rational-EFH Polynomial-BAC interpolation, it is enough for n  to be no more than 2 N d  . The Chebyshev coefficients in (12) can be obtained by [20]  is the Kronecker delta function, yielding 1 for 0 i  and zero otherwise. Exact coefficients could be obtained only if the number of summation terms, P  , approaches infinity. The angular velocity estimate is to be acquired by time derivative of the integrated angular velocity function in (12), with the aid of the derivative property of the Chebyshev polynomial of the first kind [26] where        Figure 5 illustrates the flowchart of the angular velocity polynomial/rational reconstruction and the Chebyshev polynomial approximation, which is to be used as the input of the subsequent attitude computation by functional iteration.

B. Numerical Results of Attitude Computation
The classical coning motion, widely accepted as the standard algorithm criterion, is used here to evaluate the interpolation quality and its effect on attitude computational accuracy. The angular velocity of the classical coning motion is described by In the above,  denotes the coning angle and 2 c f    denotes the angular frequency and c f is the coning frequency. The angular increment measurement is assumed and the sampling rate is nominally set to 1000 s f  Hz. The following principal angle metric is used to quantify the attitude computation error where q denotes the quaternion estimate computed by attitude algorithms,  means the quaternion multiplication, and the operator   2:4  takes the vector part of the error quaternion. The attitude computation by functional iterative approach is exemplified by QuatFIter [10] hereafter. Figure 6 compares the attitude errors of polynomial, EFH rational and polynomial-BAC angular velocity reconstruction for the first ten seconds of coning motion. The number of gyro samples is set to 8 N d   with 1   degree and 50 c f  Hz. Both EFH rational and polynomial-BAC interpolations lead to smaller attitude errors by about two orders in the middle part of the interval. Specifically, there exists accuracy degradation in the first and the last computing interval of eight samples in which the normal polynomial interpolation has to be used for the sake of being at boundary intervals. The subfigure on the right hand presents the angular error if the accuracy assessment is performed by excluding the two boundary intervals. It is interesting to see that the polynomial-BAC interpolation outperforms the EFH rational Later d samples at BAC-based interpolation (11) Integrated angular velocity Chebyshev poly approx.
Angular velocity Chebyshev poly approx.
Attitude computation by functional iteration

Block of gyro samples
interpolation even by about four orders, while the latter, as a blend of polynomials, has a similar error shape with the normal polynomial. It should be stressed that excluding the boundary intervals for accuracy assessment is not over-optimistic for those applications that the motion is benign at the start and in the end. Figure 7 plots the attitude error trend as a function of the coning frequency ratio relative to the sampling rate, namely c s f f . The coning frequency ranges from 1 Hz to 200 Hz. The EFH rational and polynomial-BAC interpolations do not bring accuracy benefit relative to the polynomial interpolation until the coning frequency is larger than 50 Hz (Fig. 7, left subfigure). In contrast, when the accuracy evaluation starts from the second computation interval, the advantage of the EFH rational and polynomial-BAC interpolations is prominent, and notably the polynomial-BAC interpolation gaining about five orders of magnitude in accuracy as compared with the polynomial interpolation (Fig. 7, right subfigure). Specifically, the attitude error of the polynomial interpolation begins to quickly rise from 10 Hz onwards, while that of the polynomial-BAC interpolation does not significantly increase until 100 Hz and after. This turning frequency could be further postponed by using larger N, the number of gyro samples in one computation interval. As shown in the right subfigure of Fig.  7, 16 N d   reduces the attitude error at 200 Hz by nearly six orders of magnitude relative to 8 N d   . In other words, the EFH rational and polynomial-BAC interpolations, especially the latter, help enhance the capability in coping with attitude computation of highly-dynamic motions. Finally, practical situations with noisy gyroscope measurement are next investigated, as the high-order/sample algorithms tend to be much more sensitive to narrow-band noises that might lead to pseudo-coning [27].    and 10 degrees, respectively. The first case simulates high vibration environment and the second case mimics highly dynamic motions like those experienced by smart munitions [28]. A common set of random gyroscope noises are generated and fed to all interpolation methods for a fair comparison. Figure 8 plots the results for errorcontaminated measurements while excluding the boundary intervals in accuracy assessment. The considered coning frequency ranges from 10 Hz to 200 Hz. The error behaviors are quite the same in both cases. For coning frequencies below 60 Hz, the gyroscope measurement errors dominate the attitude computation accuracy and the attitude errors are flat lines for all interpolations. The normal and EFH rational interpolations begin to rise at 60 Hz and 150 Hz, respectively, where the computation errors start to take the lead, while the polynomial-BAC interpolation just keeps flat in the whole range of considered coning frequencies. Compared with the normal polynomial interpolation at the coning frequency of 200 Hz, the EFH rational interpolation reduces the attitude errors by about two orders, and the polynomial-BAC interpolation reduces by about four orders, totally insensitive to the relative coning frequency in the considered range. It demonstrates the powerfulness of the BAC-based interpolations in depressing the Runge effect in attitude computation.

IV. CONCLUSIONS
The Runge effect is a frequently-encountered problem in interpolations on equispaced samples. This work reviews its underlying rationale, as well as the polynomial and rational interpolants in the barycentric form. The BAC strategy using additional measurements in the neighbouring intervals is shown effective to reduce the Lebesgue constant of the rational interpolant. Inspired by this observation, we investigate the BAC strategy in polynomial interpolation and find that the resultant polynomial-BAC interpolation is even more effective in depressing the Runge effect. Both the EFH rational and polynomial-BAC interpolations are applied to reconstruct the angular velocity from gyroscope angular incremental measurements for attitude computation using the recent functional iteration approach. Results under the classical coning motion show that the polynomial-BAC interpolation outperforms the EFH rational interpolation in the context of attitude computation, and can significantly enhance the capability in coping with attitude computation for sustained highly-dynamic motions. It would also potentially bring significant benefits to the inertial navigation computation.