Estimating Relative Angles Using Two Inertial Measurement Units Without Magnetometers

Inertial measurement units (IMUs) are used in biomechanical and clinical applications for quantifying joint kinematics. This study aimed to assist researchers new to IMUs and wanting to develop an inexpensive IMU system to estimate the relative angle between IMUs, while understanding the different algorithms for estimating angular kinematics. Thus, there were three subgoals: 1) to present a low-cost and convenient IMU system utilizing two 6-axis IMUs for computing the relative angle between the IMUs; 2) to examine seven methods for estimating the angular kinematics of an IMU; and 3) to provide an open-source code and working principles of these methods. The raw gyroscopic and accelerometer data were preprocessed. The seven methods included gyroscopic integration (GI), accelerometer inclination (AC), basic complementary filter (BCF), Kalman filter (KF), digital motion processor (DMP, a proprietary algorithm), Madgwick filter (MW), and Mahony filter (MH). An apparatus was designed to test nine conditions that computed angles for rotation about three axes (roll, pitch, yaw) and three movement speeds (50°/s, 150°/s, 300°/s). Each trial lasted 25 min. The root-mean-squared error (RMSE) between the gold-standard value measured from the apparatus’ encoder and the value calculated from each of the seven methods was determined. For roll and pitch, all methods accurately quantified angles (RMSE < 6°) at all speeds. For yaw, all methods except AC and DMP displayed RMSE < 6° at all speeds. AC could not be used for yaw angle computation, and DMP displayed RMSE >6°. Researchers can utilize appropriate methods based on their study’s application.


I. INTRODUCTION
R ESEARCHERS have utilized compact, accurate, and affordable inertial measurement unit (IMU) sensors and developed algorithms to compute angular kinematics of one or more body segments in biomechanics and clinical applications [1], [2], [3], [4], [5], [6]. An IMU provides a 3-D vector representing the spatial orientation of the IMU. Thus, two IMUs can be used to measure the relative joint angle between two adjacent body segments by placing the IMUs on the body segments. This approach can be used for computing joint kinematics in biomechanical (e.g., knee flexion/extension angle during slow gait) [1], [2], [4], [7] and clinical applications (e.g., quantifying elbow joint kinematics during a clinical assessment of abnormal muscle behavior) [5], [6]. Previous studies have investigated the accuracy of highperformance IMUs under dynamic conditions encountered in human subject testing or test apparatuses. Seel et al. developed a method for quantifying human knee and ankle angles using 6-axis (i.e., 3-axis accelerometers and gyroscopes) IMUs while exploiting kinematic constraints. They validated the accuracy of IMUs via human subject testing using a motion capture system [1]. Favre et al. [8] introduced a method for measuring 3-D knee joint motion by aligning the reference frames of the two 6-axis IMUs of the thigh and shank segments. Lebel et al. validated multiple commercially available 9-axis IMUs under various testing conditions (e.g., various speeds and rotation axes) using an instrumented gimbal table [2], [9]. Ricci et al. [10] quantified the accuracy of commercial 9-axis IMUs under various speeds and rotation axes using a robotic arm. Bergamini et al. [11] proposed an adaptive Kalman filter (KF) method with a 6-axis IMU for quantifying the trunk angles of sprint runners. Caruso et al. [12], [13] utilized a magneto-IMU and presented a method for suboptimally tuning the input parameters for the IMU filter designs. However, these studies had some limitations such as the 1) dependence on magnetometers; 2) use of costly IMUs; 3) incomprehensive study design; and 4) lack of an open-source code for computing angles from IMU data. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ First, some studies relied on IMUs with magnetometers (i.e., 9-axis IMUs) for computing the orientation of the IMUs [2], [9], [10], [12], [13]. While the magnetometers provide a measurement of absolute heading angle, the magnetometers are vulnerable to magnetic interferences caused by ferrous objects and electrical appliances (e.g., space heaters, metal filing cabinets), making magnetometers difficult to use in indoor environments [14], [15]. Although calibration procedures exist to compensate for the magnetic interferences [16], the procedures must be conducted periodically, making the device less practical in an indoor setting. Other studies have observed or noted the effect of magnetic disturbances in an indoor lab environment [9], [10], [17], [18]. However, a 6-axis IMU (i.e., 3-axis accelerometers and gyroscopes) does not use magnetometers to compute the orientation of the IMUs. While a 6-axis IMU cannot compute absolute yaw angle like the 9-axis IMU, the 6-axis IMU can still compute relative yaw angle which is still valuable in computing angular kinematics in biomechanics or clinical application [1], [2], [3], [4], [5], [6]. Thus, we utilized 6-axis IMUs in this study, since they are more robust, simpler, and equally effective choices than a 9-axis IMU for our application (i.e., indoor applications).
Second, some studies did not fully control the dynamic test conditions [1], [8]. Movement speeds, range of motion, rotation axes, and test duration are critical factors that can affect the accuracy and behavior of the IMUs [12], [13], [19], [20]. However, some studies only validated their methods using an arbitrarily defined movement speed and range of motion, or a very short test duration (<2 min) [9], [10], [11], [12].
Lastly, most of the studies did not provide an open-source code for researchers in other fields to easily implement these newly proposed methods [1], [8], [10], [11]. Researchers unfamiliar with control/estimation theory may have difficulty replicating and utilizing the newly proposed methods.
Thus, the goals of this study were to: 1) provide a low-cost and convenient IMU system utilizing two 6-axis IMUs without magnetometers for computing a relative angle between a stationary and a moving IMU; 2) validate various methods of computing 3-D absolute angles (roll, pitch, yaw) of an IMU following strict test conditions (i.e., long test duration, three movement speeds, three rotation axes); and 3) provide an open-source code and working principles of these methods to assist researchers unfamiliar with IMUs. The open-source code and data are at: https://github.com/ssong47/compute_relative_angle_between_ two_IMUs, https://ieee-dataport.org/open-access/estimatingrelative-angle-between-two-6-axis-inertial-measurement-unitsimus

II. METHODS A. Hardware and Test Apparatus
We developed an apparatus that performed programmable and repeatable rotations of IMUs and a gold standard measurement of the rotation angles (via an encoder) to validate the accuracy of our IMU system (see Fig. 1). Unlike other studies that validated their system through human subject testing by using optoelectronic motion capture data as the gold standard [1], [21], [22], our IMUs were mounted on this apparatus for testing to remove the errors due to motion artifacts from the attachment on a body segment. The apparatus consisted of two IMUs (MPU-6050; TDK-InvenSense, San Jose, CA), a stepper motor (NEMA 17, 42BYGH3401; Han Ding Motor, China) with an optical rotary incremental encoder (HEDS-9040#T00; Broadcom, USA) attached at the shaft, a stepper motor driver (A4988; Allegro Microsystems, Manchester, NH), and two microcontrollers (Teensy 3.6; PJRC, Sherwood, OR, and Arduino Uno; Arduino LLC, Italy). The MPU-6050 was chosen due to its compact size (smaller than a U.S. quarter), low power consumption, and low cost (∼$US 4). An onboard proprietary digital motion processor (DMP 1 ) algorithm was available for computing the 3-D orientation of each IMU, represented in quaternion form. Other system components are also economical and popular off-the-shelf choices (Arduino Uno ∼$US 27, Teensy 3.6 ∼$US 27, and NEMA 17 stepper motor ∼$US 12).
The test apparatus used the Teensy 3.6 microcontroller to ensure a consistent sampling rate (100 Hz) for data collection and reliable transmission of data [i.e., timestamp, encoder angle, two IMUs' output data (quaternion from DMP, 3-D linear accelerations, and 3-D angular velocities in x, y, z axes)] to a PC. The I2C protocol ensured synchronous sampling of the data from the two IMUs since the sampling of each IMU was done simultaneously by using a shared clock signal. The Arduino Uno sent the PWM signals stepper motor driver to control the motor speed and range.
The encoder measurement served as the ground truth to evaluate the performance of the IMU system. The two IMUs, fixed on their custom mounts, were placed on different parts of the apparatus: IMU 1 on the stationary motor mount and IMU 2 on the rotating motor hub. The custom 3-D printed parts were designed to allow relative rotation about three orthogonal axes (i.e., yaw, pitch, roll) between the two IMUs (see Fig. 2). Depending on the axis of rotation of interest, the IMU mounts were configured differently (see Fig. 2). Note that the yaw rotation was configured such that the rotation axis was parallel to the gravity vector [see Fig. 2(c)] to simulate a case in which gravity could no longer be a reference vector for the angle computation.

B. Testing Protocol
A total of nine test trials (three rotation axes × three movement speeds) with each trial length of 25 minutes and a range of motion of 180 • was performed. 25 min was chosen as the test length since previous studies involving clinical research on neuro-rehabilitation lasted approximately 25 min for each test subject [3], [23], [24], [25]. The range of motion was set to 180 • since most anatomical joints do not exceed 180 • . Three speeds were chosen as the following: slow (max speed, ω max = 50 • /s), medium (ω max = 150 • /s), and fast (ω max = 300 • /s) to simulate the kinematics of an elbow or wrist joint during a clinical assessment (20 • /s-180 • /s) [6], [26], [27], [28] as well as of a knee joint during the walking phases at a slow cadence (23 • /s-250 • /s) [29], [30], [31]. Similar movement speeds and rotation axes were used in other studies [9], [10]. The motor driver controlled the angular position and speed profile shown in Fig. 2(d).
Before each trial, the two IMUs were removed from the test apparatus and calibrated to remove any unwanted bias and scaling factor of the IMU readings. The calibration protocol followed the common six-position static test mentioned in [32], [33]. At the beginning of each trial, the principal axes of IMU 1 and 2 were aligned parallel to each other. Also, the encoder angle was set to zero. During each trial, the two IMU and encoder readings were sampled at 100 Hz.

C. Methods for Computing Relative Angles
Seven methods of computing relative angles about a single rotation axis were investigated. These included the DMP method, gyroscopic integration (GI), accelerometer inclination (AC), basic complementary filter (BCF), KF, Madgwick filter (MW), and Mahony filter (MH). The mathematical formulas, as well as the pros and cons of each method, are discussed here.
The testbed setup had one stationary IMU (IMU 1) and one rotating IMU (IMU 2), where two IMUs' local frames were always aligned by construction in all test configurations (see Fig. 2). The relative angle (θ m,û ) was defined as how much the moving IMU (IMU 2) rotated relative to stationary IMU 1 around the testing axis, where m ∈ {GI, AC, BCF, KF, MW, MH, DMP}, indicating the algorithm being used for calculation. For pitch and roll [see Fig. 2(a) and (b)], θ m,û was obtained by calculating the inclination angle of each IMU relative to the horizontal plane and then subtracting the two inclination angles (1). For yaw [see Fig. 2(c)], since the motion is in the horizontal plane, instead of inclination angle, the change in the angle with respect to the initial angle of each IMU was first calculated and then subtracted to obtain θ m,û where the superscript denoted the IMU of interest andû denoted the rotation axis (û ∈ {î(roll),ĵ(pitch),k(yaw)}).
In the case of pitch and roll rotations, to obtain the inclination angle of each IMU, when using quaternion-based methods (i.e., DMP, MW, and MH), θ 1 m,û and θ 2 m,û were obtained by converting the quaternion of each IMU into intrinsic Euler angles ("ZYX" order) and then selecting the Euler angle of interest (roll, pitch, or yaw angles), whereas when using axisangle-based methods (i.e., GI, AC, BCF, and KF), θ 1 m,û and θ 2 m,û were obtained directly around the rotation axis being tested.
Each IMU outputted 3-D gyroscope values, accelerometer values, and a unit quaternion vector (q) from DMP. q contained four values (a, b, c, d, where a defined the amount of rotation and b, c, d defined the axis of rotation in the 3-D Cartesian space) to quantify rotation relative to the reference frame of the IMU (î,ĵ,k) (2) [34] This article follows the vector notation recommended by Arkfen and Weber [35]: bold serif for vectors (e.g., q), bold lowercase letter with circumflex for unit vectors of a reference frame (e.g.,î), and left superscript and subscript indicating the reference frame { A} and inertial frame { A }, respectively (e.g., A A q). The frames of IMU 1 and 2 are defined Each quaternion can be converted into Euler angles using MATLAB scripting (3) [34], [36].

1) Digital Motion Processor (DMP):
The DMP method, developed commercially by TDK-InvenSense, is a proprietary algorithm that fuses the 3-D acceleration and 3-D gyroscope data and outputs the 3-D orientation of an IMU in quaternion form [37]. While the computation of a quaternion from an IMU is proprietary, the computation of relative angles of the two IMUs is not. These quaternions ( 1 1 q or 2 2 q) were converted to Euler angles (θ i z , θ i y , θ i x ) using (3). The angle of interest from the Euler angles was used to compute the absolute angle of IMU 1 or IMU 2, and (1) was used to compute the relative angle between IMU 1 and 2.

Algorithm 1 Digital Motion Processor TM (DMP)
2) Gyroscopic Integration (GI): GI integrates gyroscopic data at each sampled time [38]. Since integrating raw gyroscopic data would cause drift due to inherent noise and bias instability [38], the raw gyroscopic signals were preprocessed using a fourth-order Butterworth high-pass filter with a cutoff frequency of 0.07 Hz. A similar approach was done by other studies [39], [40], [41]. While the most optimal cut-off frequency may be unique for different movement speeds, we wanted to investigate if a single filter design was robust enough to accommodate a range of movement speeds. The cutoff frequency was chosen low enough to ensure that signals related to the motion at all speeds were collected but high enough to reduce the effect of the long-term effect of drifting of gyroscopes [42], [43]. The lowest frequency of motion was 0.139 Hz (=(1/360 • /(50 • /s))), which was determined by calculating the time it takes to complete one cycle of motion (θ c = 360 • = 2 × 180 • , since the motor rotated 180 • in one direction and came back) at the slowest rotational speed (ω = 50 • /s). Thus, a cut-off frequency of 0.07 Hz (<0.139 Hz) ensured that signals related to motions at all speeds were collected with minimal unwanted gyroscopic drift. tû , and t were the angles from GI at the tth data point about the rotation axisû, the angle from GI at t − 1th data point aboutû, the high-pass-filtered angular velocity gyro at the tth data point data aboutû, and sampling time (=0.01 s), respectively. The advantages of GI are easy implementation and low computational cost [44], [45].
3) Accelerometer Inclination (AC): AC solely uses the accelerometer data and computed the inclination angle using the gravity vector and its projection on the axes of the accelerometer [45]. The accelerometer data were preprocessed using a low-pass fourth-order Butterworth filter with a cut-off frequency of 4 Hz to fully capture the signals with a frequency (4 Hz > 0.833 Hz = (1/2(180 • )/(300 • /s))) of the motion even at the fastest speeds [46]. A single cut-off frequency was chosen for similar reasons as the high-pass filter mentioned above. The yaw angle cannot be obtained for this method since the yaw has a rotation axis parallel to the gravity vector g, making the projections of g on the accelerometer to be zero.
, and A i t,k were angles from AC aboutî (roll),ĵ (pitch) axes and low-pass-filtered acceleration readings aboutî,ĵ,k axes, respectively. The advantages of the AC method were drift-free computation of angle, simple implementation, and low computational Algorithm 3 Accelerometer Inclination (AC) [45] cost [45]. The disadvantages of AC were the high sensitivity to motion/dynamic accelerations (i.e., acceleration not due to gravity) and the inability to compute yaw angles [42], [47]. The angle from AC contained high-frequency noise, which could be mitigated using low-pass-filtered accelerometer readings [42], [43], [45].

4) Basic Complementary Filter (BCF):
The BCF method combines two filters: a low-pass filter for accelerometer signals and a high-pass filter for gyroscopic signals [42], [48], [49]. The two filters are complements of each other, so the summation of the two filtered signals results in a gain of one. The advantages of BCF were smooth, accurate, and drift-free angle computation that combined the benefits of GI and AC while maintaining low computational costs. Also, the tuning process was simple due to only one tunable parameter (γ ) that determined the cut-off frequency of the low-and high-pass filters [50]. γ was determined by the sampling rate (100 Hz) and time constant (τ ) of the system. τ was determined by the gyroscope's drift rate (ω drift ) and the tolerable error (e tol ). For our test setup and data, the measured ω drift of MPU 6050 was 0.25 • /s. For e tol , 3 • was chosen to ensure the tolerable error was at least half of the threshold for allowable error in biomechanics studies (6 • ) [1], [2], [21], [51]. Thus, τ was set to 0.083 s −1 (=0.25/3), which meant that BCF would perform GI for signals with frequencies larger than 0.083 Hz, while the accelerometer readings would be filtered out. For signals with frequencies smaller than 0.083 Hz, BCF would rely on AC angles, while the gyroscopic readings would be ignored. Thus, γ = 0.11 ∼ = 1 − (0.083/0.083 + 0.01) was chosen. For pitch and roll angles, the raw gyroscopic signals were not filtered for BCF, since the BCF algorithm already corrects for the gyroscopic drift using accelerometer values. Having extra high-pass filters is redundant and may increase unwanted filter artifacts (e.g., delays) that have an adverse effect on real-time applications. For yaw angles, a small value of γ (= 10 −7 ) was used to put less weight in the accelerometer readings since the accelerometer did not provide any relevant information regarding the yaw rotation. Also, the gyroscopic signals were high-pass filtered through the same filter settings used in GI, since the accelerometer readings could not properly correct for the gyroscopic drift. [42], [48], [49] τ =

Algorithm 4 Basic Complementary Filter (BCF)
CFû was the BCF angle about theû axis. 5) Kalman Filter (KF): KF computes angle via an optimal state estimator algorithm, consisting of prediction and update steps [16], [42], [52]. The prediction step involves estimation of the a priori (i.e., at current time step without any observation information from current time step) system states ( and θ i t,ACx .

III. Update
Step gyro was the bias of the gyroscope. Note that the initial error covariance ( P 0 ) was set to a relatively high value in our study to ensure fast convergence [53].
Q was a diagonal matrix since we assumed that the variance of accelerometer noise and gyroscope noise was independent. Q 22 was equivalent to the variance (σ 2 ω ) of raw gyroscopic noise. Assuming all zero mean gyroscopic noise (p ω = 0), σ 2 ω = p 2 ω,RMS , where p ω,RMS was the rms noise of our gyroscope. p ω,RMS = p ω,nspd √ t −1 , where t −1 was our sampling rate (=100 Hz), and p ω,nspd was the noise spectral density (=0.005 • /s/ √ Hz) of MPU6050 [37]. p ω,RMS was computed using t −1 and p ω,nspd rather than using the reported total rms noise value, since the gyroscopic noise level depends not only on the performance of the sensor, but also on the sampling rate [54]. Thus, Q 22 = 2.5 × 10 −3 . Like other studies [42], Q 11 was set to 1.0 × 10 −3 to maintain similar order of magnitude as Q 22 . R(=I 2 × 3.76 • ) was determined by empirically measuring the variance of the AC angle's noise. For pitch and roll, the raw gyroscopic signals were not filtered for KF, since the KF algorithm already tries to remove the gyroscopic bias through estimation of the bias within the state space model. For yaw rotations, R was set to a matrix with large numbers (=I 2 × 10 8 ) since the accelerometer measurements were not reliable for yaw. Like BCF, the gyroscopic signals were high-pass filtered through the same filter used in GI, since the accelerometer readings could not properly correct for the gyroscopic drift. The benefit of KF is that it optimally computes the gain that balances the inaccuracy of the measurements and the model at every time step.

6) Madgwick Filter (MW):
The MW is a sensor fusion algorithm that combines the gyroscope and accelerometer data in a quaternion representation to analytically derive and optimize a gradient-descent (GD) algorithm to compute a quaternion derivative of the IMU [47], [55]. First, the MW estimates the IMU (i ) quaternion derivative ( i i qω,t ) using the vector ( i ω t ) containing gyroscopic measurement ( i ω x,t , i ω y,t , or i ω z,t ) and previous estimate of orientation ( i i q est,t −1 ). Only one gyroscopic signal of i ω t was nonzero, since the IMU is rotated about a single axis. To minimize the drift from this gyroscopic measurement, a corrective quaternion derivative term ( i i q,t ) is computed by estimating the direction of the gravitational vector i g est using the GD-based algorithm on the accelerometer data ( i i a t ). Note that we assumed that the accelerometer mostly measures the direction of i g est . In the GD, objective function f is minimized by using the Jacobian matrix ( J) of f which computed the error direction ((∇ f / |∇ f |)) on the solution surface. A corrected quaternion derivative term ( i i qest,t ) is computed by fusing the quaternion derivative estimations from the gyroscope ( i i qω,t ) and the accelerometer ( i i q,t ) using an adjustable gain β. β represents all mean zero gyroscope measurement errors (ω β ), expressed as the quaternion derivative form.ω β (= 0.05 • /s) is dependent on the sensor's gyroscopic performance. The value of β = 0.043 was chosen by using the equation (β = √ 3/4ω β ) suggested by Madgwick et al. [47]. For yaw rotations, β = 0 since the accelerometer readings did not describe the movement of the IMU properly when the rotation axis is parallel to g. The IMU's orientation estimation ( i i q est,t ) was computed by integrating the corrected quaternion derivative ( i i qest,t ). Similar procedure from DMP involving conversion from quaternion to Euler angles was used in MW to compute the relative angles between IMU 1 and 2.
For pitch and roll, the raw gyroscopic signals ( i ω x,t , i ω y,t ) were used without a high-pass filter since the gyroscopic drifting behavior was corrected by the accelerometer data. An additional high-pass filter was redundant and added delays. However, for yaw, high-pass-filtered gyroscopic signals ( i ω hf,z,t ) were used, since MW cannot rely on the accelerometer data to compensate for gyroscopic bias.
The advantages of MW are its low computational load, low sampling rate requirement, and simple tuning method (only one adjustable parameter-β). Since MW uses a quaternion representation, the filter is not subject to gimbal lock problems associated with Euler angle representations. Note that a potential issue with MW is the unpredictable convergence behavior and a more recent version of this algorithm with enhanced robustness was reported in [56]. [47], [55] i ω t = i ω x,t 0 0 for roll

7) Mahony Filter (MH):
The MH is another variation of a BCF that corrects the gyroscopic measurement by a proportional-integral (PI) compensator [55], [57]. First, MH computes the error ( i e t ) between the estimated gravity vector ( i a est,t ) and accelerometer readings ( i a t ) represented by the cross product between i a t and i a est,t . Then, the gyroscopic vector ( i ω q t ) containing gyroscopic measurements ( i ω x,t , i ω y,t , or i ω z,t ) is combined with i e t adjusted by PI gains (K p , K i ) to obtain a corrected estimate of the gyroscopic vector ( i ω q est,t ). Afterward, the quaternion derivative ( i i qest,t ) is obtained from the previous estimate of quaternion ( i i q est,t −1 ) and i ω q est,t . Like MW, only one gyroscopic signal of i ω t is nonzero. The IMU's orientation quaternion ( i i q est,t ) is obtained by integrating the quaternion derivative ( i i qest,t ). A similar procedure as DMP, that is, conversion from quaternion to Euler angles, was used in MH to compute the relative angles between IMUs 1 and 2. The filter gains (K p , K i ) were initially set to their default values (1.0, 0.3) and then fine-tuned heuristically (5.0, 0.3) for pitch and roll rotations. For yaw rotations, K p = K i = 0 due to similar reasons as MW. Like MW, raw gyroscopic signals were used for pitch and roll, while high-pass-filtered gyroscopic signals were used for yaw. The MH has been reported to be an efficient and effective solution by other studies [55], [58].

D. Data Processing
The performance of each method was evaluated using root-mean-squared error (RMSE) between the encoder and the seven computed methods for evaluating the accuracy of Algorithm 7 Mahony Filter (MH) [55], [57] i a est, the computed angle's magnitude. For RMSE computation, a window size of one minute, or 6000 sample points (=60 s × 100 points/s) were chosen for two reasons: 1) to better visualize the trends of RMSE changes and 2) since each test trial from other research studies lasted approximately 1 min [11], [59], [60]. RMSE of 6 • was chosen as the maximum acceptable limit for biomechanical research studies [1], [2], [21], [51]. All RMSE data were truncated such that the start (0-1 min) and end (24-25 min) data were removed. The artifacts were due to the preprocessing of gyroscopic and accelerometer signals. For each speed, the average (RMSE) and standard error (σ RMSE ) RMSE were computed for the truncated data. In addition to the RMSE data, the average computational time (t comp ) required to execute a single iteration for each algorithm was calculated for all seven methods. All data were processed using MATLAB R2020a.

III. RESULT A. Roll and Pitch
All methods exhibited high accuracy (RMSE < 6 • ) at all movement speeds and no signs of drift across the entire duration (see Figs. [3][4][5]. Also, the RMSE of these methods was generally proportional to speed (see Fig. 4). While the RMSE differences between these methods were less than 3 • , the methods displayed differences in terms of consistency of RMSE trends across time and movement speeds.
The four sensor fusion algorithms (BCF, KF, MW, MH) displayed the lowest RMSE across all speeds for all rotation axes (see Figs. 4 and 5). In addition, the RMSE remained consistent across time with only minor fluctuations (∼0.5 • ). The RMSE was proportional to movement speeds, since the RMSE was the highest (3 • -5 • ) at fast speed and the lowest (1 • -3 • ) at slow speed.
AC, GI, and DMP displayed slightly larger errors than the previously mentioned four methods. The magnitude of RMSE was higher by 1 • -3 • for all speeds. While GI had a relatively large error and fluctuations compared to others, GI's RMSE did not drift and remained consistently below the 6 • threshold.

B. Yaw
Five of the methods (GI, BCF, KF, MW, MH) displayed high accuracy (RMSE < 6 • ) and consistency (σ RMSE < 3 • ) without any drifting behavior for yaw rotations. Like pitch and roll, the RMSE of these methods was generally proportional to the movement speed with a few exceptions (e.g., for yaw, RMSE CF at slow >RMSE CF at medium speed). GI exhibited the highest error among the five methods, followed by KF and BCF. MW and MH showed the lowest error. Note that MW and MH had identical RMSE values (see Figs. 4 and 5).
The other two methods (AC, DMP) either could not be computed at all (AC) or displayed high errors (DMP). DMP drifted, resulting in high RMSE and high σ RMSE .

C. Computational Execution Time
The average computational execution times required for a single update step were lowest for the GI, AC, followed by BCF, KF, MH, and MW (see Fig. 6). The execution time for  BCF was marginally higher (×8) than GI and AC. For KF, the execution time was higher than GI and AC by ∼65×. The execution times for MW and MH were significantly higher than GI, and AC by a factor of almost 160×. The computational time for DMP could not be calculated since DMP was an onboard proprietary algorithm.

IV. DISCUSSION
The purpose of this study was to systematically validate seven computational methods (GI, AC, DMP, BCF, KF, MW, MH) to the gold standard (encoder) using two low-cost 6-axis IMUs without magnetometers. An explanation of the inner workings and open-source code of the seven methods were presented to assist newcomers to the field of IMUs. A custom test apparatus and strict test protocols were used to control for critical factors affecting IMU accuracies, such as movement speeds, range of motion, rotation axes, and sufficiently long test durations.
The result of this study suggests that most methods (all methods for pitch/roll. GI, BCF, KF, MW, MH for yaw) using 6-axis IMUs can be used in a controlled indoor environment to accurately quantify angles for various movement speeds (50 • /s-300 • /s) and test durations of 25 min. The RMSE for these methods were all below the maximum acceptable RMSE for biomechanical study (=6 • ). The RMSE for our methods was generally proportional to movement speed, which was also observed in other studies [9], [10], [12], [19]. Noise introduced from the larger accelerations at faster speeds may be the cause for the slightly increased error. For pitch and roll rotations, all methods did not require retuning of filter parameters for this study's movement speeds, illustrating the method's robustness to the given movement speeds.
Adjustments in the filter parameter and preprocessing of gyroscopic signals had to be made for yaw trials. Filter parameters for BCF, KF, MW, and MH had to be changed to rely only on gyroscopic reading due to the absence of a magnetometer and unusable accelerometer readings for yaw rotations. Some methods could not be used (AC) or were inaccurate due to drift (DMP). Also, high-pass filtering of raw gyroscopic signals was necessary to remove the gyroscopic bias since accelerometer signals could not correct for the drifting of gyroscopic signals.

A. Roll and Pitch
BCF, KF, MW, and MH demonstrated the best performance (i.e., lowestRMSE) and consistency (i.e., lowest σ RMSE ), compensating for the weaknesses of each sensor to accurately compute the angles (Figs. 4 and 5, top and mid). While the sensor fusion algorithm innerworkings differed greatly for each method, the accuracy and consistency of each method were similar in our study. With carefully tuned filter parameters, relatively simple sensor fusion algorithms, such as BCF and KF performed as well as more advanced algorithms, such as MW and MH. BCF, KF, MW, and MH were easy to tune and implement since they required only one or two tuning parameters. The ability to tune the filter algorithm specifically to our test setup was critical to remove noise and correct estimates of IMUs, as others have reported [12].
AC and DMP also demonstrated sufficiently accurate and consistent but slightly poorer performance than the previously mentioned methods. AC displayed no drifting behavior and slightly higher errors due to the sole dependence on accelerometers that are sensitive to external accelerations (e.g., sudden acceleration/deceleration of the IMU) irrelevant to the gravity vector [42]. It is difficult to discuss the details of the DMP since the inner workings of DMP were proprietary. The default parameter values for DMP were not tunable and thus could not be optimized for our test setup. Others also have noted the importance of adjusting parameters for the given specific experimental scenarios [19], [20]. GI using high-pass-filtered gyroscopic signals displayed no drifting behavior and high enough accuracy and consistency for all speeds. The high-pass-filtered gyroscopic signals were useful for obtaining accurate and nondrifting angles since the high-pass filter inherently removes the gyroscopic bias [dc (=0 Hz) offset of a signal]. This bias instability of the gyroscopes is one of the main causes of the drifting behavior of gyroscopes [61]. Interestingly, a single high-pass filter could be used on different movement speeds ranging from 50 • /s to 300 • /s to provide drift-free and accurate GI angles.

B. Yaw
GI, BCF, KF, MW, and MH demonstrated stable and accurate computation of angles across all speeds during yaw rotations for the entire trial albeit in the absence of magnetometers (Figs. 4 and 5,bottom). For these methods, only high-passfiltered gyroscopic signals were used to remove the gyroscopic bias since accelerometer readings could not be used to correct the gyroscopic drift. Thus, the use of a high-pass filter was necessary since the drift could not be compensated due to the absence of magnetometers, which provide the reference for the heading. The accelerometer could not read any meaningful signals related to the motion of the IMU but the gravity vector in the z-axis.
Filter parameters for BCF, KF, MW, and MH were modified to minimize the effect of accelerometers on the angle computation. For BCF, γ was set to close to zero such that the majority of computation was based on gyroscopic signals. For KF, R was set to a larger value to minimize the dependency on accelerometers. For MW and MH, β, K p , and K i were all set to zero to compute angles solely using the gyroscopic signals. Since MW and MH computed angles from gyroscopic signals using identical equations (i.e., integrating the quaternion derivative of the IMU and converting the quaternion back to Euler angles), the computed yaw angle from MW and MH was identical. Thus, the RMSE for MW and MH was identical for yaw angles.
AC and DMP could not be used for computing yaw angles. AC could not be implemented since the rotation axis was parallel to the gravity vector, making the accelerometer reading impossible to quantify the relevant motion of the IMU. DMP yaw angles drifted significantly from the beginning, exceeding the RMSE threshold after minutes of running time. The accuracy of DMP was worse than the other fusion filters since DMP uses a default set of parameters for its algorithm which may not guarantee accurate results in our experimental setup. The ability to tune filter parameters is essential for achieving accurate angles [12].
If the z-axis of the IMUs were not parallel to gravity, the computation of the yaw angle would be similar to roll and pitch angles: AC and DMP could be viable options, but the filter parameters of BCF, KF, MW, and MH do not need to be modified, and all methods except for GI would not require high-pass filtering of GI. Conversely, had the x-or y-axes been parallel to gravity then angle calculations for roll or pitch would become difficult and require modifications in preprocessing and filter parameters. Recall that roll, pitch, and yaw are just rotational terms relative to the body, not gravity.

C. Computational Execution Time
The computational execution time for nonsensor fusion algorithms (i.e., GI, AC) was significantly faster than the other sensor fusion algorithms (i.e., BCF, KF, MW, MH). The BCF, a relatively simple sensor fusion algorithm, required the computation of GI and AC simultaneously, increasing its computational time. The KF iteratively and optimally computed the Kalman gain at every time step, which involved more computationally heavy matrix operations, increasing the computational time even more. Finally, the advanced sensor fusion algorithms (i.e., MW, MH) required the most computational time due to the heavy computation of optimization using GD algorithms (MW), more estimations and corrections of the error terms and reference vectors (MH), and conversion from quaternion to Euler angles (MW, MH). Also, the MW and MH computed the full 3-D orientation of the IMUs, whereas KF only calculated a single angle at a time, which resulted in higher computational time for MW and MH. Thus, depending on the available computational resources, one may select more advanced (e.g., KF, MW, MH) or simpler algorithms (e.g., GI, AC, BCF).

D. Notes on High-Pass Filtering of Gyroscopic Signals
A few key points need to be considered when applying a high-pass filter to the raw gyroscopic signals. First, the use of a high-pass filter should be used as a last resort to remove the gyroscopic drift behavior, and the use of sensor fusion algorithms is preferred. This is because high-pass filtered gyroscopic signals cannot be used to estimate static angles, but only the change in angle between two instants of time, therefore not suitable for applications requiring the knowledge of the IMU static angle (e.g., gravity compensation of human limb in the field of biomechanics).
Second, a proper selection of high-pass filter type, order, and cut-off frequency is critical for removing the gyroscopic bias while maintaining the integrity of the signals of interest. In practice, the cut-off frequency of the high-pass filter for raw gyroscopic signals can be set below the lowest frequency of the motion experienced by the gyroscope. However, if the cut-off frequency is set too close to zero, edge effects (i.e., large spikes at the initial and end of the signal) and delays become larger [62]. If the cut-off frequency is too large, it may compromise the integrity of the signals of interest. Certain methods can be used to reduce the magnitude of edge effects [63], such as: 1) padding the data with the reflected version of the signal in the beginning and at the end of the signal or 2) including a small amount of "warm-up" time before conducting the experimental trial.
Once a high-pass filter is properly selected, it can be used to output accurate dynamic angle estimations for a range of speeds rather than a single speed. We demonstrated that, even though we utilized a single high-pass filter, the RMSE of yaw angles was below the acceptable threshold for all movement speeds for all methods except for AC and DMP.

E. Other Limitations and Future Work
Although this study was able to systematically validate various angle computational methods, there were some limitations.
First, our study only focused on computing the relative angle of a single rotation axis at a time. Interestingly, multiaxis rotation may only have a small effect on IMU accuracy as other studies reported similar accuracy during multiaxis versus single-axis rotations [9], [10]. Also, only one IMU was rotated, while the other was kept stationary. Thus, the accuracy of the computed angle depended on the static estimates from the stationary IMU and dynamic estimates of the moving IMU. This setup was appropriate in clinical applications (e.g., spasticity assessment about elbow joints), but may not be appropriate for applications involving extreme movements (e.g., jumping or sprinting). Second, the test apparatus did not have a moment arm to understand the effect of acceleration artifacts due to the offset between the rotation axis and IMU introduced to the system. While our previous studies displayed a similar level of IMU accuracy even when a moment arm was introduced [3], more rigorous testing may be needed to fully evaluate the performances of the computational methods. In cases where the IMUs experience significant motion-related disturbances that introduce a large amount of noise to IMUs, KF, MW, and MH may perform better than other methods since KF, MW, and MH not only fuse multiple sensor readings, but also removes the effect of these external accelerations by updating the gains of the filter every time step (KF) or using a GD-based algorithm (MW) or a PI compensator (MH). The AC and BCF do not have any algorithms for reducing the effect of large external accelerations like KF, MW, and MH, so the accuracy may suffer in dynamic conditions. Third, the maximum angular speed of the moving IMU was not high or low enough to simulate all possible movements observed in biomechanics. Higher speeds may result in higher error regardless of the rotation axis, as observed in other studies [9], [10]. This may be because a greater range of speed requires a different set of optimal filter parameters [12]. In addition, low movement speeds may require different filter parameters (e.g., lower cut-off frequency for a high-pass filter of gyroscopic readings) since the frequency of motion is lower. Fourth, this study did not test the IMU system on human subjects. In real-world biomechanics, the computation of joint angles becomes more complex: 1) IMU axes may not be parallel to the body segment of interest and move separately from the body segment due to soft tissue, irregular muscle geometry, and skin motion artifacts [1]; 2) the rotation axis of human joints is not always clear and may change during motion [12]; and 3) the IMU axes may not be always parallel to the underlying anatomical axes [64]. Fifth, the filter parameters in this study were specifically tuned for our test setup. These parameters should be custom-tuned for different applications, as others have reported [12]. For different applications and models of IMUs, movement speeds, and rotation axis, the filter parameters may be completely different. Lastly, the presented methods in this study required the two 6-axis IMUs to be parallel throughout the data acquisition, which may present a practical limitation due to the irregular shapes of the human body. Other calibration methods can be explored [1]. Albeit these limitations, the proposed 6-axis IMU system with these algorithms can be valuable for certain biomechanical and clinical applications (e.g., slow gait analysis, knee/elbow flexion/extension assessments, spasticity assessment in which one body segment is moving while the other is relatively stationary. For future work, other advanced estimation algorithms ( [17], [18], [65], [66], [67]) that utilize kinematic constraints or spatial relationships between the IMUs can be investigated. These algorithms exploit the kinematic constraints (e.g., ball-socket or hinge joint) to accurately compute relative orientation between two 6-axis IMUs. In addition, the use of combinations of different estimation algorithms for the moving and stationary IMUs can be explored. For example, one may utilize AC (i.e., a simple, drift-free, computationally light, and accurate method) and use only accelerometer readings for computing the orientation of the stationary IMU while using more advanced sensor fusion algorithms that use both gyroscopic and accelerometer readings for computing the orientation of the moving IMU.

V. CONCLUSION
This study aimed to assist researchers new to the field of IMUs by providing a basic understanding of different IMU algorithms and an open-source code of the algorithms. In addition, the study systematically evaluated the performances of various angle computational methods (GI, AC, BCF, KF, DMP, MW, and MH) using 6-axis IMUs without magnetometers on a test apparatus following a strict testing protocol. Most (GI, BCF, KF, MW, and MH for yaw) if not all methods (for pitch/roll) accurately computed relative angles that were within the acceptable RMSE range of 6 • at all speeds and test duration. For yaw, filter parameters were modified, and highpass filtered gyroscopes were used. Even simple algorithms such as GI can output accurate and drift-free angles for various speed and rotation axes by properly preprocessing the raw signals. Proper filtering parameters for methods and preprocessing are essential and can be done by utilizing the information (e.g., movement speed, rotation axis) of the motion that IMU experiences.