Generally speaking, a 3D human motion is considered as a set of multiple joints and limbs movement, which can be described using mathematical expressions of coordinates and direction. As shown in Figure 1, the whole body is divided into five connected parts, namely left upper limb,

right upper limb, trunk, left lower limb and right lower limb. We use dots to represent joints and lines for limbs. The whole human body acts as a unity of these limbs and joints.

**Figure 1.** Human skeleton model, and the key distance and angle parameters related with human motion tracking. (**a**) the entire human body and (**b**–**f**) the distances and angles between joint and limbs.

As shown in Figure 1, the human body consists of 15 joints and 14 limbs. To better describe the complicated human body movements, Denavit–Hartenbe (D-H) equation [22,25] is generally used in the description of HMT features and parameters. Its measurement parameters mainly include the distances between joints (Figure 1b,c) joint and limbs (Figure 1d) and the relative angles between limbs (Figure 1e,f). To simplify the above tracking problem, we abstract the HMT serialization process as shown in Figure 1, representing the target's coordinate and measurement parameters in time sequence *T* = {*t*0, *t*1, ··· , *tn*}.

**Figure 2.** The definition of Gauss–Markov random process towards human motion tracking and its key parameters.

On the basis of random motion model [25], we assume a Gauss–Markov process to describe the movement of the joints of the body, as shown in Figure 2. Thanks to Inertial Measurement Unit (IMU), acceleration and angular data can be obtained and, furthermore, the rotation angle of the target position within the history frames can be calculated by means of Attitude and Heading Reference System (AHRS) [20]. In addition, distance-related information can be obtained using the TOA technique [19], after a few simple calculations. Compared with IMU-based dead reckoning method [25], TOA-based tracking method does not have to face the problem of drifting.

#### *3.2. Error Sources*

As mentioned above, the measurement parameters of the system are derived from IMU and TOA. The measurement error sources can be divided into two aspects: distance measurement error and angle measurement error. Thereinto, target node's distance parameter ˆ*dt* at time *t* is defined as:

$$d\_t = d\_t + n\_d \tag{1}$$

where *dt* is the actual distance value and *nd* is Gauss distributed measurement noise whose variance is *σ*2 *d* , namely *nd* ∼ N (0, *σ*2 *d* ). Vector *d*ˆ *t* = [ ˆ*d*1, ˆ*d*2, ··· , ˆ*dn*+*m*−1] *T* is introduced to collect distance information, namely the distance between target node and reference nodes, respectively, standing for one body joint.

Horizontal heading estimates from inertial-sensor-based approach *α* reads:

$$
\hbar\_t = \mathfrak{u}\_t + \mathfrak{u}\_{t\prime}\mathfrak{u}\_t \sim \mathcal{N}(0, \mathfrak{e}\_t^2) \tag{2}
$$

where *αt* is the true value of actual horizontal heading. *ut* is uncorrelated zero-mean Gaussian random variable with variance 2 *t* , which is independent from z direction and is also uncorrelated with ˆ*d*. The horizontal heading is an instantaneous measurement parameters, which may introduce cumulative error by the characteristics of IMU. Namely, *α*ˆ*t* is a time-dependent parameter. Vector *αt* = [*α*1, *α*2, ··· , *αNα* ] *T* is introduced to collect *αt*, i.e., the angular cumulative information in the past *Nα* frames.

Vertical elevation estimates from inertial-sensor-based approach *β* reads:

$$
\hat{\beta}\_t = \beta\_t + v\_{t\prime} v\_t \sim N(0, \zeta\_t^2) \tag{3}
$$

where *βk* is the actual vertical elevation. *vt* is uncorrelated zero-mean Gaussian random variable with variances *ξ*2 *t* , which is independent from z direction and is also uncorrelated with both *α*ˆ and ˆ*d*. *β*ˆ*t* is also a time-independent parameter. Vector *βt* = [*β*1, *β*2, ··· , *βNβ* ] *T* is introduced to collect *βk*, i.e., the angular cumulative information in the past *Nβ* frames. As above mentioned angles *α* and *β* are all measured from IMU, thus, for ease of expression, we take *σIMU* as the variance of both these two parameters.

#### *3.3. Model Presentation*

In continuous human motion tracking process, some body parts remain relatively steady in comparison to the end joints or limbs of the human body, such as the trunk. Thus, these positions may be set as reference nodes, defined as P*a* = {1, ..., *m*}, in which the *k*th reference node's coordinate is defined as *ak*. *m* indicates the number of reference nodes. In addition, the position state of the target is denoted by *pt* = [*pX t pY t pZ t* ] *T*, where *pX t* , *pY t* and *pZ t* are the coordinates in the 3-D positioning system, and *T* is the transpose operator. The target nodes set is represented as P*x* = {1, 2, ..., *n*}, and the *i*th target node's coordinate is defined as *pi*. *n* indicates the number of target nodes.

According to the Bayesian theorem [22,25], at time *t*, the state estimation *xt*−1 could be represented as:

$$p\_t = \mathbf{f}\_t(p\_{t-1}) + \mathbf{q}\_t \tag{4}$$

where Equation (4) is the prediction equation, and **f***t*(·) is the transition equation, and the noise of *xt* follows standard normal distribution with variance *Qt*, namely **q***t* ∼ N (0, *Qt*).

To describe the measurement parameters of human motion in a more general way, we define the measurement parameters of the system as:

$$\mathbf{z}\_l = \mathbf{h}\_l \left( \mathbf{d}(p\_l, k), \mathbf{q} \right) + \mathbf{r}\_l \tag{5}$$

where the measurement vector **z***t*, whose noise follows normal distribution, i.e., **r***t* ∼ N (0, *Rt*), contains two parts, namely TOA measurement parameter (distance) and IMU measurement parameter (acceleration, angular velocity, etc.). **ht**(·) is nonlinear observation equation, which is related with actual measurement values of target **pt**, including **d**(**pt**, **k**) (the distance vector at time *t*), and *ϕ* (the angle vector at time *t*). **k** = [*k*1, ··· , *kj*, ··· , *km*+*n*−1] is distance measurement indicator, such as line-of-sight (LOS)/non-line-of-sight (NLOS), the reference node accuracy indicator, etc. *m* + *n* − 1 is the total number of both reference and target nodes. *ϕ* = [*ϕ*1, ··· , *ki*, ··· , *ϕs*] is the vector of angle measurement values, whose amount is *s*. *k* and *ϕ* are auxiliary parameters that are not essential but may help increase the localization accuracy and could be calculated in real time. The effects of these parameters can be verified in subsequent simulation experiments.

#### **4. IMU/TOA Fusion Based HMT Method**

Our analysis considers the possible unknown random factors that may influence the human motion tracking. Considering the characteristics of sequential human motion tracking, we define the measuring variable at time slot *t* as *θt* (for simplicity, it is abbreviated as *θ*), namely

$$\boldsymbol{\theta} \stackrel{\scriptstyle \Delta}{=} \begin{bmatrix} \boldsymbol{p}\_t^T & \boldsymbol{p}\_{t-1}^T & \boldsymbol{k}^T & \boldsymbol{\varrho}^T \end{bmatrix}^T \tag{6}$$

where *pt* is the target position vector at the moment, and *pt*−1 is historical target location information at the former time slot.

Cramér–Rao lower bound (CRLB) [22] is defined as the inverse of Fisher Information Matrix (FIM) and represents the theoretical lower bound of the observed variance of unbiased estimators [22,25]. Suppose that *p*(*θ*, *zt*) represents the joint probability density function of observation vectors *zt* and *θ*, then FIM could be derived as the variance of its log likelihood function gradient, namely

$$\mathbf{J}(\theta) \stackrel{\triangle}{=} \mathbb{E}[\Delta\_{\theta}^{\theta} \ln p(\theta, \mathbf{z}\_{\mathbf{t}})] \tag{7}$$

where E{·} represent the expectation of parameters. ∇*θ* = [ *∂ ∂θ*1 , ··· , *∂ ∂θn* ], Δ*β α* = ∇*β*∇*T α* . Thus,

$$\mathbf{Cov}\_{\theta}(\boldsymbol{\theta}) \succeq \left\{ \mathbf{J}(\boldsymbol{\theta}) \right\}^{-1} \tag{8}$$

where *A* ≥ *B* for the sake of simplicity, represents that matrix *A* − *B* is non-negatively defined.

According to the Bayesian theorem, *p*(*θ*, *zt*) = *p*(*zt*|*θ*)*p*(*θ*), then **J**(*θ*) could be divided into two parts, namely

$$\mathbf{J}(\theta) = \mathbf{J}\_{\mathbf{D}}(\theta) + \mathbf{J}\_{\mathbf{P}}(\theta) \tag{9}$$

where **JD** is the information matrix obtained by measuring parameters. **JP** is the information matrix obtained by prior information. We describe the representation of these two variables separately in the following section.

#### *4.1. Measurement Parameter Information Matrix*

Due to the measurement equation, we get **h** = **h***t* (**d**(*pt*, *k*),*ϕ*), and then, based on the chain rules [48], **JD** could be represented as:

$$\mathbf{J}\_{\mathbf{D}} = \mathbf{H} \cdot \mathbf{J}\_{\mathbf{h}} \cdot \mathbf{H}^{T} \tag{10}$$

where **H** = [∇*θh*], *Jh* is the FIM of *h*, namely

$$\mathbf{J\_{h}} = \mathbb{E}\left\{ \nabla\_{\mathbf{h}} \ln p(\mathbf{z}\_{l}|\boldsymbol{\theta}) \left[ \nabla\_{\mathbf{h}} \ln p(\mathbf{z}\_{l}|\boldsymbol{\theta}) \right]^{T} \right\} \tag{11}$$

where transition matrix **H** could be further decomposed into four components

$$\mathbf{H} = \begin{bmatrix} \mathbf{H}\_t & \mathbf{H}\_{t-1} & \mathbf{K} & \boldsymbol{\Phi} \end{bmatrix}^T \tag{12}$$

where **H***t* = [∇*pth*], **H***t*−1 = [∇*pt*−1*h*], **K** = [∇*kh*], *Φ* = [∇*ϕh*]. Because the distance is a instantaneous measurement parameter, it has no error accumulation problems. Thus, **d** has no relationship with prior information *pt−***1**, namely **H***t*−1 = 0. Then,

$$\mathbf{H} = \begin{bmatrix} \mathbf{H}\_t & \mathbf{0} & \mathbf{K} & \boldsymbol{\Phi} \end{bmatrix}^T \tag{13}$$

For **J***h*, we can use diagonal matrices of order *n* to represent it,

$$\mathbf{J}\_{\mathbf{h}} = \mathbf{A} = \mathbf{diag}(\lambda\_1, \dots, \lambda\_{\mathbf{j}}, \dots, \lambda\_n) \tag{14}$$

where *λj* stands for the *j*th measurement value. Then, **JD** is represented as a (4 + *m* + *n* + *s*) × (4 + *m* + *n* + *s*) matrix,

$$\begin{aligned} \mathbf{J}\_{\mathbf{D}} &= \mathbf{H} \cdot \mathbf{A} \cdot \mathbf{H}^{T} \\ &= \begin{bmatrix} \mathbf{D}\_{11} & \mathbf{0} & \mathbf{D}\_{13} & \mathbf{D}\_{14} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{D}\_{13}^{T} & \mathbf{0} & \mathbf{D}\_{33} & \mathbf{D}\_{34} \\ \mathbf{D}\_{14}^{T} & \mathbf{0} & \mathbf{D}\_{34}^{T} & \mathbf{D}\_{44} \end{bmatrix} \end{aligned} \tag{15}$$

where

$$\begin{aligned} \mathbf{D}\_{11} &= \mathbf{H}\_{t} \boldsymbol{\Lambda} \mathbf{H}\_{t}^{T} & \mathbf{D}\_{13} &= \mathbf{H}\_{t} \boldsymbol{\Lambda} \mathbf{K}^{T} \\ \mathbf{D}\_{14} &= \mathbf{H}\_{t} \boldsymbol{\Lambda} \boldsymbol{\Phi}^{T} & \mathbf{D}\_{33} &= \mathbf{K} \boldsymbol{\Lambda} \mathbf{K}^{T} \\ \mathbf{D}\_{34} &= \mathbf{K} \boldsymbol{\Lambda} \boldsymbol{\Phi}^{T} & \mathbf{D}\_{44} &= \boldsymbol{\Phi} \boldsymbol{\Lambda} \boldsymbol{\Phi}^{T} \end{aligned} \tag{16}$$

#### *4.2. State Parameter Information Matrix*

It could be known by the whole probability formula, for the variable *θ*, the probability *p*(*θ*) = *p*(*pt|pt−***1**)*p*(*k*)*p*(*ϕ*). Then, log-likelihood function could be represented as:

$$\ln p(\theta) = \left[ \ln p(p\_t | p\_{t-1}) \right] + \ln p(k) + \ln p(\mathcal{g}) \tag{17}$$

where *p*(*k*) and *p*(*ϕ*) are corresponding independent information of *pt* and *pt*−1. We decompose the vector *θ* into two components, namely state vector [*pt pt−***1**] *T* and observation vector [*k ϕ*] *T*. Thus, *JP* could be represented as:

$$\begin{aligned} \mathbf{J}\_{\mathbf{P}} &= \mathbb{E}[\boldsymbol{\Delta}\_{\boldsymbol{\theta}}^{\boldsymbol{\theta}} \ln p(\boldsymbol{\theta})] \\ &= \begin{bmatrix} \mathbf{J}\_{\mathbf{P}\_{11}} & \mathbf{J}\_{\mathbf{P}\_{12}} \\ \mathbf{J}\_{\mathbf{P}\_{12}}^{T} & \mathbf{J}\_{\mathbf{P}\_{22}} \end{bmatrix} \end{aligned} \tag{18}$$

where *JP***11** could be obtained with the recursion of vectors *pt* and *pt−***1**, which, according to Tichavsky et al. [43], could be described as:

$$\mathbf{J}\_{\mathbf{P}\_{11}} = \begin{bmatrix} \mathbf{M}\_{11} & \mathbf{M}\_{12} \\ \mathbf{M}\_{12}^T & \mathbf{M}\_{22} + \mathbf{J}(p\_{t-1}) \end{bmatrix} \tag{19}$$

*Electronics* **2019**, *8*, 142

where

$$\begin{aligned} \mathbf{M}\_{11} &= \mathbf{Q}\_{t}^{-1} \\ \mathbf{M}\_{12} &= \nabla\_{p\_{t-1}} \mathbf{f}\_{t}(p\_{t-1}) \mathbf{Q}\_{t}^{-1} \\ \mathbf{M}\_{22} &= \nabla\_{p\_{t-1}} \mathbf{f}\_{t}(p\_{t-1}) \mathbf{Q}\_{t}^{-1} \left[\nabla\_{p\_{t-1}} \mathbf{f}\_{t}(p\_{t-1})\right]^{T} \end{aligned} \tag{20}$$

**J**(*pt*−1) is the FIM of *pt*−1. As *p***(***k***)** and *p***(***ϕ***)** are independent from *pt* and *pt*−1, *JP***12** is **0**, namely

$$\mathbf{J}\_{\mathbf{P}\_{12}} = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \tag{21}$$

Similarly, the prior probability *p*(*pt*|*pt*−1) are independent from *k* and *ϕ*. Thus, *JP***12 =***JT P***21 = 0**. Finally,

$$\mathbf{J} \mathbf{p}\_{22} = \begin{bmatrix} \mathbf{J}\_{\mathbf{K}} & \mathbf{0} \\ \mathbf{0} & \mathbf{J}\_{\Phi} \end{bmatrix} \tag{22}$$

where **JK** and **J***Φ* are, respectively, the FIM of vectors *k* and *ϕ*:

$$\begin{aligned} \mathbf{J}\_{\mathbf{K}} &= \mathbb{E}[\boldsymbol{\Delta}\_{\mathbf{k}}^{\mathbf{k}} \ln p(\boldsymbol{k})] \\ \mathbf{J}\_{\boldsymbol{\Phi}} &= \mathbb{E}[\boldsymbol{\Delta}\_{\boldsymbol{\Phi}}^{\boldsymbol{\Phi}} \ln p(\boldsymbol{\varrho})] \end{aligned} \tag{23}$$

#### *4.3. Integral Information Matrix*

Due to the analysis above, substituting Equations (15) and (18) into Equation (9), the FIM of variable *θ* could be represented as:

$$\mathbf{J}(\boldsymbol{\theta}) = \begin{bmatrix} \mathbf{M}\_{11} + \mathbf{D}\_{11} & \mathbf{M}\_{12} & \mathbf{D}\_{13} & \mathbf{D}\_{14} \\ \mathbf{M}\_{12}^T & \mathbf{M}\_{22} + \mathbf{J}(\boldsymbol{x}\_{t-1}) & \mathbf{0} & \mathbf{0} \\ \mathbf{D}\_{13}^T & \mathbf{0} & \mathbf{D}\_{33} + \mathbf{J}\mathbf{x} & \mathbf{D}\_{34} \\ \mathbf{D}\_{14}^T & \mathbf{0} & \mathbf{D}\_{34}^T & \mathbf{D}\_{44} + \mathbf{J}\boldsymbol{\phi} \end{bmatrix}. \tag{24}$$

However, only the lower bound of *pt* is of interest, namely a small submatrix of **J**(*θ*), **J**−1(*pt*) = [**J**−1(*θ*)]2×2.

Furthermore, due to the Schur complement theorem [49], **J**(*pt*) could be divided into two components, namely state matrix *JS* and measurement matrix *JC*, then

$$\mathbf{J}(p\_t) = \mathbf{J}\_S - \mathbf{J}\_C \tag{25}$$

where

$$\begin{aligned} \mathbf{J\_S = M\_{11} + D\_{11} - M\_{12} \left(\mathbf{M\_{22}} + \mathbf{J}(\boldsymbol{\mu}\_{t-1})\right)^{-1} \mathbf{M\_{12}^T} \\ \mathbf{J\_A = \begin{bmatrix} \mathbf{D\_{13}} & \mathbf{D\_{14}} \end{bmatrix} \begin{bmatrix} \mathbf{D\_{33} + \mathbf{J\_K} & \mathbf{D\_{34}} \\ \mathbf{D\_{34}^T} & \mathbf{D\_{44} + \mathbf{J\_C}} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{D\_{13} & \mathbf{D\_{14}} \end{bmatrix}^T \end{aligned} \tag{26}$$

As mentioned above, we demonstrate a comprehensive analysis of how distance and angle measurements are related with human motion tracking accuracy in IMU/TOA fusion systems. The detailed calculation process is shown in Algorithm 1. In the next two sections, we further evaluate the overall performance of the fusion system from two aspects, namely spatial and temporal performance, in theory.

#### **Algorithm 1** Proposed IMU/TOA fusion-based human motion tracking method.

0: **Initialization**:

1: Define measuring variable at time slot *t* as *θt* - *pT t pT t*−1 *kT ϕT T* . Set the time slot index

*t* = 0. Choose an initial value for *pt*−1 of variable *θt*, generally as 0.

#### **Part 1: measurement parameter information**

#### **Part 2: state parameter information**

9: Update the information matrix obtained by prior information, i.e., **JP** = ⎣ **JP**11 **JP**12 **J***T* **P**12 **JP**22 ⎦.

#### **Part 3: integral information**

#### **Part 4: recursive computation**

13: *t* + 1 → *t*, go to step 2.

#### **5. Spatial Performance Analysis**

When it comes to spatial performance, we suppose that the tracking system has time independence, namely *J***(***pt***)** and *J***(***pt−***1)** are independent from each other. Thus, *M***12**, *M***22**, *Jpt−***1** , *JK*, *D***13**, *D***33** and *D***34** are all **0**, namely

$$\mathbf{J}(p\_l) = \mathbf{M}\_{11} + \mathbf{D}\_{11} - \mathbf{D}\_{14} [\mathbf{D}\_{44} + \mathbf{J}\_{\Phi}]^{-1} \mathbf{D}\_{14}^T \tag{27}$$

⎡

⎤

Movement characteristics of each body part should be taken into consideration. Because of the spatial features of human motion sensing, it could be classified into two categories: 2-D and 3-D conditions. 2D motions refer to the body movements that can be captured by a plane, such as raising one's arms into Y pose or T pose. 3D motions refer to stereoscopic movements such as walking, running or swing hands. The different relative positions may contribute to various capture accuracy. Therefore, both the CRLB of human body in 2D and 3D scenarios are taken into consideration. CRLB of each location is derived as following demonstration and performance of proposed fusion method is verified. See Appendix A for the detailed derivation process of CRLB for spatial performance analysis.

#### *5.1. Scenario Setup*

As mentioned above, the trunk remains relatively steady in comparison to the end joints or limbs of the human body. In consideration of end-effector problem, we set the joints of neck, chest, left and right hip as reference nodes for error bounds estimation. A constrained scenario of 2 m × 2 m is set up for further experiments. The entire human body is supposed to located in the center of this region. In contrast, two sequences of reference nodes combination are considered, namely *Case*1 = {*Neck*, *Chest*, *LHip*, *RHip*} and *Case*2 = {*Neck*, *LShoulder*, *RShoulder*, *Chest*, *LHip*, *RHip*}. The body is allowed to move freely in this space.

#### *5.2. Performance Analysis*

To demonstrate the typical calculation results, CRLB, when the variance of distance measurement *σd* is 0.2 and *σIMU* is 0.1, is selected upon common commercial motion sensing systems, such as Xsens. From the experiment results, as shown in Figure 3, the following conclusions could be obtained:

**Figure 3.** *Cont.*

**Figure 3.** CRLB comparison in two-dimensional and three-dimensional condition: (**a**,**b**,**e**) CRLB in condition that both distance and IMU information are considered using proposed model; and (**c**,**d**,**f**) CRLB considered with distance information solely. *σd* is set as 0.2 and *σIMU* is 0.1.

To clearly show the different capture capability of distance measurement and IMU, Figure 4 is presented. As shown, under different *σ*2 *d* , CRLB increases along with the increase of *σ*2 *IMU*. When *σ*2 *IMU* is fixed, CRLB also increases along with the increase of *σ*2 *d* . However, the change is more apparent when *σ*2 *d* is fixed. It means that the accuracy of IMU measurement matters more in the fusion process, but the distance measurement helps restrain the bounds of CRLB. It is worth mentioning that, for a given *σ*2 *IMU*, CRLB shows little variation with the increase of *σ*2 *d* if it is larger than 0.2. It indicates that the distance measurement may be not the main limiting factor of CRLB, but *σ*2 *IMU* is.

**Figure 4.** CRLB distribution under different *σ*2 *d* when *σ*2 *IMU* changes.

#### **6. Temporal Performance Analysis**

We suppose that the tracking system has cumulative recursion characteristics, namely *J***(***pt***)** and *J***(***pt−***1)** are not independent from each other. Based on [22], Equation (26) could be simplified into a recursive form, namely Posterior Cramér–Rao Lower Bound (PCRLB) [22,25].

$$\mathbf{J}(\mathbf{x}\_{t}) = \mathbf{M}\_{11} + \mathbf{D}\_{11} - \mathbf{M}\_{12} \left[\mathbf{M}\_{22} + \mathbf{J}(\mathbf{x}\_{t-1})\right]^{-1} \mathbf{M}\_{12}^{T} \tag{28}$$

See Appendix B for the detailed derivation process of PCRLB for temporal performance analysis.

#### *6.1. Scenario Setup*

Even though there are elegant expressions to recursively calculate the FIM, Equation (27) usually does not have analytical close-form solution. To deal with that, we employ the Monte Carlo approach [50] to convert those continuous integrals into discrete summations, and finally work out the PCRLB. The root-mean-square of PCRLB is given by 1 *L* ∑*L i*=1 *Pi k*, where *Pi k* is the PCRLB on the Root-Mean-Square Error (RMSE) of joint at step *k* in the *i*th Monte Carlo trial. *L* is the total Monte Carlo trial number (*L* = 1000 was used in this study). Note that, for each Monte Carlo trial, we randomly selected the initial location for PCRLB calculation to get a fair average of the entire human body moving process.

#### *6.2. Performance Analysis*

We mainly analyzed the influence of different factors on PCRLB from the following two perspectives: fusion method factor and topology factor.

**(1) Fusion method factor:** Considering different *Qt* and *Rt* in IMU/TOA fusion methods, as well as using sole TOA or IMU method, may lead to different results of human motion tracking accuracy.

Figure 5 shows in IMU/TOA fusion methods, the lower bounds under different *Qt* and *Rt* as well as in the use of sole TOA or IMU method. As shown in the figure:

1. When only inertial sensors are adopted in the tracking system, as indicated by the black solid line in Figure 5, the accumulative errors may tend to be diverging. Theoretically, this confirms that IMU based HMT system faces the problem of accumulative errors. However, IMU/TOA fusion method can avoid this divergence. The performance curves of proposed approaches achieve stability after certain steps, i.e., their errors converge.

**Figure 5.** Temporal performance as a function of step index.

**(2) Topology factor:** As already indicated in Section IV-B, topologies have much influence on the performance of certain tracking systems. We work out the PCRLB of different topologies shown in Figure 6, where *Qt* is 0.1 and *Rt* is 0.05. As shown in Figure 7, we could draw the following conclusions:

**Figure 6.** Topologies configuration. ◦ represents free joints, and represents selected anchors.

**Figure 7.** Root mean square of PCRLB as a function of index of steps. Different receiver topologies are considered.

#### **7. Practical Use Case**

In above sections, the theoretical effectiveness of our proposed human motion tracking method is verified in terms of both spatial and temporal performances. The introduced IMU/TOA fusion method utilizes both distance and IMU information to lower the expected error bounds in theory. In this section, we take advantage of this method in 3D practical application scenarios to verify its performance.

#### *7.1. Platform Overview and Experiment Setup*

A minimized wearable sensing platform was specially designed for the tracking of human motion in 3D scenario. It aimed to collect spatial and temporal information during the human motion process. Sensing nodes were designed and intended to be put on joints to capture the movement conditions. Target data included accelerated velocity, angular velocity and the distances between joints and body parts. Among these, distance information was especially different when compared with other platforms, such as Xsens and Invensense. We sampled data using above-mentioned platform at 10 Hz. All processing work were performed online in Matlab with PC (Intel Core i7-6700M CPU, 16GB RAM). The communication between wearable nodes and the PC was via Bluetooth.

Our integrated wearable sensor system was composed of two parts: one control unit and several data acquisition units. The control unit worked as a gateway to control the whole operation process via Bluetooth communication. Data acquisition unit was mainly responsible for data sensing. Each data acquisition unit had a six-axis sensor (MPU6050, which integrated a triaxial accelerometer and a triaxial gyroscope), a barometer sensor (MS5611) and a UWB-TOA ranging module (DWM1000) [19]. The MEMS sensors were connected to a micro-controller (STM32F103) for the sake of sampling efficiency at a rate of 10 Hz. Data were transferred to the control unit in real-time and also written into its SD card as backups for other offline analysis and applications. Three sensing nodes were mounted onto the ankle, knee and hip joint. Distance measurements were conducted between these nodes. The whole system architecture is demonstrated in Figure 8.

**Figure 8.** Experimental Platform Settings. Three sensing nodes were mounted onto the ankle, knee and hip joint. Each node mainly consisted of micro-controller unit (STM32), IMU unit (MPU6050 and MS5611) and TOA unit (UWB module).

We proposed an IMU/TOA fusion method for human motion tracking applications, introducing distance information into inertial systems. However, one of the most obvious difference between the proposed method and that in [23] is: the algorithm in [23] needs to take advantage of external UWB stations as reference nodes, namely anchors. However, in our method, we simply utilized internal tags (mounted onto body joints) to range and conduct the fusion process, which can significantly reduce the complexity of tracking systems.

Based on above considerations, human lower limb movement analysis in 3D scenario was chosen as a typical verification experiment, for better demonstrate the performance of proposed method when it is applied to human motion tracking.

#### *7.2. Practical Use Case in 3D Scenario*

A spiral-stairs (between two floors) scenario Was selected as a 3D use case and four UWB anchors were located at each floor (totally eight anchors were used), as shown in Figure 9. The location of anchors in each floor are referred to those in [23] and only used in comparative experiments. Three comparative algorithms were also tested:

Experiment results are shown in Figures 9 and 10.

For better comparison, ground truth and tracking trajectories applying all mentioned methods (the proposed method and three comparative ones) are drawn in Figure 9. It is clearly seen that, when applying the proposed method, the experiment results were far closer to the ground truth, while the results when applying methods of IMU [51], IMU/ex-TOA [23] and IMU-OEKF [52] were drifting away as time accumulates. The result remained similar at the very beginning; however, the gap became larger over time.

**Figure 9.** Walking trajectory when climbing a spiral-stairs with our proposed method and two comparative methods.

**Figure 10.** The accumulative distribution function (CDF) curve in 3D scenario when algorithms were applied.

Detailed accumulated errors are shown in Figure 10. IMU, IMU/ex-TOA and IMU-OEKF faced the problem of drift error, especially when the experimenter turned around at stair corners. Time accumulated errors were more likely to be caused when moving direction changes. With IMU/ex-TOA and IMU-OEKF methods applied into the algorithm, localization errors to some extent could be fixed, but still existed and were critical at turns. On the contrary, our proposed human motion tracking method showed a good stability and no clear drift errors were observed in our experiment. This is because, with the use of proposed method, body connections are considered and temporal transitions are estimated. Especially at turns, body constraints could to some extent fix time accumulated errors. The results suggest that our proposed method performed well while applied in human lower limb analysis. This also confirmed the conclusion obtained in theoretical simulation (Sections 4 and 5). Methods with only inertial sensors faced the problem of drift errors, but our proposed model could fix time accumulative errors.

From above analysis, our proposed method had significantly higher accuracy, as well as little drift problem. Besides, compared with IMU/ex-TOA method [23] and IMU-OEKF method [52], our method did not need external anchors and had higher tracking accuracy, thus is more suitable for wearable motion tracking applications.

#### **8. Conclusions**

A IMU/TOA fusion-based human motion tracking method is proposed in this paper. IMU based HMT technique faces the tough problem of accumulative errors and drift. Time of Arrival (TOA) is considered to compensate the drift and accumulation caused by IMU. On this basis, Cramér–Rao lower bound (CRLB) is derived in detail with consideration of both temporal and spatial related factors. Simulation results show that the IMU/TOA fusion-based HMT method proposed in this paper has significantly higher accuracy, as well as little drift problem, compared with the traditional independent IMU or TOA tracking methods. The comprehensive fundamental limits analysis conducted in this paper can provide a theoretical basis for further system design and algorithm analysis. In the practical use case, proposed IMU/TOA fusion method also shows obvious performance advantages when compared with state-of-the-art methods. Besides, it does not need external anchors, thus it is more suitable for wearable motion tracking applications.

**Author Contributions:** C.X. and J.H. conceived and designed the experiments; C.X. and X.Z. (Xinghang Zhou) analyzed the data and contributed to the writing the manuscript. X.Z. (Xiaotong Zhang) and S.D. reviewed the analyzed data and contributed to the writing. C.X. and J.H. are co-first authors and contributed equally to this work.

**Funding:** This work was supported by The National Key R & D Program of China, No. 2016YFC0901303; National Natural Science Foundation of China (NSFC) project No. 61671056, No. 61302065, No. 61304257, and No. 61402033; Beijing Natural Science Foundation project No. 4152036; and Tianjin Special Program for Science and Technology No. 16ZXCXSF00150.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Derivation of the CRLB**

CRLB is represented as a theoretical lower limit for any unbiased estimation and is widely used to assess localization performance. Thus, we comprehensively derive the CRLB for 3D localization of IMU/TOA fusion method in WSNs to evaluate its spatial performance. The following are some definitions. If *pk* is an unbiased estimate of *pk*, then

$$E(\mathfrak{p}\_k - p\_k)^2 \ge \mathcal{C}RLB = \text{tr}\left\{ I(p\_k)^{-1} \right\} \tag{A1}$$

where *J*(*pk*) is the Fisher information matrix [43,46]. Before solving the Fischer information matrix, we need to first define the joint probability density function as

$$\begin{split} &p(\hat{d}\_{k},\hat{\alpha}\_{k},\hat{\beta}\_{k},\hat{p}\_{k})\\ &= \left\{ \prod\_{n=1}^{N} p(\hat{d}\_{n\boldsymbol{k}} | p\_{\boldsymbol{k}}) \right\} p(\hat{\alpha}\_{k} | p\_{\boldsymbol{k}-1\boldsymbol{\nu}},p\_{\boldsymbol{k}}) p(\hat{\beta}\_{\boldsymbol{k}} | p\_{\boldsymbol{k}-1\boldsymbol{\nu}},p\_{\boldsymbol{k}}) \end{split} \tag{A2}$$

where *p*(*α*ˆ *k*|*pk*−1, *pk*), *p*(*β*ˆ *k*|*pk*−1, *pk*) and *p*( ˆ*dn*,*k*|*pk*) can be obtained according to Equations (1)–(3).

According to the joint probability density function, we can define the Fischer information matrix as

$$f(p\_k)\_{i,j} = -E\left[\frac{\partial^2 \ln p(\hat{d}\_{k\prime}\hat{u}\_{k\prime}\hat{\rho}\_{k\prime}\hat{\rho}\_k)}{\partial p\_{k\prime}\partial p\_{k,j}}\right], i, j = 1, 2, 3\tag{A3}$$

Bringing Equation (A2) into Equation (A3) can calculate each element of the Fischer Information Matrix, and further obtain CRLB.

#### **Appendix B. Derivation of the PCRLB**

Sequential tracking problem is a temporal as well as spatial problem. This continuous information could be used to evaluate the performance of given algorithms. Thus, we extend the above CRLB to PCRLB by considering posterior information.

Before the derivation, we redefine the joint probability density function as

$$\begin{aligned} &p(\hat{d}, \pounds, \hat{\mathfrak{p}}, \hat{\mathfrak{p}})\\ &= p(\hat{d}\_0 | \hat{\mathfrak{p}}\_0) \prod\_{k=1}^{K} p(\pounds\_k | p\_{k-1}, p\_k) p(\hat{\mathfrak{p}}\_k | p\_{k-1}, p\_k) p(\hat{d}\_k | p\_k) \end{aligned} \tag{A4}$$

To calculate the Fischer information matrix at state *k*, we define

$$p\_k = p(\hat{d}\_{0:k\prime} \land\_{0:k\prime} \hat{\rho}\_{0:k\prime} \not p\_{0:k\prime}) \tag{A5}$$

where ˆ*d*0:*k*, *α*ˆ 0:*k*, *β*ˆ 0:*k*, and *p*ˆ0:*k* represent ranging, vertical angle, horizontal angle and target coordinate vector from the start state to state *k*, respectively.

Therefore,

$$\begin{aligned} &= \begin{bmatrix} J(p\_{0:k}) \\ E \end{bmatrix} = \begin{bmatrix} E \begin{Bmatrix} -\Delta\_{p\_{0:k-1}}^{p\_{0:k-1}} \ln p\_k \\ E \end{Bmatrix} &= \begin{Bmatrix} -\Delta\_{p\_{0:k-1}}^{p\_{0:k}} \ln p\_k \\ -\Delta\_{p\_{0:k}}^{p\_{0:k}} \ln p\_k \end{Bmatrix} \\ &= \begin{bmatrix} A\_k & B\_k \\ B\_k^T & C\_k \end{bmatrix} \end{aligned} \tag{A6}$$

According to Tichavsky et al. [43], the sub-matrix *Jk* can be obtained by pseudo-inverse of the matrix *J*(*p*0:*k*), i.e.,

$$J\_k = \mathbb{C}\_k - B\_k^T A\_k^{-1} B\_k \tag{A7}$$

According to Equations (A5) and (A6), the joint probability density for the *k* + 1 state is

$$p\_{k+1} = p\_k p(\pounds\_{k+1} | p\_{k\prime} p\_{k+1}) p(\hat{\beta}\_{k+1} | p\_{k\prime} p\_{k+1}) p(\hat{d}\_{k+1} | p\_{k\prime} p\_{k+1}) \tag{A8}$$

According to the joint probability density of state *k* + 1, we can find that

$$J(p\_{0:k+1}) = \begin{bmatrix} A\_k & B\_k & 0 \\ B\_k^T & \mathbb{C}\_k + H\_k^{11} & H\_k^{12} \\ 0 & H\_k^{12} & \gamma\_{k+1} + H\_k^{12} \end{bmatrix} \tag{A9}$$

where *H*11 *k* , *H*12 *k* , and *H*22 *k* reflect the posterior information from state *k* to state *k* + 1, and *γk*+1 reflects the location information based on TOA ranging [19].

From *J*(*p*0:*k*+1) and *Jk*, we can get the Fischer information matrix for state *k* + 1, i.e.,

$$\begin{aligned} &f\_{k+1} \\ &= \gamma\_{k+1} + H\_k^{22} - \begin{bmatrix} 0 & H\_k^{12} \end{bmatrix} \begin{bmatrix} A\_k & B\_k \\ B\_k^T & C\_k + H\_k^{11} \end{bmatrix} \begin{bmatrix} 0 \\ H\_k^{12} \end{bmatrix} \\ &= \gamma\_{k+1} + H\_k^{22} - H\_k^{12} (f\_k + H\_k^{11})^{-1} H\_k^{12} \end{aligned} \tag{A10}$$

Due to the step error and the directional error obey Gaussian distribution, *H*11 *k* = *H*12 *k* = *H*22 *k* = *Hk* can be calculated.

The solution of *Hk* can be found in [53]. In summary, the posterior Fisher information matrix is

$$J\_{k+1} = \gamma\_{k+1} + H\_k - H\_k (f\_k + H\_k)^{-1} H\_k \tag{A11}$$

According to the SMW (Sherman–Morrison–Woodbury) formula [54], it can be further simplified as

$$J\_{k+1} = \gamma\_{k+1} + (H\_k^{-1} + J\_k^{-1})^{-1} \tag{A12}$$

where *γk*+1 reflects the information based on TOA, *Hk* reflects information based on IMU.

Then, sequential simulation can be conducted by Monte Carlo methods [50].

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Implementation and Assessment of an Intelligent Motor Tele-Rehabilitation Platform**

**Yves Rybarczyk 1,2,3,\*, Jorge Luis Pérez Medina 1, Louis Leconte 4, Karina Jimenes 1, Mario González 1 and Danilo Esparza 1**

Received: 18 October 2018; Accepted: 24 December 2018; Published: 4 January 2019

**Abstract:** Over the past few years, software applications for medical assistance, including tele-rehabilitation, have known an increasing presence in the health arena. Despite the several therapeutic and economic advantages of this new paradigm, it is important to follow certain guidelines, in order to build a safe, useful, scalable, and ergonomic tool. This work proposes to address all these points, through the case study of a physical tele-rehabilitation platform for patients after hip replacement surgery. The scalability and versatility of the system is handled by the implementation of a modular architecture. The safeness and effectiveness of the tool is ensured by an artificial intelligence module that assesses the quality of the movements performed by the user. The usability of the application is evaluated by a cognitive walkthrough method. Results show that the system (i) is able to properly assess the correctness of the human's motion through two possible methods (Dynamic Time Warping and Hidden Markov Model), and (ii) provides a good user experience. The discussion addresses (i) the advantages and disadvantages of the main approaches for a gesture recognition of therapeutic movements, and (ii) critical aspects to provide the patient with the best usability of a tele-rehabilitation platform.

**Keywords:** eHealth; software engineering; gesture recognition; Dynamic Time Warping; Hidden Markov Model; usability

#### **1. Introduction**

Over the past few years, alternative healthcare deliveries have been developed. For instance, the advances in the telecommunications have permitted the emergence of telehealth systems. One specialized field of telehealth is tele-rehabilitation, which allows for the implementation of a therapeutic program via an interactive multimedia web-based platform [1,2]. A remarkable increase in the number of patients treated by tele-rehabilitation has been noticed since the beginning of the 21st century, especially in physiotherapy [3]. Despite the several medical and economic advantages of this new paradigm, the development of a tele-rehabilitation platform has to follow specific rules to be safe and efficient. Two fundamental rules are (i) to build a platform based on an affordable system to capture the human movements and (ii) to make sure the patient is performing the therapeutic exercises correctly [4].

The present work proposes a web-based platform for motor tele-rehabilitation applied to patients after hip arthroplasty surgery. This orthopedic procedure is an excellent case study, because it involves people who need a postoperative functional rehabilitation program to recover strength and joint mobility. The development of a tele-rehabilitation system is justified by the condition of these individuals that makes difficult their transportation to and from the physiotherapist's office. The proposed approach considers two fundamental conditions for the development of a suitable tele-rehabilitation platform. First, the system is based on a modular architecture composed of a low-cost motion capture device, in order to ensure the viability and scalability of the tool. Second, the platform detects automatically the correctness of the executed movement to provide the patient with real-time feedback [5]. Since exercises are carried out in front of a machine instead of a therapist, the subjects may produce incorrect movements that could be harmful, especially after a surgery. In addition, the lack of human presence diminishes the motivation, which may slow down the recovery process. Here, two possible approaches are tested to assess the movements and display an appropriate feedback: (i) Dynamic Time Warping (DTW) [6], and (ii) Hidden Markov Models (HMMs) [5].

The remainder of the manuscript is organized into six parts. Section 2 is a presentation of related work. Section 3 is a general description of the web-based platform. Section 4 focuses on a first study to assess the movements and it is based on DTW. Section 5 presents a second study on movement assessments, which uses HMMs. The reliability of both approaches to discriminate between correctly and incorrectly performed movements is tested through laboratory experiments that compare the evaluation of the computer to the evaluation of the therapists. Section 6 consists of an analysis of the usability of the system, through an experiment that involves end users. Finally, the last section (i) discusses the results of the automatic assessments and the usability tests, and (ii) draws conclusions on the most suitable solutions to develop an efficient tele-rehabilitation system.

#### **2. Related Work**

Different approaches are proposed to provide the patient with a tele-rehabilitation system. The motion capture can be based on inertial wearable sensors or visual sensors. [7] propose a cloud-assisted wearable system (Rehab-aaService) that enables a general motor rehabilitation, even if the application is optimized for the upper limbs. The platform is scalable and can be integrated into a body sensor network (BodyCloud) [8] for the monitoring of different physiological parameters. However, the system does not have an artificial intelligent module that allows for a rigorous assessment of the rehabilitation exercises. In addition, the wearable devices require the inertial sensors to be precisely placed on the body and/or to involve a calibration stage. Thus, another approach consists of using a vision-based motion capture. This system also presents certain limitations, such as the occultation problem, but it has the advantage to provide an easy setup. The occultation refers to a configuration in which a body joint disappears temporarily behind a part of the individual's body. This situation can be overcome (at least partially) by asking the individuals to change the orientation of their body according to the plane in which the movement is performed. Most of the systems use the Kinect, because it is an affordable piece of equipment and its accuracy is good enough for functional assessment activities [9]. An experiment that consists of assessing the accuracy of the Kinect by comparison to an accelerometer shows a percentage of correlation between the two sensors equal to 96% [4]. An avatar evolving in a serious game is usually used to motivate the users to regularly practice their therapeutic exercises [10]. This gamification approach is also applied for physical exercises, such as Tai Chi [11]. The application can indicate how well the player imitates the instructor, even if the system is not designed to classify gestures. Another type of Tai Chi prototype platform is built to rehabilitate patients with movement disorders [12]. The patient's movements are compared with pre-recorded movements of a coach and further evaluated by using a fuzzy logic algorithm [13]. Additionally, a component-based application framework is proposed by [14] for the development of 3D games by combining already existing 3D visual components. The preliminary results indicate that the prototypes can be used as serious game for physical therapy. However, this framework is not open source and the matching between the avatar and the user is not robust, especially when occultations occur.

In general, it is very rare to find a system that provides both an attractive virtual environment and an algorithm to assess the quality of the rehabilitation movements executed by the patients. Instead of evaluating a spontaneous movement, other studies propose to guide the gesture through visual [15] or haptic [16] feedback, in order to avoid wrong motions. The disadvantage of these approaches is to induce a too stereotyped movement, which reduces the functional benefit of the rehabilitation. Thus, the technological implementation of an appropriate program of motor rehabilitation must involve an expert system that could substitute the therapist and provide the patient with feedback. Recent studies applied algorithms based on rules [17] or Dynamic Time Warping [18] for the recognition of therapeutic movements, and fuzzy logic [19] for the diagnosis of physical impairments. Nevertheless, the current systems are essentially able to discriminate between a correct and an incorrect gesture, but they do not give a targeted feedback on the type of error or they do not consider compensatory movements when an exercise is wrongly executed. To get such a feedback, it seems necessary to build a model of the therapeutic exercises as suggested by [20], who propose a theoretical modeling in UML for the reeducation of the upper limbs. Our work proposes a more advanced approach, since we developed and implemented a statistical model to assess the rehabilitation movements, which can be applied on any part of the body (upper and lower limbs) and can precisely identify the cause of a bad performance to provide the users of the platform with comprehensive feedback regarding the corrections to be made to the gesture.

#### **3. Functionalities and Architecture of the Platform**

As represented in Figure 1, the platform is composed of two main blocks dedicated to (i) the Core and (ii) the Real Time Evaluation.

**Figure 1.** Architecture of the tele-rehabilitation platform.

#### *3.1. Core*

The core block allows for empowering patients, who can perform autonomously a set of rehabilitation activities from a program elaborated by health professionals (e.g., physiotherapists) [21]. Several functionalities are implemented in the platform, in order to support patients in their self-recovery process. They have access to a set of multi-modal learning resources on the appropriate use of the platform and the correct procedure to complete the therapeutic program. The patient can initially visualize a list of learning resources, select one of them and display an additional list of associated multimedia files. Then, it is possible to choose the preferred learning resources among different modalities, review the information and return to the general learning menu.

In addition, this block is connected to the database that contains a set of rehabilitation exercises. The computational implementation is based on Django 2.0.5. This web framework uses the Model-View-Controller (MVC) design pattern and works under the two central principles: (i) maintaining loose coupling between the layers of the framework, and (ii) Don't Repeat Yourself (DRY) [22]. According to this pattern, the application functionality is divided into three kinds of components. The model represents the knowledge and has the logic to update the controller. The view represents the visualization of the data. The controller determines the data flow into a model object and updates the view whenever data change. It keeps view and model separated and acts on both.

#### *3.2. Real Time Evaluation*

#### 3.2.1. Functionalities

The rehabilitation sessions are planned by the physiotherapist. They are composed of a set of exercises. An exercise is based on the composition pattern [23], which means that it is a composition of simple or compound exercises. Before performing the exercises, patients have to answer a questionnaire regarding their health status. The responses are processed by the application to authorize or disavow them to start a rehabilitation session. If they are authorized to initiate a session, the system activates the execution interface of the exercises. An exercise consists of a set of repetitions predefined in the recovery plan. There is a short break between each repetition to provide the patient with a feedback regarding the quality of the movement. The trials are assessed in terms of range of motion (ROM) and compensatory movements. Before starting an exercise, the patients can consult a set of instructions to execute the movements correctly. During the completion of the exercise, an avatar that maps in real time the patient's movement is displayed in a game-based virtual environment (Figure 2). The movement quality and facial expressions are also assessed in real time by artificial intelligence algorithms. Once a trial is completed, the performance is graphically presented to the user. Finally, the patients have the possibility to suspend a rehabilitation session at any time. In this case, they must provide the physiotherapist with a reason to abort the session.

**Figure 2.** User's interface of the game-based exercises.

#### 3.2.2. Implementation

The real time evaluation block is composed of two modules, the first one for processing skeleton data with a RESTful API, and the second one for acquiring movement data as a client application. They communicate to each other for processing the movement data and for giving the results using JSON data format. The RESTful API was developed using Django Rest Framework 3.8.2. It is a

powerful and flexible toolkit of Django to build RESTful Web Services [24]. The client application is built with the Kinectron 1.4.2 and Processing 1.4.8. Kinectron is an open source tool to capture the movement data. It has two components: a server that broadcasts Kinect data, and an API that receives Kinect data [25]. A Kinect v2 camera is used to capture the movements of the patient and extracts the coordinates of the main body joints. A cross-browser JavaScript Library is incorporated (Three.js) to display in a web browser the game-based exercises, which are composed of 3D animated objects and the user's avatar. Additional technologies, such as Boostrap 4.1.13 and JQuery 3.1.14 are used to build the interface. Boostrap is a CSS framework that provides a characteristic of responsive web design. It allows developing interfaces compatible with different kinds of browsers and terminal platforms [26]. JQuery is a Javascript library which has CSS selector, flexible animation system and event system, rich plugins, and solutions for browser compatibility issues.

The API REST package works with the exercises packages and contains the intelligent modules to assess both the quality of the movements and the emotions of the patients when performing an exercise. It has a couple of services to process the skeleton and the images frames through different artificial intelligence algorithms. The first service processes the skeleton to evaluate if the patient is exercising correctly and the result is transferred as a JSON structure with a set of angles. Each angle has a type, a DTW metric that compares the performed movement to the reference movement. A graphical-based feedback is displayed on the user's interface to inform the patients on the quality of their performance (Figure 2). A red color represents bad movement amplitude for a determined limb or body part. This visual feedback is reinforced by a textual description of potential errors in the motion. Other information, such as a progress bar and a score are also displayed. The second service processes the image frame every second to determine if the patient feels pain. A machine-learning algorithm based on Support Vector Machine (SVM) is implemented to operate a Boolean classification of facial expressions in pain (1) vs. no pain (0). A sequence of 1 and 0 is obtained at the end of each repetition. This information is used by the therapist to identify if the intensity of the exercises is appropriate or not for a determined patient.

#### **4. Movement Assessment: DTW Approach**

#### *4.1. Dynamic Time Warping*

Dynamic Time Warping (DTW) is an algorithm for measuring similarity between two temporal sequences, which may vary in speed. It is perfectly adapted for the assessment of therapeutic movements for which the velocity of the execution is not relevant. This characteristic explains why we have chosen this technique, because it enables us to assess the quality of the movement by a single comparison to a reference. Any distance (Euclidean, Manhattan, ... ) which aligns the *i*-th point on one-time series with the *i*-th point of another produces a poor similarity score. On the contrary, a non-linear alignment produces a more intuitive similarity measure, allowing similar shapes to match even if they are out of phase in the time axis. The implementation of such an algorithm is obtained through a matrix sized with the length of each signal (Figure 3). On each point of the matrix, the distance between the points associated with each of the two signals is calculated. The best alignment of two signals is given by the path through the grid that minimizes the total distance between these two signals. Our algorithm implements a quadratic distance (*D*), as defined in Equation (1).

$$D\_{\langle i,j \rangle} = \left(\mathbf{x}\_{\langle j \rangle} - \mathbf{y}\_{\langle i \rangle}\right)^2 \tag{1}$$

Thus, a distance matrix is obtained. Then, a cumulated distance matrix (*C*) is created, based on the distance matrix (*D*). The latter is built thanks to a definition of distance very close to the 8-connectivity in a square tiling (mathematic topology) for each point of the matrix, as presented in Equation (2).

$$\mathbb{C}\_{\{i,j\}} = \min \left\{ \mathbb{C}\_{\{i-1,\,j-1\}}, \mathbb{C}\_{\{i-1,j\}}, \mathbb{C}\_{\{i,j-1\}} \right\} + D\_{\{i,j\}} \tag{2}$$

Finally, we have to look for the path that minimizes the cumulated distance. To do so, we start from the last point and look for the point with the lower cumulated distance between all its three inferior neighbors (within the meaning of Moore neighborhood). We reiterate the process until we reach a side or the originate point *C*(0, 0). Then, we just have to calculate the length of this path. The result is small when signals are close, and becomes larger when signals differ. An optimization can be provided to this algorithm, using a warping window. It is based on the fact that a good alignment path is unlikely to wander too far from the diagonal and, thus, limits the window length. Since the Kinect frequency is about 33 Hz and the movements are short, we do not need to implement this optimization, because only few points (between 100 and 200) have to be analyzed.

**Figure 3.** Cumulated distance matrix and its optimal path (in red), between a reference movement and a movement to be tested (blue signals).

#### *4.2. Trigonometric Parametrization*

As we expect our algorithm to be invariant to 3D translational motion and to different bodies, we have to focus not on coordinates of each joint, but on angles. This section describes the mathematical calculation of the (i) working angles (ROM) and (ii) compensation angles for the four main therapeutic movements: hip abduction, slow flexion hip and knee, hip extension, and forward-sideway-backward sequence. During a rehabilitation session it is important to identify both the movement amplitude the patient should reach (ROM) and the compensatory movements that should be avoid for an appropriate re-education of the injured limbs.

For hip abduction, we identified one working angle (Figure 4, picture on the left) and three compensation angles (Figure 4, three pictures on the right). To obtain the angle of the hip abduction (θ*aw*), we have to project −−−→ 16\_18 on the frontal plane, and proceed with the calculation of Equations (3)–(5).

$$m\_1 = \left[\sum\_{k=x,y,z} (k\_{18} - k\_{16})(k\_{16} - k\_0)\right]^2 \tag{3}$$

$$d\_1 = (\mathbf{x}\_{16} - \mathbf{x}\_0)^2 + (y\_{16} - y\_0)^2 + (z\_{16} - z\_0)^2 \tag{4}$$

*Electronics* **2019**, *8*, 58

$$\cos\Theta\_{4w} = \frac{|y\_{16} - y\_{18}|}{\sqrt{\frac{\eta\_1}{d\_1} + \left(y\_{16} - y\_{18}\right)^2}}\tag{5}$$

where *n*1 and *d*1 are intermediate variables used to calculate the angle θ*aw*, and k is successively equal to *x*, *y*, and *z* in the calculation of the sum. '0' stands for spine base, '16' stands for left hip, and '18' stands for left ankle. The leg compensation angle (θ*lc*) is obtained by the calculations of Equations (6)–(8).

$$m\_2 = \left[\sum\_{k=x,z} (k\_{18} - k\_{16})(k\_{15} - k\_{14})\right]^2 \tag{6}$$

$$d\_2 = \left(\mathbf{x}\_{15} - \mathbf{x}\_{14}\right)^2 + \left(z\_{15} - z\_{14}\right)^2\tag{7}$$

$$\cos\theta\_{l\mathfrak{c}} = \frac{|y\_{16} - y\_{18}|}{\sqrt{\frac{n\_2}{d\_2} + \left(y\_{16} - y\_{18}\right)^2}}\tag{8}$$

where *n*2 and *d*2 are intermediate variables used to calculate the angle θ*lc*. '14' stands for right ankle, and '15' stands for left foot. Equations (9)–(11) are used to calculate the frontal torso compensation angle (θ*ftc*).

$$m\_3 = \left[\sum\_{k=x,z} (k\_0 - k\_{20})(k\_{15} - k\_{14})\right]^2 \tag{9}$$

$$d\_3 = \left(\mathbf{x}\_{15} - \mathbf{x}\_{14}\right)^2 + \left(z\_{15} - z\_{14}\right)^2\tag{10}$$

$$\cos\theta\_{ftc} = \frac{|y\_0 - y\_{20}|}{\sqrt{\frac{y\_3}{d\_3} + \left(y\_0 - y\_{20}\right)^2}}\tag{11}$$

where *n*3 and *d*3 are intermediate variables used to calculate the angle θ*ftc*. '20' stands for spine shoulder. Equations (12)–(14) are used to calculate the lateral torso compensation angle (θ*ltc*).

$$m\_4 = \left[\sum\_{k=x,y,z} (k\_{20} - k\_0)(k\_{16} - k\_0)\right]^2 \tag{12}$$

$$d\_4 = (x\_{16} - x\_0)^2 + (y\_{16} - y\_0)^2 + (z\_{16} - z\_0)^2 \tag{13}$$

$$\cos\theta\_{ltc} = \frac{|y\_0 - y\_{20}|}{\sqrt{\frac{\mathbf{u}\_4}{d\_4} + \left(y\_0 - y\_{20}\right)^2}}\tag{14}$$

**Figure 4.** Expected and possible movements involved in hip abduction (from left to right): working angle, leg compensation, front torso compensation, and lateral torso compensation.

*Electronics* **2019**, *8*, 58

For slow flexion of hip and knee, we identified two working angles (Figure 5, picture on the left) and two typical compensation angles (Figure 5, two pictures on the right). Equations (15)–(17) are used to calculate the hip-working angle (θ*hw*).

$$n\_5 = \left[\sum\_{k=x,y,z} (k\_{17} - k\_{16})(k\_{15} - k\_{14})\right]^2 \tag{15}$$

$$d\_5 = \left(\mathbf{x}\_{15} - \mathbf{x}\_{14}\right)^2 + \left(y\_{15} - y\_{14}\right)^2 + \left(z\_{15} - z\_{14}\right)^2\tag{16}$$

$$\cos\theta\_{hw} = \frac{|y\_{17} - y\_{16}|}{\sqrt{\frac{y\_5}{d\_5} + \left(y\_{17} - y\_{16}\right)^2}}\tag{17}$$

where *n*5 and *d*5 are intermediate variables used to calculate the angle θ*hw*. '14' stands for right ankle, '15' stands for left foot, '16' stands for right hip, and '17' stands for left knee. Equations (18)–(20) are used to calculate the knee-working angle (θ*kw*).

$$m\_6 = \left[\sum\_{k=x,y,z} (k\_{18} - k\_{17})(k\_{15} - k\_{14})\right]^2 \tag{18}$$

$$d\_6 = (x\_{15} - x\_{14})^2 + (y\_{15} - y\_{14})^2 + (z\_{15} - z\_{14})^2 \tag{19}$$

$$\cos\theta\_{kw} = \frac{|y\_{18} - y\_{17}|}{\sqrt{\frac{n\_6}{d\_6} + (y\_{18} - y\_{17})^2}}\tag{20}$$

where *n*6 and *d*6 are intermediate variables used to calculate the angle θ*kw*. '14' stands for right ankle, '15' stands for left foot, '17' stands for left knee, and '18' stands for left ankle. Equations (21)–(23) are used to calculate the thigh compensation angle (θ*thc*).

$$m\_7 = \left[\sum\_{k=x,y,z} (k\_{17} - k\_{16})(k\_{12} - k\_{16})\right]^2 \tag{21}$$

$$d\tau = \left(\mathbf{x}\_{12} - \mathbf{x}\_{16}\right)^2 + \left(y\_{12} - y\_{16}\right)^2 + \left(z\_{12} - z\_{16}\right)^2\tag{22}$$

$$\cos\theta\_{thc} = \frac{|y\_{17} - y\_{16}|}{\sqrt{\frac{\eta\_7}{d\_7} + \left(y\_{17} - y\_{16}\right)^2}}\tag{23}$$

where *n*7 and *d*7 are intermediate variables used to calculate the angle θ*thc*. '12' stands for left hip, '16' stands for right hip, and '17' stands for left knee. Equations (24)–(26) are used to calculate the tibia compensation angle (θ*tic*).

$$m\_8 = \left[\sum\_{k=x,y,z} (k\_{18} - k\_{17})(k\_{12} - k\_{16})\right]^2 \tag{24}$$

$$d\_8 = (\mathbf{x}\_{12} - \mathbf{x}\_{16})^2 + (y\_{12} - y\_{16})^2 + (z\_{12} - z\_{16})^2 \tag{25}$$

$$\cos\theta\_{\rm{fic}} = \frac{|y\_{18} - y\_{17}|}{\sqrt{\frac{n\_8}{d\_8} + \left(y\_{18} - y\_{17}\right)^2}}\tag{26}$$

where *n*8 and *d*8 are intermediate variables used to calculate the angle θ*tic*. '12' stands for left hip, '16' stands for right hip, '17' stands for left knee, and '18' stands for left ankle. The same frontal and lateral compensation movements of the trunk as described in hip abduction are also calculated for this exercise (Equations (9)–(14)).

**Figure 5.** Expected and possible movements involved in slow flexion of hip and knee (from left to right): working angle for hip (in green) and knee (in blue), thigh compensation, and tibia compensation.

For hip extension, the working angle is the same as the compensation angle of hip abduction, and the compensation angle is equal to the working angle of hip abduction (Figure 6). The same frontal and lateral compensation movements of the trunk as described in hip abduction are also calculated for this exercise (Equations (9)–(14)).

**Figure 6.** Working angle for hip extension.

There are two working angles for the forward-sideway-backward sequence, which are the working and leg compensation angles of the hip abduction (Figure 7). Again, the frontal and lateral compensation movements of the trunk described in hip abduction are considered for this exercise (Equations (9)–(14)).

**Figure 7.** Working angle for forward–sideway–backward sequence.

#### *4.3. Filtering and Implementation*

Since the raw input data provided by the Kinect are too noisy, a low pass filter is applied to the signal (Figure 8). The Butterworth Filter available in the python library 'scipy.signal' is used to smooth the signal. It is necessary to create two different filters according to the noisiness of the data (a moderate and a strong one). The best ratio signal/noise for the moderate filter is obtained empirically with the parameters as follows:

The settings of the strong filter are:

**Figure 8.** Examples of the filtering effect for the moderate (**left panel**) and the strong (**right panel**) filter.

The assessment module is implemented to compare the DTW cost of each described angle (ROM and compensations) to an empirically defined threshold (the optimum value is defined in the next section). If the distance calculated by the DTW is higher than the threshold, the angle of the movement is classified as wrong. Since different DTW run in parallel (one for each specific working and compensation movements), it is possible to provide the patients with a detailed feedback of their errors for helping them to correct their performance in the next trial. Whatever the correctness of the movement, a graphical message is always displayed to the users at the end of the trials (see Figure 2). Once a session is finished, a JSON file is created with every angle data cost, pictographic feedback (good, bad, or neutral smileys), and six possible messages that are displayed on the patient's interface (evaluation of the ROM and potential compensation movements). Then, these data are stored in the platform database.

#### *4.4. Validation*

#### 4.4.1. Experimental Protocol

The experiment consisted of comparing the assessment performed by four physiotherapists to the assessment made by the algorithm. Seven healthy subjects participated in the study. They were informed about the purpose of the study and signed a consent form to take part to the experiment. They were asked to take place at approximately two meters distance from a single Kinect camera. The motion capture device was placed at the height of the subject's xiphoid apophysis. Each of the four movements

described in Section 4.2 was repeated eleven times (one of them was used as reference). The participant had to introduce a different error in the execution of the movement at each repetition (e.g., wrong range of motion or different kinds of compensation). To do so, they followed a script unknown to the therapists and different from one participant to another. The only similarity between the scripts is the fact that they were composed of six accurate executions, two movements with an incorrect ROM, and three movements with compensation errors. For each trial, the physiotherapists were asked to evaluate all the angles (ROM and compensations) on a 4-point Likert scale (0–3). The data were processed with the python library *Pandas.* Due to the inter-individual differences between the assessments made by the physiotherapists, the scale of the marks was reduced to two levels. All the assessments of 0 and 1 were grouped in a single category that corresponded to a bad execution of the exercise. And the assessments of 2 and 3 were all classified as a good movement. The final therapist's evaluation was obtained by averaging the score of the four professionals. In case of draw, the trial was discarded. Then, these Boolean variables (bad vs. good) were compared with the assessment performed by the algorithm. The outcome was expected to be different depending on the decision threshold.

#### 4.4.2. Results

Overall, the outcome of the algorithm matches at 88% with the assessment made by the physiotherapists. The accuracy of this matching depends on the threshold of tolerated difference between the reference movement and the analyzed movement (Figure 9). A certain variability in the performance of the classification is observed from one exercise to another. For instance, the best results are obtained for the assessment of the ROM for hip abduction and forward-sideway-backward sequence (percentage of accuracy >90%).

**Figure 9.** Overall variation of accuracy depending on the threshold.

A confusion matrix permits to identify the percentage of movements assessed as correct by the therapists and predicted as correct by the algorithm (True Positives—TP), the incorrect movements predicted as incorrect (True Negatives—TN), the correct movements predicted as incorrect (False Negatives—FN), and the incorrect movements predicted as correct (False Positives—FP). The accuracy of the algorithm to predict the TP and TN is almost the same, even if it is slight higher for the former than the latter (Table 1). This small difference can be explained by the fact that there is more diversity in the execution of the bad (errors of ROM and compensations) than the good movements.

**Table 1.** Confusion matrix providing a detailed comparison between the movement assessments of the physiotherapists and the algorithm.

The misclassified movements (FP and FN) can be explained by two main reasons. First, there is only a partial consensus between the assessments carried out by the four therapists. Nevertheless, the discrepancy between the health professionals is less than 8% and, consequently, the human judgment is considered as an acceptable reference to evaluate the algorithm. Second, some angles suffer from occultation and cannot be evaluated correctly by the computer. This situation happens occasionally when a body joint (e.g., the hip) disappears temporarily behind a part of the participant's body (e.g., the hand).

Table 2 provides a quantification of the missing data in terms of the percentage of sampling errors and occultations. Sampling errors are defined as technical issues caused by faults in the capture sampling rate of the Kinect. The occultations refer to situations in which the participants cross their arms, cross their arms in front of the torso or cross the legs while performing the exercises. They are computed considering the pose arrangement [27] that identifies intersections between the torso, the arms, and/or the legs. The results show that the sampling errors are almost insignificant (<0.01%) and the occultations represent less than 6% of the missing data. Some movements are more affected by occultation than others, which suggests defining an optimal orientation of the body with respect to the vision-based motion capture sensor, for each exercise.

**Table 2.** Missing data defined in terms of the percentage of sampling errors and joint occultations, for each exercise (slow flexion hip and knee—SFHK; hip abduction—HA; hip extension—HE; forward-sideway-backward—FSB).

#### **5. Movement Assessment: HMM Approach**

#### *5.1. Hidden Markov Model*

Hidden Markov Model is a probabilistic approach that aims to model a given action into hidden states. It is composed of initial, transitional and emission probabilities. Initial probabilities are the distribution of probabilities of 'being in a state' before a sequence is observed. Transitional probabilities are represented by a matrix, in which the probabilities indicate the possible changes from one state to another. Finally, the emission probabilities model the variance of each state's associated values (mostly Gaussian Probability Density Functions—PDFs) obtained from continuous variable observations. These model parameters can be learned with the use of the Expectation-Maximization (EM) algorithm. Signals can be classified by looking at the probability that a signal is generated from a trained HMM. This is done by the use of the forward algorithm.

Learning the model parameters (states and transitions) by optimizing the likelihood is essential to make meaningful use of the HMM in classification. The distribution function defined by a Gaussian, Mixed Gaussian or multinomial density function, as well as the covariance type, needs to be characterized prior to this process. An observation is merely a noisy and variable representation of a related state. A state is a clustering of observations that relates to a distribution with a specific mean in the parameter space. A likely state is retrieved by finding the cluster that the observation is member of. Also, the transition probabilities between states creates a sequence of the most likely temporal succession of states. Estimating the model parameters is done by utilizing the Baum-Welch Expectation-Maximization algorithm, which is based on a forward-backward algorithm used in classifying Hidden Markov Chains [28,29]. The probabilities are calculated at any point of a sequence by inspecting previous observations, to find out how well the model describes the data, and following observations, to conclude how well the model predicts the rest of the sequence. This is an iterative

process, in which the objective is to find an optimal solution (state sequence) for the HMM. This optimal sequence of states is inferred using the Viterbi algorithm. In addition, the forward algorithm can be used to calculate the probability that a sequence is generated by a specific trained HMM, making it applicable for classification.

This classification is based on training an individual HMM per subclass of an exercise. For instance, one HMM could be trained on 'running' while another one would be trained on 'walking' (both subclasses of the human locomotion class). When calculating the forward probabilities of a sequence of observations and comparing the probabilities of all the HMMs, the sequence is classified as the category that provides the highest probability, as described in Equation (27).

$$\text{Class} = \arg\_{i=1}^{n} \max \left[ \Pr(O|\lambda\_i) \* \Pr(\lambda\_i) \right] \tag{27}$$

where *λi* represents a determined model and *O* is a sequence of observations. The amount of states in a HMM is a free parameter. The Bayesian Information Criteria (BIC) is a technique that aids to define a determined number of parameters by taking into account the possibility of overfitting the data when the number of states increases. BIC penalizes HMMs that have a high number of states, as described in Equation (28).

$$\text{BIC} = \ln(n)\text{s} - 2\ln(MLE) \tag{28}$$

where *n* is the data size and *s* the amount of states. Therefore, the optimal amount of states is retrieved by selecting the model with the lowest BIC score. HMMs trained with multiple states are evaluated by cross-validation on their Maximum Likelihood Estimation (MLE) and the previously mentioned penalizing term.

#### *5.2. Gesture Representation*

A skeletonized 3D image from a Kinect camera provides Cartesian *x*, *y* and *z* coordinates of twenty body joints. The gesture representation is chosen to be a skeletonized image as this has been shown to improve the model accuracy [30]. This representation depends on the position of the subject in relation to the camera and the roll, yaw and pitch angles of the device. The causal relationships between different joints are not captured by this representation. This means that physical constrains, such as a movement of the ankle that could be influenced by bending the knee, are not accounted for. To overcome these limitations, the joints are used to create a new representation that contains angles of multiple joints in respect to the frontal and sagittal planes, as well as multiple angles between relevant limbs. Figure 10 shows a graphical representation of the features in relation to the skeleton image. Table 3 describes the feature vector of the joint movements according to the anatomical terminology.

**Figure 10.** Graphical representation of the features used in the study. The movements in the egocentric frontal plane and sagittal plane are represented in green and red, respectively. The purple arrow represents the angle of the knee (independent from any plane).

**Table 3.** Feature vector describing the joint movements. For the purpose of the assessment these features are also transformed into speed and acceleration.

In this study, the motion is defined from the following joints: ankles, knees, hips, and spine. The angles of the knees are obtained by calculating the angle between ankle, knee and hip. The orientation of the knees induced by hip activity are expressed in four angular representations, following the two opposite directions for both sagittal and frontal planes. The same method is applied to describe the orientation of the torso, by finding the displacement of the center between the two shoulders in relation to the hips. This leads to a description of the movement into fourteen features. It is the principal representation followed by a first order and second order derivatives of these features that provide speed and acceleration of the movement. Overall, a total amount of 42 features is used. Figure 11 represents a diagram of the HMM's implemented in this study, in which *State*(*t*) and *Observation*(*t*) are state id and associated feature values at t time, respectively.

**Figure 11.** Graphical representation of the HMM for the exercise assessment. stands for: takes value out of; and ∀ stands for: out of all. *F* and *S* are the collection of features (42) and states (13), respectively. The definition of the optimal number of states is explained in Section 5.3.4. Each state is dependent on its previous state and observations are samples of the associated current state.

#### *5.3. Experiment*

#### 5.3.1. Protocol

Four healthy subjects participated in the experiment. They were informed about the purpose of the study and signed a consent form to take part in the experiment. They were asked to stand at approximately two meters distance front of a Kinect camera. The motion capture device was placed at the height of the subject's xiphoid apophysis. Each participant executed 70 movements leading to a total of 280 records. The rehabilitation exercise was a sequence, in which the subjects had to do one step forward, one step sideways and one step backward, with variations. These variations are staged executions of errors or compensatory movements that can occur during the rehabilitation in practice. The exercise was performed in batches of ten and the experiment was divided into two independent parts. In the first part, the subjects had to execute the movements as follows: (I) correct execution, (II) steps too short, (III) execution without moving the center of mass, (IV) steps too large, (V) steps with bended knee, (VI) steps with bended knee and flexed torso. The objective of this first experimental part was to assess the accuracy of the models to classify the movements. In the second part, the subjects had to perform ten trials of a partially wrong executions of the exercise (VII), in which the faults II to VI are only occurring in the beginning, middle or end of the sequence. This second experiment is used to evaluate the real-time applicability of the HMM technique.

An application was created to capture the skeletonized image of the subjects performing an exercise. Python 2.7 was used to create a graphical user interface with the option to name, start and stop a recording. In addition, Python was used for the later processing steps, which were feature transformations and classifications. The application communicated with the Kinect SDK and wrote the data into a CSV file, with a frequency of 60 Hz, during the recording mode. The developed programs and raw data are available at http://docentes.fct.unl.pt/y-rybarczyk/files/programs.rar and http: //docentes.fct.unl.pt/y-rybarczyk/files/data.rar, respectively. The HMMLearn package for python 2.7 was used for the training and application of the HMMs (https://github.com/hmmlearn/hmmlearn).

#### 5.3.2. Evaluation Method

Using the BIC score to select the appropriate amount of states was done for each type of trained HMMs (I–VI). It provided insight on the semantic variation within the exercises. For instance, less states assigned to a faulty movement relative to the good execution implies that 'there was something missing in the execution', whereas the detection of extra states implied that 'there was something added to the movement'.

For each type of execution (I–VI) an HMM was trained, leading to a total of six distinct HMMs. In order to build a general model that could assess the movements of any subject, models that classified the executions of a subject were exclusively trained on the recordings of the other remaining subjects. This led to a total of 24 trained models (6 per subject). The initial model parameters were set with a Gaussian density function and a full covariance matrix type (the initializations of the model parameters were done randomly, HmmLearn standard initialization was used, and EM iterated a 1000 times unless log-likelihood gain was less than 0.01). The HMM topological structure was fully connected, because priory knowledge about the expected state sequence could not be estimated with sufficient certainty. In addition, as the outcome of six classifiers determined the most likely model that was associated with the sequence, unpredicted variance in a signal (or noise) should not drastically influence the likelihood of the signal. The outcome of the six classifiers was calculated by means of the forward probabilities. Then, these probabilities were ranked from 1 to 6, where 1 and 6 were assigned to the highest and the lowest probability, respectively. A confusion matrix was used to map these values in terms of average prediction rank of each type of execution. In addition, an indication of the similarity of a type of execution in relation to the combination of all the other executions was provided. Finally, a range of sliding temporal windows was used to evaluate the real-time suitability of the approach. It was applied to assess the correct detection of the present types of faults in executions VII. These windows classified a subsequence in a fixed number of samples, which partially overlapped over time.

#### 5.3.3. Validation Method

To get a reliable result that validated the models, it was important that the test data were different from the training data. Both training and test sets had to be produced by independent sampling from an infinite population, which avoided a misleading result that would not reflect what it was expected when the classifier is deployed. This section describes the used method that enabled us to apply this mandatory rule in machine learning. A 10 times repeated random sub-sampling, Monte Carlo Cross-Validation (MCCV) was used to evaluate the performance of each model (6 per subject and 24 in total). Results with Gaussian mixtures on real and simulated data suggested that MCCV provided genuine insight into cluster structure [31]. This method selects the most appropriate model for classification. To assure the ability of the model to generalize well, the validation was executed by applying, for each subject, the other three subject's recordings. This means that each trained HMM was used as classifier of the data of an unrelated subject. Each fold contained a trial of the three subjects. The split (80% train, 20% validation) was newly created during every validation (10 times) with a random assignment of the trials in the training and test sets. This led to a model trained with 24 exercises (8 of each subject). To perform the random assignment, the python built-in random function that implements the Mersenne Twister regenerator method was used. During each validation, 6 HMMs (models I–VI) were trained such as the best performing (based on the validation score) set of HMMs was selected as models set for classification. Forward probabilities were calculated for each HMM. When the correct HMM output the highest probability, the classification value became 1 and contrary 0. Per fold, each HMM classified the remaining 6 exercises where the performance per fold was the fraction correctly classified exercises (sum of classification results) of the total classifications (36) of the 6 HMMs combined. The model's parameters differed slightly between the sets as the random data selection altered the learned state Probability Density Functions (PDFs) per fold. The best performing model set out of the 10 validations was then selected to perform the classification for the test subject.

#### 5.3.4. Results

#### 1. State assignment and classification

The MLE for each HMM up to twenty states was used to define the BIC scores (see Equation (28)) against the amount of states (Figure 12). The profile of the BIC score against the amount of states was similar between the HMMs. The consensual lowest BIC score was obtained for an amount of states equal to thirteen (value that corresponds to the minimum of the average curve—the blue bold broken line in Figure 12). Thus, thirteen states were used to model the six movements. This makes intuitively sense as the exercise was constructed out of three distinctive parts (a multiple of three is expected), plus an initial/ending part (inactive state). Hence, each part in the exercise was described by four states.

**Figure 12.** BIC scores for each type of execution (I–VI) and an averaged BIC score over these executions (blue bold broken line). The black vertical line indicates the optimal amount of states.

The classification performance shows a high level of accuracy (Table 4) in classifying a whole sequence into the classes (I–VI). A value of 1 means the model always gave the highest probability, with respect to the other models and for any sequence of the related movement, whereas a value of 6 indicates the lowest probability. The values in this table are averaged prediction ranks for each model of each movement (I–VI). The average prediction rank of HMM I is the highest (2.78), which means that the execution type I (correct movement) is most closely related to all the other types. The overall performance of the classification for each class (I–VI) is shown in Table 5. It is to note that the execution type III (i) is more likely to be classified as type I, and (ii) has the lowest prediction accuracy compared with the other classes. Movement III is the one that presents the subtlest difference from movement I. The only distinction between these two executions is the fact that the participants do not move their center of mass in the former case. In other words, the trunk does not go with the rest of the movement. This fine difference could explain the limitation of the algorithm to discriminate between the two types of movement. A solution to overcome this issue could be to implement a better staging of the execution III and/or getting a higher descriptive power in the gesture representation (e.g., including additional features related to the upper-part of the body).

**Table 4.** Confusion matrix of executions (I–VI). Each column represents the types of movement and each row the output prediction ranks of the HMMs (I–VI). The closer is the value to 1 (green cells) the better is the prediction.

**Table 5.** Performance of the classification of movements I–VI.

#### 2. Real-time testing

The aim for the platform is not only to provide an overall classification of the movement (correct vs. types of fault), but also to identify a potential real-time transition from one classification to another during a single exercise. This can aid the patient to be aware on the phases of the movement in which certain errors tend to occur. The result of this instantaneous classification could be displayed as a real-time feedback when the patient is executing any exercise.

The samples of execution type VII are used to evaluate the ability for the trained models I–VI to classify partially incorrect movement (i.e., only certain phases of the movement are wrong). This execution is a rehabilitation movement in which the individual performs one step to the front, one step to the side and one step to the back, in a continuous sequence. The classification takes place over a selection of frames within the movement. The forward algorithm, as described previously, is used to perform this classification. Three different sizes of windows are used for the classification: 100, 60 and 20 frames. These different samplings are made to study the effect of the window size on the consistency and accuracy of the assessment. After classifying the frames of a determined window size, the window shifts half the number of frames in the total sequence and the classification is repeated until the end of the sequence is reached. This so-called overlapping window is used to obtain a smoother classification path over time. There are multiple classification values during the full exercise. At each newly created classification moment in the exercise the values of the six classifiers are normalized in a fashion that the highest value becomes 1 and the other values are expressed as a fraction of this value. Detection is considered accurate if the majority of the movement's phase where the error occurred assigns the value of 1 to the expected error type. In the case of the execution type VII, there are three phases: step forward, step sideways, and step backward.

Figure 13 shows that whatever the window size is the performance maintains quite accurate, even though the result is a little bit noisier when the sizes get smaller. For instance, there is a very high detection rate (21/24) when errors of types IV to VI are present in the sequence of the movements, for any sampling size. Detecting execution types II and III are less successful (9/16). It can be explained by the fact that these two types share high similarities with execution I (see Table 4).

**Figure 13.** Execution of type VII where the middle part (i.e., step to the side) is performed as type V. Three different window sizes are represented: 100, 60 and 20 samples. The red dotted line represents the prediction of HMM V and the blue line indicates the prediction of HMM I (correct movement).

There is a certain trade-off for choosing the window size. A smaller window can provide a frequent feedback, but a slightly noisier prediction. Figure 13 presents the example of a correct sequence, except the middle part is performed as execution V (step with bended knee). In this figure the amount of feedback moments is represented on the *x*-axis and a normalized classification value on the *y*-axis. The sampling rate is 50 Hz (20 ms per sample). Window sizes of 100, 60 and 20 represent approximately every second, twice a second and five times a second feedback, respectively. This example shows that the accuracy of the prediction (identification of correct vs. incorrect executions) is not significantly altered by the window sizes, which confirms the pertinence of an HMM approach for real-time applications.

#### **6. Usability of the System**

#### *6.1. Participants and Procedure*

In order to evaluate if the platform can be easily used by the end users, a usability test was carried out. A total of 41 right handed participants (distributed in 32 males and 9 females) volunteered for the experiment. No specific skill was necessary to take part to the study. 92.3% of the subjects was between 18 and 30 years old. All participants reported that they had never used a tele-rehabilitation system before.

Participants were informed about the purpose of the study and signed a consent form to participate. Then, subjects were asked to take an accommodation phase, in which they followed a tutorial on the use of the platform. Once the scenario of tele-rehabilitation was clarified and after the participants were confident with the application, they were asked to complete five tasks described in the next section. The total duration for the completion of the tasks was around 15 min. Once the whole tasks were completed, the subjects were requested to fill in the IBM Computer System Usability Questionnaire (CSUQ) described in Table 6 [32]. CSUQ is a questionnaire based on a Likert scale, in which a statement is made and the respondent then indicates the degree of agreement or disagreement with the statement on a 7-point scale, from 'strongly disagree' (1) to 'strongly agree' (7).

#### *6.2. Tasks*

The rehabilitation program is implemented in the platform according to 3 chronological stages. The intensity of the therapeutic exercises increases progressively from stage 1 (the week following the surgery) to stage 3 (about 3 weeks after the surgery). During the usability test, participants had to perform five tasks related to activities implemented in these different stages.

#### 6.2.1. Task 1: Consult Therapeutic Instructions of Stage 1

In this task, subjects were asked to access and understand the documentation provided by the core component regarding 'how to use the platform'. They could select among different learning modalities (text, audio, and video).

#### 6.2.2. Task 2: Perform Rehabilitation Exercises of Stage 1

This task involved several subtasks. First, participants had to answer a questionnaire to self-assess their health status and unlock the exercise interface (meaning that they are in physical condition to pursue). Second, they had to select an exercise of stage 1 (low intensity exercises) and consult recommendations to perform it properly.  Task 3: Consult Therapeutic Instructions of Stage 2

Again, participants had to access and consult general recommendation documents developed for the therapeutic education of the patients in stage 2 (medium intensity exercises). At the end, and before continuing to the rehabilitation program, they had to accept the conditions by clicking a check-box.

#### 6.2.4. Task 4: Perform Rehabilitation Exercises of Stage 2

As in task 2, subjects were requested to perform an exercise available in stage 2. An addition requirement for this task was to abort the series of repetition by using the functionality 'cancel' and, then, to fill in the questionnaire to justify the cause of the suspension.

#### 6.2.5. Task 5: Consult and Send a Message to the Therapist

In this task, participants had to access the mailbox integrated in the platform, in order to consult new messages and send a text-based or audio-based email to their physiotherapists.

#### *6.3. Results*

Overall, the participants mostly agreed with the statements. In most cases, the averages of the responses were greater than 5 with a standard deviation (SD) less than 1.19. The global result to the questions (#09, #10 and #11) that belong to the information quality or INFOQUAL (Mean = 5.49, Median = 6.00, SD = 1.15) suggests the platform provides a good online aid and documentation. Nevertheless, question #09 is assessed less positively (Mean = 4.77, Median = 5.00, SD = 1.51), which means that users tend to consider the error messages as not helpful to solve problems. Finally, the result to question #15 (Mean = 5.92, Median = 6.00, SD = 0.84) clearly indicates that the participants agree on the organization of the information of the user interface.

**Table 6.** Structure and results of the IBM Computer System Usability Questionnaire.

The average values for the IBM CSUQ measures are presented in Table 6. Concerning the category SYSUSE (Mean = 5.84, Median = 6.00, SD = 0.94), the results show that the participants tend to consider the system as useful. The results to the INFOQUAL category (Mean = 5.49, Median = 6.00, SD = 1.15) also suggest a good quality of the information presented on the user interface. The participants have indicated that the information provided by the system was easily understood and helpful to complete the tasks and scenarios. Nevertheless, future work will concentrate on improving the feedback provided by the system, especially when errors occur. Interface quality or INTERQUAL (Mean = 5.50, Median = 6.00, SD = 1.11) shows that the user interfaces of the platform are adapted to the end user. However, it seems that the tool requires some additional improvements, such as: (i) a simplification of the visual interface for the execution of the rehabilitation exercises, and (ii) the incorporation of a speech-based remote interaction with the application. Overall satisfaction or OVERALL (Mean = 5.74, Median = 6.00, SD = 0.88) suggests that the system is positively evaluated by the users. It can be concluded that the system usability of the current version of the platform is assessed positively by the

participants, who seem to like the general organization and functionalities of the application, even if some upgrades are required.

#### **7. Discussion**

This study describes the development of a tele-rehabilitation platform to support motor recovery. Exercises designed for patients after hip arthroplasty are presented as a case study. Nevertheless, this platform is built on a modular architecture that facilitates its adaptation for any other physical rehabilitation programs (e.g., hemiparetic patients). In addition, an artificial intelligence module is integrated in the platform, in order to assess the quality of the therapeutic movements performed by the patients. Two possible methods are tested: Dynamic Time Warping and Hidden Markov Models.

First, a DTW approach was used to evaluate the motion by comparison to a reference. The experimental results show a high correlation between the assessment performed by the computer and by the health professionals. The slight discrepancy between the two evaluations has both technical and human explanations. The accuracy is affected by partial occultations of some joints inherent to the motion capture device. In addition, we observed certain inconsistences in the judgment of the therapists, which limit the reliability of the human reference.

Second, a HMM approach is tested to detect variance within movement, caused by errors or compensatory movements that may occur during the completion of the therapeutic exercise. HMMs are trained on these errors and compensatory actions. Although the setting of the experiment was controlled, the classification also included intrapersonal and interpersonal variances, as a model that classified a determined subject was merely trained with the data of the other participants. It suggests that the proposed assessment algorithm has a fair capability of generalization. A high classification accuracy of the movements is obtained by building a general model that can be applied to any subject. A real-time analysis enables us to detect four out of five faulty movements, when these errors briefly occur in the beginning, middle or end of a correct execution of the exercise. The same level of accuracy is maintained whatever the detection rate. These findings demonstrate that the HMM is an appropriate method to provide real-time feedback regarding the correctness of the rehabilitation movement performed by a patient. This approach is successfully applied on a real-time assessment of components of the movement, which are discriminated in several classes that differ on extremely subtle aspects.

Both methods, DTW and HMM, have advantages and disadvantages. It is not necessary to collect lots of data to perform an assessment based on DTW, because just a single reference is required to apply the technique. On the contrary, HMM is a probabilistic approach, in which the accuracy of the classification depends on the quantity of available data. The bigger is the dataset, the higher is the precision. This aspect is particularly relevant if you presume a difficult access to collect enough data to train the model. On the other hand, HMM is more adapted than DTW to perform a real-time evaluation, as it creates a model that can be dynamically generated and adjusted. Another advantage of utilizing HMM is the fact that it is more robust than DTW when the quantity of features to discriminate increases.

Finally, the usability of the platform is positively assessed by end users. The empirical results based on subjective perception and self-reported feedback show that the application is useful, effective, efficient, and easy to use. In addition, the evaluation of the user experience enables us to identify usability aspects that should be implemented, in order to improve the visual interface. The experiment indicates that the error messages should be as detailed as possible. Our results suggest that providing the patients with a systematic and non-ambiguous feedback is a critical aspect to ensure a perfect acceptance of a tele-rehabilitation system.

The main contribution of this manuscript is to provide a holistic description of a solution to implement a tele-rehabilitation platform that presents three advantages rarely available in other similar systems. The first is the implementation of a scalable and modular architecture to facilitate (i) the adaptation to other physical rehabilitation programs, and (ii) the evolution of the system for integrating new software (e.g., game engines) and hardware (e.g., motion capture sensor) technologies. The second is the integration of an artificial intelligence (AI) module that assesses the quality of the movement through two possible solutions: DTW or HMM. This aspect, which is absent in comparable applications [12,33,34], is however crucial to ensure the efficiency and safety of a remote motor rehabilitation. The third is the gamification of the therapeutic exercises, in order to enhance the motivation of the patients to complete the rehabilitation program at home.

Nevertheless, the current version of the platform presents certain limitations that will be overcome in the future. The vision-based motion capture implies occultation issues that alter the accuracy of the assessment algorithms. A possible solution would be to implement a hybrid system that integrates both Kinect and inertial sensors [35]. In addition, it will be necessary to include additional features to get a more advanced discrimination of the movements assessed by a HMM approach. In particular, these features should stand for kinematic variables related to the upper part of the body, which will improve the identification of compensatory movements. In terms of evaluation of the AI algorithms, the future experimental protocols will simplify the assessment tasks of the therapists, in order to reduce the discrepancies. For instance, this assessment could be based on videos and not performed on the fly, and/or each expert could have to gauge fewer parameters simultaneously. Finally, we are aware that it is required to carry out usability tests on real patients and during the whole completion period of the rehabilitation program, as indicated in [36]. For this purpose, we plan to apply usability questionnaires specifically designed to measure the user experience with e-Health platforms, such as USEQ [36] and TUQ [37].

To conclude, we would like to highlight the fact that the application is built considering both the health professional requirements and the patient limitations. The architecture of the system is modular and flexible enough to be integrated in a global smart home. In addition, the proposed approach intends to deliver a therapeutic program as a service over the Internet. It means that the therapists can implement a personalized rehabilitation program that their patients will be able to execute anywhere through the platform. In that sense, the system promotes a smart medical environment, in which both patients and health professionals can interact remotely from home, for the former, and the medical office, for the latter; transforming the application into a ubiquitous tool to support the patient's recovery.

**Author Contributions:** Conceptualization, Y.R., J.L.P.M. and D.E.; methodology, Y.R., J.L.P.M. and L.L., software, K.J., L.L., Y.R. and J.L.P.M.; validation, Y.R., L.L. and J.L.P.M.; data curation, L.L., J.L.P.M. and D.E.; data analysis, Y.R., L.L., J.L.P.M. and M.G.; writing and editing, Y.R., J.L.P.M., L.L. and K.J.; project administration, Y.R.; funding acquisition, Y.R.

**Funding:** This research was funded by REDCEDIA, grant number CEPRA XI-2017-15.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Spline Function Simulation Data Generation for Walking Motion Using Foot-Mounted Inertial Sensors**

#### **Thanh Tuan Pham and Young Soo Suh \***

Electrical Engineering Department, University of Ulsan, Ulsan 44610, Korea<EMAIL_ADDRESS>**\*** Correspondence<EMAIL_ADDRESS>Tel.: +82-52-259-2196

Received: 22 October 2018; Accepted: 20 December 2018; Published: 23 December 2018

**Abstract:** This paper investigates the generation of simulation data for motion estimation using inertial sensors. The smoothing algorithm with waypoint-based map matching is proposed using foot-mounted inertial sensors to estimate position and attitude. The simulation data are generated using spline functions, where the estimated position and attitude are used as control points. The attitude is represented using B-spline quaternion and the position is represented by eighth-order algebraic splines. The simulation data can be generated using inertial sensors (accelerometer and gyroscope) without using any additional sensors. Through indoor experiments, two scenarios were examined include 2D walking path (rectangular) and 3D walking path (corridor and stairs) for simulation data generation. The proposed simulation data is used to evaluate the estimation performance with different parameters such as different noise levels and sampling periods.

**Keywords:** motion estimation; inertial sensors; simulation; spline function; Kalman filter

#### **1. Introduction**

Motion estimation using inertial sensors is one of the most important research topics that is increasingly applied in many application areas such as medical applications, sports, and entertainment [1]. Inertial measurement unit (IMU) sensors are commonly used to estimate human motion [2–5]. Inertial sensors can be used alone or combined with other sensors such as cameras [6,7] or magnetic sensors [8]. In personal navigation and healthcare applications [9–11], foot-mounted inertial sensors without using any additional sensors is a key enabling technology since it does not require any additional infrastructure. If additional sensors (other than inertial sensors) are also used, sensor fusion algorithms [12,13] can be used to obtain more accurate motion estimation.

To evaluate the accuracy of any motion estimation algorithm, it is necessary to compare the estimated value with the true value. The estimated value is usually verified through experiments with both IMU and optical motion tracker [14,15]. In [2], the Vicon optical motion capture system is used to evaluate the effectiveness of human pose estimation method and its associated sensor calibration procedure. An optical motion tracking is used as a reference to compare with motion estimation using inertial sensors in [15].

Although the experimental validation gives the ultimate proof, it is not only more convenient but also more important to work with simulation data, which is statistically much more relevant for testing and the development of algorithms. For example, the effect of gyroscope bias on the performance of an algorithm cannot be easily identified in experiments. In this case, simulation using synthesized data is more convenient since gyroscope bias can be changed arbitrarily.

In [16–18], IMU simulation data are generated using optical motion capture data for human motion estimation. In [19], IMU simulation data is generated for walking motion by approximating walking trajectory as simple sinusoidal functions. In [20–22], IMU simulation data are also generated for flight motion estimation using an artificially generated trajectory. There are many different methods to represent attitude and position simulation data. Among them, the spline function is the most popular approach [23]. The B-spline quaternion algorithm provides a general curve construction scheme that extends the spline curves [24]. Reference [25] presents an approach for attitude estimation using cumulative B-splines unit quaternion curves. Position can be represented by any spline function [23]. Usually, position data are represented using cubic spline functions. The cubic spline function gives a continuous position, velocity, and acceleration. However, the acceleration is not sufficiently smooth. This is a problem for the inertial sensor simulation data generation since the smooth accelerometer data cannot be generated. Thus, we use the eighth-order spline function for position data since it gives position data with smooth acceleration, whose jerk (third derivative) can be controlled in the optimization problem.

In this paper, we use cumulative B-spline for attitude representation and eighth-order algebraic spline for position representation. The eighth-order algebraic spline is a slightly modified version of the spline function from a previous study [26].

The aim of this paper is to obtain simulation data combining inertial sensors (accelerometer and gyroscope) and waypoint data without using any additional sensors. A computationally efficient waypoint-based map matching smoothing algorithm is proposed for position and attitude estimation from foot-mounted inertial sensors. The simulation data are generated using spline functions, where the estimated position and attitude are used as control points. This paper is an extended version of a work published in the conference paper [27], where a basic algorithm with a simple experimental result is given.

The remainder of the paper is organized as follows. Section 2 describes the smoothing algorithm with waypoint-based map matching to estimate motion for foot-mounted inertial sensors. Section 3 describes the computation of the spline function to generate simulation data. Section 4 provides the experiment results to verify the proposed method. Conclusions are given in Section 5.

#### **2. Smoothing Algorithm with Waypoint-Based Map Matching**

In this section, a smoothing algorithm is proposed to estimate the attitude, position, velocity and acceleration, which will be used to generate simulation data in Section 3.

Two coordinate frames are used in this paper: the body coordinate frame and the navigation coordinate frame. The three axes of the body coordinate frame coincide with the three axes of the IMU. The *z* axis of the navigation coordinate frame coincides with the local gravitational direction. The choice of the *x* and *y* axes of the navigation coordinate frame can be arbitrarily chosen. The notation [*p*]*n* ([*p*]*b*) is used to denote that a vector *p* is represented in the navigation (body) coordinate frame.

The position is defined by [*r*]*n* ∈ *R*3, which is the origin of the body coordinate frame expressed in the navigation coordinate frame. Similarly, the velocity and the acceleration are denoted by [*v*]*n* ∈ *R*3 and [*a*]*n* ∈ *R*3, respectively. The attitude is represented using a quaternion *q* ∈ *R*4, which represents the rotation relationship between the navigation coordinate frame and the body coordinate frame. The directional cosine matrix corresponding to quaternion *q* is denoted by *C*(*q*) ∈ *SO*(3).

The accelerometer output *ya* ∈ *R*3 and the gyroscope output *yg* ∈ *R*3 are given by

$$\begin{aligned} y\_{\mathfrak{a}} &= \mathbb{C}(q)\mathfrak{F} + a\_{\mathfrak{b}} + n\_{\mathfrak{a}\prime} \\ y\_{\mathfrak{k}} &= \omega + n\_{\mathfrak{k}\prime} \end{aligned} \tag{1}$$

where *ω* ∈ *R*3 is the angular velocity, *ab* ∈ *R*3 is the external acceleration (acceleration related to the movement, excluding the gravitational acceleration) expressed in the body coordinate frame, and *na* ∈ *R*3 and *ng* ∈ *R*3 are sensor noises. The vector *g*˜ is the local gravitational acceleration vector. It is assumed that *g*˜ is known, which can be computed using the formula in [28]. The sensor biases are assumed to be estimated separately using calibration algorithms [29,30].

Let *T* denote the sampling period of a sensor. For a continuous time signal *y*(*t*), the discrete value is denoted by *yk* = *y*(*kT*). The discrete sensor noise *na*,*k* and *ng*,*k* are assumed to be white Gaussian sensor noises, whose covariances are given by

$$\begin{aligned} R\_d &= \mathbb{E}\left\{ n\_{d,k} n\_{d,k}^T \right\} = r\_d I\_3 \in R^{3 \times 3}, \\ R\_{\mathfrak{F}} &= \mathbb{E}\left\{ n\_{\mathfrak{g},k} n\_{\mathfrak{g},k}^T \right\} = r\_{\mathfrak{F}} I\_3 \in R^{3 \times 3}. \end{aligned} \tag{2}$$

where *ra* > 0 and *rg* > 0 are scalar constants.

We assume that a walking trajectory consists of straight line paths, where the angle between two adjacent paths can only take one of the following angle {90◦, −90◦, 180◦}. An example of a walking trajectory is given in Figure 1, where there is 90◦ turn at *P*2 and −90◦ turn at *P*3. Waypoints are denoted by *Pm* ∈ *R*3 (1 ≤ *m* ≤ *M*), which include positions with turn events (*P*2 and *P*3 in Figure 1), the initial position (*P*1) and the final position (*P*4).

**Figure 1.** Example of waypoints (*P*1 ∼ *P*4).

#### *2.1. Standard Inertial Navigation Using an Indirect Kalman Filter*

In this subsection, *q*, *r*, and *v* are estimated using a standard inertial navigation algorithm with an indirect Kalman filter [31].

The basic equations for inertial navigation are given as follows:

$$\begin{aligned} \eta\_{\ \prime} &= \frac{1}{2} \Omega(\omega) q\_{\prime} \\ \psi\_{\ n} &= a\_{\ n} = \mathbb{C}^{T}(q) a\_{\mathbb{W}} \\ \dot{v}\_{\ n} &= v\_{\mathbb{W}} \end{aligned} \tag{3}$$

where symbol Ω is defined by

$$
\Omega(\omega) = \begin{bmatrix} 0 & \omega^T \\ \omega & -[\omega \times ] \end{bmatrix},
$$

and [*ω*×] ∈ *R*3×3 denotes the skew symmetric matrix of *ω*:

$$\begin{aligned} \left[\omega \times \right] = \begin{bmatrix} 0 & -\omega\_{\mathbb{Z}} & \omega\_{\mathbb{Y}} \\ & \omega\_{\mathbb{Z}} & 0 & -\omega\_{\mathbb{X}} \\ -\omega\_{\mathbb{Y}} & \omega\_{\mathbb{X}} & 0 \end{bmatrix} . \end{aligned} $$

Let *q*ˆ*k*, *r*ˆ*k* and *v*ˆ*k* be the estimated values of *q*, *r* and *v* using (3), where *ω* and *ab* are replaced by *yg* and *ya* − *C*(*q*ˆ)*g*˜:

$$\begin{aligned} \dot{\mathfrak{q}} &= \frac{1}{2} \Omega(y\_{\mathcal{J}}) \dot{\mathfrak{q}},\\ \dot{\mathfrak{d}} &= \mathbb{C}^{T}(\dot{\mathfrak{q}}) y\_{\mathcal{A}} - \overline{\mathfrak{g}},\\ \dot{\mathfrak{r}} &= \mathfrak{d}. \end{aligned} \tag{4}$$

Coriolis and Earth curvature effects are ignored since we use a consumer grade IMU for limited walking distances where these effects are below the IMU noise level [32].

When the initial attitude (*q*ˆ0) is computed, the pitch and roll angles can be computed from *ya*,0. However, the yaw angle is not determined since there is no heading reference sensor such as a magnetic

sensor. In the proposed algorithm, the initial yaw angle can be arbitrarily chosen since the yaw angle is automatically adjusted later (see Section 2.3).

*q*¯*e*, *re* and *ve* denote the estimation errors in *q*ˆ, *r*ˆ, and *v*ˆ, which are defined by

$$\begin{aligned} \eta\_{\varepsilon} & \stackrel{\triangle}{=} \eta \otimes \mathfrak{H}\_{k'}^{\*}\\ \tau\_{\varepsilon} & \stackrel{\triangle}{=} r\_{\text{fl}} - \hat{r}\_{k'}\\ \upsilon\_{\varepsilon} & \stackrel{\triangle}{=} \upsilon\_{\text{h}} - \mathfrak{H}\_{k'} \end{aligned} \tag{5}$$

where ⊗ is the quaternion multiplication and *q*∗ denotes the quaternion conjugate of a quaternion *q*. Assuming that *q*¯*e* is small, we can approximate *q*¯*e* as follows:

$$
\eta\_e \approx \begin{bmatrix} & \mathbf{1} \\ & q\_e \end{bmatrix} \tag{6}
$$

The state of an indirect Kalman filter is defined by

$$\mathbf{x} \triangleq \begin{bmatrix} q\_{\ell} \\ r\_{\ell} \\ \upsilon\_{\ell} \\ \upsilon\_{\ell} \end{bmatrix} \in \mathbb{R}^{9 \times 1}. \tag{7}$$

The state equation for the Kalman filter is given by [31]:

$$\mathbf{x} = \begin{bmatrix} [-y\_{\mathcal{R}} \times] & 0 & 0 \\ 0 & 0 & I\_3 \\ -2\mathcal{C}(\hat{q})^T [y\_{\mathcal{d}} \times] & 0 & 0 \end{bmatrix} \mathbf{x} + \mathbf{w} \tag{8}$$

where the covariance of the process noise *w* is given by

$$\mathbb{E}\left\{ww^{\mathsf{T}}\right\} = \begin{bmatrix} 0.25R\_{\mathcal{S}} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & C(\mathfrak{q})^{T}R\_{\mathcal{A}}C(\mathfrak{q}) \end{bmatrix}.$$

The discretized version of (8) is used in the Kalman filter as in [31].

Two measurement equations are used in the Kalman filter: the zero velocity updating (ZUPT) equation and the map matching equation.

The zero velocity updating uses the fact that there is almost periodic zero velocity intervals (when a foot is on the ground) during walking. If the following conditions are satisfied, the discrete time index *k* is assumed to belong to zero velocity intervals [31]:

$$\begin{array}{ll}||y\_{\mathcal{S},i}|| \le B\_{\mathcal{S}'} & k - \frac{N\_{\mathcal{S}}}{2} \le i \le k + \frac{N\_{\mathcal{S}}}{2} \\ ||y\_{a,i} - y\_{a,i-1}|| \le B\_{a\prime} & k - \frac{N\_{\mathcal{S}}}{2} \le i \le k + \frac{N\_{\mathcal{S}}}{2} \end{array} \tag{9}$$

where *Bg*, *Ba*, *Ng* and *Na* are parameters for zero velocity interval detection. This zero velocity updating algorithm is not valid when a person is on a moving transportation such as an elevator. Thus, the proposed algorithm cannot be used when a person is on a moving transportation.

During the zero velocity intervals, the following measurement equation is used:

$$\begin{bmatrix} 0\_{3 \times 1} - \vartheta\_k = \begin{bmatrix} 0\_{3 \times 3} & 0\_{3 \times 3} & I\_3 \end{bmatrix} \mathbf{x}\_k + \mathbf{n}\_{\mathbb{P}\_\ell} \tag{10}$$

where *nv* is a fictitious measurement noise representing a Gaussian white noise with the noise covariance *Rv*.

The second measurement equation is from the map matching, which is used during the zero velocity intervals. From the assumption, a straight line path is parallel to either *x* axis or *y* axis. If a path is parallel to *x*(*y*) axis, *y*(*x*) position is constant and this can be used in the measurement equation.

Let *ei* ∈ *R*3×1(1 ≤ *i* ≤ 3) be the unit vector whose *i*-th element is 1 and the remaining elements are 0. Suppose a person is on the path [*Pm*, *Pm*+1] and

$$e\_j^T(P\_{m+1} - P\_m) = 0, \quad j = 1 \text{ or } 2. \tag{11}$$

For example, (11) is satisfied with *j* = 1 if the path is parallel to the *y* axis.

When (11) is satisfied with *j* (*j* = 1 or 2), then the map matching measurement equation is given by

$$\begin{array}{cc} \mathbf{e}\_{\mathbf{j}}^T (P\_m - \mathfrak{k}\_k) = \begin{bmatrix} \mathbf{0}\_{1 \times 3} & \mathbf{e}\_{\mathbf{j}}^T & \mathbf{0}\_{1 \times 3} \end{bmatrix} \mathbf{x}\_k + n\_{r, 12 \text{'}} \end{array} \tag{12}$$

where *nr*,12 is the horizontal position measurement noise whose noise covariance is *Rr*,12.

If the path is level (that is, *eT* 3 (*Pm*+1 − *Pm*) = 0), the *z* axis value of *rk* is almost the same in the zero velocity intervals when the foot is on the ground. In this case, the following *z* axis measurement equation is also used:

$$e\_3^T (P\_m - \mathfrak{k}\_k) = \begin{bmatrix} \ 0\_{1 \times 3} & e\_3^T & 0\_{1 \times 3} \end{bmatrix} \mathbf{x}\_k + n\_{r, 3 \prime} \tag{13}$$

where *nr*,3 is the vertical position measurement noise whose noise covariance is *Rr*,3.

The proposed filtering algorithm is illustrated in Figure 2. To further reduce the estimation errors of the Kalman filter, the smoothing algorithm in [31] is applied to the Kalman filter estimated values.

**Figure 2.** Proposed measurement update flowchart.

#### *2.2. Path Identification*

To apply map matching measurement Equation (12), a current path must be identified. Using the zero velocity intervals, each walking step can be easily determined. Let *Sl* be the discrete time index of the end of *l*-th zero velocity interval (see Figure 3).

**Figure 3.** Walking step segmentation.

Let *ψSl* be the yaw angle at the discrete time *Sl* and *ψl* be defined by

$$
\triangle \,\Psi \rangle \stackrel{\triangle}{=} \Psi\_{\mathcal{S}\_{l+1}} - \Psi\_{\mathcal{S}\_l}.\tag{14}
$$

Please note that *ψl* denotes the yaw angle change during *l*-th walking step.

The turning (that is, the change of a path) is detected using *ψl*. Suppose a current path is [*Pm*, *Pm*+1]. Then the next turning angle is determined by two vectors *Pm*+1 − *Pm* and *Pm*+2 − *Pm*+1. For example in Figure 1, if a current path is [*P*1, *P*2], the next turning angle is 90◦. Let *δ* be the next turning angle. Turning is determined to occur at *l*-th walking step if the following condition is satisfied:

$$|\sum\_{j=l-M\_l}^{l+M\_l} \triangle \psi\_j - \delta \mid < \psi\_{\text{thv}} \tag{15}$$

where *ψ*th is a threshold parameter and *Ml* is a positive integer parameter. Equation (15) detects the turning if the summation of the yaw angle changes during 2*Ml* + 1 walking steps is similar to the next expected turning angle *δ*. The use of parameter *Ml* is due to the fact that turning does not occur during a single step. Once turning is detected in (15), a current path is modified to [*Pm*+1, *Pm*+2] from [*Pm*, *Pm*+1].

#### *2.3. Initial Yaw Angle Adjustment*

The initial yaw angle is arbitrarily chosen and adjusted as follows. Suppose the first *N*first walking steps belong to the first path [*P*1, *P*2]. Let *xl* and *yl* (1 ≤ *l* ≤ *Nl*) be defined by

$$
\begin{bmatrix} \mathbf{x}\_{l} \\ \mathbf{y}\_{l} \end{bmatrix} = \begin{bmatrix} e\_{1}^{T} \\ e\_{2}^{T} \end{bmatrix} \mathbf{f}\_{\mathbf{S}\_{l'}} \quad \mathbf{1} \le l \le N\_{l}. \tag{16}
$$

The line equation representing (*xl*, *yl*) is computed using the following least squares optimization:

$$\min\_{a,b,c} \sum\_{l=1}^{N\_l} \left\| ax\_l + by\_l + c \right\|\_2^2. \tag{17}$$

The angle between the line (defined by the line equation *a*, *b*, *c*) and [ *e*1 *e*2 ] *T*(*P*2 − *P*1) is computed. Using this angle, the initial yaw angle is adjusted so that the walking direction coincides with the direction dictated by the waypoints.

#### **3. Spline Function Computation**

In this section, spline function *q*¯(*t*) (quaternion spline) and *r*¯(*t*) (position spline) are computed using *q*ˆ*k*, *r*ˆ*k* and *v*ˆ*k* as control points.

#### *3.1. Cumulative B-Splines Quaternion Curve*

Since the quaternion spline *q*¯(*t*) is computed for each interval [*kT*,(*k* + 1)*T*), it is convenient to introduce *q*¯*k*(*u*) notation, which is defined by

$$\overline{q}(t) = \overline{q}\_k(u) \triangleq \overline{q}(kT + u), \ t \in [kT, (k+1)T), \ 0 \le u < T. \tag{18}$$

To define *q*¯*k*(*u*), *ω*ˆ *k* is defined as follows:

$$
\hat{\omega}\_k \triangleq \frac{1}{T} \log(\mathbb{C}(\mathfrak{q}\_{k+1})\mathbb{C}^T(\mathfrak{q}\_k)),
\tag{19}
$$

where the logarithm of an orthogonal matrix is defined in [33]. Please note that *w*ˆ *k* is a constant angular velocity vector, which transform from *q*ˆ*k* to *q*ˆ*k*+1.

Cumulative B-spline basis function is used to represent the quaternion spline, which is first proposed in [24] and also used in [34].

$$\mathfrak{q}\_k(\boldsymbol{\mu}) = (\prod\_{i=1}^3 \mathcal{S}\_i(\boldsymbol{\mu})) \mathfrak{q}\_{k-1} \tag{20}$$

where

$$S\_i(\mu) \stackrel{\triangle}{=} \exp\left(\frac{1}{2}\vec{B}\_i(\mu)\Omega(\hat{\omega}\_{k-2+i})\right).$$

The cumulative basis function *B*˜ *i*(*u*) ∈ *R*3 (1 ≤ *i* ≤ 3) is defined by

$$\begin{bmatrix} \mathcal{B}\_1(u) \\ \mathcal{B}\_2(u) \\ \mathcal{B}\_3(u) \end{bmatrix} = D \begin{bmatrix} 1 \\ \frac{u}{T} \\ \frac{u^2}{T^2} \\ \frac{u^3}{T^3} \end{bmatrix} \tag{21}$$

where *D* is computed using the matrix representation of the De Boor–Cox formula [34]:

$$D = \frac{1}{6} \left[ \begin{array}{cccc} 5 & 3 & -3 & 1 \\ 1 & 3 & 3 & -2 \\ 0 & 0 & 0 & 1 \end{array} \right].$$

To generate the gyroscope simulation data, the angular velocity spline *ω*¯(*t*) is required. From the relationship between *q*¯(*t*) and *ω*¯(*t*) [35], *ω*¯(*t*) is given by

$$
\omega(t) = 2\Xi(\vartheta(t))\vartheta(t) \tag{22}
$$

where

$$
\Xi(q) \triangleq \begin{bmatrix} -q\_1 & -q\_2 & -q\_3 \\ q\_0 & -q\_3 & q\_2 \\ q\_3 & q\_0 & -q\_1 \\ -q\_2 & q\_1 & q\_0 \end{bmatrix}
$$

The derivative of *q*¯(*t*) (*kT* ≤ *t* < (*k* + 1)*T*) is given by

$$
\hat{q}(t) = \frac{d\vec{q}\_k(u)}{du} = (\dot{S}\_3 S\_2 S\_1 + S\_3 \dot{S}\_2 S\_1 + S\_3 S\_2 \dot{S}\_1) \hat{q}\_{k-1} \tag{23}
$$

.

where

$$
\dot{S}\_i = \frac{dS\_i(\mu)}{du} = \frac{1}{2} \dot{B}\_i(\mu) \Omega(\hat{\omega}\_{k-2+i}) S\_i
$$

$$
\dot{B}(\mu) = \begin{bmatrix}
\dot{B}\_1(\mu) \\
\dot{B}\_2(\mu) \\
\dot{B}\_3(\mu)
\end{bmatrix} = D \begin{bmatrix}
0 \\
\frac{1}{T} \\
\frac{2u}{T^2} \\
\frac{3u^2}{T^3}
\end{bmatrix}.
$$

#### *3.2. Eighth-Order Algebraic Splines*

The position *r*¯(*t*) is represented by eighth-order algebraic spline [26], where *k*-th spline segment is given by

$$\vec{r}(t) = \vec{r}\_{k,m}(\boldsymbol{u}) = \sum\_{j=0}^{7} p\_{kj,m} \boldsymbol{u}^{j}, \ t \in [kT, (k+1)T). \tag{24}$$

where *pkj*,*m*(0 ≤ *k* ≤ *N*, 0 ≤ *j* ≤ 7, 1 ≤ *m* ≤ 3) are the coefficients of the *k*-th spline segment of the *m*-th element of *r*¯*k* (for example, *r*¯*k*,1 is the *x* position spline). Continuity constraints up to third derivative of *r*¯*k*,*m*(*u*) are imposed on the coefficients.

Let *a*ˆ*k* be the estimated value of the external acceleration in the navigation coordinate frame, which can be computed as follows (see (3))

$$\mathfrak{d}\_{k} = \mathbb{C}^{T}(\mathfrak{d}\_{k})\mathfrak{y}\_{a,k} - \mathfrak{g}.\tag{25}$$

The coefficients *pkj*,*m* are chosen so that *r*¯(*kT*), *r*¯˙(*kT*) and *r*¨¯(*kT*) are close to *r*ˆ*k*, *v*ˆ*k* and *a*ˆ*k*. An additional constraint is also imposed to reduce the jerk term (the third derivative of *r*¯(*t*)), which makes *r*¯(*t*) smooth rather than jerky. Thus, the performance index for *m*-th element (1 ≤ *m* ≤ 3) of *r*¯(*t*) is defined as follows:

$$\begin{split} J\_{m} &= \mathfrak{a}\_{1} \sum\_{k=0}^{N} (\mathfrak{r}\_{k,m}(T) - \mathfrak{r}\_{k,m})^{2} + \mathfrak{a}\_{2} \sum\_{k=0}^{N} (\mathfrak{r}\_{k,m}(T) - \mathfrak{d}\_{k,m})^{2} \\ &+ \mathfrak{a}\_{3} \sum\_{k=0}^{N} (\mathfrak{r}\_{k,m}(T) - \mathfrak{a}\_{k,m})^{2} + \mathfrak{f} \sum\_{k=0}^{N} \int\_{0}^{T} (\overleftarrow{\ } \mathfrak{r}\_{k,m}(t))^{2} dt. \end{split} \tag{26}$$

This minimization of *Jm* can be formulated as a quadratic minimization problem of the coefficients *pkj*,*m* as follows:

$$\min\_{X} J\_{\mathcal{W}} = X^T M\_1 X + 2M\_2 X + M\_3 \tag{27}$$

where *M*1 ∈ *R*8(*N*−1)×8(*N*−1), *M*2 ∈ *R*1×8(*N*−1) and *M*3 ∈ *R* can be computed from (26).

Although the matrix size of *M*1 could be large, it is a banded matrix with a bandwidth of 7. An efficient Cholesky decomposition algorithm is available for matrices with a small bandwidth [36] and (27) can be solved efficiently.

Let *p*¯*kj*,*m* be the minimization solution of (26). Then *r*¯ can be computed by inserting *p*¯*kj*,*m* into (24) and *v*¯, *a*¯ can be computed by taking the derivatives.

The weighting factors *αi*(1 ≤ *i* ≤ 3) and *β* determine the smoothness of *r*¯, *v*¯ and *a*¯. If *β* is large, we obtain smoother *r*¯, *v*¯ and *a*¯ curves at the sacrifice of nearness to the control points *r*ˆ, *v*ˆ and *a*ˆ.

Let *y*¯*a* and *y*¯*g* be accelerometer and gyroscope output generated using the spline functions. These values can be generated using (1), where *q*, *ab*, *ω* are replaced by *q*¯, *C*(*q*¯)*a*¯, *ω*¯ with appropriate sensor noises. An overview of the simulation data generation system is shown in Figure 4.

**Figure 4.** System overview for simulation data generation.