# Chaotic Properties of Gait Kinematic Data

Michal Piorek

► To cite this version:

Michal Piorek. Chaotic Properties of Gait Kinematic Data. 14th Computer Information Systems and Industrial Management (CISIM), Sep 2015, Warsaw, Poland. pp.111-119, 10.1007/978-3-319-24369-6_9. hal-01444509

HAL Id: hal-01444509

https://inria.hal.science/hal-01444509

Submitted on 24 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Chaotic properties of gait kinematic data

Michal Piorek

Department of Computer Engineering
Wroclaw University of Technology
Janiszewskiego 11/17, 50-372 Wroclaw, Poland
michal.piorek@pwr.edu.pl

**Abstract.** Time delay reconstruction for real systems is a widely explored area of nonlinear time series analysis. However, the majority of related work relates only to univariate time series, while multivariate time series data are common too. One such example is human gait kinematic data. The main goal of this article is to present a method of nonlinear analysis for kinematic time series. This nonlinear analysis is designed for detection of chaotic behavior. The presented approach also allows for the largest Lyapunov's exponent estimation for kinematic time series. This factor helps in judging the stability of the examined system and its chaotic properties.

**Keywords:** deterministic chaos, nonlinear time series analysis, quaternions, human motion analysis, human gait data

## 1 Introduction

Multivariate time series data are common in real systems. Many of those systems are results of the evolution of nonlinear systems dynamics. To assess the chaotic origin of time series, time delay reconstruction of a phase space is required. This step provides a view of the dynamics of the underlying system and allows the estimation of other properties of the investigated system (e.g. Lyapunov's largest exponent or one of the fractal dimensions of the reconstructed attractor). Based on the estimations of several parameters characterizing nonlinear processes, a decision about chaotic origin of time series may be taken.

According to the embedding theorem ([27]), for recovering dynamics only a univariate time series is needed, but often measurements of more than one quantityes related to the same dynamical system are available. One such case is the gait kinematic data for patients suffering from various diseases affecting walking ability.

In related work, there is a large number of applications of embedding theorem including univariate time series embedding ([24], [25], [28]), multivariate time series embedding ([3], [10], [18]), modelling ([1], [9], [11]), chaos controll ([2], [5], [22]), noise reduction ([4], [12], [19]) and signal classification ([26]).

The method presented in this paper is designed for a nonlinear analysis procedure for gait kinematics time series (Euler Angles or Quaternion) aimed at
computing the Largest Lyapunov Exponent for reconstructed dynamics for the purpose of identifying of deterministic chaos in gait kinematic data.

The presented method does not assume that investigated system is a deterministic chaotic system in the first place. This empirical approach is aimed to investigate the origin of the examined time series.

In this approach the Largest Lyapunov Exponent is also treated as a measure of human locomotion stability. In this case stability is defined as the sensitivity of a dynamic system to perturbations. This is a very popular factor in biomedical applications ([7], [6]).

The second section describes gait kinematic data time series - a subject of further investigations. It also includes information about conversion to quaternion's angle time series which allows for further efficient computations. The third section includes information about methods used in a nonlinear analysis procedure. The fourth section presents the numerical results. The conclusions are presented in section five.

## 2 Gait Kinematics Data Time series

There are three kinds of parameterisation of orientation space. Euler angles are good for human understanding of angular position. Matrices are able to do the calculations. However, quaternions are computationally efficient whilst also avoiding the singularities of Euler angles ([20]).

Quaternion parameterisation describes human motion which requires significant amounts of information about rotational displacement of selected segments. The method presented in this paper is designed for data captured from physical systems. We assumed that, to make a proper estimate of finite-time Lyapunov exponents experimentally, it is necessary to collect time series data captured from a large number of consecutive strides of gait ([14]).

The data are recorded as time series formed by Euler Angles. It is the most common parameterisation especially for biomedical applications.

$$s(n) = [\beta_1(t_0 + n\Delta t), \beta_2(t_0 + n\Delta t), \beta_3(t_0 + n\Delta t)] \quad (1)$$

where $s(n)$ is a n-th sample of Euler angles measured in interval time $\Delta t$ from initial time $t_0$ and $\beta_1, \beta_2, \beta_3$ are Euler angles in a X-Y-Z sequence. Axes X, Y and Z are defined as unit vectors:

$$X = (1, 0, 0), Y = (0, 1, 0), Z = (0, 0, 1) \quad (2)$$

Due to the fact that an Euler Angles time series is a multivariate time series, the following procedure was used to obtain better efficiency of computations.

Based on the assumption that axes X, Y and Z are unit vectors, rotation coding in Euler's Angles and in quaternions are identical. We can define quaternions
for each base rotation in a Euler sequence:

$$ 
\begin{aligned}
q_x(\beta_1) &= \cos(\frac{\beta_1}{2}) + \underbrace{(\sin(\frac{\beta_1}{2}), 0, 0)}_{\rightarrow} \\
q_y(\beta_2) &= \cos(\frac{\beta_2}{2}) + \underbrace{(0, \sin(\frac{\beta_2}{2}), 0)}_{\rightarrow} \\
q_z(\beta_3) &= \cos(\frac{\beta_3}{2}) + \underbrace{(0, 0, \sin(\frac{\beta_3}{2}))}_{\rightarrow}
\end{aligned} 
\quad (3) $$

Finally a quaternion equivalent to Euler angles representation can be calculated from three consecutive rotations described by the following quaternions:

$$ 
\begin{aligned}
q_{xyz}(\beta_1, \beta_2, \beta_3) &= q_x(\beta_1)q_y(\beta_2)q_z(\beta_3) = \\
&\quad (\cos(\frac{\beta_1}{2})\cos(\frac{\beta_2}{2})\cos(\frac{\beta_3}{2}) - \sin(\frac{\beta_1}{2})\sin(\frac{\beta_2}{2})\sin(\frac{\beta_3}{2})) + \\
&\quad i(\sin(\frac{\beta_1}{2})\cos(\frac{\beta_2}{2})\cos(\frac{\beta_3}{2}) - \cos(\frac{\beta_1}{2})\sin(\frac{\beta_2}{2})\sin(\frac{\beta_3}{2})) + \\
&\quad j(\cos(\frac{\beta_1}{2})\sin(\frac{\beta_2}{2})\cos(\frac{\beta_3}{2}) - \sin(\frac{\beta_1}{2})\cos(\frac{\beta_2}{2})\sin(\frac{\beta_3}{2})) + \\
&\quad k(\cos(\frac{\beta_1}{2})\cos(\frac{\beta_2}{2})\sin(\frac{\beta_3}{2}) - \sin(\frac{\beta_1}{2})\sin(\frac{\beta_2}{2})\cos(\frac{\beta_3}{2}))
\end{aligned} 
\quad (4) $$

Quaternion's time series is formed by conversion of each value of a Euler's angles time series to unit quaternion

$$ q(n) = q_{xyz}(s(n)) = q_{xyz}(\beta_1(t_0 + n\Delta t), \beta_2(t_0 + n\Delta t), \beta_3(t_0 + n\Delta t)) \quad (5) $$

Based on the assumption that there is greater variability in the quaternion's angle than its axis, nonlinear analysis directed at identifying the presence of deterministic chaos and local stability investigation is performed on the time series formed by angles of quaternion $q(n)$

$$ \alpha(n) = 2\arccos(\operatorname{real}(q(t_0 + n\Delta t))) \quad (6) $$

## 3 Nonlinear analysis procedure for quaternion angle time series

The nonlinear analysis procedure consists of two steps: time delay reconstruction and the largest Lyapunov exponent estimation.

### 3.1 Time delay reconstruction

According to the Takens embedding theorem ([27]), its possible to reconstruct the state trajectory from a single time series using the algorithm below:

$$ y(n) = [s(n), s(n + T), \dots, s(n + (d - 1))T], \quad (7) $$
where $T$ is a time delay and $d$ is an embedding dimension, which estimates a real dimension of the observed system. The main point of the state space reconstruction method is $T$ and $d$ estimation. To estimate time delay $T$, the average mutual information $I$ has been used, while for the embedding dimension the false nearest neighbor method ([1]).

The mutual information approach is based on information theory and transformation of linear autocorrelation to non-linear systems. More precisely, this method consists of 2-dimensional adaptive histogram ([9]).

Let's assume that there are two nonlinear systems: A and B. The outputs of these systems are denoted as $a$ and $b$, while the values of these outputs are represented by $a_i$ and $b_k$. The mutual information factor describes how many bits of $b_k$ could be predicted where $a_i$ is known.

$$I_{AB}(a_i, b_k) = \log_2 \left( \frac{P_{AB}(a_i, b_k)}{P_A(a_i)P_B(b_k)} \right), \qquad (8)$$

where $P_A(a_i)$ is the probability that $a = a_i$ and $P_B(b_k)$ is the probability that $b = b_k$ and $P_{AB}(a_i, b_k)$ is the joint probability that $a = a_i$ and $b = b_k$.

The average mutual information factor can be described by:

$$I_{AB}(T) = \sum_{a_i, b_k} P_{AB}(a_i, b_k)I_{AB}(a_i, b_k). \qquad (9)$$

In order to use this method to assess the correlation between different samples in the same time series, the Average mutual information factor is finally described by the equation:

$$I(T) = \sum_{n=1}^{N} P(S(n), S(n+T)) \log_2 \left( \frac{P(S(n),S(n+T))}{P(S(n))P(S(n+T))} \right). \qquad (10)$$

Fraser and Swinney ([9]) propose that $T_m$ where the first minimum of $I(T)$ occurs as a useful selection of time lag $T_d$. This selection guarantees that the measurements are somewhat independent, but not statistically independent. In case of absence of the average mutual information clear minimum, this criterion needs to be replaced by choosing $T_d$ as the time for which the average mutual information reaches four-fifths of its initial value:

$$\frac{I(T_d)}{I(0)} \approx \frac{4}{5}. \qquad (11)$$

The false nearest neighbours method is based on determining an acceptable minimum embedding dimension by looking at the behaviour of near neighbours under changes in the embedding dimension from $d$ to $d+1$. The most important assumption is that all points in the attractor that are close in $\mathbb{R}^m$ should be also close in $\mathbb{R}^{m+1}$. The false nearest neighbour is a point that appears to be a nearest neighbor because the embedding space is too small. When the number of false nearest neighbors arising through projection is zero in dimension $d_E$, the attractor has been unfolded in this dimension ([17]).
Assume that the dimension of space is *d*. The *r*-th nearest neighbor of *y*(*n*) is denoted by *y*<sup>*r*</sup>(*n*). The distance between point *y*(*n*) and it's *r*-th nearest neighbor is a square of the Euclidean distance.

$$R_d^2(n, r) = \sum_{k=0}^{d-1} [x(n + kT) - x^r(n + kT)]^2 \quad (12)$$

In going from dimension *d* to *d* + 1 by time delay embedding new coordinate *x*(*n* + *Td*) is added onto each delayed vectors *y*(*n*). The distance between points before and after adding new coordinates is now compared. A point is designated as a false neighbour when increase of distance is too large and the criteria below are fulfilled.

$$\left[ \frac{R_{d+1}^2(n, r) - R_d^2(n, r)}{R_d^2(n, r)} \right]^{\frac{1}{2}} = \left| \frac{x(n + Td) - x^r(n + Td)}{R_d(n, r)} \right| > R_{Tol} \quad (13)$$

where $R_{Tol}$ is some threshold.

Authors of the method ([17]), in numerical investigations, proved that for $R_{Tol} \ge 10$ the false nearest neighbors are clearly identified.

An acceptable minimum embedding dimension is chosen by looking at the percentage of false nearest neighbors during the addition of *d* + 1 components to the delayed vectors. When the percentage of false nearest neighbors drops to zero, a proper embedding dimension may be obtained.

## 3.2 The largest Lyapunov's exponent

Lyapunov's exponents examine the action of the dynamics defining the evolution of trajectories. The largest Lyapunov exponent describes the mean divergence between neighboring trajectories in the phase space by the following formula

$$d(t) = D e^{\lambda_1 t} \quad (14)$$

where *D* is the initial separation between neighboring points and $\lambda_1$ is the largest Lyapunov exponent. There are a few algorithms designed for the Largest Lyapunov's exponent estimation ([29], [23] and [16]). In this investigation Rosten's algorithm was used.

The first stage of the Rosten algorithm is time delay reconstruction. After reconstruction, for each point on the trajectory, the nearest neighbor $x_{\hat{j}}$ is found. This point minimizes the distance to the particular reference point, $x_j$ as described below

$$d_j(0) \approx \min_{x_{\hat{j}}} \|x_j - x_{\hat{j}}\| \quad (15)$$

where $d_j(0)$ is the initial distance from the *j*-th point to its nearest neighbor, and $\|\dots\|$ is the Euclidean norm.

From the definition of $\lambda_1$ in eq. (14) authors ([23]) assumed the *j*-th pair of nearest neighbors diverge approximately at a rate given by the largest Lyapunov exponent.

$$d_j(i) \approx D_j e^{\lambda_1 (i \Delta t)} \quad (16)$$
where $D_j$ is the initial separation. We can take the logarithm of both sides eq. (17)

$$ \ln(d_j(i)) \approx \ln(D_j) + \lambda_1(i\Delta t) \quad (17) $$

Eq. (17) represents a set of approximately parallel lines (for j = 1,2,...,M), each with a slope roughly proportional to $\lambda_1$. The largest Lyapunov exponent is calculated then by linear regression

$$ y(i) = \frac{1}{\Delta t} \langle \ln(d_j(i)) \rangle \quad (18) $$

where $\langle...\rangle$ denotes the average over all values of j.

## 4 Numerical results

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Stages of the nonlinear analysis procedure: mutual information, percentage of false nearest neighbors, logarithmic divergence reconstructed attractor</figcaption></figure>
The experiment involves the analysis of gait sequences which were recorded in the Human Motion Laboratory (HML) of the Polish-Japanese Institute of Information Technology by means of the Vicon Motion Kinematics Acquisition and Analysis System ([21], [15], [13], [8]).

Gait sequences were recorded in Euler angles. Six kinds of time series were recorded - movements of femurs, tibias and feet (left and right). The Experiment's aim was the investigation of local (associated with body parts) chaotic behaviour occurring in a human's gait. There were six investigations associated with movement of all the mentioned parts of the body's skeleton.

Fig. 1 illustrates results of successive stages of the nonlinear analysis procedure for the left femur: a) mutual information, b) percentage of false nearest neighbors, c) divergence, d) reconstructed attractor. Based on the above figure one can see that from the mutual information chart time delay embedding $T = 21$ could be obtained when the first minimum of $I(T)$ occurs. From the false nearest neighbours chart embedding dimension $d = 4$ is obtained when the percentage of nearest neighbours drops to zero. From the linear regression of the logarithmic divergence the largest Lyapunov exponent is estimated $\lambda_1 = 0.3157$.

Table 1 includes results of the nonlinear analysis procedure for all mentioned body parts. For each of them the time delay $T$, embedding dimension $d$ and the largest Lyapunov exponent have been estimated. For right foot and left tibia negative LLEs values have been observed. It could be caused by too large value of the estimated embedding dimension. However, negative value of the largest Lyapunov exponent does not undermine the possible chaotic nature of the investigated data. The positive value of the Largest Lyapunov exponent is only one indicator of the presence of deterministic chaos. The others are fractal structure of the reconstructed attractor or limited prediction of the investigated system.

<table><caption>**Table 1.** Results from nonlinear analysis procedure for each body part</caption><thead><tr><td>Body part</td><td>Time delay</td><td>Embedding dimension</td><td>LLE</td></tr></thead><tbody><tr><td>Left femur</td><td>21</td><td>4</td><td>0.3157</td></tr><tr><td>Right femur</td><td>20</td><td>4</td><td>0.0084</td></tr><tr><td>Left foot</td><td>15</td><td>4</td><td>0.0476</td></tr><tr><td>Right foot</td><td>16</td><td>5</td><td>-0.0206</td></tr><tr><td>Left tibia</td><td>33</td><td>7</td><td>-0.0450</td></tr><tr><td>Right tibia</td><td>23</td><td>5</td><td>0.1600</td></tr></tbody></table>
# 5 Conclusion

In this article, six time series captured from a large number of consecutive strides of gait were examined. Measurements came from three parts of the body: femur, foot and tibia. For each time series time delay reconstruction and largest Lyapunov's exponent estimation has been carried out. Values of all of the computed parameters for all examined time series are gathered in table 1.

Time delay embedding for time series has already been widely explored by various contributors. However, most of the work in the published literature concerns only scalar time series. In this paper the method for multivariate kinematic's time series is presented. The results are promising for practical applications in human gait analysis.

Analyzing time delays estimated for these time series, it can be stated that proper time delay is different for each time series. On the other hand all of the values are in the range [15; 33]. Embedding dimension values are from the range [4; 7]. However the most frequent values are 4 and 5. It is probable that the embedding dimension value for left tibia is a numerical error. One factor to consider is the averaging time delay and embedding dimension values in further investigations.

Most of the largest Lyapunov's exponents calculated values are positive. Due to this fact human gait kinematics data exhibit the properties of chaotic behavior. It could be also stated that during human locomotion, all mentioned body parts are highly sensitive to initial conditions.

**Acknowledgments.** The author is grateful for the reviews, which certainly helped raise the quality of the work. The author would like to acknowledge the Polish-Japanese Institute of Information Technology for the gait sequences recorded in the Human Motion Laboratory (HML).
Bibliography

[1] H. Abarbanel. Analysis of observed chaotic data. Springer-Verlag New York, 1996.

[2] B. R. Andrievskii and A. L. Fradkov. Control of chaos: Methods and applications. *Automation and Remote Control*, 64(5):673–713, 2003.

[3] L. Cao, A Mees, and K. Judd. Dynamics from multivariate time series. *Physica D*, 121:75–88, 1998.

[4] Martin Casdagli, Stephen Eubank, J Doyne Farmer, and John Gibson. State space reconstruction in the presence of noise. *Physica D: Nonlinear Phenomena*, 51(1):52–98, 1991.

[5] A. Cordoba, M. C. Lemos, and F. Jimnez-Morales. Periodical forcing for the control of chaos in a chemical reaction. *Journal of Chemical Physics*, 124(1):1–6, 2006.

[6] JB Dingwell, JP Cusumano, PR Cavanagh, and D Sternad. Local dynamic stability versus kinematic variability of continuous overground and treadmill walking. *Journal of biomechanical engineering*, 123(1):27–32, 2001.

[7] Jonathan B Dingwell and Laura C Marin. Kinematic variability and local dynamic stability of upper body motions when walking at different speeds. *Journal of biomechanics*, 39(3):444–452, 2006.

[8] Wiktor Filipowicz, Piotr Habela, Krzysztof Kaczmarski, and Marek Kulbacki. A generic approach to design and querying of multi-purpose human motion database. In *Computer Vision and Graphics*, pages 105–113. Springer, 2010.

[9] A.M. Fraser and H.L. Swinney. Independent coordinates for strange attractors from mutual information. *Physical Review A*, 1986.

[10] Sara P. Garcia and Jonas S. Almeida. Multivariate phase space reconstruction by nearest neighbor embedding with different time delays. *Phys. Rev. E*, 72, Aug 2005.

[11] P. Grassberger and I. Procaccia. Measuring the strangeness of strange attractors. *Physica D9*, 9, 1983.

[12] Peter Grassberger, Rainer Hegger, Holger Kantz, Carsten Schaffrath, and Thomas Schreiber. On noise reduction methods for chaotic data. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 3(2):127–141, 1993.

[13] Bartosz Jabłoński. Application of quaternion scale space approach for motion processing. In *Image Processing and Communications Challenges 3*, pages 141–148. Springer, 2011.

[14] Bartosz Jablonski. Quaternion dynamic time warping. *Signal Processing, IEEE Transactions on*, 60(3):1174–1183, 2012.

[15] Henryk Josiński, Agnieszka Michalczuk, Adam Świtoński, Romualda Mucha, and Konrad Wojciechowski. Quantifying chaotic behavior in treadmill walking. In *Intelligent Information and Database Systems*, pages 317–326. Springer, 2015.
[16] H. Kantz. A robust method to estimate the maximal lyapunov exponent of a time series. *Physics Letters A*, 185:77-87, 1994.

[17] Matthew B Kennel, Reggie Brown, and Henry DI Abarbanel. Determining embedding dimension for phase-space reconstruction using a geometrical construction. *Physical review A*, 45(6):3403, 1992.

[18] K Kocak, L Saylan, and J Eitzinger. Nonlinear prediction of near-surface temperature via univariate and multivariate time series embedding. *Eco-logical Modelling*, 173:1-7, 2004.

[19] Eric J Kostelich and Thomas Schreiber. Noise reduction in chaotic time-series data: a survey of common methods. *Physical Review E*, 48(3):1752, 1993.

[20] Jack B Kuipers. *Quaternions and rotation sequences*, volume 66. Princeton university press Princeton, 1999.

[21] Bogdan Kwolek, Tomasz Krzeszowski, and Konrad Wojciechowski. Real-time multi-view human motion tracking using 3d model and latency tolerant parallel particle swarm optimization. In *Computer Vision/Computer Graphics Collaboration Techniques*, pages 169-180. Springer, 2011.

[22] P. Parmananda. Controlling turbulence in coupled map lattice systems using feedback techniques. *Physical Review E*, 56(1):239-244, 1997.

[23] Colins J. J. Rossenstein M. T. and de Luca C. J. A practical method for calculating largest lyapunov exponents from small data sets. *Physica D*, 65:117-134, 1993.

[24] T. Sauer, J.A Yorke, and M. Casdagli. *Embedology*. Journal Of Statistical Physics, 1991.

[25] J. Stark. *Delay Embeddings for Forced Systems. I. Deterministic Forcing*. Journal of Nonlinear Science, New York, 1999.

[26] George Sugihara. Nonlinear forecasting for the classification of natural time series. *Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences*, 348(1688):477-495, 1994.

[27] F. Takens. *Detecting strange attractors in turbulence*. Springer-Verlag, Berlin, 1981.

[28] H. Whitney. *Differentiable manifolds*. Ann.Math., 1936.

[29] A. Wolf, J.B. Swift, H.L. Swinney, and J.A Vastano. Determining lyapunov exponents from a time series. *Physica 16D*, 16:285-317, 1985.