2.2.1. Calibration

When an IMU is at rest, the accelerometer should only measure gravity, and the gyroscope should show zero rotational velocity. In contrast, they always show some data while still in real-world circumstances. On the other hand, the magnetometer has issues with magnetic interference [21]. To minimize these issues, calibration was therefore required for the magnetometer, gyroscope, and accelerometer. With the aid of WitMotion's built-in software, each of the three sensors—accelerometer, gyroscope, and magnetometer—were calibrated. Gyroscope calibration was sufficient using WitMotion software, but the accelerometer and magnetometer required more than this procedure. In light of this, further calibration of these two sensors was essential, and an additional simple calibration was performed as an early step in the algorithm to reduce or eliminate the bias in data after performing calibration using the manufacturer software. If the sensor was not calibrated using this two-stage algorithm, the results would be different from the actual value.

This calibration process and the calculation are simple, less time-consuming, and easy to follow. Thus, this extra process did not increase the calculation burden. Table 2 depicts the shortcomings of various calibration techniques described in other research articles and the advantages of this proposed calibration process over the previous research.

**Table 2.** Inertial sensor calibration shortcomings mentioned in the previous literature and advantages of the proposed calibration process.


The acceleration (*a<sup>i</sup>* ), angular velocity (*g<sup>i</sup>* ), and magnetic field strength (*m<sup>i</sup>* ) of the IMU about the *x*, *y*, and *z* axes were measured with a tri-axis accelerometer, gyroscope, and magnetometer, and they were expressed as

$$a^i = \begin{bmatrix} a\_{xi} & a\_{yi} & a\_{zi} \end{bmatrix} \tag{1}$$

$$\mathbf{g}^i = \begin{bmatrix} \omega\_{xi} & \omega\_{yi} & \omega\_{zi} \end{bmatrix} \tag{2}$$

$$m^i = \begin{bmatrix} m\_{xi} & m\_{yi} & m\_{zi} \end{bmatrix} \tag{3}$$

where *i* = 1 (for forearm sensor), 2 (for upper arm sensor), 3 (for reference sensor).

In this system, the accelerometer was calibrated using the bias correction vector. During the measurement, the sensors were in static condition for at least 10 s. The first hundred raw accelerometer data were used for determining the bias correction vector for the calibration of the accelerometer. The bias correction vector and calibrated value of the accelerometer were calculated from Equations (4) and (5), respectively.

$$b\_i = \begin{bmatrix} \frac{\sum\_{0}^{100} \left( a\_{xi\_{bias}} \right)}{100} & \frac{\sum\_{0}^{100} \left( a\_{yi\_{bias}} \right)}{100} & \frac{\sum\_{0}^{100} \left( a\_{zi\_{bias}} \right)}{100} \end{bmatrix} \tag{4}$$

$$a\_c^k = (a^k - b\_i) \tag{5}$$

where

*bi* = bias in accelerometer readings;

*a<sup>k</sup>* = accelerometer reading after calibration by WitMotion software;

*ak <sup>c</sup>* = calibrated accelerometer data.

When there was a value of roll and pitch, the magnetometer data in the local coordinate system were transferred to the global coordinate system. This transformation

was performed using Equations (6)–(8), where *mGx*, *mGy*, and *mGz* are the magnetometer measurements in the global coordinate frame in X, Y, and Z axes.

$$m\_{\rm Gx} = m\_{\rm xi} \cos \theta + m\_{\rm yi} \sin \theta \sin \varphi + m\_{\rm zi} \sin \theta \cos \varphi \tag{6}$$

$$m\_{Gy} = m\_{yi} \cos \varphi - m\_{zi} \sin \varphi \tag{7}$$

$$m\_{Gz} = -m\_{\rm xi} \sin \theta + m\_{\rm yi} \sin \varphi \cos \theta + m\_{\rm zi} \cos \varphi \cos \theta \tag{8}$$

The yaw angle, calculated using Equations (6)–(8), is the heading of the sensor with respect to the magnetic north. However, the magnetometer data were converted to get the initial zero yaw angle using Equations (9)–(14).

$$m\_{i(avg)} = \begin{bmatrix} m\_{X\_{iavg}} \ m\_{y i\_{avg}} m\_{z i\_{avg}} \end{bmatrix} = \begin{bmatrix} \frac{\sum\_{0}^{100} \binom{m}{x\_{i}}}{100} & \frac{\sum\_{0}^{100} \binom{m}{y\_{i}}}{100} & \frac{\sum\_{0}^{100} \binom{m}{z i\_{j}}}{100} \end{bmatrix} \tag{9}$$

$$
heading = \begin{bmatrix} 0 & 0 & A \tan 2 \left( m\_{y\_{\text{avg}} \prime} \, m\_{x i\_{\text{avg}}} \right) \end{bmatrix} \tag{10}
$$

$$q\_h = \begin{bmatrix} q\_{h0} & q\_{h1} & q\_{h2} & q\_{h3} \end{bmatrix} \tag{11}$$

$$q\_{m\_k} = \begin{bmatrix} 0 & m\_{Gx} & m\_{Gy} & m\_{Gz} \end{bmatrix} \tag{12}$$

$$q\_{m\_c^k} = \begin{pmatrix} q\_h \odot & q\_{m\_k} \end{pmatrix} \odot \ q\_h^\* = \begin{bmatrix} q\_a & q\_b & q\_c & q\_d \end{bmatrix} \tag{13}$$

$$m\_c^k = \begin{bmatrix} q\_b & q\_c & q\_d \end{bmatrix} \tag{14}$$

where

*mi*(*avg*) = average of first hundred data of magnetometer along X, Y, and Z axes;

*qh* = quaternion corresponding to the heading angle;

*qmk <sup>c</sup>* = quaternion corresponding to the corrected magnetometer data;

*mk <sup>c</sup>* = corrected magnetometer data.

#### 2.2.2. First Stage: Estimation of Orientation from Sensor Fusion Technique

The Madgwick algorithm (MAD) was implemented in this study as an orientation filter. The MAD is a sensor fusion algorithm that combines the data from the accelerometer, gyroscope, and magnetometer in a quaternion form. This algorithm computes a quaternion derivative of the IMU, which can be converted to Euler angles. Although the sample rate restricts the range of the motion or joint angle measurement accuracy, satisfactory behavior was recorded even when utilizing low sampling rates, close to 10 Hz, and this level of accuracy may be suitable for applications involving human motion [32]. The Madgwick filter is an algorithm based on gradient descent that can compensate gyroscope drift [32]. The advantages of the Madgwick filter are its low computational load and complexity, low sampling rate requirement, that it can compensate for magnetic distortion, and that it has a simple tuning method where there is only one adjustable parameter. Since MAD uses a quaternion representation, the filter is not subjected to the gimbal lock problem associated with Euler angle representation. The Madgwick algorithm is also very convenient for small-size, low-cost, and low-energy-consumption microprocessors.

The implementation of the proposed algorithm was carried out with the assistance of Matlab R2022a. After going through the calibration steps, the data from three IMUs (forearm, upper arm, and reference) acted as an input for the Madgwick filter. The filter's output was three different sets of quaternions for the upper arm, forearm, and reference sensors. These three sets of quaternions were utilized to compute the joint angle of the rigid body. The Madgwick filter's precise details are stated in Algorithm 2.

#### **Algorithm 2** Madgwick Filter [32]

First Step: Computation of the orientation from the gyroscope. Gyroscope measurement: *S<sup>ω</sup>* = **[0** *ω<sup>x</sup> ω<sup>y</sup> ωz***]** Quaternion derivative: *<sup>S</sup> E* **.** *qw***,***<sup>t</sup>* = **<sup>1</sup> 2** *S <sup>E</sup>*6*qest***,***t*−**1**<sup>⊗</sup> *<sup>S</sup>ω<sup>t</sup>* Orientation from Gyroscope, *<sup>S</sup> <sup>E</sup>qω***,***<sup>t</sup>* <sup>=</sup> **<sup>1</sup> 2** *S <sup>E</sup>*6*qest***,***t*−**<sup>1</sup>** <sup>+</sup> *<sup>S</sup> E* **.** *qw***,***t***Δ***t* Second Step: Use of accelerometer data to get the orientation quaternion. Sensor Orientation: *<sup>S</sup> <sup>E</sup>*6*<sup>q</sup>* <sup>=</sup> **[***q***<sup>1</sup>** *<sup>q</sup>***<sup>2</sup>** *<sup>q</sup>***<sup>3</sup>** *<sup>q</sup>***4]** Predefined reference direction of the field in the earth frame:

$$E\_{\hat{d}} = \begin{bmatrix} \mathbf{0} & d\_{\mathbf{x}} & d\_{y} & d\_{z} \end{bmatrix}$$

Measurement of the field in the sensor frame, *SS*<sup>6</sup> <sup>=</sup> **[0** *sx sy sz***]** The sensor orientation *<sup>S</sup> <sup>E</sup>*6*<sup>q</sup>* can be formulated as an optimization problem by

**min** *S <sup>E</sup>*6*<sup>q</sup>* <sup>∈</sup>*R***<sup>4</sup>** *f***(** *S <sup>E</sup>*6*q***,***Ed*6**,***SS*6**)** where the objective function *f***(** *S <sup>E</sup>*6*q***,***Ed*6**,***SS*6**)** can be calculated by

$$(\hat{\mathfrak{q}}\_E^\* \hat{\mathfrak{q}}^\* \otimes E\_{\hat{\mathfrak{q}}} \otimes {}^{\hat{\mathfrak{q}}}\_E \hat{\mathfrak{q}} - \mathcal{S}\_{\hat{\mathfrak{q}}})^\ast$$

Using Gradient Descent algorithm, estimated orientation based on previous one and *μ* step size:

$${}^{S}\_{E}\widehat{q}\_{k+1} = {}^{S}\_{E}\widehat{q}\_{k} - \mu \frac{\mathsf{PV}\_{E}^{S}\widehat{q}\_{\prime}E\_{\widehat{d}^{\prime}}S\_{\widehat{S}}}{||\mathsf{PV}\_{E}^{S}\widehat{q}\_{\prime}E\_{\widehat{d}^{\prime}}S\_{\widehat{S}}\||}$$

where ࢺ*f***(** *S <sup>E</sup>*6*q***,***Ed*6**,***SS*6**)** <sup>=</sup> *JT***(** *S <sup>E</sup>*6*qk***,***Ed*6**)** *<sup>f</sup>***(** *S <sup>E</sup>*6*q***,***Ed*6**,***SS*6**)** For accelerometer, the Objective function and the Jacobian matrix are *fg***(** *S <sup>E</sup>*6*q***,***Sa*6**)** and *Jg***(** *S <sup>E</sup>*6*q***)**. Third Step: Use of magnetometer data. Earth's magnetic field: *<sup>E</sup>*6*<sup>b</sup>* <sup>=</sup> **[0** *bx* **<sup>0</sup>** *bz***]** Magnetometer measurement: *Sm*<sup>6</sup> <sup>=</sup> **[0** *mx my mz***]** For Magnetometer the Objective function and the Jacobian matrix are *fb***(** *S <sup>E</sup>*6*q***,***E*6*b***,***Sm*<sup>6</sup> **)** and *Jb***(** *S <sup>E</sup>*6*q***,***E*6*b***)** respectively. Complete solution considering accelerometer and magnetometer: Objective function: *fg***,***b***(** *S <sup>E</sup>*6*q***,***Sa***,** *<sup>E</sup>*6*b***,***Sm*<sup>6</sup> **)** <sup>=</sup> *fg***(** *S <sup>E</sup>*6*q***,***Sa*6**)** *fb***(** *S <sup>E</sup>*6*q***,***E*6*b***,***Sm*<sup>6</sup> **)** Jacobian matrix: *Jg***,***b***(** *S <sup>E</sup>*6*q***,***E*6*b***)** <sup>=</sup> *Jg***(** *S E*6*q***)** *Jb***(** *S <sup>E</sup>*6*q***,***E*6*b***)** Estimated Orientation, *<sup>S</sup> <sup>E</sup>q∇***,***<sup>t</sup>* <sup>=</sup> **<sup>1</sup> 2** *S <sup>E</sup>qest***,***t*−**<sup>1</sup>** <sup>−</sup> *<sup>μ</sup> <sup>∇</sup><sup>f</sup> <sup>f</sup>*ࢺ where, *∇<sup>f</sup>* <sup>=</sup> *JT <sup>g</sup>***,***b***(** *S <sup>E</sup>*6*q***,***E*6*b***)***fg***,***b***(** *S <sup>E</sup>*6*q***,***Sa***,** *<sup>E</sup>*6*b***,***Sm*<sup>6</sup> **)** and *<sup>μ</sup>* <sup>=</sup> *<sup>α</sup><sup>S</sup> E* **.** *qw***,***t***Δ***t* **,** *α***>1** Final Step: Final estimation of the orientation: *S Eqest***,***<sup>t</sup>* = *γt S Eq∇***,***<sup>t</sup>* + **(1** − *γ***)** *S Eq<sup>ω</sup>***,***<sup>t</sup>* The simple expression after some simplifications and assumptions:

$$\,\_{\underline{E}}^{\underline{\rm S}}q\_{est,t} = -\beta \, \frac{\nabla f}{||\nabla f||} \Delta t + \,\_{\underline{E}}^{\underline{\rm S}} \hat{q}\_{est,t-1} + \,\_{\underline{E}}^{\underline{\rm S}} \dot{q}\_{\rm W,t} \Delta t$$
