Received November 6, 2020, accepted November 21, 2020, date of publication November 27, 2020,
date of current version December 9, 2020.
Digital Object Identifier 10.1109/ACCESS.2020.3041158

Principal Motion Ellipsoids: Gait Variability
Index Invariant With Gait Speed
TOMOYUKI IWASAKI, SHOGO OKAMOTO , (Member, IEEE), YASUHIRO AKIYAMA ,(Member, IEEE), AND YOJI YAMADA, (Member, IEEE)
Department of Mechanical Systems Engineering, Nagoya University, Nagoya 464-8603, Japan

Corresponding author: Shogo Okamoto (shogo.okamoto@mae.nagoya-u.ac.jp)
This work was supported in part by the Japan Agency for Medical Research and Development (AMED) under Grant 19he2002003h0302,
and in part by the JSPS KAKENHI under Grant 19K21584.

ABSTRACT The walking motion of an individual involves considerable variability. We develop a gaitvariability index that determines how and to what degree repeated human gait motions vary, based on
generalized principal motion analysis (GPMA). Principal motion analysis (PMA) is an extension of principal
component analysis and decomposes multivariate time-series data, such as human joint angles during
walking, into a linear combination of several principal motion bases. The developed gait variability index
is defined by the size of an ellipsoid referred to as a PM ellipsoid and approximates the distribution of
repeated gait motions on the motion base space. We expand the method to compute PM ellipsoids using
GPMA, which enables us to mask the effects of a certain factor affecting gait motion observed under multiple
factors. We compute the principal motion bases invariant with the gait speed from the gait motions of nine
participants at two speed levels, 3.5 and 4.0 km/h. The sizes of the PM ellipsoids computed using GPMA
do not depend on the gait speed and exhibit good agreement with the MeanSD, i.e., a typical gait variability
index with correlation coefficients greater than 0.87.
INDEX TERMS Gait speed, modal analysis, walking motion, principal motion analysis.

I. INTRODUCTION

There is large variability in the walking motion of an individual. For example, the gait of a human varies dependingon the gait speed [1], [2] and step length [3], [4] during
walking. Some studies have discussed the gait variability inindividuals during curving motion [5], [6]. The MeanSD [7]and maximum Lyapunov exponent [8] are the most popular
indices for evaluating gait variability. MeanSD is the meanvalue of the standard deviation of parameters, such as joint
angles. The maximum Lyapunov exponent, which is also used
as a gait stability index, quantifies the average logarithmic
rate of the divergence of a system after perturbation. If thevalues of these indices are large, the gait variability in an individual is considerable. In addition, multiscale entropy [9] andrecurrence quantification analysis [10] have been reported ashighly reproducible indices [11]. Thus, gait variability has
been extensively studied. Indices to represent gait variability
are invaluable in enhancing commercial applications, such as
The associate editor coordinating the review of this manuscript andapproving it for publication was Zhong Wu

213330.

walking aids, and in the clinical scenario [12], [13] where
gait variability is an important characteristic of neurologicaldisorders [14]. However, as these indices are represented bysingle scalar quantities, they are unsuitable for discussing thevariation in walking motion, although they are suitable fordiscussing the level of gait variability.
Human motion has been analyzed through principal
motion analysis (PMA) [15], [16], also called principal component analysis for time-series [17]–[24] or spectral analysis [25]. This method decomposes human motion into bases
called principal motions or synergies without losing information on the interlocked multiple degrees of freedom of
redundant systems. It reduces the entire body motion toa few-dimensional parameter space, and has mainly been
applied for understanding motion synergies [17]–[25] andgenerating human and robot motions [15], [16], [26]–[28].
PMA has also been adapted to classify individual gait [29],[30] and sit-to-stand motion [31].
We had previously proposed the principal motion ellipsoid (PM ellipsoid) as an index for the gait variability [32].
This ellipsoid is an error ellipsoid calculated from theThis work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

VOLUME 8, 2020

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

distribution of the principal motion scores. Conventional
indices aim to represent the variability as a scalar quantity.
In contrast, the PM ellipsoid can analyze it as a multivariate
time-series. The variation in the walking motion of an individual can be determined by interpreting the obtained principal
motions [32]. Furthermore, the size of the PM ellipsoid represents the level of variability of individual motion and exhibitshigh correlation with the MeanSD.
Gait speed has considerable effect on the walking motion;
hence, it may not be possible to investigate individual walking
motion, if the gait speed is not suitably controlled. Someof the gait characteristic quantities vary significantly under
different gait speeds [4], [33]–[38]. PMA explores base functions that maximize the variability among all the samples.
Hence, when the gait speed is not controlled among all thesamples, the base functions are sensitive to not only theindividual characteristic walking motion but also the gait
speed. Therefore, PMA cannot be applied for identifying theprincipal motions in experiments involving uncontrolled gait
speed.
Generalized principal motion analysis (GPMA) is amethod for extracting the principal motions related to specificexperimental factors [31]. Gait motion comprises multiple
factors such as individual differences, gait speed, and repetition errors. GPMA can mask the effects of some of thesemotion factors, and extract the corresponding invariant principal motions. For example, the masking of repetitive errors
in sit-to-stand motion has been demonstrated, and individuals
within a group have been classified [31].
In this study, we develop a new PM ellipsoid based on
GPMA, combining [31] and [32]. As mentioned above,
the PM ellipsoid is advantageous because it can analyze thevariation in the walking motion and the phase at which it
occurs within the gait cycle. Furthermore, a GPMA-based
PM ellipsoid can mask specific factors among the motion
samples, such as the gait speed, whereas the afore-mentioned
indices [7]–[10], [32] do not aim to mask the effects ofcertain factors. For example, the maximum Lyapunov components [4], [36], variabilities in the step length, interval [34],and stride time [37], and MeanSD [37], [38], all of which areused as gait variability indices, are sensitive to the gait speed,
whereas the step-length variability may not be affected by
the gait speed [37], [39]. The existing popular gait variability
indices are sensitive to the gait speed. In contrast, PM ellipsoids based on GPMA can be employed as a gait variability
index invariant with the gait speed. However, when GPMAis applied, the principal motions may not be easy to interpretbecause of its strong masking ability. Therefore, in this study,
we extend GPMA to introduce a method for adjusting its
masking ability.

II. MOTION ANALYSIS METHODS
A. PRINCIPAL MOTION ANALYSIS (PMA)

Principal motion analysis (PMA) is an extension of principal component analysis, which deals with multivariate
VOLUME 8, 2020

time-series data [15]–[29]. It extracts covariation sets in thetemporal dimension, among multiple variables. Each set iscalled principal motion and the sets are linearly independentof each other. PMA is a linear analysis method, and it is relatively easy to interpret the meaning of the obtained principal
motions.
Suppose that a human motion at a certain instant is represented by p variables. For variable i (i = 1, . . . , p) at thek-th trial (k = 1, . . . , k 0 ), the time-series data vector θ i,kcomprising discretized u moments isθ i,k = (θi,k,1 , θi,k,2 , . . . , θi,k,u ).

(1)

Using this, we create an extended column vector xk (∈ Rpu×1 )including p variables:
xk = (θ 1,k , . . . , θ i,k , . . . , θ p,k )T .

(2)

(s)

Note that xc,k represents the time-series column vector ofexperimental participant s (s = 1, . . . , s0 ) at the k-th trial
under condition c (c = 1, . . . , c0 ). If we obtain k 0 sets ofwalking motion from each of the s0 participants under c0
conditions as the motion samples, the time-series data of all
0 0 0
the motions are represented by matrix X (∈ Rs c k ×pu ) as
follows:
(1)

(s0 )

(s)

X = (x1,1 , . . . , xc,k , . . . , xc0 ,k 0 )T .

(3)

From the eigenvector expansion of covariance matrix X,
we obtain
X ∼ YV T ,
(1)

(s0 )

(4)
0 0 0

where Y = (y1,1 , . . . , yc0 k 0 )T (∈ Rs c k ×r ) is a score matrix
that represents the quantity of information of each principal motion included by each observed motion and V =(v1 , . . . , vq , . . . , vr ) (∈ Rpu×r ) is a principal motion matrix
composed of r eigenvectors corresponding to the r largest
eigenvalues; vq (∈ Rpu×1 ) is the q-th principal motion vector,
whose eigenvalue is the q-th largest. Any motion can beapproximated by a linear combination of principal motion
vectors weighted by their corresponding scores:
(s)

(s)

xc,k ∼ V yc,k .

(5)

B. GENERALIZED PRINCIPAL MOTION ANALYSIS (GPMA)

We extended linear discriminant analysis to render it applicable to time-series data and propose generalized principal motion analysis (GPMA) [31]. GPMA includes basis
functions that maximize the variation between factors, while
reducing the variation within the factors. Therefore, thismethod can determine principal motions invariant with specific conditions such as the gait speed.
0 0
A matrix B (∈ Rs c ×pu ) which exhibits the effects ofthe masked conditions within individuals is represented asfollows:
(1)

(s0 )

0

(s)
(s ) T
B = (x̄1 − x̄(1) , . . . , x̄(s)c − x̄ , . . . , x̄c0 − x̄ ) .

(6)

(s)

Here, x̄c represents the mean motion of participant s undercondition c. Let the covariance matrices of X and B be
213331

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

S (= X T X) and D (= BT B), respectively. We obtain aprincipal motion matrix and score matrix from the eigenvector expansion of matrix D−1 S. However, some of theprincipal motions obtained through this calculation cannot
be interpreted because they are considerably affected by themasking effect of D. Therefore, we extend GPMA to improve
the interpretability of its principal motions. We obtain Dr bymodifying D as
Dr = D + 10α I.

(7)Here, I (∈ Rpu×pu ) is a unit matrix, and coefficient α isan arbitrary real number. We can determine the intensityof the masking functions through the value of α. From theeigenvector expansion of matrix D−1
r S, we obtain principal
motion matrix V g and score matrix Y g . When α is large, V g issimilar to V , and GPMA and PMA are practically the same.
On the other hand, when α is small, we can mask the effect(s)of specific factors. As well as PMA, xc,k is approximated bya linear combination of the principal motion vectors:
(s)

(s)

xc,k ∼ V g yg,c,k .

(8)III. GAIT VARIABILITY INDICES
A. PRINCIPAL MOTION ELLIPSOID (PM ELLIPSOID)

We use the size of the error (standard deviation) ellipsoid
calculated individually from the distribution of the principal
motion scores as the gait variability index (Fig. 1). Let thiserror ellipsoid be the principal motion ellipsoid (PM ellipsoid) [32]. The size of the PM ellipsoid represents the level ofmotion variability, whereas the direction represents how the(s)motion varies. Suppose qm represents the length of the major
axis m (m = 1, . . . , r) of the PM ellipsoid for participant s,
with r being the number of principal motions. We define thesize of the PM ellipsoid for participant s as e(s) using thelength of each axis, as follows:
(s)

e

=

r
X

q(s)
m.

(9)

i=1

Fig. 1 depicts the PM ellipsoid for participant s on thefirst–second principal motion plane as an example. Because
the direction of the PM ellipsoid is nearly along the firstprincipal motion, we can determine that the walking motionof participant s mainly varies for the first principal motion.
Note that this ellipsoid is defined in an r-dimensional space.
B. MeanSD

The afore-mentioned gait variability indices evaluate different aspects of the walking motion. Hence, the correlation
coefficients among them are small or null [7], [11], [40].
For example, the MeanSD evaluates the variability based on
the standard deviation of the motion, whereas the maximum
Lyapunov exponent evaluates the motion stability; hence,
they are regarded independent. As there is high correlation
between the PM ellipsoid and MeanSD [32], we use MeanSD
as the reference index to investigate whether the PM ellipsoidcomputed using GPMA is valid as a gait variability index.
213332

FIGURE 1. Definition of the size of the PM ellipsoid for participant s. The
error ellipsoid (red) is calculated using the distributed samples(black dots).

MeanSD, which is the mean value of the standard deviation
based on the gait cycle, has been used for the joint angleof the lower limbs [7], [37], [41], and for the velocity [2],[36], [37], [42] and acceleration [7], [40] obtained throughan accelerometer attached to the waist.
For a participant, SD(l) represents the standard deviation at
moment l on the gait cycle of a variable. Then, the MeanSD,
which represents the gait variability for the entire gait cycleof the variable, is defined as
MeanSD = hSD(l) i, l ∈ {0, 1, . . . , 100% gait cycle},(10)

where h i is the mean value of all the temporal points [7].
We calculated (10) for the six joint angles mentioned inSec. 4.3, and defined their mean value as the MeanSD of thejoint angle.
IV. WALKING MOTION EXPERIMENT
A. EXPERIMENTIn this study, we analyze the walking data measured by
Ullauri et al. [43]; the relevant details can be found in thementioned reference. A walking experiment was conductedon nine adult males (height: 172.5 ± 0.05 cm, weight:
63.6 ± 7.2 kg) without any neurological and musculoskeletal
abnormalities.
B. METHOD

The gait motion was recorded at a frequency of 100 Hz using
a motion capture system (MAC 3D system, Motion Analysis
Corporation, U.S.) with ten cameras (nominal position resolution of submillimeter). Thirty-six reflective markers wereattached on the entire body of the subject, as recommended
by the motion module of the software for interactive musculoskeletal modeling (SIMM, Musculographics Inc., U.S.)
used for calculating the joint angles. The participants walked
on a treadmill with a 1.4×0.5 m2 walking surface (OMEGA3,
Johnson Health Tech Co., Taiwan).
Two walking trials with the different walking speed wererecorded for each subject. To avoid exhaustion, the durationof each trial was set to 6 min. A speed of 1.11 m/s (4 km/h),
VOLUME 8, 2020

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

regarded as the natural walking speed of young adults, wasselected as the regular walking speed, and a speed of 0.97 m/s(3.5 km/h) was selected as the slow speed.

TABLE 1. Contribution ratio (%) of the first–third principal motions of
PMA, GPMA2 , and GPMA−7 . PM: principal motion.

C. MOTION DATA ANALYSIS

We visually inspected the motion samples used in [43] toconfirm that there were no apparent erroneous gait motions.
The first minute of all the gait motions in each trial was
excluded from the analysis. All the data were smoothed using
a 6-Hz Butterworth filter. From a single trial, we extracted
100 motions, defined as a gait cycle from the left-heel contact
to the next left-heel contact in continuous walking. The gait
cycle was normalized to 0–100% (u = 101) for aligning
the data length of each trial (Fig. 2). The values of the joint
angles were positive in the direction of flexing and negativein the direction of extending from the basic standing posture.
We analyzed the hip flexion, knee flexion, and ankle dorsiflexion angles for the left and right legs. We defined the hip
flexion angle as the angle between the vertical line passing
through the hip joint and the straight line connecting the hipand knee joints.

FIGURE 2. Gait cycle (defined from the left-heel contact (0%) to the next
left-heel contact (100%)) normalized to 0–100%.

PMA (left in Fig. 3), whereas they were almost matched for
GPMA−7 . Those of GPMA2 were moderately apart, between
two gait-speed levels. These figures demonstrate that GPMA
masks the gait-speed factor within individuals and α value
adjusts the level of the masking function.
We investigated whether each principal motion was related
to the gait speed through a t-test for the centroids of the score
distribution between two gait speeds. The test statistics andp values are shown in Table 2. The first principal motion
indicated no significant differences between the gait speeds
in any analysis methods. For the second principal motion,
the PMA scores when the participants walked at 4.0 km/h
were greater than those at 3.5 km/h, which was not the casefor GPMA. For the third principal motion, PMA exhibited
marginal differences in the scores between the two gait-speed
levels, i.e., 3.5 km/h > 4.0 km/h, whereas the GPMA did
not. Hence, the second and third principal motions of PMA
depend on the gait speeds; however, those of GPMA donot, suggesting that GPMA extracts the principal motions
invariant with speed.
TABLE 2. Test statistics and p values computed from the centroids of thescore distribution between two gait speeds for each principal motion. If t
is positive, the scores for 4.0 km/h are greater than those for 3.5 km/h.

V. RESULTS
A. EFFECT OF GAIT SPEED ON THE PRINCIPAL MOTION
SCORE

We confirmed the expression of the gait speed factor fromthe distribution of the principal motion score, within individuals. Table 1 presents the contribution ratio of the first–thirdprincipal motions in each motion analysis method. This ratio
indicates the extent of variation in the samples for eachprincipal motion. Even in earlier studies, gait motions havebeen expressed up to a few principal motions [19], [20], [24].
Fig. 3 displays the distribution of the principal motion scoreobtained using PMA (left), GPMA2 (middle), and GPMA−7(right). GPMAα denotes GPMA with coefficient α. The
reason for determining the value of α will be explainedin Sec. V-C. The figures on top depict the first–second principal motion planes, whereas those on the bottom depict
the second–third principal motion planes. The red and blue
ellipses are the PM ellipsoids for gait speeds of 3.5 km/h and
4.0 km/h, respectively. The centroids of the sample distribution are represented as a marker for each participant.
The distances among the centroids in the score distribution between two gait-speed levels were maximum for
VOLUME 8, 2020

To determine the speed dependence of MeanSD, we compared the MeanSD values between two speed levels. The
mean and standard deviation of MeanSD among all the participants were 1.73 ± 0.23 deg and 1.80 ± 0.32 deg, for 3.5 and4.0 km/h, respectively. There were no significant differences(t(8) = 0.55, p = 0.30) indicating that the MeanSD values
were insensitive to the gait speed level in our data. This is
not fully consistent with previous studies [37], [38], where
MeanSD is sensitive to the gait speed. We discuss the reason
for this in the first paragraph of Sec. VI.
B. CONCURRENT VALIDITY AS THE GAIT VARIABILITY
INDEX

We compared the size of the PM ellipsoid obtained using
each analysis method with the MeanSD, which is a representative gait variability index, and determined the concurrent validity of the PM ellipsoid as a gait variability index.
213333

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

FIGURE 3. Scatter plots of the first–third principal motion scores of PMA (left), GPMA2 (middle), and GPMA−7 (right). Top: first–second principal motion
planes; Bottom: second–third principal motion planes.

Table 3 represents the correlation coefficients between these
indices. High correlation can be observed among all the cases,
indicating that the MeanSD and size of the PM ellipsoid aresimilar indices for the degree of gait variability.

TABLE 4. Relative error of the principal motion ellipsoids (d ) between
two gait speeds. The values for GPMA2 are comparable to those for PMA(p = 0.154), whereas those for GPMA−7 are slightly smaller than thosefor PMA (p = 0.043).

TABLE 3. Correlation coefficients between the MeanSD and size of the
PM ellipsoid for PMA, GPMA2 , and GPMA−7 at two gait speeds.

C. COMPARISON OF THE PMA- AND GPMA-ELLIPSOIDS

We verified the level to which the effect of the gait speed factor was masked for the gait variability level. We determinedthat the size of the PM ellipsoid, when participant s walked at(s)a gait speed of c was ec . Using the relative error, we defined
an index representing the level to which the gait speed factor
affects the size of the PM ellipsoid within individuals for eachmethod as follows:(s)d (s) =

e3.5
(s)

e4.0

−1 .

(11)

When d (s) is close to zero, the effect of the gait speed ismasked. Table 4 lists d computed for an individual, for eachmethod. The mean values and standard errors of d of the nine
213334

participants for PMA, GPMA2 , and GPMA−7 are 0.17 ±0.02, 0.14 ± 0.03 and 0.10 ± 0.03, respectively. There areno significant differences in ds between PMA and GPMA2(t(8) = −1.58, p > 0.05), whereas d for GPMA−7 tends
to be lesser than that for PMA (t(8) = −2.40, p < 0.05,
Bonferroni correction). Therefore, the PM ellipsoid based on
GPMA−7 is less sensitive to the gait speed factor for the gaitvariability level.
Next, we verified the level to which the effect of the gait
speed factor was masked in the direction of the PM ellipsoidby comparing the angles between the major axes of the
VOLUME 8, 2020

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

PM ellipsoids within individuals. We compared them on thefirst–second principal motion plane (φ1−2 ) and second–thirdprincipal motion plane (φ2−3 ). Table 5 presents φ1−2 and
φ2−3 computed through each method. The mean values andstandard errors of φ1−2 of the nine participants for PMA,
GPMA2 , and GPMA−7 are 20.0 ± 6.77, 27.1 ± 7.37,, and12.2 ± 2.24, respectively, and those of φ2−3 are 6.58 ± 1.94,
12.2 ± 6.87, and 26.3 ± 6.53, respectively. There are no
significant differences among these values except for φ2−3
between PMA and GPMA−7 (t(8) = −3.10, p < 0.05).
The masking ability of GPMA did not affect the PM ellipsoid
direction.
TABLE 5. Differences in the degrees of the principal motion ellipsoids
between two gait speeds, calculated on the first–second principal motion
plane (φ1−2 ) and second–third principal motion plane (φ2−3 ).

the GPMA-based ellipsoid do not depend on the gait speed,in contrast to those based on PMA. These results indicate that
the PM ellipsoids computed using GPMA function as a gaitvariability index invariant with the gait speed, whereas theother popular gait variability indices are influenced by the
gait speed, as previously mentioned in Sec. I. Nonetheless,the MeanSD values did not differ among the two speed levels
for our gait samples. The possible reason is as follows: The
gait variability index can be a quadratic U-shaped functionof the gait speed with minimum values at the potentially
preferred speed [38], [44]–[46]. The gait speeds in this study
were within the range of the participants’ preferred speed,
and the speed-dependency of the MeanSD may have been
concealed.
A. INTERPRETATION OF THE PRINCIPAL MOTIONS

We may understand the variation in the walking motion
within individuals by interpreting each principal motion.
Fig. 5 shows the time-series change of the load of thefirst–third principal motions obtained using PMA, GPMA2 ,and GPMA−7 . As the load represents the flexion of the joint
angles, if the load is zero, the flexion angle is zero, and if theload is positive or negative, the joint angle flexes or extends.
Further, Fig. 6 shows the average joint motions during a gait
cycle, which helps in interpreting the principal motions.
1) PRINCIPAL MOTION ANALYSIS

Fig. 4 depicts the mean values of d for the nine participants,
when α is varied by 0.1. It can be observed that the meanvalues of d vary considerably for 0 < α < 4. The mean valueof d for GPMA7 is 0.165, indicating that GPMA7 has a low
masking ability similar to PMA with d = 0.167. The meanvalue of d for GPMA−7 is 0.10. The mean value of d is 0.13,
when α is approximately two, indicating that the masking
ability of GPMA2 is between those of PMA and GPMA−7 .

FIGURE 4. Relationship between the mean value of d (relative error of
the ellipsoidal sizes between 3.5 km/h and 4.0 km/h) and masking
coefficient α.

VI. DISCUSSION

In Sec. V, we demonstrated that the PM ellipsoids based on
GPMA and the MeanSD, a popular gait variability index,
exhibited correlation coefficients greater than 0.87. Furthermore, the principal motion scores of GPMA and the size of
VOLUME 8, 2020

Fig. 5 (a), which is the first principal motion of PMA,
indicates that the knee joints flex (positive), hip joints flex
(positive), and ankle joints plantarflex (negative) during thefirst half of the swing phase of both legs. These represent
typical walking motion because they resemble the average
joint motion in Fig. 6. Therefore, we can interpret the firstprincipal motion as average walking motion or the degree ofentire gait motion. The variance of the gait samples along thefirst principal motion is greater than those along the otherprincipal motions suggesting that the variance of walking
motion is mainly the degree of the entire joint motion.
Fig. 5 (b), which is the second principal motion, shows thatthe hip and knee joints flex (positive), and the ankle joints
plantarflex (negative) during the first half of the swing phase,whereas the knee joints extend (negative) during the second
half of the swing phase for both legs. This motion can beconsidered to represent the advance of the gait phase because
it appears as the differential motion of the average walking
motion in Fig. 6. Fig. 7 depicts the knee flexion angles of gait
samples with the maximum and minimum scores for the second principal motion, for a participant. Gait samples withthe maximum scores for the second principal motion shift tothe swing phase earlier than the minimum ones. Therefore,
the second principal motion also represents walking motion
with long swing phases. As shown in Table 2 in Sec. 5.1,
the second principal motion is related to the gait speed, andgait samples at 4.0 km/h have larger scores for this principal
motion than those at 3.5 km/h. People, who walk faster, have
longer swing phases [19], [43].
213335

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

FIGURE 5. Load of each principal motion of PMA (left), GPMA2 (middle), and GPMA−7 (right) computed using the joint angles. Top: first principal motion;
Middle: second principal motion; Bottom: third principal motion.

FIGURE 6. Average joint motion during a gait cycle.

Fig. 5 (c), which is the third principal motion, shows thatthe knee joint of the stance leg flexes (positive) and the swing
leg extends (negative) at heel contact, after which the hip
joint of the stance leg flexes (positive) and the swing leg
extends (negative) for both legs. As in Table 2, this principal
motion is related to the gait speed. Similar to [43], our data
indicates that the gait-speed factor affects the length of the213336

FIGURE 7. Gait samples with maximum (blue line) and minimum (blackline) scores for the second principal motion of PMA. Solid line: left-knee
angle; Dotted line: right-knee angle.

swing phase and step length. In PMA, as mentioned in theprevious paragraph, the second principal motion represents
the length of the swing phase. Hence, it is considered thatthe third principal motion can be related to the step length.
People, who walk fast, have large step lengths as per our
VOLUME 8, 2020

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed

data [43]. Participants with large scores for this principal
motion tend to walk slowly, as depicted in Table 2, and thescores of this principal motion are moderately correlated withthe step length (r = −0.317, t = 14.96, p < 0.001). Samples
with the greater scores tend to involve slower motions withsmaller step lengths.
2) GENERALIZED PRINCIPAL MOTION ANALYSIS (α = 2)

Here, we interpret the principal motions obtained using
GPMA2 . Although it is occasionally difficult to interpret the principal motions of GPMA because they contain high-frequency components, we may be able to infer
their meanings from similar principal motions of PMA.
We smoothened the principal motions using a Gaussian filter
with a length of ten for visual clarity in Fig. 5.
Fig. 5 (d), which is the first principal motion, resembles thefirst principal motion of PMA. Therefore, it is considered torepresent average walking motion.
Fig. 5 (e), which is the second principal motion, showsthat the left knee joint flexes (positive) during the first halfof the stance phase, becomes zero at approximately 50% of
the gait cycle, and extends in the swing phase. The rightknee joint is around zero during the first half of the swing
phase, extends (negative) during 20%–50% of the gait cycle,
and flexes (positive) during the stance phase. The hip joints
flex (positive) in the stance phase and extend (negative orzero) in the swing phase of both legs. This shows the left-rightasymmetry of walking, suggesting that the double stancephase of right-heel contact—left-toe off is longer than that ofleft-heel contact—right-toe off. Fig. 8 shows the knee flexion
angles of two gait samples with the maximum and minimumscores, respectively, of the second principal motion for acertain participant. In this example, the gait sample with the
maximum score for the second principal motion (blue curves)
exhibits peak left and right knee flexion angles at 18% and72%, respectively. The difference in phase between these two
peaks is 54%. In contrast, the gait sample with the least score(black curves) exhibits peak left and right knee angles at 21%

FIGURE 8. Gait samples with maximum (blue line) and minimum (blackline) scores for the second principal motion of GPMA2 . Solid line:
left-knee angle; Dotted line: right-knee angle.
VOLUME 8, 2020

and 67%, respectively, resulting in a phase difference of 46%.
As such, the second principal motion pertains to the left-rightasymmetry of motion.
Fig. 5 (f), which is the third principal motion, is difficult tointerpret because of its complexity, but it resembles the thirdprincipal motion of PMA, relatively. However, as in Table 2,
the third principal motion of GPMA2 is invariant with the gait
speed, whereas that of PMA is sensitive to the gait speed.
3) GENERALIZED PRINCIPAL MOTION ANALYSIS (α = −7)

The principal motions obtained through GPMA−7 are considerably different from those of PMA and GPMA2 , and are
difficult to interpret.
The first principal motion (Fig. 5 (g)) resembles that of
PMA and GPMA2 , relatively. Therefore, it may represent
average walking motion.
The second principal motion of GPMA−7 (Fig. 5 (h))
resembles that of GPMA2 , partially. From this point of view,it may represent the left-right asymmetry of the double stancephases.
The third principal motion of GPMA−7 (Fig. 5 (i)) does notresemble any of the principal motions of PMA and GPMA2 .
Further, the joint profiles differ among the first and second halves of the gait cycle, indicating laterally asymmetricmotion. However, this principal motion could not be directlyinterpreted.
B. LIMITATIONS OF PM ELLIPSOIDS AND THISPRESENT STUDY

As indicated in Sec. V-B, the MeanSD and size of the
PM ellipsoid are similar indices for the degree of gait variability. As mentioned in Sec. III-B, some of the representative
gait variability indices do not exhibit high correlation coefficients with each other because they are based on inherently
different aspects of walking motions [7], [11], [40]. Hence,
it is considered that the correlation coefficient between thesize of the PM ellipsoid and the other gait variability indices,
including the maximum Lyapunov component, is small. Thisis a common problem in gait variability indices.
PM ellipsoids consider the motion repetition error as variability. When PMA or GPMA is applied to motions, except
walking, the motion start and finish points need to be clearly
defined, in continuously recorded motions. Thus, aperiodic
motion or motion with no apparent cycle is unsuitable foranalysis. Further, although this study targeted the gait variability of healthy participants, it is important to investigate
whether PM ellipsoids are available for gait with abnormality
because gait fluctuation is the key for analyzing such gait
patterns [12]–[14]. This remains to be studied in future. Similarly, although we expect that our method is applicable toextremely slow or fast walking motions, this study considers
only the typical gait speed range; this aspect also needs to be
demonstrated in future.
As shown in Fig. 3, GPMA can mask the effect of the
gait speed factor within individuals because the distancesamong the centroids of the scores between two gait speed
213337

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speedlevels are smaller than those of PMA. Furthermore, as Table 4shows, GPMA can mask the effect of gait speed on the degreeof individual gait variability. On the other hand, there wereno significant differences between PMA and GPMA for the
direction of the PM ellipsoid, as indicated by Table 5. In thisstudy, as in (6), we defined the effect of the gait speed as thedifference between the mean motions at 3.5 and 4.0 km/h.
Hence, the difference between the mean motions for the two
gait speed levels can be masked by GPMA, but it does notmask the effect of the gait speed level on the directions of the
PM ellipsoids.
As previously mentioned, GPMA-based PM ellipsoids
employ the difference between the mean motions of the two
gait speed levels to mask the effect of the gait speed factor.
If the gait speeds are discrete, as in 3.5 km/h and 4.0 km/h,
we can analyze the variability between them; however, thiscannot be done, if the gait speeds vary continuously.
VII. CONCLUSION

In this study, we developed a new index for gait variability
using GPMA, which is an extended PMA method, that canmask the effect of a certain factor for computing the principal
motions. This index is based on the variability ellipsoidsof motion samples in the principal motion space. Hence,
it provides the magnitude and properties of motion variability, unlike previous gait variability indices. Furthermore,
we modified the previously proposed GPMA by introducing afunction to adjust the masking ability of GPMA. Using PMA
and GPMA, we computed the PM ellipsoids of the walking
motions recorded at two gait-speed levels. The GPMA ellipsoids were less sensitive to the gait speed than the PMA ellipsoids, as indicated by the distance among the centroids of thescore distribution between the gait speeds and the size of the
PM ellipsoid within individuals. Ellipsoids computed using
PMA and GPMA exhibited high correlation coefficients with
the MeanSD, which is a representative gait variability index.
Furthermore, we demonstrated that the interpretation of theprincipal motions can be easier by adjusting the masking
capability of GPMA. A PM ellipsoid based on GPMA is a
useful variability index because it indicates the degree of gaitvariability and the variation in individual gait motion, withthe ability to mask the effect of experimental factors such as
the gait speed.
REFERENCES[1] S. A. England and K. P. Granata, ‘‘The influence of gait speed on local
dynamic stability of walking,’’ Gait Posture, vol. 25, no. 2, pp. 172–178,
Feb. 2007.
[2] S. M. Bruijn, J. H. van Dieën, O. G. Meijer, and P. J. Beek, ‘‘Is slow walking
more stable?’’ J. Biomech., vol. 42, no. 10, pp. 1506–1512, Jul. 2009.
[3] T. M. Owings and M. D. Grabiner, ‘‘Variability of step kinematics in young
and older adults,’’ Gait Posture, vol. 20, no. 1, pp. 26–29, Aug. 2004.
[4] L. Hak, H. Houdijk, P. J. Beek, and J. H. van Dieën, ‘‘Steps to take toenhance gait stability: The effect of stride frequency, stride length, andwalking speed on local dynamic stability and margins of stability,’’ PLoS
ONE, vol. 8, no. 12, Dec. 2013, Art. no. e82842.
[5] Y. Akiyama, H. Toda, T. Ogura, S. Okamoto, and Y. Yamada, ‘‘Classification and analysis of the natural corner curving motion of humans based on
gait motion,’’ Gait Posture, vol. 60, pp. 15–21, Feb. 2018.
213338[6] Y. Akiyama, S. Okamoto, H. Toda, T. Ogura, and Y. Yamada, ‘‘Gait motion
for naturally curving variously shaped corners,’’ Adv. Robot., vol. 32, no. 2,
pp. 77–88, Jan. 2018.
[7] J. B. Dingwell, J. P. Cusumano, P. R. Cavanagh, and D. Sternad, ‘‘Local
dynamic stability versus kinematic variability of continuous overgroundand treadmill walking,’’ J. Biomechanical Eng., vol. 123, no. 1, pp. 27–32,
Feb. 2001.
[8] J. B. Dingwell and J. P. Cusumano, ‘‘Nonlinear time series analysis ofnormal and pathological human walking,’’ Interdiscipl. J. Nonlinear Sci.,
vol. 10, no. 4, pp. 848–863, Dec. 2000.
[9] M. Costa, C.-K. Peng, A. L. Goldberger, and J. M. Hausdorff, ‘‘Multiscale
entropy analysis of human gait dynamics,’’ Phys. A, Stat. Mech. Appl.,
vol. 330, nos. 1–2, pp. 53–60, Dec. 2003.
[10] F. Sylos Labini, A. Meli, Y. P. Ivanenko, and D. Tufarelli, ‘‘Recurrence
quantification analysis of gait in normal and hypovestibular subjects,’’ Gait
Posture, vol. 35, no. 1, pp. 48–55, Jan. 2012.
[11] F. Riva, M. C. Bisi, and R. Stagni, ‘‘Gait variability and stability measures:
Minimum number of strides and within-session reliability,’’ Comput. Biol.
Med., vol. 50, pp. 9–13, Jul. 2014.
[12] S. Qiu, L. Liu, Z. Wang, S. Li, H. Zhao, J. Wang, J. Li, and K. Tang, ‘‘Body
sensor network-based gait quality assessment for clinical decision-support
via multi-sensor fusion,’’ IEEE Access, vol. 7, pp. 59884–59894, 2019.
[13] Y. Yan, O. M. Omisore, Y. Xue, H. Li, Q. Liu, Z. Nie, J. Fan, andL. Wang, ‘‘Classification of neurodegenerative diseases via topologicalmotion analysis—A comparison study for multiple gait fluctuations,’’
IEEE Access, vol. 8, pp. 96363–96377, 2020.
[14] Y. Moon, J. Sung, R. An, M. E. Hernandez, and J. J. Sosnoff, ‘‘Gait
variability in people with neurological disorders: A systematic review andmeta-analysis,’’ Hum. Movement Sci., vol. 47, pp. 197–208, Jun. 2016.
[15] F. C. Park and K. Jo, ‘‘Movement primitives and principal component analysis,’’ in On Advances in Robot Kinematics. Dordrecht, The Netherlands:
Springer, 2004, pp. 97–106.
[16] B. Lim, S. Ra, and F. C. Park, ‘‘Movement primitives, principal component
analysis, and the efficient generation of natural motions,’’ in Proc. IEEE
Int. Conf. Robot. Autom., Apr. 2005, pp. 4630–4635.
[17] Z.-M. Li and J. Tang, ‘‘Coordination of thumb joints during opposition,’’J. Biomechan., vol. 40, no. 3, pp. 502–510, Jan. 2007.
[18] N. St-Onge and A. G. Feldman, ‘‘Interjoint coordination in lower limbs
during different movements in humans,’’ Exp. Brain Res., vol. 148, no. 2,
pp. 139–149, Jan. 2003.
[19] N. A. Borghese, L. Bianchi, and F. Lacquaniti, ‘‘Kinematic determinantsof human locomotion,’’ J. Physiol., vol. 494, no. 3, pp. 863–879, 1996.
[20] C. D. Mah, M. Hulliger, R. G. Lee, and I. S. O’Callaghan, ‘‘Quantitative
analysis of human movement synergies: Constructive pattern analysis for
gait,’’ J. Motor Behav., vol. 26, no. 2, pp. 83–102, Jun. 1994.
[21] T. Funato, S. Aoi, H. Oshima, and K. Tsuchiya, ‘‘Variant and invariant
patterns embedded in human locomotion through whole body kinematiccoordination,’’ Exp. Brain Res., vol. 205, no. 4, pp. 497–511, Sep. 2010.
[22] S. Vernazza-Martin, N. Martin, and J. Massion, ‘‘Kinematic synergies andequilibrium control during trunk movement under loaded and unloaded
conditions,’’ Exp. Brain Res., vol. 128, no. 4, pp. 517–526, Oct. 1999.
[23] T. Funato, S. Aoi, N. Tomita, and K. Tsuchiya, ‘‘Human gait control suggested by the evaluation of the fluctuation of synergy,’’ in Proc. IEEE/SICE
Int. Symp. Syst. Integr. (SII), Dec. 2011, pp. 267–272.
[24] H. Oshima, S. Aoi, T. Funato, N. Tsujiuchi, and K. Tsuchiya, ‘‘Variant andinvariant spatiotemporal structures in kinematic coordination to regulate
speed during walking and running,’’ Frontiers Comput. Neurosci., vol. 13,
p. 63, Sep. 2019.
[25] K. Fujii, N. Takeishi, B. Kibushi, M. Kouzaki, and Y. Kawahara, ‘‘Datadriven spectral analysis for coordinative structures in periodic humanlocomotion,’’ Sci. Rep., vol. 9, no. 1, pp. 1–14, Dec. 2019.
[26] K. Mishima, S. Kanata, H. Nakanishi, T. Sawaragi, and Y. Horiguchi,‘‘Extraction of similarities and differences in human behavior using singular value decomposition,’’ IFAC Proc. Vols., vol. 43, no. 13, pp. 436–441,
2010.
[27] H. Kang and F. C. Park, ‘‘Humanoid motion optimization via nonlinear
dimension reduction,’’ in Proc. IEEE Int. Conf. Robot. Autom., May 2012,
pp. 1444–1449.
[28] S. Morishima, K. Ayusawa, E. Yoshida, and G. Venture, ‘‘Whole-body
motion retargeting using constrained smoothing and functional principle
component analysis,’’ in Proc. IEEE-RAS 16th Int. Conf. Humanoid Robots(Humanoids), Nov. 2016, pp. 294–299.
VOLUME 8, 2020

T. Iwasaki et al.: PM Ellipsoids: Gait Variability Index Invariant With Gait Speed[29] T. Akiduki, K. Kawamura, Z. Zhang, and H. Takahashi, ‘‘Extraction andclassification of human gait features from acceleration data,’’ Int. J. Innov.
Comput., Inf. Control, vol. 14, no. 4, pp. 1361–1370, 2018.
[30] T. Akiduki, A. Uchida, Z. Zhang, T. Imamura, and H. Takahashi, ‘‘Extraction of human gait feature from acceleration data,’’ ICIC Exp. Lett., B,
Appl., vol. 7, no. 3, pp. 649–656, 2016.
[31] T. Iwasaki, S. Okamoto, Y. Akiyama, and Y. Yamada, ‘‘Generalized principal motion analysis: Classification of Sit-to-Stand motions,’’ in Proc. IEEE
8th Global Conf. Consum. Electron. (GCCE), Oct. 2019, pp. 679–681.
[32] T. Iwasaki, S. Okamoto, Y. Akiyama, and Y. Yamada, ‘‘Principal motion
ellipsoids: Gait variability index based on principal motion analysis,’’ inProc. IEEE/SICE Int. Symp. Syst. Integr. (SII), Jan. 2020, pp. 489–494.
[33] J. M. Hausdorff, M. E. Cudkowicz, R. Firtion, J. Y. Wei, and
A. L. Goldberger, ‘‘Gait variability and basal ganglia disorders: Strideto-stride variations of gait cycle timing in parkinson’s disease and Huntington’s disease,’’ Movement Disorders, vol. 13, no. 3, pp. 428–437,
May 1998.
[34] K. Jordan, J. H. Challis, and K. M. Newell, ‘‘Walking speed influences on
gait cycle variability,’’ Gait Posture, vol. 26, no. 1, pp. 128–134, Jun. 2007.
[35] M. D. Latt, H. B. Menz, V. S. Fung, and S. R. Lord, ‘‘Walking speed,
cadence and step length are selected to optimize the stability of headand pelvis accelerations,’’ Exp. Brain Res., vol. 184, no. 2, pp. 201–209,
Nov. 2007.
[36] J. B. Dingwell and L. C. Marin, ‘‘Kinematic variability and local dynamic
stability of upper body motions when walking at different speeds,’’
J. Biomech., vol. 39, no. 3, pp. 444–452, Jan. 2006.
[37] H. G. Kang and J. B. Dingwell, ‘‘Separating the effects of age and walking
speed on gait variability,’’ Gait Posture, vol. 27, no. 4, pp. 572–577,
May 2008.
[38] M. Yamasaki, T. Sasaki, and M. Torii, ‘‘Sex difference in the pattern oflower limb movement during treadmill walking,’’ Eur. J. Appl. Physiol.
Occupational Physiol., vol. 62, no. 2, pp. 99–103, 1991.
[39] J. L. Helbostad and R. Moe-Nilssen, ‘‘The effect of gait speed on lateralbalance control during walking in healthy elderly,’’ Gait Posture, vol. 18,
no. 2, pp. 27–36, Oct. 2003.
[40] P. Terrier and O. Dériaz, ‘‘Kinematic variability, fractal dynamics and local
dynamic stability of treadmill walking,’’ J. NeuroEngineering Rehabil.,
vol. 8, no. 1, p. 12, Dec. 2011.
[41] D. H. Gates, B. J. Darter, J. B. Dingwell, and J. M. Wilken, ‘‘Comparisonof walking overground and in a computer assisted rehabilitation environment (CAREN) in individuals with and without transtibial amputation,’’J. NeuroEngineering Rehabil., vol. 9, no. 1, p. 81, 2012.
[42] M. Punt, S. M. Bruijn, H. Wittink, and J. H. van Dieën, ‘‘Effect of arm
swing strategy on local dynamic stability of human gait,’’ Gait Posture,
vol. 41, no. 2, pp. 504–509, Feb. 2015.
[43] J. B. Ullauri, Y. Akiyama, S. Okamoto, and Y. Yamada, ‘‘Technique toreduce the minimum toe clearance of young adults during walking tosimulate the risk of tripping of the elderly,’’ PLoS ONE, vol. 14, no. 6,
Jun. 2019, Art. no. e0217336.
[44] F. Danion, E. Varraine, M. Bonnard, and J. Pailhous, ‘‘Stride variability inhuman gait: The effect of stride frequency and stride length,’’ Gait Posture,
vol. 18, no. 1, pp. 69–77, Aug. 2003.
[45] J. M. Hausdorff, ‘‘Stride variability: Beyond length and frequency,’’ Gait
Posture, vol. 20, no. 3, p. 304, Dec. 2004.
[46] J. M. Hausdorff, ‘‘Gait variability: Methods, modeling and meaning,’’
J. NeuroEng. Rehabil., vol. 2, p. 19, Dec. 2005.

VOLUME 8, 2020

TOMOYUKI IWASAKI received the B.S. degreefrom Nagoya University, in 2019, where he is currently pursuing the degree with the Department of
Mechanical Systems Engineering.

SHOGO OKAMOTO (Member, IEEE) received
the Ph.D. degree in information sciences fromTohoku University, in 2010. Since 2010, he has
been with Nagoya University. He is currently
an Associate Professor with the Department of
Mechanical Systems Engineering, Nagoya University. His research interests include haptics, assistive robotics, and affective engineering.

YASUHIRO AKIYAMA (Member, IEEE) received
the B.E. degree in engineering from the Tokyo
Institute of Technology, Tokyo, Japan, in 2006, andthe M.S. and Ph.D. degrees in engineering fromthe University of Tokyo, Tokyo, in 2008 and 2011,
respectively. Since 2016, he has been an Assistant Professor with Nagoya University, Japan. His
main research interests include mechanical safety,
human–robot interaction, and biomechanics.

YOJI YAMADA (Member, IEEE) received the
Ph.D. degree from the Tokyo Institute of Technology. Since 1983, he has been with the Toyota
Technological Institute, Nagoya, Japan. In 2004,
he joined the National Institute of Advanced
Industrial and Science Technology (AIST), as a
Leader of the Safety Intelligence Research Group.
In 2009, he moved to the Department of Mechanical Science and Engineering, Graduate School ofEngineering, Nagoya University, as a Professor.

213339