Research Article  Open Access
Attitude Estimation Using Kalman Filtering: External Acceleration Compensation Considerations
Abstract
Attitude estimation is often inaccurate during highly dynamic motion due to the external acceleration. This paper proposes extended Kalman filterbased 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, oneaxis and multiple axes sensor movement. Five approaches were used to test the estimation of the attitude: (1) the KFbased model without compensating for external acceleration, (2) the proposed KFbased model which employs the external acceleration compensation model, (3) the twostep KF using weightedbased switching approach, (4) the KFbased model which uses the thresholdbased approach, and (5) the KFbased model which uses the thresholdbased approach combined with a softened part approach. The proposed algorithm showed high effectiveness during the oneaxis 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, sittostand [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 randomwalk errors. In static or slow movement, the accelerometer measures roll and pitch by leveling to correct the gyrounbounded 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 weightedswitching approach. In contrast, Sabatini in [11] used a thresholdbased approach, while Harada et al. [12] proposed a thresholdbased 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.
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 wellknown 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 threeaxis 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 threeaxis 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 steadystate 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 lefthalf plane. There are two possibilities for poles locations on the lefthalf side, such as complexconjugate 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 lefthalf 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 steadystate response for and the steadystate error as inThe time response for ramp input () is . The steadystate output and error are calculated as follows:A rectangular pulse of amplitude and of duration , , has a zero steadystate output and also a zero steadystate error.
For a parabolic function input, , the time response is calculated as follows, with an infinity steadystate output and infinity steadystate 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 steadystate 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 twostep EKF which adjusts the measurement noise covariance using the weightedswitching approach. Mode 4 is the quaternionbased EKF, which uses a thresholdbased 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 thresholdbased 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 rad^{2}. And for the gyroscope along the , , and axes, the values are , , and in rad/s^{2}. 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} g^{2}. 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 twostep EKF [1] used separate measurements between the gyro and accelerometer. The external acceleration compensation technique is based on weightedswitching 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 twostep 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 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 weightedswitching 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 quaternionbased 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 thresholdbased 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 quaternionbased UKF. The method is thresholdbased 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 singleaxis 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 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 shoetype measurement device. An IMU sensor was placed on the top of one shoe. The subject was asked to walk straight, 34 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 singleaxis 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 RSManager 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 RS232 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 crosscorrelation 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 degree^{2} 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 g^{2}, = 0.4, and = 80. The experimental results of Mode 1 and Mode 3 are presented in Tables 1(a) and 1(b).
(a)  
 
(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 deg^{2} 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 oneaxis movement test is as follows. The MSE trend decreases as the value of λ increases or decreases. The best combination for oneaxis movement in this experiment is λ = 200 and = 0.01; the MSE is 8.57 and 1.48 deg^{2} 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/s^{2}) 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 deg^{2}.
(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 dashdot 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 dashdot 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 g^{2}, = 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 g^{2}, = 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 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 oneaxis 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 deg^{2}, 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 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), 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 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 thresholdbased 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 firstorder approximation error. Mode 2 uses EKFbased filtering that employs linearizing the nonlinear model. The firstorder 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 firstorder 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 filterbased 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 freehand 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 weightedswitching method (Mode 3); the fourth is a quaternionbased EKF using a thresholdbased method (Mode 4); and the fifth (Mode 5) uses an unscented Kalman filter and is thresholdbased 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 firstorder 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 ( 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.
References
 Y.S. Suh, S.K. Park, H.J. Kang, and Y.S. Ro, “Attitude estimation adaptively compensating external acceleration,” JSME International Journal, Series C: Mechanical Systems, Machine Elements and Manufacturing, vol. 49, no. 1, pp. 172–179, 2006. View at: Publisher Site  Google Scholar
 P. D. Groves, Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, Artech House (Horizon House), London, UK, 2nd edition, 2013.
 C. Wada, D. Takigawa, F. Wada, K. Hachisuka, T. Ienaga, and Y. Kimuro, “Improvement research of shoetype measurement device for a walking rehabilitation support system,” in CrossCultural Design. Cultural Differences in Everyday Life: 5th International Conference, CCD 2013, Held as Part of HCI International 2013, Las Vegas, NV, USA, July 21–26, 2013, Proceedings, Part II, P. L. P. Rau, Ed., vol. 8024 of Lecture Notes in Computer Science, pp. 157–164, Springer, Berlin, Germany, 2013. View at: Publisher Site  Google Scholar
 B. Mariani, C. Hoskovec, S. Rochat, C. Büla, J. Penders, and K. Aminian, “3D gait assessment in young and elderly subjects using footworn inertial sensors,” Journal of Biomechanics, vol. 43, no. 15, pp. 2999–3006, 2010. View at: Publisher Site  Google Scholar
 F. Dadashi, B. Mariani, S. Rochat, C. J. Büla, B. SantosEggimann, and K. Aminian, “Gait and foot clearance parameters obtained using shoeworn inertial sensors in a largepopulation sample of older adults,” Sensors, vol. 14, no. 1, pp. 443–457, 2014. View at: Publisher Site  Google Scholar
 D. M. Karantonis, M. R. Narayanan, M. Mathie, N. H. Lovell, and B. G. Celler, “Implementation of a realtime human movement classifier using a triaxial accelerometer for ambulatory monitoring,” IEEE Transactions on Information Technology in Biomedicine, vol. 10, no. 1, pp. 156–167, 2006. View at: Publisher Site  Google Scholar
 J. Klenk, C. Becker, F. Lieken et al., “Comparison of acceleration signals of simulated and realworld backward falls,” Medical Engineering & Physics, vol. 33, no. 3, pp. 368–373, 2011. View at: Publisher Site  Google Scholar
 C. Wada, Y. Tang, and T. Arima, “Development of a standingup motion guidance system using an inertial sensor,” in Advances in Intelligent Systems and Computing: Advanced Intelligent Systems, Y. S. Kim, Y. J. Ryoo, M. S. Jang, and Y. C. Bae, Eds., vol. 268, pp. 179–187, Springer, Berlin, Germany, 2014. View at: Google Scholar
 J. Musić, R. Kamnik, and M. Munih, “Model based inertial sensing of human body motion kinematics in sittostand movement,” Simulation Modelling Practice and Theory, vol. 16, no. 8, pp. 933–944, 2008. View at: Publisher Site  Google Scholar
 S. Romaniuk and Z. Gosiewski, “Kalman filter realization for orientation and position estimation on dedicated processor,” Acta Mechanica et Automatica, vol. 8, no. 2, pp. 88–94, 2014. View at: Publisher Site  Google Scholar
 A. M. Sabatini, “Quaternionbased extended Kalman filter for determining orientation by inertial and magnetic sensing,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 7, pp. 1346–1356, 2006. View at: Publisher Site  Google Scholar
 T. Harada, T. Mori, and T. Sato, “Development of a tiny orientation estimation device to operate under motion and magnetic disturbance,” The International Journal of Robotics Research, vol. 26, no. 6, pp. 547–559, 2007. View at: Publisher Site  Google Scholar
 R. B. Widodo, H. Edayoshi, and C. Wada, “Complementary filter for orientation estimation: adaptive gain based on dynamic acceleration and its change,” in Proceedings of the 15th International Symposium on Soft Computing and Intelligent Systems (SCIS '14), Joint 7th International Conference on and Advanced Intelligent Systems (ISIS '14), pp. 906–909, Kitakyushu, Japan, December 2014. View at: Publisher Site  Google Scholar
 J. Calusdian, X. Yun, and E. R. Bachmann, “Adaptivegain complementary filter of inertial and magnetic data for orientation estimation,” in Proceedings of the IEEE International Conference on Robotics and Automation (ICRA '11), pp. 1916–1922, Shanghai, China, May 2011. View at: Publisher Site  Google Scholar
 J. J. D'Azzo, C. H. Houpis, and S. N. Sheldon, Linear Control System Analysis and Design with MATLAB, Marcel Dekker, 5th edition, 2003.
 T. Ozyagcilar, “Implementing a tiltcompensated ecompass using accelerometer and magnetometer sensors,” Application Note AN4248, Freescale Semiconductor, 2012. View at: Google Scholar
 M. S. Grewal and A. P. Andrews, Kalman Filtering Theory and Practice Using MATLAB, John Wiley & Sons, Hoboken, NJ, USA, 3rd edition, 2008.
 J. C. F. de Winter, “Using the student's ttest with extremely small sample sizes,” Practical Assessment, Research and Evaluation, vol. 18, no. 10, pp. 1–12, 2013. View at: Google Scholar
 N. H. Q. Phuong, H.J. Kang, Y.S. Suh, and Y.S. Ro, “A DCM based orientation estimation algorithm with an inertial measurement unit and a magnetic compass,” Journal of Universal Computer Science, vol. 15, no. 4, pp. 859–876, 2009. View at: Google Scholar
 J. J. LaViola Jr., “A comparison of unscented and extended kalman filtering for estimating quaternion motion,” in Proceedings of the American Control Conference, vol. 3, pp. 2435–2440, June 2003. View at: Google Scholar
Copyright
Copyright © 2016 Romy Budhi Widodo and Chikamune Wada. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.