Attitude Estimation Using Kalman Filtering : External Acceleration Compensation Considerations

Attitude estimation is often inaccurate during highly dynamic motion due to the external acceleration. This paper proposes extended Kalman filter-based attitude estimation using a new algorithm to overcome the external acceleration. This algorithm is based on an external acceleration compensation model to be used as a modifying parameter in adjusting the measurement noise covariance matrix of the extended Kalman filter. The experiment was conducted to verify the estimation accuracy, that is, one-axis andmultiple axes sensormovement. Five approacheswere used to test the estimation of the attitude: (1) theKF-basedmodel without compensating for external acceleration, (2) the proposed KF-based model which employs the external acceleration compensation model, (3) the two-step KF using weighted-based switching approach, (4) the KF-based model which uses the threshold-based approach, and (5) the KF-based model which uses the threshold-based approach combined with a softened part approach. The proposed algorithm showed high effectiveness during the one-axis test. When the testing conditions employed multiple axes, the estimation accuracy increased using the proposed approach and exhibited external acceleration rejection at the right timing. The proposed algorithm has fewer parameters that need to be set at the expense of the sharpness of signal edge transition.


Introduction
The inertial measurement unit (IMU) is typically used to determine attitude, that is, roll and pitch, by fusing accelerometer and gyro data.The most notable disturbance of attitude determination is external acceleration [1], which is caused by a change of velocity in magnitude or direction [2].Recently, the IMU has been widely used in many applications, such as gait assessment systems [3][4][5]; classification of human movement and falls [6,7]; and specific modeling on certain body movements-for example, sit-to-stand [8,9].Attitude estimation using IMU is an important theme in the literature, specifically on how to improve its accuracy and reliability.This leads us to the aim of this study, which is to present a new algorithm to overcome the existence of external acceleration in order to improve attitude estimation accuracy.
The attitude solution provided by the gyro is prone to being unbounded, to bias, and to random-walk errors.In static or slow movement, the accelerometer measures roll and pitch by leveling to correct the gyro-unbounded error.This is due to the trustworthiness of the gravitational measurement.Therefore, a proper fusion of IMU data and the algorithm to compensate for external acceleration is needed to overcome the shortcomings of each sensor and the effect of external acceleration.The fusion technique evolved along two major paths: one approach incorporates the use of a Kalman filter [10][11][12] while the other algorithm consists of a complementary filter [13,14].
Suh et al. [1], Sabatini [11], and Harada et al. [12] proposed the adapted measurement noise covariance matrix to overcome some disturbances, including external acceleration.The adapted algorithm in [1] is based on the weighted-switching approach.In contrast, Sabatini in [11] used a threshold-based approach, while Harada et al. [12] proposed a threshold-based approach combined with softened part.
Inspired by Suh et al. [1], Sabatini [11], and Harada et al. [12], this paper proposed the extended Kalman filter-(EKF-) based attitude estimation using a new algorithm to overcome the existence of external acceleration.We extracted external acceleration from the accelerometer signal with a model and used it to adjust the measurement noise covariance matrix.The aforementioned algorithm is verified by conducting oneaxis and multiple axes movement experiments.
The rest of the paper is organized as follows: Section 2 describes attitude determination and the proposed algorithm.Section 3 discusses the experimental setup and results.Section 4 then presents a discussion of the algorithm approaches and a comparison regarding the experimental results.In Section 5, we present our conclusion.

Method
2.1.Sensor Model.The relative attitude (roll and pitch) of two frames of reference used a direction cosine matrix (rotation matrix) for coordinate transformation or rotation of object. (roll),  (pitch), and  (yaw) are the rotation angles about the -, -, and -axes on the sensor body, as shown in Figure 1. / indicates a rotation matrix that represents the relative attitude of frame  (inertial frame) with respect to frame  (body of sensor), thus allowing the relation as follows: where x  is the coordinate vector in frame  and x  is the coordinate in frame .
The well-known rotation matrix, with the function of  →  →  Euler angles, is as follows: where  and  are the short forms of sine and cosine, respectively.Measurements model for attitude determination using accelerometers is constructed by the third row of  / that contains  and .It is independent of rotation along the axis ().
Sensor signal from the accelerometer is modeled as follows: where   is a vector of   ,   , and   ; the acceleration is measured by the three-axis accelerometer on the body frame.
The unwanted external acceleration information () in the accelerometer will be estimated.This acceleration is caused by the change of velocity in magnitude and/or direction.The external acceleration is the main source of errors in attitude estimation.

External Acceleration Compensation Method.
In this section, we introduce the method for compensating for the external acceleration.Our method is based on the view that external acceleration is a disturbance in the original signal and therefore should be removed from the original signal.The disturbance rejection model was derived from a common control ratio model [15].From this point on, the term "compensation model" will be used to replace the term "disturbance rejection model."The three-axis accelerometer data is input for the compensation model.Each axis will be treated in the same step; therefore, for the sake of convenience, only the -axis will be discussed.In (4), ideally the compensated -axis accelerometer signal ( *   ) was determined by subtracting the compensation model in the axis (  ) from the original signal (  ): The external acceleration (  ) is extracted from   by the compensation model.Figure 2 shows the ideal model diagram of (4).However, due to the dynamism of the system and many uncertainties, there is no ideal model.Consequently, we will use the result of the compensation model as one parameter for adjusting the noise covariance matrix in sensor fusion in the next section.The design approach for a compensation model is based on the specification that the model should not be affected by the disturbance input ().In other words, the steady-state output should be zero if possible or a small constant; () ss = 0, () ̸ = 0.A control ratio in the  domain [15] is commonly written as follows: Assuming a step disturbance signal () =  0 /, the output will be where  is gain constant,  > 0,  and  are constant coefficients, and  0 is the final value of step disturbance.The condition () ss = 0 will be achieved if  0 = 0 and requires the numerator of ()/() to have at least one zero at the origin of the -plane.On the other side, the location of the poles on the -plane determined its corresponding response.
To produce the damped response, the poles should be on the left-half -plane.There are two possibilities for poles locations on the left-half side, such as complex-conjugate poles and real poles, which have the output characteristics of exponentially damped sinusoid and damped exponential response, respectively.In our model, we chose damped exponential response because the disturbance should be suppressed without a sinusoidally oscillating component.Finally, we determined the transfer function of our model: there is one zero at the origin and one real pole on the left-half of the -plane which is based on (5) and becomes as follows: The compensating model in (7) attenuates low frequencies and as a result a high frequency signal is obtained.This is in line with the fact that the external acceleration occurs when the change of velocity in magnitude and/or direction of each axis exceeds a certain limit.
The analysis of the stability of model ( 7) will be introduced in the time domain as follows.In the time domain, the stability will be achieved based on the principle that there is a limited response to the limited input.We discuss the time response of singularity input functions to this model, such as step, ramp, rectangular pulse, and parabolic functions.The time response model (7) for a unit step input, (), is () =  − 0  .We get steady-state response for () and the steadystate error as in err The time response for ramp input (() = ) is () = / 0 (1−  − 0  ).The steady-state output and error are calculated as follows: A rectangular pulse of amplitude 1/ and of duration , () = (() − ( − ))/, has a zero steady-state output and also a zero steady-state error.For a parabolic function input, () = 0.5 2 , the time response is calculated as follows, with an infinity steady-state output and infinity steady-state error: These four kinds of singularity inputs and their combination represent the possibility of disturbance waveforms.All of them have the tendency to be rejected and damped by model (7).The response of the parabolic input has not shown constant steady-state output; however, it represented the exponentially damped factor.Figures 3(a), 3(b), 3(c), and 3(d) illustrate all time responses for each basic input signal, using model ( 7) within  = 1.Each figure has two outputs, that is,  1 () and  2 (), where  0 in  1 () is greater than  0 in  2 ().

Attitude Determination through Sensor Fusion Method.
Attitude determination from gyro had errors accumulating along the process of numerical integration, but the gyro was sensitive to the change of rotation.While the attitude determination by the accelerometer does not show accumulation errors or divergence, the error in the roll and pitch axis is too large because those two axes are mutually affected.When there is a movement in the roll axis, the pitch axis also moves up and down with respect to the horizontal plane, and vice versa.In the proposed method, both sensors are complementary to each other to achieve performance in the sensor fusion technique using the extended Kalman filter (EKF).In this paper, the work will be divided into five modes.Mode 1 is the standard EKF using the proposed system model but without compensating the external acceleration.Mode 2 is the proposed system model, also containing the proposed external acceleration compensation, with a preprocessing technique for accelerometer data using model (7) before the EKF process.The difference between Mode 2 and Mode 1 is the addition of the measurement error covariance based on the external acceleration compensation model.A detailed explanation of Mode 1 and Mode 2 is in Section 2.4.Mode 3 is derived from [1]: the two-step EKF which adjusts the measurement noise covariance using the weighted-switching approach.Mode 4 is the quaternion-based EKF, which uses a threshold-based approach to set the measurement error covariance.This mode is derived from [11].Mode 5 is derived from [12], which uses an unscented Kalman filter (UKF).This method works as threshold-based approach combined with a softened part of measurement noise covariance matrix adjustment.

Algorithm Description of Mode 1 and Mode 2
2.4.1.Process and Measurement Models.We initially defined the states and the observation variables for the system model.We set attitude and gyro bias as state variables, since the bias errors are a highly complex function to the ambient temperature.The Euler angle was the angle representation.The state variable () and the measurement variable () are defined as follows: where  (roll) and  (pitch) are the rotation angles about the and -axes and  is yaw angle, but it is not of concern in this study.These come from the integration of the rate of change from gyros, while   ,   , and   are biases from gyro in -, -, and -axis, respectively.We use the measurement from the accelerometer in order to calculate   ,   , and   as measurement variables.The system equation is given by ẋ () =  ( ()) +  () , where (()) is a nonlinear function representing the relation between gyros data and kinematic equation for the Euler angles ( →  → );  is state to measurement matrix as shown in (14); and () and V() are process noise and measurement noise, respectively, which are assumed to be uncorrelated Gaussian distributed white noise.Consider Measurement update: In (14), , , and  are angular velocities measured from the gyro along the -, -, and -axes, respectively.Symbols of , , , and  are short forms of sine, cosine, tangent, and secant, respectively.This process model is nonlinear since it contains trigonometric functions.Therefore, we use Jacobian values of the nonlinear model of () to replace state transition matrix (), where the values in that matrix are obtained by applying the previous estimated states in the EKF process.In contrast, matrix  is the linear expression.
The measurement from the accelerometer which is used to calculate   and   as measurement variables is a common initialized leveling equation [2,16]: Figure 4 shows the structure of the proposed algorithm combined with the extended Kalman filter algorithm.The notations and meanings of x , x−  ,   , and  −  are estimates of state, prediction of state, estimate of the error covariance, and the prediction of the error covariance, respectively.

Error Covariance.
In Figure 4, the estimated error covariance is computed by   =  −  −    −  .Error covariance indicates the difference between state estimation (x  ) and the unknown true state value (  ).It also can be defined as the following [17]: The error covariance of attitude states was initialized by calculating the angle variance using (15).The rest of the error covariances of the bias states were approximated from gyro data.All the variances were calculated using ( 17) and (18).Equation ( 17) was used to calculate the estimated mean value and ( 18) was used to calculate its covariance as follows: where  is a number of samples and ⃗   is a vector of accelerometer data on the -, -, and -axes.
All calculations were taken from IMU and placed on the table without movement for approximately one minute with sampling period  = 0.05 sec.Figure 5 shows the accelerometer and gyroscope signals in static conditions on all three axes.The experimentally determined variance values for attitude states are  2  = 6.05 − 6 and  2  = 4.98 − 6 in rad 2 .And for the gyroscope along the -, -, and -axes, the values are Figure 5: Three-axis gyroscope and accelerometer in static condition.

Initial Values and Covariance of the Noise.
The initial value of the bias gyro was determined heuristically through static IMU data in rad/s unit.The initial value for attitude states was chosen as one degree for each state.The initial states value in EKF model is as follows: Since the noise is assumed to be in normal distribution and independent on each axis, we use the variance of noise () and V().Noise process matrix () was selected heuristically as  11 ,  22 , and  33 are 1 and  44 ,  55 , and  66 are 0.3.
The measurement noise covariance matrix ( nom ) is defined by (20) and the measurement comes from the attitude in (15); however, the measurement noise covariance matrix is obtained by calculating the variance of accelerometer signal through (17) and (18).These signals have a direct relation to the attitude measurement: To determine the diagonal element of  nom , the variances of   ,   , and   are calculated:  11 = 9.04 × 10 −6 ,  22 = 1.04 × 10 −5 , and  33 = 1.80 × 10 −5 g 2 .Since the accelerometer is sensitive to disturbances, for the sake of practicality we tuned the value of  nom based on the application.Related to the proposed method,   from the external acceleration model is uncorrelated to noise V  , and the discrete form of  is now modified to be as follows: where  is a constant coefficient that is heuristically chosen based on the application; in Mode 1,  is set to zero.Particularly in Mode 2,  is tuned based on the external acceleration model and  2  is the element of the diagonal matrix, which consists of  2  ,  2  , and  2  .The diagonal matrix of  2  in (21) plays a key role in compensating for the external acceleration by increasing the   value.During static conditions or slow movement (i.e.,  2  ≈ 0), the covariance matrix will be the same as  nom .
The above process and measurement models construct the procedure of the EKF as follows: (1) Set the initial values for states and error covariance x0 ,  0 . ( (2) Predict states and error covariance (3) Compute the Kalman gain (4) Compute the states estimate (5) Compute the error covariance where In Mode 3, there will be a setting for the measurement error covariance matrix, in which  1 and  2 belong to the accelerometer and  3 belongs to the gyroscope, as in The two-step measurement updated Kalman filter was used to update the Kalman gain, state estimation, and covariance matrix.First, the measurement was updated only for variance of the error measurement in the gyroscope ( 3 ).Second, the measurement was updated for  1 and  2 using a weightedswitching rule as in (32), based on the threshold rule of (31) to detect external acceleration.The threshold rule is derived from a necessary condition for acceleration free movement, as in (30).The weighted-switching rule works on the principle that when there is external acceleration, the gyroscope outputs should be trusted more.This can be done by making  1 and  2 larger.For a more detailed explanation, refer to [1].Consider the following: The mechanism of adaptation of the measurement noise covariance matrix of acceleration is implemented using a threshold-based approach: The measured acceleration magnitude is tested in advance to recognize deviations from gravity.If there is deviation greater than the threshold, the observation variance   2  is set to a high value.accelerometer is detected as reliable if it satisfies the following condition: |‖ ⃗   ‖ − ‖ ⃗ ‖| <  th , where ⃗   is the output of the accelerometer on the body frame.The mechanism of the adaptation of the measurement covariance matrix is as follows: where  0 is the offset and ] . (35)

Various Test Conditions.
We have five modes of attitude estimation, as stated in Section 2.3.There will be three tests in each mode, aimed at studying the performance under various dynamic conditions.Test A was performed by placing the IMU sensor on the slider table of a MISUMI RSH3 single-axis robot.By moving the slider back and forth with the robot controller, external acceleration is applied on the -axis of the IMU sensor.Test A intended to test proposed method on the lateral movement as is done in [1], which will affect the pitch angle.In Test A, we conducted the trial by using such conditions as follows: (1) four acceleration coefficients in ascending order from the MISUMI software settings, that is, 0.1, 0.75, 1.5, and 2.5 m/s 2 ; (2) four  0 parameter values, that is, 0.5, 0.1, 0.05, and 0.01; (3) six  values, that is, 35, 80, 100, 150, 170, and 200.
Test B involved movement using free hand in the and -axes of the IMU as performed in our previous study [13] by mimicking the pendulum swing movement.The movement along the -axis was considered as roll and the movement along the -axis was considered as pitch angle.In contrast with Test A, Test B was performed using external acceleration on a combination of axes.
Test C was executed based on the application of attitude estimation on a shoe-type measurement device.An IMU sensor was placed on the top of one shoe.The subject was asked to walk straight, 3-4 strides, at a normal speed.

Test Setup.
To verify and validate the proposed method, all tests were conducted using an IMU sensor consisting of an accelerometer (±16 g), a gyro (±1500 deg/s), and a magnetometer (±0.9 Ga) (from Logical Product, Japan).Sensor data were transmitted to the PC wirelessly.A combination of NI LabVIEW5 and Simulink5 was used for data acquisition and the MATLAB5 program was used to execute the proposed method.In Test A, the MISUMI RSH3 single-axis robot (from MISUMI Group Inc., Japan) had a max speed of 300 mm/s, had an effective stroke of 500 mm, and was controlled by the RS-Manager support software used as a testbed (see Figures 6(a) and 6(b)).Note that the attitude is not changed because the IMU sensor was moving in a lateral direction without rotation; that is, in order to validate the proposed method, the reference values are  = 0 ∘ and  = 0 ∘ .
In Test B, the IMU sensor was placed in plastic jar surrounded by styrofoam to avoid magnetic interference between the IMU and the electromagnetic motion tracking system receiver (Fastrak from Polhemus Inc., USA).We used Fastrak as the reference by recording the roll and pitch angles of the receiver using C# data acquisition program.Fastrak attitude data were transmitted to the PC via a cable using an RS-232 protocol (see Figure 7).
In Test C, we mounted the IMU sensor on the shoe.For a reference measurement, the experiment setup contained six OptiTrack5 cameras, and four reflective markers were placed at the forefoot and heel.Figures 8(a) and 8(b) show the measurement setup for validating Test C. In Tests B and C, we calculated cross-correlation in order to compare references and IMU data from different data acquisition programs.The time when the correlation was maximal will be used to synchronize both measurements.To accomplish time normalization between the reference and the IMU sensor, cubic spline data interpolation was employed.For the purpose of looking at the effect of the reduction of external acceleration, the Mode 2 experiment was performed using a different compensating parameter value.Table 2 shows the average MSE of the attitude estimation from four different accelerations, with some variation in the values of  0 and .The variation value of  0 is 0.5, 0.1, 0.05, and 0.01, all of them using  = 1.While the variation value of  is 35, 80, 100, 150, 170, and 200, in ascending order, the value settings of the EKF parameters are  nom = 1.5.The relationship between  and  0 in one-axis movement test is as follows.The MSE trend decreases as the value of  increases or  0 decreases.The best combination for one-axis movement in this experiment is  = 200 and  0 = 0.01; the MSE is 8.57 and 1.48 deg 2 for pitch and roll, respectively.
The influence of the presented  0 and  for Mode 2 in Table 2 was also displayed as a column bar in Figures 9 and  10 for the pitch and roll, respectively.In the meantime, the MSE difference -test [18] between Mode 2 and Mode 1, Mode 2 and Mode 3, and Mode 2 and Mode 4 as well as that between Mode 2 and Mode 5 was evolved in those figures.In Figures 9  and 10, "a," "b," "c," and "d

Test B.
Test B was executed by rotating the IMU sensor in a random manner by hand along the and -axis.In order to determine the optimal  0 value for Mode 2, we conduct the variability of the  0 test to calculate the MSE and maximum error along the timeline of the whole experiment.Table 4 presents the MSE and maximum errors using different values of  0 .It was found that  0 = 0.05 has an optimum performance among all the possible choices.Figures 12(a Figures 13(a), 13(b), and 13(c) present the accelerometer signal and Figures 13(d), 13(e), and 13(f) present the compensation model signal along the -, -, and -axis of the accelerometer, using model parameters  = 1 and  0 = 0.05.When the sensor is not experiencing external acceleration, the compensating signal tends to be around zero.While the sensor is being moved quickly, it experiences two major external acceleration periods along the and -axes as shown in Figures 13(a) and 13(b).In Figure 13(a), we call the periods Phase I and Phase II (indicated by the dashdot line); Phase I is for 2.4 s during 5.8-8.q 0 = 0.5 q 0 = 0.1 q 0 = 0.05 q 0 = 0.01 Mode 4 Attitude reference signal from Fastrak is shown in Figures 15(a) and 15(g) for pitch and roll, respectively.The quantitative evaluation using MSE between each mode at the time the external acceleration occurred, that is, during Phase I to Phase IV, is shown in Table 5. Figures 16(a 5 for the MSE of pitch and roll, respectively.In Table 6 and Figure 17 we conclude the MSE result of Test B, the average MSE of Phases I and II as pitch estimation error, and the average of Phases III and IV as roll estimation error.
Some parameters values are used in Test C for each mode.Mode 1 uses  nom = 1.5, while Mode 2 uses  0 = 20,  nom = 1.5, and  = 34.Mode 3 uses thresholds  = 0.2 g 2 ,  1 = 0.8,  2 = 0.3, and  1 = 0.093,  2 = 0.034, and  3 = 1.5.Mode 4 uses the thresholds  th = 0.8 and  2  = 60.Mode 5 works using  2 0 = 1.5 and the factor of softened part   = 50.3 shows that the MSE of pitch and roll estimation is 25.25 and 3.65 deg 2 , respectively.The attitude estimation accuracy of Mode 1 is lower than Mode 2. This result confirms the findings presented in the literature [1], where the standard EKF has lower accuracy in the one-axis test type.

Discussion
The result of Mode 1 in Test B, as shown in Figures 15(b) and 15(h), demonstrates that the effect of external acceleration during Phases I to IV on the pitch and roll estimation is still dominant.The final result of Mode 1 in Test B as shown in Table 6 and Figure 17 indicates that Mode 1 estimation accuracy is also lower than Mode 2, Mode 3, and Mode 5. Kalman filtering employs a compensation mechanism that surpasses Mode 1's accuracy, except in Mode 4. Mode 1 has a MSE of pitch 36% greater than that of Mode 2, 26% greater than that of Mode 3, and 43% greater than that of Mode 5.The result of Mode 1 in Test C is also in line with the result in Test B, which is that Mode 1 estimation accuracy is lower than that in Modes 2, 3, and 5.

Mode 2.
The proposed algorithm in Mode 2 involved the external acceleration compensation models (7) and covariance matrix updating process (21), which plays an important x (q 0 = 0.1) x (q 0 = 0.3)    role in improving the attitude estimation accuracy during fast movement.If  = 1 was imposed for all experiments, then the only setting parameters are  0 and .
The two parameters  0 and  are chosen heuristically.However, there is some consideration when choosing these parameters.As shown in (7), a smaller  0 reduces the cutoff frequency.The value of  0 that determines the cutoff frequency affects the value of the compensating model output,   , where the effect continues to the value of the measurement noise covariance matrix () in attitude estimation.The impact of the  0 variability to the square output of   along the -axis was presented in Figures 12(a) and 12(b).It appears that the smaller value of  0 has an increasing amplitude of  2  , which will simultaneously increase the value of the measurement error covariance matrix.In the meantime, the effect of the  is as an amplified factor for the softened part ( 2  ).The increasing value of  simultaneously increases the value matrix , which reduces the Kalman gain.Therefore, the contribution of the measurement to the estimation process also decreases; thus, the estimate is less affected by measurement from an accelerometer.
The result of Mode 2 in Test A is shown in Figures 9 and 10 and Table 3 as a concluding table.As can be seen, the Mode 2 MSE of pitch is 8.57 and 1.48 deg 2 , which is the lowest among all those tested.The significant difference with Mode 1 is generally started at  = 150 and with  0 up to 0.1.
For Test A, the variability of  0 influences the estimation error; however, the attitude reference is always zero degrees, and therefore we cannot observe the effect of the  0 variability to the attenuation of some important signal by model (7).Therefore, we need a test like Test B and Test C to observe the variability effect of  0 on the appearance of the signal.In Test B, quantitative assessment of MSE on pitch and roll is done in Phases I to IV, as shown in Table 6 and Figure 17.Mode 2 has a MSE value of 34.36 and 30.24 degree 2 for pitch and roll estimation, respectively.The MSE of pitch is lower than that in Modes 1, 3, and 4, and Mode 5's accuracy surpasses the accuracy of Mode 2. In roll estimation, the accuracy of Mode 2 outperformed the accuracy of Modes 1 and 4. It was observed that Mode 2 is able to reduce the effect of external acceleration on the correct timing when it occurs, that is, during Phases I, II, III, and IV.The results of real world application, such as in Test C, are in Table 7 and Figures 18(a 6.It is important that we do not make the criticism that Mode 2 is superior to Mode 3 in Test B, since different settings on the combinations of ,  1 , and  2 as in (31) and (32) might result in a better performance.
The result in Table 7 for Test C indicates that Mode 3 outperformed all modes.The number of setting parameters in Mode 3 is higher than in other modes (i.e., three setting parameters), which provides a more softened setting.However, this also requires more effort than in the other modes.Our proposed algorithm in Mode 2 has fewer parameter settings, that is,  0 and .Furthermore, the execution of the compensation algorithm is fully dependent on the existence and magnitude of the square of external acceleration model, as in ( 7) and ( 21), rather than threshold-based approach, as in (31).However, we observed that the measurement model in Mode 3 (as in (28)) has an advantage over Mode 2. This model incorporates the data from gyro as well as accelerometer data.Therefore, the setting mechanism of matrix  consists of two parts: accelerometer and gyro.When the level of trust in the accelerometer lowers then it is possible to set the level of trust higher in gyro.

Mode 4.
Of all the tests, the estimation accuracy of Mode 4 is the lowest.We suspect that this is caused by the presence of magnetometer measurement vectors in the measurement model of Mode 4. The earth's magnetic vector,  , has a magnitude that is always changing over a large range of time [19].Furthermore, the experiment room we used is not guaranteed to be free of magnetic interference and soft    iron distortion.In [19] it is suggested that, to overcome this problem, initialization must be done carefully to find out the exact magnitude and orientation of the magnetic field.This vector can be used during the experiment.for pitch and roll, respectively (as shown in Table 6).This also happened in Test C. As shown in Table 7, the estimation accuracy of Mode 5 is 4.5% and 6% over Mode 2's in pitch and roll estimation, respectively.In Test A, Mode 2 outperforms Mode 5 in all  and  0 combinations.The difference between Mode 2 and Mode 5 is in the softened part of the measurement error covariance matrix.Mode 5 in (34) uses the absolute difference between measurement and predicted acceleration.However, in Mode 2 we use the model of external acceleration as the softened part.

Experimental Limitation and Future Work.
In all tests, the proposed Mode 2 outperformed Mode 1. From these results, we were able to ensure that a mechanism of external acceleration compensation has the influence to improve estimation accuracy.Even though there is an advantage to improving the estimation accuracy, some of the major limitations to the experiments will be described.
First is the limitation of measurement model.Modes 1, 2, and 3 used Euler representation.The model in Mode 1 and Mode 2 in (12) did not include the measurement from a gyroscope.With a slow motion sensor this measurement model might be not a problem, because it is not necessary to compensate the external acceleration.The setting of the measurement noise covariance matrix () in (21) primarily relies on the roll and pitch from accelerometer data.When the value of  becomes larger due to the presence of external acceleration, the estimation process in Kalman filtering is less affected by accelerometer, but at the same time we cannot increase the level of trust in the gyroscope.Furthermore, it is important to consider a modification of the measurement model in future work.
The second limitation is that the proposed external acceleration compensation in (7) did not include automatic calculation of parameter  0 .The model in (7) is used as an additional part for the measurement noise covariance in (21).This model relies on the frequency of the application that the model will use.Before using the model, the determination of the  0 value is important.One prospective improvement of this model is the additional step of calculating the optimum value of  0 from IMU data before the Kalman filtering process.
The third limitation is considering the application of Mode 2 whenever the pitch angle () reached the /2 radian, even though there is an advantage of Euler representation over the quaternion representation.The limitation of Euler representation in this experiment is that whenever pitch angle () reaches /2, the state in Mode 2 will be singular.This limitation need not be considered as long as the application is still able to accept the range of angle.
The fourth limitation is related to the linearization process in the process model of Mode 2, which leads to firstorder approximation error.Mode 2 uses EKF-based filtering that employs linearizing the nonlinear model.The first-order linearization might be the cause of degraded accuracy in all modes that employed EKF, that is, Mode 1, Mode 2, Mode 3, and Mode 4.However, Mode 5 employs an unscented Kalman filter (UKF) that is free from linearizing through Jacobian matrix.As a result, Mode 5 outperformed almost all modes in Test B and outperformed Mode 2 in Test C. One consideration to avoid the first-order approximation error besides using UKF is using a direction cosine matrix (DCM) method.

Conclusion
In this paper, the main contribution is the algorithm for external acceleration compensation, which aims to improve the attitude estimation accuracy.A Kalman filter-based attitude estimation using a compensating algorithm has been discussed.The experiment was performed using three types of tests: movement on one axis, dynamic movement using free-hand movement, and walking.We employed five different approaches to deal with the dynamic test and the proposed method is placed on Mode 2. The first approach is the standard KF model, without using external acceleration compensation (Mode 1).The second approach is the modified KF model, using the proposed compensating procedure (Mode 2); the third is a weighted-switching method (Mode 3); the fourth is a quaternion-based EKF using a threshold-based method (Mode 4); and the fifth (Mode 5) uses an unscented Kalman filter and is threshold-based combined with a softened part.
The experiment results showed that, by using the external acceleration compensation process, the estimation accuracy of the proposed algorithm is improved when compared with the standard EKF procedure in Mode 1 in all tests.Mode 2 also outperformed all modes in Test A by using the optimal parameter setting.In dynamic Test B, Mode 5 outperformed other modes; we suspect that this is caused by using UKF in Mode 5.The UKF was free from first-order approximation of a nonlinear system.The advantage of Mode 3 over other modes is presented in Test C. Mode 3 used a measurement model that included accelerometer and gyroscope data, while the measurement model in Mode 2 was related to accelerometer data.
There is a lack of efficiency comparison to some modes in the experiments.Compared to the other modes, the advantage of Mode 2 over Mode 3 is the number of parameters set; Mode 2 has fewer parameters.Mode 3 takes a twostep EKF, which leads to additional computational overhead.The advantage of Mode 2 over Mode 4 and Mode 5 is the parameterized spatial rotation; quaternion as used in Mode 4 and Mode 5 is hardly used because it is a burden to update its four variables [19].However, Euler needs to update two variables.Specifically, comparing Mode 2 and Mode 5, their effectiveness has been unexplored; the computational time of the extended Kalman filter is much lower than in the unscented Kalman filter [20].
Nevertheless, as a future problem to be solved, in order to increase the estimation accuracy potential for other applications, it needs the addition of a step that can perform adaptive parameter settings ( 0 and ) based on the present input from the IMU.Using UKF and DCM is also one consideration to improve accuracy in order to be free from linearization approximation error.

Figure 2 :
Figure 2: The input-output diagram of external acceleration compensation.

Figure 4 :
Figure 4: Structure of the proposed algorithm.

2. 7 .
System Model of Mode 5. Mode 5 is proposed by Harada et al.[12] using quaternion-based UKF.The method is threshold-based combined with a softened part of the measurement covariance matrix adaptation.The output of the

Figure 6 :
Figure 6: Test A setup for validating attitude estimation.The RSH3 table slider moves along the -axis.

Figure 7 : 2 , 1 .
Figure 7: Test B setup consisted of Fastrak5 attitude reference system and IMU inside a plastic jar.

Figure 8 :
Figure 8: Measurement setup for Test C. (a) Room setup for recording the attitude; (b) four reflective markers were placed on the forefoot and heel.

Figure 9 :
Figure 9: MSE pitch of Mode 2 in Test A: influence of the variability of  0 and ."a, b, c, d" denote significant difference from Mode 1, Mode 3, Mode 4, and Mode 5, respectively.

Figure 10 :
Figure 10: MSE roll of Mode 2 in Test A: influence of the variability of  0 and ."a, b, c, d" denote significant difference from Mode 1, Mode 3, Mode 4, and Mode 5, respectively.
) and 16(b) present the column bar of Table

4. 1 .
Mode 1. External acceleration compensation does not play a role in Mode 1. Mode 1 used  nom in the measurement noise covariance matrix, as in (20), and did not use the external acceleration compensation model.An experiment result of Mode 1 in Test A was presented in Tables 1(a) and 1(b), as well as in Figures 9 and 10 . The comparison in Table

Figure 13 :
Figure 13: External acceleration compensation: original accelerometer signal is shown in (a), (b), and (c); compensation model output is shown in (d), (e), and (f) along the -, -, and -axis, respectively.In (a), the dash-dot line indicates Phases I and II, respectively, while in (b), the dash-dot line indicates Phases III and IV, respectively.

Figure 14 :
Figure 14: (a), (b), and (c): graphs of the square of compensation model output for -, -, and -axis of Mode 2. As a comparison is shown in (d): necessary condition of Mode 3.
), 18(b), 18(c), 18(d), 18(e), 18(f), 18(g), 18(h), 18(i), 18(j), 18(k), and 18(l); the MSE of Mode 2 is 175.97 deg 2 and 245.8 deg 2 in pitch and roll estimation, respectively.This means that the average estimation error in every point calculation is around 13.3 deg and 15.7 deg.In Test C, Mode 2 accuracy outperformed the accuracy of Mode 1 and Mode 4 in pitch estimation and outperformed Mode 1 in roll estimation.Overall, for Test C, Mode 3 outperformed the accuracy of other modes on both pitch and roll estimation.

4. 3 .
Mode 3.  In our Test A and Test B experiment, Mode 3 has lower estimation accuracy in pitch estimation (a bigger MSE) as compared to Mode 2 and Mode 5, as shown in Figure9and Table

Figure 15 :
Figure 15: Test B attitude estimation result: five modes in comparison to the references (a and g) for pitch (left column) and roll (right column).

)Figure 16 : 5 Figure 17 :
Figure 16: Test B: (a) MSE pitch of five modes, (b) MSE roll of five modes during external acceleration, all in deg 2 .

4. 5 .
Mode 5.In our experiments, the accuracy of Mode 5 outperformed Mode 2 in Test B by around 8.1% and 15.3% the softened part of the adaptation mechanism.

Table 3 :
Test A: MSE of all modes, presented as mean (standard deviation).

Table 4 :
Test B: MSE () (in deg 2 ) and maximum error () (in deg) of attitude estimation of Mode 2 at various values of  0 during all periods of time.

Table 5 :
Test B: MSE () (in deg 2 ) and maximum error () (in deg) for all modes during major external acceleration periods (in detail).

Table 6 :
Test B: MSE of all modes in degree 2 .

Table 7 :
Test C: MSE in degree 2 of all modes.