Abstract
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 and multiple axes sensor movement. Five approaches were used to test the estimation of the attitude: (1) the KF-based model 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.
1. 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–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–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 one-axis 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.
2. 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 is the coordinate vector in frame and 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.
2.2. 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; , . A control ratio in the domain [15] is commonly written as follows:Assuming a step disturbance signal , the output will bewhere is gain constant, , and are constant coefficients, and is the final value of step disturbance. The condition will be achieved if 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 . We get steady-state response for and the steady-state error as inThe time response for ramp input () is . The steady-state output and error are calculated as follows:A rectangular pulse of amplitude and of duration , , has a zero steady-state output and also a zero steady-state error.
For a parabolic function input, , 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 . Each figure has two outputs, that is, and , where in is greater than in .
(a)
(b)
(c)
(d)
2.3. 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.
2.4. 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 bywhere 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 are process noise and measurement noise, respectively, which are assumed to be uncorrelated Gaussian distributed white noise. Consider
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 , , , and are estimates of state, prediction of state, estimate of the error covariance, and the prediction of the error covariance, respectively.
2.4.2. Error Covariance
In Figure 4, the estimated error covariance is computed by . Error covariance indicates the difference between state estimation () 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 and in rad2. And for the gyroscope along the -, -, and -axes, the values are , , and in rad/s2. These values will be the diagonal elements of the initial .
2.4.3. 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 . Noise process matrix () was selected heuristically as , , and are 1 and , , and are 0.3.
The measurement noise covariance matrix () 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 , the variances of , , and are calculated: = 9.04 × 10−6, = 1.04 × 10−5, and = 1.80 × 10−5 g2. Since the accelerometer is sensitive to disturbances, for the sake of practicality we tuned the value of based on the application. Related to the proposed method, from the external acceleration model is uncorrelated to noise , 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 is the element of the diagonal matrix, which consists of , , and . The diagonal matrix of in (21) plays a key role in compensating for the external acceleration by increasing the value. During static conditions or slow movement (i.e., ), the covariance matrix will be the same as .
The above process and measurement models construct the procedure of the EKF as follows:(1)Set the initial values for states and error covariance(2)Predict states and error covariance(3)Compute the Kalman gain(4)Compute the states estimate (5)Compute the error covariance
2.5. System Model of Mode 3
The two-step EKF [1] used separate measurements between the gyro and accelerometer. The external acceleration compensation technique is based on weighted-switching in setting the noise covariance of those measurements and was proposed by Suh et al.
The state for Mode 3 is and the measurement variable is . The system equation is given bywhere In Mode 3, there will be a setting for the measurement error covariance matrix, in which and belong to the accelerometer and 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 (). Second, the measurement was updated for and using a weighted-switching 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 and larger. For a more detailed explanation, refer to [1]. Consider the following:
2.6. System Model of Mode 4
Mode 4 is proposed by Sabatini [11] using quaternion-based EKF. The discrete state for Mode 4 is composed of the rotation quaternion, accelerometer, and magnetometer bias vectors . The measurement model is constructed by accelerometer and magnetometer measurement vectors . 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 is set to a high value.
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 accelerometer is detected as reliable if it satisfies the following condition: , 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 is the offset and is the softened part of the adaptation mechanism. is the predicted acceleration that comes from . If the orientation quaternion , where is the vector part and is the scalar part of the quaternion, then
3. Experimental Setup and Results
3.1. 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/s2; (2) four 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.
3.2. 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 LabVIEW® and Simulink® was used for data acquisition and the MATLAB® 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°.
(a)
(b)
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 OptiTrack® 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.
(a)
(b)
The quantitative performance assessment for Test A, Test B, and Test C is using mean squared error (MSE or ) in degree2 and maximum error () in degree between predicted attitude and reference attitude as described in
3.3. Experiment Results
3.3.1. Test A
Mode 1 is not related to the compensation model in (7), while the measurement error covariance in (21) was set to = 1.5 and λ = 0. As with Mode 1, neither the compensation model nor updating covariance were used in Mode 3, but the constant parameters were heuristic choices based on trial and error; that is, = 0.1 g2, = 0.4, and = 80. The experimental results of Mode 1 and Mode 3 are presented in Tables 1(a) and 1(b).
To follow the adaptation of measurement noise covariance matrix (33) in Mode 4, the threshold used is 0.02 and = 1.5 as used in Mode 1 and Mode 2. To meet the mechanism of external acceleration (34) in Mode 5, is 1.5, which is the same value as we used in Mode 1 and Mode 2. is determined to be 5, which comes from experiments using = 2 to 10. The MSE will decrease in value when becomes greater, although the MSE decrease is smooth; it is in line with the function of as the softened part of adaptation mechanism. The results of Mode 4 and Mode 5 are also presented in Tables 1(a) and 1(b). Overall, the results of Modes 1, 3, 4, and 5 are shown in Tables 1(a) and 1(b). The mean MSE of Mode 5 is the lowest among Modes 1, 3, and 4, that is, 23.63 and 2.82 deg2 for pitch and roll, respectively.
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 and λ. The variation value of is 0.5, 0.1, 0.05, and 0.01, all of them using . While the variation value of λ is 35, 80, 100, 150, 170, and 200, in ascending order, the value settings of the EKF parameters are = 1.5. The relationship between λ and in one-axis movement test is as follows. The MSE trend decreases as the value of λ increases or decreases. The best combination for one-axis movement in this experiment is λ = 200 and = 0.01; the MSE is 8.57 and 1.48 deg2 for pitch and roll, respectively.
The influence of the presented 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” denote the significant difference between MSE of Mode 2 and Mode 1, Mode 3, Mode 4, and Mode 5 (), respectively. Figures 11(a), 11(b), 11(c), 11(d), 11(e), 11(f), 11(g), 11(h), 11(i), and 11(j) show an example of the attitude estimation. We choose the fastest movement of the MISUMI table slider (2.5 m/s2) in this test as a figure of attitude example. As can be seen, Mode 2 (as compared to Modes 1, 3, and 5) has the smallest amplitude of MSE, that is, 13.47 deg2.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
3.3.2. 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 value for Mode 2, we conduct the variability of the 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 . It was found that = 0.05 has an optimum performance among all the possible choices. Figures 12(a), 12(b), and 12(c) show the influence of variability to the square output of the proposed compensation model , as in the proposed model of Mode 2 in (7) and (21). Figures 13(a), 13(b), 13(c), 13(d), 13(e), 13(f), 14(a), 14(b), 14(c), 14(d), 15(a), 15(b), 15(c), 15(d), 15(e), 15(f), 15(g), 15(h), 15(i), 15(j), 15(k), and 15(l) are produced by using this parameter value to present the analysis and the result of attitude estimation.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
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.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 dash-dot line); Phase I is for 2.4 s during 5.8–8.2 s and followed by Phase II, 2 s during 12.8–14.8 s. Observing accelerometer signal along the -axis in Figure 13(b), we use Phase III and Phase IV (indicated by the dash-dot line) to refer to the external acceleration experienced for 2 s during 9.2–11.2 s and 1.9 s during 16.2–18.1 s, respectively. The external acceleration during Phases I and II is closely related to the pitch angle while Phases III and IV are related to the roll angle.
Mode 2 works based on an automatically adjusted measurement noise covariance matrix, as (21). The square of compensation model output was needed to adjust the measurement error covariance matrix and the squared signal output shown in Figures 14(a), 14(b), and 14(c). The timing of the occurrence of the enlarged amplitude of the square of the signal is aligned with the timing of the external acceleration in Phases I to IV. This means that in (21) will increase during the dynamic condition of this phase. However, Figure 14(d) presents the necessary condition of Mode 3 as a comparison.
Parameter values are used in each mode in Test B. Mode 1 uses = 1.5, while Mode 2 uses = 0.05, = 1.5, and λ = 3.5. Mode 3 uses thresholds δ = 0.1 g2, = 0.4, = 0.3, and = 0.093, = 0.034, and = 1.5. Mode 4 uses the thresholds = 0.1 and = 0.3. We did experiments with some value of , that is, 0.2, 0.3, and 0.4; = 0.3 produced a minimum MSE. Mode 5 works using = 1.5 and the factor of softened part = 3. One of the results of Test B using five modes is presented in Figures 15(b), 15(c), 15(d), 15(e), and 15(f) for pitch and Figures 15(h), 15(i), 15(j), 15(k), and 15(l) for roll estimation.
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) and 16(b) present the column bar of Table 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.
(a)
(b)
3.3.3. Test C
Test C was executed by walking straight forward three strides along the -axis of the measurement room coordinate frame. Figures 18(a), 18(b), 18(c), 18(d), 18(e), 18(f), 18(g), 18(h), 18(i), 18(j), 18(k), and 18(l) show the result of attitude estimation using five modes and the MSE result is presented in Table 7.
Some parameters values are used in Test C for each mode. Mode 1 uses = 1.5, while Mode 2 uses = 20, = 1.5, and λ = . Mode 3 uses thresholds δ = 0.2 g2, = 0.8, = 0.3, and = 0.093, = 0.034, and = 1.5. Mode 4 uses the thresholds = 0.8 and = 60. Mode 5 works using = 1.5 and the factor of softened part = 50.
4. Discussion
4.1. Mode 1
External acceleration compensation does not play a role in Mode 1. Mode 1 used 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 3 shows that the MSE of pitch and roll estimation is 25.25 and 3.65 deg2, 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.
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.
4.2. 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 role in improving the attitude estimation accuracy during fast movement. If was imposed for all experiments, then the only setting parameters are and λ.
The two parameters and λ are chosen heuristically. However, there is some consideration when choosing these parameters. As shown in (7), a smaller reduces the cutoff frequency. The value of 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 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 has an increasing amplitude of , 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 (). 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 deg2, which is the lowest among all those tested. The significant difference with Mode 1 is generally started at λ = 150 and with up to 0.1.
For Test A, the variability of influences the estimation error; however, the attitude reference is always zero degrees, and therefore we cannot observe the effect of the 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 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 degree2 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), 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 deg2 and 245.8 deg2 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 Figure 9 and Table 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 δ, , and 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, 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.
4.4. 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.
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% 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 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.
4.6. 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 . 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 value is important. One prospective improvement of this model is the additional step of calculating the optimum value of from IMU data before the Kalman filtering process.
The third limitation is considering the application of Mode 2 whenever the pitch angle (θ) reached the 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 , 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 first-order 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.
5. 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 two-step 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 ( 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.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.