Nonlinear methods to quantify Movement Variability in Human-Humanoid Interaction Activities
SStrengths and Weaknesses of Recurrent QuantificationAnalysis in the context of Human-Humanoid Interaction
Miguel Xochicale , Chirs Baber , School of Engineering, University of Birmingham, Birmingham, UK* [email protected]
Abstract
Movement variability is defined as the variations that occur in motor performanceacross multiple repetitions of a task and such behaviour is an inherent feature withinand between each persons’ movement. Quantifying movement variability is still an openproblem, particularly when traditional methods in time domain and frequency domainfail to detect tiny modulations in frequency or phase for time series. We thereforeconsider methodologies from nonlinear dynamics such as reconstructed state space(RSS), uniform time-delay embedding (UTDE), recurrence plots (RPs) and recurrencequantification analysis (RQA) metrics. These algorithms require time series measuredwith sensors that provide well sampled data with little noise. However, for this study,we are interested in the analysis of data collected through wearable inertial sensors(IMUs) and its post-processing effects in each of the nonlinear dynamics methodologies.Twenty right-handed healthy participants (mean and standard deviation (SD) age ofmean=19.8 (SD=1.39)) performed an experiment in human-humanoid imitationactivities (H-H-I). Then IMUs were attached to the participants and the humanoidrobot performed simple vertical and horizontal arm movements in normal and fasterspeed. With four window lengths and three levels of smoothed time series, we foundvisual differences in the patterns of RSSs and RPs. Then using metrics of RQA, we findout that the type of movements and the level of smoothness affects those metrics. Inparticular, we found that entropy values from RQA were well distributed and presentingvariation in all the conditions for time series. Hence, we demonstrated the potential ofnonlinear techniques to quantify human movement variability in the context of H-H-Ican enhance the development of better diagnostic tools for various applicationsrehabilitation, sport science or for new forms of human-robot interaction.Movement variability is defined as the variations that occur in motor performanceacross multiple repetitions of a task and such behaviour is an inherent feature withinand between each persons’ movement. Quantifying such movement variability is still anopen problem, particularly when traditional methods in time domain and frequencydomain fail to detect tiny modulations in frequency or phase for time series. For thiswork, we therefore consider methodologies from nonlinear dynamics such asreconstructed state space (RSS), uniform time-delay embedding (UTDE), recurrenceplots (RPs) and recurrence quantification analysis (RQA) metrics (REC, DET, RATIO,and ENT). These algorithms generally require time series measured with sensors thatprovide well sampled data with little noise. However, for this study, we are interested inthe analysis of data collected through wearable inertial sensors and its post-processingeffects in each of the nonlinear dynamics methodologies. With that in mind, twentyright-handed healthy participants (mean and standard deviation (SD) age of mean=19.8(SD=1.39)) performed an experiment in the context of human-humanoid imitationOctober 23, 2018 1/40 a r X i v : . [ ee ss . SP ] O c t ctivities. Then inertial sensors were attached to the participants and the humanoidrobot that performed simple vertical and horizontal arm movements in normal and fasterspeed. With four window lengths and three levels of normalised smoothed of time series,we evidently found visual differences in the patterns of RSSs and RPs. Then usingmetrics of RQA, we find out that the type of movement (vertical or horizontal) and thelevel of smoothness affects those metrics. In particular, we found that ENT values werewell distributed and presenting variation in all the conditions for time series whichhappends partially in remained metrics. We demonstrated the potential of nonlineartechniques to quantify human movement variability in the context of human-humanoidimitation activities can enhance the development of better diagnostic tools for variousapplications rehabilitation, sport science or for new forms of human-robot interaction. Introduction
Human movement is a complex system where not only multiple joints and limbs areinvolved for a specific task in a determined environment but also how we process theexternal information with all of our available senses and use our experiences, which playa crucial role in the way each person moves. Recent studies in human motionrecognition have revealed the possibility to estimate features from lower dimensionssignals to distinguish differences between styles of pedalling motion [1, 2], to performgait identification [3, 4] or to do pattern recognition from biological signals [5]. The lowerdimension signals from biological signals are time series of 1 − dimension in R whichcommonly have high nonlinearity, complexity, and non-stationarity [5], where methodsof nonlinear dynamics can objectively quantify such human movement variability [1–7].With this in mind, [8] reviewed methods for nonlinear time series analysis based on theappropriate estimation of the embedding parameters ( m embedding dimension and τ embedding delay) to reconstruct the state space, where an n -dimensional reconstructedstate space using 1 − dimensional time series, can preserve the topological properties ofan unknown M -dimensional state space [9]. Also [8] reviewed the use of RecurrencePlots (RP), a graphical representation of a two-dimensional map which showblack/white dots as recurrences in a given n -dimensional system, and RecurrenceQuantification Analysis (RQA), metrics that can pick out important directions andstatistics in RPs. RPs and RQA help to have a more intuitive meaning of the timeseries, for instance, RQA is quantitatively and qualitatively independent of embeddingdimension which is also verified experimentally [10]. However, the estimation ofembedding parameters and finding the right parameters to perform RQA is still an openproblem. For instance, [8] pointed out that there is no general technique that can beused to compute the embedding parameters since the time series are system-dependentwhich means that once computing the values for embedding parameters, these may onlywork for one purpose (e.g., prediction) and may not work well for another purpose (e.g.,computing dynamical invariants). Additionally, the methodologies of nonlineardynamics for computing the embedding parameters e.g., autocorrelation, mutualinformation, and nearest neighbour require data which is well sampled and with littlenoise [11] or require purely deterministic signals [12]. Similarly, these methodologies forcomputing the embedding parameters can break down with real-world datasets whichhave generally different length, different values of accuracy and precision (roundingerrors due to finite precision of the measurement apparatus which include frequencyacquisition [4]), and data may be contaminated with different or unknown sources ofnoise [11]. It is surprising that even with the previous constraints with regard to thequality of data, and the problem with the estimation of embedding parameters, theresults of analysis using nonlinear dynamics have proven to be helpful to understandand characterise time series [1–5, 5–8]. Another point to consider with time seriesOctober 23, 2018 2/40nalysis using nonlinear dynamics is the appropriate use of different post-processingtechniques such as interpolation, filtering or normalisation.With that in mind, it can be said that there is little research and understanding onthe effects for post-processing techniques in the interpretation of reconstructed statespaces, RPs and RQA. For this work, we are interested in investigating and exploringthe effects of different features of time series (e.g. levels of smoothness, types ofmovements with different velocities and window length) to estimate appropriateembedded parameters for uniform time-delay embedding and for the computation forrecurrence plots and recurrence quantification analysis. To explore those questions, weconducted an experiment with twenty right-handed healthy participants in the contextof human-humanoid imitation activities where participants were asked to imitate simplearm movements performed by a humanoid robot. Therefore, the primary aim of thepaper is to explore the following questions: • What are the effects on the three nonlinear methods (uniform time-delayembedding, RPs, and RQA), for different embedding parameters, differentrecurrence thresholds and different characteristics of time series (window length,smoothness of the signal, structure)?
Materials and methods
State Space Reconstruction
The method of state space reconstruction was originally proposed by [13] and formalisedby [9]. Since then, different investigations and disciplines have benefited from the use ofthe method of state space reconstruction [2–4, 7, 14]. The method of state spacereconstruction is based on uniform time-delay embedding methodology which is a simplematrix implementation that can reconstruct an unknown d − dimensional manifold M from a scalar time series (e.g. one-dimensional time series in R ). A manifold, in thiscontext, is a multidimensional curved surface within a space (e.g. a saddle) [15].The use of a scalar time series is the main advantage of the method of state spacereconstruction which in essence preserve dynamic invariants such as correlationdimension, fractal dimension, Lyaponov exponents, Kolmogorov-Sinai entropy anddetrended fluctuation analysis [1, 2, 8, 16, 17]. However, selecting appropriate embeddingparameters which are sued to apply the state space reconstruction is still an openchallenge for which we present introductions for the methodologies to compute suchembedding parameters With that in mind, in the following subsections, we describe inmore detail the state space reconstruction theorem (RSSs), uniform time-delayembedding theorem (UTDE), false nearest neighbours (FNN), average mutualinformation (AMI) and other methodologies for state space reconstruction. State Space Reconstruction Theorem
Following the notation employed in [9, 11, 18–21], state space reconstruction is definedby: s ( t ) = f t [ s (0)] , (1)where s , s : A → M given that A ⊆ R and M ⊆ R d , represents a trajectory whichevolves in an unknown d − dimensional manifold M , f : M → M is an evolution functionand f t , with time evolution t ∈ N , is the t -th iteration of f that corresponds to aninitial position s (0) ∈ M [9]. Then, a point of a one-dimensional time series x ( t ) in R ,can be obtained with x ( t ) = h [ s ( t )] , (2)October 23, 2018 3/40here h is a function, h : M → R , defined on the trajectory s ( t ). Reconstructed statespace can then be described as an n − dimensional state space defined by y ( t ) = Ψ[ X ( t )]where X ( t ) = { x ( t ) , x ( t − τ ) , ..., x ( t − ( m − τ ) } is the uniform time-delay embeddingwith a dimension embedding m and delay embedding τ and Ψ : R m → R n is a furthertransformation of dimensionality (e.g. Principal Component Analysis, Singular ValueDecomposition, etc) being n ≤ m . With that in mind, uniform time-delay embedding, X ( t ), defines a map Φ : M → R m such that X ( t ) = Φ( s ( t )), where Φ is a diffeomorphicmap [9] whenever τ > m > d box and d box is the box-counting dimension of M [11]. Then, if Φ is an embedding of evolving trajectories in the reconstructed spacethen a composition of functions represented with F t is induced on the reconstructedstate space determined: X ( t ) = F t [ X (0)] = Φ ◦ f t ◦ Φ − [ X (0)] . (3)With this in mind, an embedding is defined as ”a smooth one-to-one coordinatetransformation with a smooth inverse” and the total reconstruction map is defined asΞ = Ψ ◦ Φ [18]. Fig 1 illustrates the state space reconstruction.
Uniform Time-Delay Embedding (UTDE)
Frank et al. and Sama et al. refer to the state space reconstruction as ”time-delayembeddings” or ”delay coordinates” [3, 4]. However, we consider the term ”uniformtime-delay embedding” as more descriptive and appropriate terminology for our work.The uniform time-delay embedding is represented as a matrix of uniform delayedcopies of the time series { x n } Nn =1 where N is the sample length of { x n } and n is indexfor the samples of { x n } . { x n } Nn =1 has a sample rate of T . The delayed copies of { x n } are uniformly separated by τ and represented as { ˜ x n − iτ } where i goes from0 , , . . . , ( m −
1) (Fig 2). Generally speaking, { ˜ x n − iτ } contains information ofunobserved state variables and encapsulates the information of the delayed copies of theavailable time series in the uniform time-delay embedding matrix X mτ , X mτ ∈ R m ,defined as X mτ = ˜ x n ˜ x n − τ ˜ x n − τ ... ˜ x n − ( m − τ (cid:124) , (4)where m is the embedding dimension, τ is the embedding delay and (cid:124) denotes thetranspose. m and τ are known as embedding parameters. The matrix dimension of X mτ is defined by N − ( m − τ rows and m columns and N − ( m − τ defines the length ofeach delayed copy of { ˜ x n } in X mτ . For further details and explicit examples of uniformtime-delay embedding methodology, we refer the reader to the S1 Appendix A. Estimation of Embedding Parameters
The estimation of the embedding parameters ( m and τ ) is a fundamental step for thestate space reconstruction with the use of uniform time-delay embedding method. Withthis in mind, we review two of the most common algorithms, which will be used in ourwork, to compute the embedding parameters: the false nearest neighbour (FNN) andthe average mutual information (AMI). False Nearest Neighbours
To select the minimum embedding dimension m , Kennel et al. [23] used the method offalse neighbours which can be understood as follows: on the one hand, when theOctober 23, 2018 4/40 ig 1. State space reconstruction methodology. State space reconstruction isbased on x ( t ) = h [ s ( t )] = h [ f t [ s (0)]] where f t is the true dynamical system, s ( t )indicates the state, s , at time, t , and h [] the measurement function. The time-delayembedding represented as the Φ, maps the original d − dimensional state s ( t ) into the m − dimensional uniform time-delay embedding matrix X ( t ). The transformation mapΨ then maps X ( t ) into a new state y ( t ) of dimensions n < m . (A) M − dimensionalmanifold representing the state space (e.g. Lorenz system); (B) Delayed copies of1 − dimensional x ( t ) from the Lorenz system; (C) m − dimensional reconstructed statespace with m and τ , and (D) y ( t ) is the n − dimensional reconstructed state space. Thetotal reconstruction map is represented as Ξ = Ψ ◦ Φ where Φ is the delayreconstruction map and Ψ is the coordinate transformation map. This figure is adaptedfrom [1, 18, 20] and R code to reproduce it is available from [22].embedding dimension is too small to unfold the attractor ”not all points that lie closeeach other will be neighbours and some points appear as neighbours because of theattractor has been projected down into an smaller space”, on the other hand, whenincreasing the embedding dimension ”points that are near to each other in the sufficientembedding dimension should remain close as the dimension increase from m to m + 1 [17]”. From a mathematical point of view, the state space reconstruction theoremis done when the attractor is unfolded with either the minimum embedding dimension,October 23, 2018 5/40 ig 2. Uniform time-delay embedding. UTDE is illustrated as m − { x n } , uniformly separated by τ and represented as { ˜ x n , . . . , ˜ x n − ( m − τ } (Eq. 4). R code to reproduce the figure is available from [22]. m , or any other embedding dimension value where m ≥ m [23]. On the contrary, anylarge value of m leads to excessive computations [8]. With this in mind, Cao [24]proposed an algorithm based on the false neighbour method where only the time-seriesand one delay embedding value are necessary to select the minimum embeddingdimension. Cao’s algorithm is based on E ( m ) which is the mean value of all a ( i, m ),October 23, 2018 6/40oth defined as follows E ( m ) = 1 N − mτ N − mτ (cid:88) i =1 a ( i, m )= 1 N − mτ N − mτ (cid:88) i =1 || X i ( m + 1) − X n ( i,m ) ( m + 1) |||| X i ( m ) − X n ( i,m ) ( m ) || (5)where X i ( m ) and X n ( i,m ) ( m ) are the time-delay embeddings with i = 1 , , . . . , N − ( m − τ and n ( i, m ) = 1 ≤ n ( i, m ) ≤ N − mτ . From Eq. 5, it can beseen that E ( m ) is only dependent on m and τ for which E ( m ) is defined as E ( m ) = E ( m + 1) E ( m ) . (6) E ( m ) is therefore considered to investigate the variation from m to m + 1 in order tofind the minimum embedding dimension m (Eq. 6). As described in [24]: ” E ( m ) stopschanging when m is greater than some m , if the time series comes from amultidimensional state space then m + 1 is the minimum dimension”. Additionally,Cao proposed E ( m ) to distinguish deterministic signals from stochastic signals. E ( m )is defined as E ( m ) = E ∗ ( m + 1) E ∗ ( m ) , (7)where E ∗ ( m ) = 1 N − mτ N − mτ (cid:88) i =1 | X i ( m + 1) − X n ( i,m ) ( m + 1) | . (8)For instance, when the signal comes from random noise (values that are independentfrom each other), all E ( m ) values are approximately equal to 1 (e.g. E ( m ) ≈ E ( m ) is not constant for all m (e.g. E ( m ) (cid:54) = 1).As an example of the use of E ( m ) and E ( m ) values, we consider two time series:the solution for the x variable of the Lorenz system (Fig 3E), and a Gaussian noise timeseries with zero mean and a variance of one (Fig 3F). We then compute E ( m ) and E ( m ) values for each time series. The E ( m ) values for the chaotic time series appearto be constant after the dimension is equal to six. The determination of six is given thatany value of m can be used as they are within the threshold of 1 ± .
05 (Fig 3A). E ( m )values, for chaotic time series, are different to one (Fig 3C), for which, it can beconcluded that for the chaotic time series the minimum embedding dimension the timeseries comes from a deterministic signal. With regard to the noise time series, E ( m )values appeared to be constant when m is close to thirteen, which is defined by thethreshold of 1 ± .
05 (Fig 3B). E ( m ) values then indicate the minimum embeddingdimension of the noisy time series is thirteen, however all of the E ( m ) values areapproximately equal to one (Figure 3D) for which it can be concluded that noise timeseries is a stochastic signal.It is important to note that for this work not only E ( m ) and E ( m ) are computedbut also a variation of τ from 1 to 20 is presented. The purpose of such variation for τ is to show its independence with regard to E ( m ) and E ( m ) values as τ is increasing(Fig 3A,B,C, and D) However, one negative of the Cao’s algorithm [24] is the definitionof a new threshold where m values appear to be constant in E ( m ). In the case of thegiven examples and reported results, we defined such threshold as 0.05. Furtherinvestigation is required for the selection of the threshold in the E ( m ), as the selectionof the threshold in this work is base on no particular method but visual inspection.October 23, 2018 7/40 B CD EF
Fig 3. Minimum dimension embedding values with Cao’s method. (A, B) E ( m ) values and (C, D) E ( m ) values with variations of τ values from one to twentyfor (E) chaotic and (F) random time series. R code to reproduce the figure is availablefrom [22]. Average Mutual Information
When selecting the delay dimension parameter, τ , one can consider the following twocases: (i) when τ is too small, the elements of time-delay embedding will be along thebisectrix of the phase space and the reconstruction is generally not satisfactory, (ii) onthe contraty, when τ is too large the elements of the uniform time-delay embedding willbecome spread and uncorrelated which makes recovering the underlying attractor moredifficult if not impossible [18, 25, 26]. With regard to the algorithms to compute τ ,Emrani et al. [25], for instance, used the autocorrelation function in which the first zerocrossing is considered as the minimum delay embedding parameter. However, theautocorrelation function is a linear statistic over which the Average Mutual Information(AMI) algorithm is preferred because the AMI takes into account the nonlineardynamical correlations [17, 27]. With this in mind, the AMI algorithm is describedbelow to estimate the minimum delay embedding parameter, τ .To compute the AMI, an histogram of x ( n ) using n bins is calculated and then aprobability distribution of data is computed [12]. AMI is therefore denoted by I ( τ )which is the average mutual information between the original time series, x ( n ), and thedelayed time series, x ( n − τ ), delayed by τ [28]. AMI is defined by I ( τ ) = N (cid:88) i,j p ij log p ij p i p j . (9)Probabilities are defined as follows: p i is the probability that x ( n ) has a value inside the i -th bin of the histogram, p j is the probability that x ( n + τ ) has a value inside the j -thbin of the histogram and p ij ( τ ) the probability that x ( n ) is in bin i and x ( n + τ ) is inbin j . The AMI is measured in bits (base 2, also called shannons) [12, 29]. For small τ ,AMI will be large and it will then decrease more or less rapidly. As τ increase and goesto a large limit, x ( n ) and x ( n + τ ) have nothing to do with each other and p ( ij ) isfactorised as p i p j for which AMI is close to zero. Then, in order to obtain τ , ”it has tobe found the first minimum of I ( τ ) where x ( n + τ ) adds maximal information to theknowledge from x ( n ), or, where the redundancy is the least” [12].For example, we compute the AMI for two time series: A) the x solution of theOctober 23, 2018 8/40orenz system, and B) a noise time series using a normal distribution with mean zeroand standard deviation equal to one. From Fig 4, it can then be concluded that theamount of knowledge for any noise time series is zero for which the first minimumembedding parameter is τ = 1. On the contrary, the first minimum of the AMI for thechaotic time series is τ = 17 which is the value that maximize the independencebetween x ( n ) and x ( n + τ ) in the reconstructed state space [8]. Similarly as Cao’s AB CD
Fig 4. Minimum delay embedding values with AMI’s method. (A, B) AMIvalues where its first minimum value in the curve is the minimum time delay embedding( τ ), for (C) a chaotic and (D) noise time series. R code to reproduce the figure isavailable from [22].algorithm negatives, AMI’s algorithm is not an exception for negatives, which areworthwhile to mention for further investigations. For instance, (i) is not clear why thechoose of the first minimum of the AMI is the minimum delay embeddingparameter [12] and (ii) the probability distribution of the AMI function is computedwith the use of histograms which depends on a heuristic choice of number of bins forwhich AMI depends on partitioning [26]. Other methodologies for state space reconstruction.
It is important to note that other methods for state space reconstruction have beenrecently investigated, for instance: (i) the nonuniform time-delay embeddingmethodology where the consecutive delayed copies of { x n } are not equidistant. Suchmethod has been proben to create better representations of the dynamics of the statespace to analyse, for instance, quasiperiodic and multiple time-scale time series over theconventional uniform time-delay embedding algorithm [1, 2, 16, 20, 30], and (ii) Uniform 2time-delay embedding method which takes advantage of finding an embedding windowinstead of the traditional method of finding the embedding parameters separately [5]. Ingeneral, Uniform 2 time-delay embedding method computes m with False NearestOctober 23, 2018 9/40eighbour (FNN) algorithm and τ is computed as τ = d w / ( m − d w is givenby the minimisation of the Minimum Description Length [31]. However, the previousmethods are out of the scope of the paper but we think is important to refer readers tothose references for further investigations in this regard. Recurrence Quantification
Recurrence Plots
Originally, Henri Poincar´e in 1890 introduced the concept of recurrences in conservativesystems, however such discovery was not put into practice until the development offaster computers [32], for which Eckmann et al. [33] in 1987 introduced a method whererecurrences in the dynamics of a system can be visualised using Recurrence Plots. Theintention of Eckmann et al. [33] was to propose a tool, called Recurrence Plot (RP),that provides insights into high-dimensional dynamical systems where trajectories arevery difficult to visualise. Therefore, ”RP helps us to investigate the m − dimensionalphase space trajectories through a two-dimensional representation of itsrecurrences” [34]. Similarly, Marwan et al. [34] pointed out that additionally to themethodologies of the state space reconstruction and other dynamic invariants such asLyapunov exponent, Kolmogorov-Sinai entropy, the recurrences of the trajectories in thephase space can provide important clues to characterise natural process that present, forinstance, periodicities (as Milankovitch cycles) or irregular cycles (as El Ni˜no SouthernOscillation). Such recurrences can not only be presented visually using Recurrence Plots(RP) but also be quantified with Recurrence Quantification metrics, which leads toapplications of these in various areas such as economy, physiology, neuroscience, earthscience, astrophysics and engineering [32].For the creation of a recurrence plot based on time series { x n } , it is first computedthe state space reconstruction with uniform time-delay embedding X ( i ) = { ˜ x n , . . . , ˜ x n − ( m − τ } where i = 1 , . . . , N , N is the number of considered statesof X ( i ) and X ( i ) ∈ R m [33]. The recurrence plot is therefore a two-dimensional N × N square matrix, R , where a black dot is placed at ( i, j ) whenever X ( i ) is sufficientlyclose to X ( j ): R mi,j ( (cid:15) ) = Θ( (cid:15) i − || X ( i ) − X ( j ) || (10)where i, j = 1 , . . . , N , (cid:15) is a threshold distance, ||· || a norm, and Θ( · ) is the Heavisidefunction (i.e. Θ( x ) = 0, if x <
0, and Θ( x ) = 1 otherwise) (Fig 5) [32–34]. RP is alsocharacterised with a line of identity (LOI) which is a diagonal line due to R i,j = 1( i, j = 1 , . . . , N ). Structures of Recurrence Plots
Pattern formations in the RPs can be designated either as topology for large-scalepatterns or texture for small-scale patterns. In the case of topology, the followingpattern formations are commonly presented: (i) homogeneous where uniform recurrencepoints are spread in the RP e.g., uniformly distributed noise (Fig 6A), (ii) periodic andquasi-periodic systems where diagonal lines and checkerboard structures representoscillating systems, e.g., sinusoidal signals (Fig 6B), (iii) drift where paling or darkeningrecurrence points away from the LOI is caused by drifting systems, e.g., logistic map(Fig 6C), and (iv) disrupted where recurrence points are presented white areas or bandsthat indicate abrupt changes in the dynamics, e.g. Brownian motion (Fig 6D) [33, 34].Texture patterns in RPs can be categorised as: (i) single or isolated recurrence pointsthat represent rare occurring states, and do not persist for any time or fluctuate heavily,(ii) dots forming diagonal lines where the length of the small-scale parallel lines in theOctober 23, 2018 10/40 B Fig 5. Recurrence Plots. (A) State space of the Lorenz system with controllingparameters ( ρ = 28 , σ = 10 , β = 8 / j , in trajectory X () which falls into theneighborhood (black circle) of a given point at i is a recurrent point and is representedas a black dot in the recurrence plot at location ( i, j ) or white otherwise. (B)Recurrence plot using the three components of the Lorenz system and the RP with noembeddings and threshold (cid:15) = 5. This figure is adapted from [34] and R code toreproduce it is available from [22].diagonal are related to the ratio of determinism or predictability in the dynamics of thesystem, and (iii) dots forming vertical and horizontal lines where the length of the linesrepresent a time length where a state does not change or change very slowly and thesepatterns formation represent discontinuities in the signal, and (iv) dots clustering toinscribe rectangular regions which are related to laminar states or singularities [34]. A B C D
Fig 6. Patterns in Recurrence Plots.
Time-series with its respective recurrenceplot for: (A) uniformly distributed noise, (B) super-positionet harmonic oscillation( sin ( ∗ t ) ∗ sin ( ∗ t )), (C) drift logistic map ( x i +1 = 4 x i (1 − x i )) corrupted with alinearly increase term (0 . ∗ i ), and (D) disrupted brownian motion( x i +1 = x i + 2 ∗ rnorm (1)). This figure is adapted from [34] and R code to reproduce itis available from [22].Although, each of the previous pattern descriptions of the structures in the RP offeran idea of the characteristics of dynamical systems, these might be misinterpreted andconclusions might tend to be subjective as these require the interpretation of aparticular researcher(s). Because of that, recurrence quantification analyis (RQA) offerOctober 23, 2018 11/40bjective methodologies to quantify such visual characteristics of previous recurrentpattern structures in the RP [35]. Recurrence Quantifications Analysis (RQA)
Originally, Zbilut et al. [35] proposed metrics to investigate the density of recurrencepoints in RPs, then histograms of lengths for diagonal lines in RPs were studied by [36]which were the introduction to the term recurrence quantification analysis (RQA) [37].RQA has been applied in many fields such as life science, engineering, physics, andothers [37] Particularly in human movement to investigate noise and complexity ofpostural control [40], postural control [38] or interpersonal coordination [39]. Thesuccess of RQA is not only due to its simple algorithmic implementation but also to itscapacity to detect tiny modulations in frequency or phase which are not detectableusing standard methods e.g. spectral or wavelet analysis [6], and that RQA’s metricsare quantitatively and qualitatively independe]nt of embedding dimension which isverified experimentally by [10]. RQA metrics comprehend percentage of recurrence,percentage of determinism, ratio, Shannon entropy of the frequency distributions of theline lengths, maximal line length and divergence, trend and laminariy [32, 34]. For thiswork, we considered only four RQA metrics, due to its consistency with our preliminaryexperiments, which are described below. Such metrics are computed thenonlinearTseries R package [29].
REC values
The percentage of recurrence (REC) is defined as
REC ( (cid:15), N ) = 1 N − N N (cid:88) i (cid:54) = j =1 R mi,j ( (cid:15) ) , (11)which enumerates the black dots in the RP excluding the line of identity. REC is ameasure of the relative density of recurrence points in the sparse matrix [34]. DET values
The percent determinism (DET) is defined as the fraction of recurrence points thatform diagonal lines and it is determined by
DET = (cid:80) Nl = d min lH D l (cid:80) Ni,j =1 R i,j ( (cid:15) ) , (12)where H D ( l ) = N (cid:88) i,j =1 (1 − R i − ,j − ( (cid:15) ))(1 − R i + l,j + l ( (cid:15) )) l − (cid:89) k =0 R i + k,j + k ( (cid:15) ) (13)is the histogram of the lengths of the diagonal structures in the RP. DET can beinterpreted as the predictability of the system for periodic signals which, in essence,have longer diagonal lines than the short diagonals lines for chaotic signals or absentdiagonal lines for stochastic signals [32, 34]. Similarly, DET is considered as ameasurement for the organisation of points in RPs [10]. RATIO values
RATIO is defined as the ratio between DET and REC and it is calculated from thefrequency distributions of the lengths of the diagonal lines. RATIO is useful to discoverdynamic transitions [34].October 23, 2018 12/40
NT values
ENT is the Shannon entropy of the frequency distribution of the diagonal line lengthsand it is defined as
EN T = − N (cid:88) l = d min p ( l ) lnp ( l ) with p ( l ) = H D ( l ) (cid:80) Nl = d min H D ( l ) . (14)ENT reflects the complexity of the deterministic structure in the system. For instance,for uncorrelated noise or oscillations, the value of ENT is rather small and indicates lowcomplexity of the system, therefore ”the higher the ENT is the more complex thedynamics are” [32, 34]. Sensitivity and robustness of RPs and RQA.
RP and RQA are a very young field in nonlinear dynamics and many questions are stillopen, for instance, different parameters for window length size of the time series,embedding parameters or recurrence threshold can generate different results in RQA’smetrics [6, 33].The selection of recurrence threshold, (cid:15) , can depend on the system that is analysed.For instance, when studying dynamical invariants (cid:15) require to be very small, fortrajectory reconstruction (cid:15) requires to have a large thresholds or when studyingdynamical transition there is little importance about the selection of the threshold [6].Other criteria for the selection of (cid:15) is that the recurrence threshold should be five timeslarger than the standard deviation of the observational noise or the use of diagonalstructures within the RP is suggested in order to find the optimal recurrence thresholdfor (quasi-)periodic process [6].Similarly, Iwanski et al. [10] highlighted the importance of choosing the rightembedding parameters to perform RQA for which many experiments have to beperformed using different parameters in order to have a better intuition of the nature ofthe time series and how this is represented by using RQA.With that in mind, this work explores the sensitivity and robustness of the windowsize of time series, embedding parameters for RSS with UTDE and recurrence thresholdfor RP and RQA in order to gain a better insight into the underlying time seriescollected from inertial sensors in the context of human-humanoid imitation activities.
Experiment
We conducted an experiment in the context of human-humanoid imitation (HHI)activities where participants were asked to imitate simple horizontal and vertical armmovements performed by NAO, a humanoid robot [41]. Such simple movements wererepeated ten times for the participant who copied NAO’s arm movements in aface-to-face imitation activity. Also, wearable inertial measurement unit (IMU) sensorswere attached to the right hand of the participant and to the left hand of the robot(Figure 7 A,C). Data were then collected with four NeMEMSi IMU sensors withsampling rate of 50Hz provinding tri-axial data of the accelerometer, gyroscope andmagnetometer sensors and quaternions [42]. A further description of the NeMEMSiIMU sensors is given in S1 Appendix B.
Participants
Twenty-three participants, from now on defined as pN where N is the number ofparticipant, were invited to do the experiment. However, data for three participantsOctober 23, 2018 13/40 eMEMSi C NeMEMSi
Fig 7. Human-humanoid imitation activities.
Face-to-face human-humanoidimitation (HHI) activities for (A) HHI of horizontal arm movement, (B) Humanoidhorizontal arm movement, (C) HHI of vertical arm movement, and (D) Humanoidvertical arm movement.were not used because the instructions for p
01, who was the only left-handed, weremistakenly given in a way that movements were performed different from what hadbeen planned, and for participants p
13 and p
16 data were corrupted because bluetoothcommunications problems with the sensors. With that in mind, data for twentyparticipants were analysed in this work.Of the twenty participants, all of them are right-handed healthy participants ofwhom four are females and sixteen are males with a mean and standard deviation (SD)age of mean=19.8 (SD=1.39). All participants provided informed consent forms prior toparticipation in the experiment.
Human-humanoid imitation activities
For human-humanoid imitation (HHI) activities four neMEMSi sensors were used, twoof which were attached to the right hand of the participant and the other two to the lefthand of the humanoid robot. Then, each participant was asked to imitate repetitions ofsimple horizontal and vertical arm movements performed by the humanoid robot in thefollowing conditions: (i) ten repetitions of horizontal arm movement at normal (HN)and faster (HF) speed (Figure 7 A), and (ii) ten repetitions of vertical arm movement atnormal (VN) and faster (VF) speed (Figure 7 C). The normal and faster speed of armmovements is defined by the duration in number of samples of one repetition of NAO’sarm movements. We select NAO’s arm movements duration to distinguish betweennormal and faster arm movements as NAO’s movements have less variation betweenrepetition to repetition. The duration for one repetition of the horizontal armmovement at normal speed, HN, is about 5 seconds considering that each repetition lastaround 250 samples. For horizontal arm movement at faster speed, HF, each repetitionwere performed in around 2 seconds which correspond to 90 samples of data. Thevertical arm movement at normal speed, VN, were performed in 6 seconds which isaround 300 samples of data. For vertical arm movement at faster speed, VF, eachrepetition lasts about 2.4 seconds which correspond to 120 samples of data. To visualisethe distinction between normal and faster speed for horizontal and vertical armmovements, Fig 8 shows smoothed time series for axes Z and Y of the gyroscope sensorswith four window lengths: 2-sec (100-samples), 5-sec (250-samples), 10-sec (500-samples)and 15-sec (750-samples).October 23, 2018 14/40
BC D
Fig 8. Time series duration of horizontal and vertical arm movements.
Time series of smoothed data from gyroscope sensor for different speed arm movementsperformed by NAO: (A) Horizontal Normal arm movement (HN), (B) Horizontal Fasterarm movement (HF), (C) Vertical Normal arm movement (VN) and (D) Vertical Fasterarm movement (VF). Additionally, (A) shows window sizes for 2-seconds (100 samples),5-seconds (250 samples), 10-seconds (500 samples) and 15-seconds (750 samples) whichare also presented in (B), (C) and (D). R code to reproduce figure is available [22].
Data from Inertial Measurement Units
Raw data
Considering the work of [43] which provided evidence of an improvement in recognitionactivities when combining data from accelerometer and gyroscope. We focus ouranalysis from data of the accelerometer and gyroscope of the NeMEMsi sensors [42] andleave the data of the magnetometer and quaternions for further investigation because oftheir possible variations with regard to magnetic disturbances.Data from the accelerometer is defined by triaxial time series A x ( n ), A y ( n ), A z ( n )which forms the matrix A (Eq. 15), and the same for data from the gyroscope which isdefined by triaxial time-series of G x ( n ), G y ( n ), G z ( n ) representing the matrix G (Eq. 16). Both triaxial time series of each sensor, a and g , are denoted with its respectiveaxes subscripts x, y, z , where n is the sample index and N is the same maximum lengthof all axes for the time series. Matrices A and G are represented as follows A = A x ( n ) A y ( n ) A z ( n ) = a x (1) , a x (2) , . . . , a x ( N ) a y (1) , a y (2) , . . . , a y ( N ) a z (1) , a z (2) , . . . , a z ( N ) , (15)October 23, 2018 15/40nd G = G x ( n ) G y ( n ) G z ( n ) = g x (1) , g x (2) , . . . , g x ( N ) g y (1) , g y (2) , . . . , g y ( N ) g z (1) , g z (2) , . . . , g z ( N ) . (16) Postprocessing data
After the collection of raw data from four NeMEMsi sensors, time synchronisationalignment and interpolation were performed in order to create time series with the samelength and synchronised time. We refer the reader to [42] for further details about thetime synchronisation process.
Data normalization
Data is normalised to have zero mean and unit variance using sample mean and samplestandard deviation. The sample mean and sample standard deviation using x ( n ) isgiven by µ x ( n ) = 1 N ( N (cid:88) i =1 x ( i )) , σ x ( n ) = (cid:115) (cid:80) N ( x ( i ) − µ x ( n ) ) N − , (17)and the normalised data, ˆ x ( n ), is computed as followsˆ x ( n ) = x ( n ) − µ x ( n ) σ x ( n ) . (18) Smoothing data
Commonly, a low-pass filter is the method either to capture the low frequencies thatrepresent %99 of the human body energy or to get the gravitational and body motioncomponents of accelerations [44]. However, for this work the elimination of certainrange of frequencies is not the main focus but the conservation of the structure in thetime series in terms of the width and heights where, for instance, Savitzky-Golay filtercan help to accomplish such task [45]. Savitzky-Golay filter is based on the principle ofmoving window averaging which preserves the area under the curve (the zeroth moment)and its mean position in time (the first moment) but the line width (the secondmoment) is violated and that results, for example, in the case of spectrometric datawhere a narrow spectral line is presented with reduced height and width. With that inmind, the aim of Savitzky-Golay filtering is to find filter coefficients c n that preservehigher momentums which are based on local least-square polynomialapproximations [45–47]. Therefore, Savitzky-Golay coefficients are therefore computedusing an R function sgolay(p,n,m) where p is the filter order, n is the filter length(must be odd) and m is the m -th derivative of the filter coefficients [48]. Smoothedsignal is represented with a tilde over the original signal: ˜ x ( n ). Window size data
With regard to the window size, [43] investigated its effects using seven window lengths(2, 5, 10, 15, 20, 25, 30 seconds) and combination of inertial sensors (accelerometer,gyroscope and linear acceleration sensor) to improve the activity recognitionperformance for repetitive activities (walking, jogging and biking) and less repetitiveactivities (smoking, eating, giving a talk or drinking a coffee). With that in mind,Shoaib et al. [43] concluded that the increase of window length improve the recognitionof complex activities because these requires a large window length to learn the repetitivemotion patterns. Particularly, one of the recommendations is to use large window sizeOctober 23, 2018 16/40o recognise less repetitive activities which mainly involve random hand gestures.Therefore, for the four activities (HN, HF, VN, and VF) in this work, which are mainlyrepetitive, we select only four window sizes for analysis: 2-s window (100 samples), 5-swindow (250 samples), 10-s (500 samples) and 15-s window (750 samples) (Fig 8).
Results
Once time series were collected, we investigated the robustness and weaknesses of thereconstructed state spaces (RRSs) using the uniform time-delay embedding technique(UTDE) and recurrence plots (RPs) for recurrent quantification analysis (RQA)methodologies in the following conditions: • Three levels of smoothness for the normalised data (sg0zmuv, sg1zmuv andsg2zmuv), computed from two different filter lengths (29 and 159) with the samepolynomial degree of 5 using the function sgolay(p,n,m) [48], • Four velocities arm movement activities: horizontal normal (HN), horizontal faster(HF), vertical normal (VN), and vertical faster (VF), and • Four window lengths: { } . Time series
After the data collection, raw time series were windowed, normalised and smoothed. Weonly present 10-sec (500 samples) window length time series, due to space limitations,for three participants (p01, p01 and p03) performing horizontal arm movements (axisGyroZ) and vertical arm movements (axis GyroY) (Figs 9 and 10).
A B C
Fig 9. Time series for horizontal arm movements. (A) raw-normalised(sg0zmuvGyroZ), (B) normalised-smoothed 1 (sg1zmuvGyroZ) and (C)normalised-smoothed 2 (sg2zmuvGyroZ). Time series are only for three participants(p01, p02, and p03) for horizontal movements in normal and faster velocity (HN, HF)with the normalised GyroZ axis (zmuvGyroZ) and with one sensor attached to theparticipant (HS01) and other sensor attached to the robot (RS01). R code to reproducethe figure is available from [22].October 23, 2018 17/40
B C
Fig 10. Time series for vertical arm movements. (A) raw-normalised(sg0zmuvGyroY), (B) normalised-smoothed 1 (sg1zmuvGyroY) and (C)normalised-smoothed 2 (sg2zmuvGyroY). Time series are only for three participants(p01, p02, and p03) for vertical movements in normal and faster velocity (VN, VF) withthe normalised GyroY axis (zmuvGyroY) and with one sensor attached to theparticipant (HS01) and other sensor attached to the robot (RS01). R code to reproducethe figure is available from [22].
Minimum embedding parameters
Minimum embedding parameters were firstly computed to create the reconstructedstate spaces with UTDE and RQAs with RPs.
False Nearest Neighbour
Considering the time series for twenty participants, minimum embedding dimensionswere computed using False Nearest Neighbour for horizontal and vertical armmovements. Figs 11 and 12 show that minimum embedding values appear to be moreconstant for sensor RS01 than the slightly variations of such values for sensor HS01. Itcan also be seen that there is a minor decrease of minimum embedding values assmoothness of time series increase. To have an overall minimum dimension value thatrepresent participants, sensors and activities, a sample mean were computed over all theminimum values in Figs 11 and 12 which results in m = 6. Average Mutual Information
Similarly, considering the time series for twenty participants, minimum delay valueswere computed as the first minimum values of the Average Mutual Information (AMI)for horizontal and vertical arm movements (Figs 13 and 14).For horizontal arm movements, Fig 13A shows that values tend to be more spread asthe smoothness is increased which is different for Fig 13C where values show no effect asthe smoothness of time series increase. In contrast, Fig 13B shows the values are lessspread as smoothness is increased which we believe the reason for that is due to thehigh frequencies on robots movements in the horizontal normal movement. However,values in Fig 13D tend to be spread as smoothness is increasing which are due to verydifferent curves in the AMI. With regard to vertical arm movements, values in Figs 14Aand 14C show an slightly increase of the spread values as the smoothness increase andvalues in Fig 14B appear to have less variation as the smoothness of the signals isincreasing. However, that do not happen for the second smoothed valuesOctober 23, 2018 18/40
CB D
HN HF H S R S Fig 11. Minimum embedding dimensions for horizontal arm movements. (A, B) Horizontal Normal (HN), (C, D) Horizontal Faster (HF) movements, (A, C)sensor attached to participants (HS01), and (B, D) sensor attached to robot (RS01).Minimum embedding dimensions are for twenty participants (p01 to p20) with threesmoothed signals (sg0zmuvGyroZ, sg1zmuvGyroZ and sg2zmuvGyroZ) and windowlenght of 10-sec (500 samples). R code to reproduce the figure is available from [22].
A CB D
VN VF H S R S Fig 12. Minimum embedding dimensions for vertical arm movements. (A,B) Vertical Normal (VN), (C, D) Vertical Faster (VF) movements, (A, C) sensorattached to participants (HS01), and (B, D) sensor attached to robot (RS01). Minimumembedding dimensions are for twenty participants (p01 to p20) with three smoothedsignals (sg0zmuvGyroY, sg1zmuvGyroY and sg2zmuvGyroY) and window length of10-sec (500 samples). R code to reproduce the figure is available from [22].(sg2zmuvGyroY) in Fig 14D. We also computed an overall minimum delay value thatrepresent participants, sensors and activities, using a sample mean of all values inFigs 13 and 14 which results in τ = 8. A CB D
HN HF H S R S Fig 13. First minimum AMI values for horizontal arm movements. (A, B)Horizontal Normal (HN), (C, D) Horizontal Faster (HF) movements, (A, C) sensorattached to participants (HS01), and (B, D) sensor attached to robot (RS01). Firstminimum AMI values are for twenty participants (p01 to p20) with three smoothedsignals (sg0zmuvGyroZ, sg1zmuvGyroZ and sg2zmuvGyroZ) and window lenght of10-sec (500 samples). R code to reproduce the figure is available from [22].October 23, 2018 19/40
CB D
VN VF H S R S Fig 14. First minimum AMI values for vertical arm movements. (A, B)Vertical Normal (VN), (C, D) Vertical Faster (VF) movements, (A, C) sensor attachedto participants (HS01), and (B, D) sensor attached to robot (RS01). First minimumAMI values are for twenty participants (p01 to p20) with three smoothed signals(sg0zmuvGyroZ, sg1zmuvGyroZ and sg2zmuvGyroZ) and window lenght of 10-sec (500samples). R code to reproduce the figure is available from [22].
Reconstructed State Spaces using Uniform Time-DelayEmbedding
Although the implementation of Uniform Time-Delay Embedding for the reconstructedstate space, one of the main challenges of the latter is the selection of embeddingparameters because each time series is unique in terms of its structure (modulation ofamplitude, frequency and phase) [3, 4, 8]. With that in mind, the problem is not tocompute individual embedding parameters for each of the time series but to deal with aselecting of two parameters that can represent all the time series. Our solution for thatwas to compute a sample mean over all values in each of the conditions of the timeseries of Figs 11, 12 for minimum dimension values and Figs 13 and 14 for minimumdelay values, resulting in an average minimum embedding parameters of ( m = 6, τ = 8). Then, the reconstructed state spaces were computed with ( m = 6, τ = 8) andthe first three axis of the rotated data of the PCA are shown in Figs 15 for horizontalarm movements and Figs 16 for vertical arm movements.Evidently, it is easy to observe by eye the differences in each of the trajectories inthe reconstructed state spaces (Figs 15, 16), however one might be not objective whenquantifying those differences since those observation might vary from person to person.With that in mind, we tried to objectively quantify those differences using euclideandistances between the origin to each of the points in the trajectories, however thesecreated suspicious metric, specially for trajectories which looked very messy. With thatin mind, we computed Recurrence Quantification Analysis to objectively quantify thedifferences in each of the cases of the time series. Recurrences Plots
Considering the time series of Figs 9 and 10, we computed its Recurrence Plots forhorizontal arm movements (Fig 17) and for vertical arm movements (Fig 18) using theaverage embedding parameters ( m = 6, τ = 8) and an recurrence threshold of (cid:15) = 1.With regard to the selection of recurrence threshold, Marwan et al. [6] pointed out thatchoosing an appropriate recurrence threshold is crucial to get meaningfulrepresentations in RPs, however, for our work where quantifying movement variability isour aim, we give little importance to the selection of the recurrence threshold as as longas it is able to represent the dynamical transitions in each of the time series.October 23, 2018 20/40 N H S R S p H S R S p H S R S p H S R S p H S R S p H S R S p sg0zmuvGyroZ sg1zmuvGyroZ sg2zmuvGyroZ sg0zmuvGyroZ sg1zmuvGyroZ sg2zmuvGyroZHF Fig 15. RSSs for horizontal arm movements.
Reconstructed state spaces fortime series of Figure 9. Reconstructed state spaces were computed with embeddingparameters m = 6, τ = 8. R code to reproduce the figure is available from [22].As similar as with the Reconstructed State Spaces, the differences in the RPs can beeasily noticed by eye for different conditions of the time series (Figs 18, Fig 17), whichlead us to apply Recurrence Quantification Analysis to have an objective quantificationof each of the time series. Recurrence Quantification Analysis
RPs provide pattern formations for each of the time series conditions (Figs 17 and 18)which are quantified with four metrics of RQA: REC, DET, RATIO and ENTR.
REC values
In Figs 19 and 20 can be seen that REC values are more spread for HN than HFmovements with data coming from HS01 sensor. In contrast, REC values appear to beconstant and present little variation for both HN and HF movements with data fromthe sensor attached to the humanoid robot RS01. With regard to the increase ofsmoothness of data (sg0zmuvGyroZ, sg1zmuvGyroZ and sg2zmuvGyroZ), REC valuespresent little variation as the smoothness is increasing for data from HS01 and RECvalues more similar as the smoothness is increasing for data from RS01.October 23, 2018 21/40 N H S R S p H S R S p H S R S p H S R S p H S R S p H S R S p sg0zmuvGyroY sg1zmuvGyroY sg2zmuvGyroY sg0zmuvGyroY sg1zmuvGyroY sg2zmuvGyroYVF Fig 16. RSSs for vertical arm movements.
Reconstructed state spaces for timeseries of Figure 10. Reconstructed state spaces were computed with embeddingparameters m = 6, τ = 8. R code to reproduce the figure is available from [22]. DET values
Little can be said with regard to the variation of DET values as these change very littleeven for type of movement or type of sensor (Figs 21 and 22). With regard to thesmoothness of time series, DET values appear to be more similar as the smoothness ofthe data is increasing.
RATIO values
RATIO values for HN movements vary less than HF movements for HS01 sensor whichis similar behaviour of RATIO values for RS01 sensors in both vertical and horizontalmovements (Figs 23 and 24). It can also noticed a decrease of variation in RATIOvalues as the smoothness of the signal is increasing.
ENTR values
ENTR values show more variation for HS01 sensor than ENTR values for RS01 sensorwhich appear to be more constant and the smoothness of data affects little to thevariation of ENTR values (Figs 25 and 26).October 23, 2018 22/40 S R S p H S R S p H S R S p H S R S p H S R S p H S R S p HNsg0zmuvGyroZ sg1zmuvGyroZ sg2zmuvGyroZ sg0zmuvGyroZ sg1zmuvGyroZ sg2zmuvGyroZHF
Fig 17. RPs for horizontal arm movements.
Recurrence plots were computedwith embedding parameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure isavailable from [22]. RQA metrics with different embedding parameters, recurrencethresholds, window lengths, levels of smoothness, and timeseries structures.
Zbilut et al. [35] established RQA metrics with the aim of determining embeddingparameters, their method consisted on creating 3D surfaces with RQA metrics with anincrease of embedding parameters ( m and τ ), then Zbilut et al. [35] exploredfluctuations and gradual changes in the 3D surfaces that provide information about theembeddings. Much recently, Marwan et al. [34] created 3D surfaces for visual selectionof not only embedding parameters but also recurrence thresholds. Following samemethodologies, we explored the stability and robustness of RQA metrics (REC, DET,RATIO and ENTR) using 3D surfaces by an unitary increase of the pair embeddingparameters (0 ≥ m ≤
10, 0 ≥ τ ≤
10) and a decimal increase of 0.1 for recurrencethresholds (0 . ≥ (cid:15) ≤
3) (Fig. 27). We also computed 3D surfaces of RQA metrics fordifferent sensors and different activities (Fig. 29). RQA metrics are also affected by thewindow length where for example four window lengths of 100, 250, 500 and 750 samples(Fig. 29). Three level of smoothness were computed for RQA metrics showing smoothed3D surfaces ad the level of smoothness increase (Fig. 30). Similarly, 3D surfaces of RQAmetrics were also computed for three participants (Fig. 31).October 23, 2018 23/40 S R S p H S R S p H S R S p H S R S p H S R S p H S R S p VNsg0zmuvGyroZ sg1zmuvGyroZ sg2zmuvGyroZ sg0zmuvGyroZ sg1zmuvGyroZ sg2zmuvGyroZVF
Fig 18. RPs for vertical arm movements.
Recurrence plots were computed withembedding parameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure isavailable from [22]. HS01 RS01
Fig 19. REC values for horizontal arm movements.
REC values (representing% of black dots in the RPs) for 20 participants performing HN and HF movements withsensors HS01, RS01 and three smoothed-normalised axis of GyroZ (sg0zmuvGyroZ,sg1zmuvGyroZ and sg2zmuvGyroZ). REC values were computed with embeddingparameters m = 6, τ = 8 and (cid:15) = 1 R code to reproduce the figure is available from [22]. Discussion
It is evidently that time series from different sources (participants, movements, axistype, window length or levels of smoothness) presents visual differences for embeddingOctober 23, 2018 24/40
S01 RS01
Fig 20. REC values for vertical arm movements.
REC values (representing % ofblack dots in the RPs) for 20 participants performing VN and VF movements withsensors HS01, RS01 and three smoothed-normalised axis of GyroY (sg0zmuvGyroY,sg1zmuvGyroY and sg2zmuvGyroY). REC values were computed with embeddingparameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure is availablefrom [22]. HS01 RS01
Fig 21. DET values for horizontal arm movements.
DET values (representingpredictability and organisation of the RPs) for 20 participants performing HN and HFmovements with sensors HS01, RS01 and three smoothed-normalised axis of GyroZ(sg0zmuvGyroZ, sg1zmuvGyroZ and sg2zmuvGyroZ). DET values were computed withembedding parameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure isavailable from [22].parameters and therefore for RRSs. For which, the selection of embedding parameterswas our first challenge where we computed embedding parameters for each time seriesand then computed a sample mean over all time series in order to get two embeddingparameters to compute all RSSs with its corresponded type of movement. Then wefound that the quantification of variability with regard to the shape of the trajectoriesin RSSs requires more investigation since our original proposed method base oneuclidean metric failed to quantify those trajectories. Specially, for trajectories whichwere not well unfolded. With that in mind, we proceed to take advantage of four RQAmetrics (REC, DET, RATIO and ENTR) in order to avoid any subjectiveinterpretations or personal bias with regard to the evolution of the trajectories in RSSs.October 23, 2018 25/40 S01 RS01
Fig 22. DET values for vertical arm movements.
DET values (representingpredictability and organisation of the RPs) for 20 participants performing VN and VFmovements with sensors HS01, RS01 and three smoothed-normalised axis of GyroY(sg0zmuvGyroY, sg1zmuvGyroY and sg2zmuvGyroY). DET values were computed withembedding parameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure isavailable from [22]. HS01 RS01
Fig 23. RATIO values for horizontal arm movements.
RATIO (representingdynamic transitions) for 20 participants performing HN and HF movements withsensors HS01, RS01 and three smoothed-normalised axis of GyroZ (sg0zmuvGyroZ,sg1zmuvGyroZ and sg2zmuvGyroZ). RATIO values were computed with embeddingparameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure is availablefrom [22]. RQA metrics with fixed parameters
Considering that RQA metrics were computed with fixed embedding parameters ( m = 6and τ = 8) and recurrence thresholds ( (cid:15) = 1), we found the following. REC values,which represents the % of black points in the RPs, were more affected with and increasein normal speed movements (HN and VN) than faster movements (HF and VF) for thesensor attached to the participants (HS01). Such decrease of REC values from normalspeed to faster speed movements is also presented in data from sensor attached to therobot (RS01), and little can be said with regard to the dynamics of the time seriescoming from RS01. Similarly, DET values, representing predictability and organisationin the RPs, present little variation in the any of the time series where little can be said.In contrast, RATIO values, which represent dynamic transitions, were more variable forfaster movements (HF and VF) than normal speed movements (HN and VN) withsensors attached to the participants (HS01). For data coming from sensors attached toOctober 23, 2018 26/40 S01 RS01
Fig 24. RATIO values for vertical arm movements.
RATIO (representingdynamic transitions) for 20 participants performing VN and VF movements withsensors HS01, RS01 and three smoothed-normalised axis of GyroY (sg0zmuvGyroY,sg1zmuvGyroY and sg2zmuvGyroY). RATIO values were computed with embeddingparameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure is availablefrom [22]. HS01 RS01
Fig 25. ENTR values for horizontal arm movements.
ENTR values(representing the complexity of the deterministic structure in time series) for 20participants performing HN and HF movements with sensors HS01, RS01 and threesmoothed-normalised axis of GyroZ (sg0zmuvGyroZ, sg1zmuvGyroZ andsg2zmuvGyroZ). ENTR values were computed with embedding parameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure is available from [22].the robot (RS01), RATIO values from horizontal movements (HN, HF) appear to varymore than values coming from vertical movmentes (VN, VF). With that, it can be saidthat RATIO values can represent a bit better than REC or DET metrics for thevariability of imitation activities in each of the conditions for time series. In the sameway, ENTR values for HN were higher than values for HF and ENTR values variedmore for sensor attached to participants than ENTR values for sensors of the robot. Itis evidently that the higher the entropy the more complex the dynamics are, however,ENTR values for HN appear a bit higher than HF values, for which we believe thishappens because of the structure the time series which appear more complex for HNthan HF movements which presented a more consistence repetition.We observed that some RQA metrics are affected by the smoothness of data. Forwhich, we also explored the effect of smoothness of raw-normalised data where, forexample, REC and DET values were not affected by the smoothness of data since theseseemed to be constants. However, for RATIO values, the effect of smoothness can beOctober 23, 2018 27/40 S01 RS01
Fig 26. ENTR values for vertical arm movements.
ENTR values (representingthe complexity of the deterministic structure in time series) for 20 participantsperforming VN and VF movements with sensors HS01, RS01 and threesmoothed-normalised axis of GyroY (sg0zmuvGyroY, sg1zmuvGyroY andsg2zmuvGyroY). ENTR values were computed with embedding parameters m = 6, τ = 8 and (cid:15) = 1. R code to reproduce the figure is available from [22]. REC DET RATIO ENTR
Fig 27. 3D surfaces for RQA metrics.
3D surfaces for REC, DET, RATIO andENTR values with increasing pair embedding parameters (0 ≥ m ≤
10, 0 ≥ τ ≤
10) andrecurrence thresholds ( 0 . ≥ (cid:15) ≤ RQA metrics with different parameters
Patterns in RPs and metrics for RQA are independent of embedding dimensionparameters [10], however, that is not the case when using different recurrence thresholds.Such changes of recurrence threshold values can modify the patterns in RPs andtherefore the values of RQA metrics. We therefore computed 3D surfaces to explore thesensibility and robustness of embedding parameters and recurrence threshold in RQAmetrics. Following the same methodology of computing 3D surfaces, we also consideredvariation of window length size to present RQA metrics dependencies with embeddingparameters, recurrence thresholds and window length size.October 23, 2018 28/40 onclusions
We conclude that using a different level of smoothness for time series help us tovisualise and to quantify the variation of movements between participants using RSSs,RPs and RQA. Also, it is important to mention that some RQA’s metrics (e.g. DETand ENTR) are more robust to the effect of smoothness of time series.Althought our our experiment is limited to twenty healthy right-handed participantsof a range age of mean 19.8 SD=1.39, RQA metrics can quantify human movementvariability. With that in mind, we conclude that qunatification of human-humanoidimitation activities is possible for participants of different ages, state of health andanthropomorphic features.In general, activity type, window length and structure of the time series affects thevalues of the metrics of RQA for which certain RQA metrics are better to describedetermined type of movement. Using determined RQA metrics depends on what onewant to quantify, for instance, one can find predicability, organisation of the RPs,dynamics transitions, or complexity and determinism.Similarly, such differences in time series created differences in each of the RQAmetrics, for instance, RATIO and ENTR are helpful to distinguish differences in any ofthe categories of the time series (sensor, activity, level of smoothness and number ofparticipant), however for certain time series (data from the sensor attached to therobot) seemed to have little variations between each of the participants. The latterphenomena was in a way evidently as robot degrees of freedom did not allow it to movewith a wide range of variability.
Future Work
Inertial Sensors
To have more fundamental understating of nature of signals collected through inertialsensors in the context of human-robot interaction, we are considering to apply derivatesto the acceleration data. We can then explore the jerkiness of movements and thereforethe nature of arm movements which typically have minimum jerk [51], its relationshipwith different body parts [49, 50] or the application of higher derivatives of displacementwith respect time such as snap, crackle and pop [52].
RQA
Having presented our results with RQA metrics, we believe that further investigation isrequired to have more robust metrics. For example, Marwan et al. [32, 34] revieweddifferent aspects to compute RPs using different criteria for neighbours, different norms( L − norm , L − norm , or L ∞− norm ) or different methods to select the recurrencethreshold (cid:15) , such as using certain percentage of the signal [53], the amount of noise orusing a factor based on the standard deviation of the observational noise among manyothers [32]. Supporting information
S1 Appendix A. Examples of Uniform Time-Delay Embedding
Twoexamples are presented in this appendix: (A.1) using a 20 sample length vector, and(A.2) using a time series from an horizontal movement of a triaxial accelerometer.October 23, 2018 29/40 .1. Vector of sample length of N = 20 . Given { x n } n =1 with a sample length N = 20, we implement an uniform time-delayembedding matrix, X mτ , with embedding dimension of m = 5 and delay dimension of τ = 3 (Eq. (4)).The representation of the uniform time-delay embedding matrix X isas follows X = ˜ x n ˜ x n − ˜ x n − ˜ x n − ˜ x n − (cid:124) (19)The dimension of the uniform time-delay embedding matrix is defined by N − ( m − τ rows and m columns. N − ( m − τ is also the sample length of the delayed copies of x n which is equal to eight (20 − ((5 − ∗
3) = 8). Therefore, X can be explicitlyrepresented as X = x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x (cid:124) . (20)After transposing X , one can see that the ranges of values of the uniform time-delayembedded matrix are between (( m − τ ) + 1 to N (for this example from 13 to 20): X = x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x = X [13] X [14] X [15] X [16] X [17] X [18] X [19] X [20] . (21) A.2 Time series from a triaxial accelerometer.
For this example, it is considered a time series of a triaxial accelerometer, (Fig 32(C)),captured from repetitions of a horizontal movements trajectory (Fig 32(A, B)). FromFig 32(C)) is evidently that the A y ( n ) is the most affected axis of the accelerometer dueto horizontal movement trajectory. With that in mind, we select A y ( n ) as the inputtime series for the uniform time-delay embedding theorem. Therefore, considering thatthe sample rate of the data is 50 Hz, we have a sample length of N = 1000. We thenselected m = 7 and τ = 11 as the minimum embedding parameters for A y ( n ). Theuniform time-delay embedding matrix, A y , has a dimension of N − ( m − τ (934)rows and m (7) columns and is represented as follows: A y = A y ( n ) A y ( n − A y ( n − A y ( n − A y ( n − A y ( n − A y ( n − (cid:124) = a y (1) . . . a y (934) a y (12) . . . a y (945) a y (23) . . . a y (956) a y (34) . . . a y (967) a y (45) . . . a y (978) a y (56) . . . a y (989) a y (67) . . . a y (1000) (cid:124) , (22)October 23, 2018 30/40 y = a y (1) a y (12) a y (23) a y (34) a y (45) a y (56) a y (67)... ... ... ... ... ... ... a y (934) a y (945) a y (956) a y (967) a y (978) a y (989) a y (1000) , (23) A y = A y [67]... A y [1000] . (24) S1 Appendix B. Inertial Sensors
For this work, data were collected usingNeMEMsi sensors, which provide triaxial accelerometer, triaxial magnetometer, triaxialgyroscope and quaternions (Fig 33). NeMEMsi sensors were tested against thestate-of-the-art device MTi-30 IMU from Xsense. The comparison values betweenNeMEMsi and MTi-30 in was in terms of standard deviation of noise of each componentof Euler angles at a static state were lower than 0.1 degrees. Additionally, the NeMEMsiprovides not only a lower-power consumption but also the smaller dimensions againstother state-of-the-art brands of IMUs [42]. We refer the reader to check [42] for furtherdetails with regard to the IMU such as sample rate, power consumption, orientation,accelerometer, gyroscope, magnetometer, microprocessor, connectivity and form factor.
Acknowledgments
This work was funded by the Mexican National Council of Science and Technology(CONACyT) as part of Miguel P Xochicale’s PhD degree at University of BirminghamUK. I would like to thank to Dr. Dolores Columba Perez Flores and Prof. Martin JRussell for their helpful comments to polish the use of the language of mathematics andto Constantino Antonio Garcia Martinez et al. for developing the R packagenonlinearTseries that significantly help to accelerate the analysis of the nonlinear timeseries in this work.
Author Contributions
Conceptualisation
MX, CB
Data Curation MX Formal Analysis MX Funding Acquisition
MX, CB
Investigation MX Methodology MX Project Administration MX Resources CB Software MX Supervision CB Validation
MXOctober 23, 2018 31/40 erification MX Writing - Original Draft Preparation MX Writing - Review CB Writing - Editing MX References
1. Quintana Duque JC. Non-linear Dynamic Invariants based on EmbeddingReconstruction of Systems for Pedaling Motion. In: Byshko R, editor.Sportinformatik 2012 : 9. vom 12.- 14. Sept. 2012 in Konstanz : Beitr¨age.Universit¨at; 2012. p. 28–38.2. Quintana Duque JC. Non-linear methods for the quantification of cyclic motion[PhD dissertation]. University of Konstanz. Konstanz, Germany; 2016.3. Sam`a A, Ruiz FJ, Agell N, P´erez-L´opez C, Catal`a A, Cabestany J. Gaitidentification by means of box approximation geometry of reconstructedattractors in latent space. Neurocomputing. 2013;121:79–88.doi:10.1016/j.neucom.2012.12.042.4. Frank J, Mannor S, Precup D. Activity and Gait Recognition with Time-DelayEmbeddings. AAAI Conference on Artificial Intelligence. 2010; p. 407–408.5. G´omez-Garc´ıa JA, Godino-Llorente JI, Castellanos-Dominguez G. Non uniformEmbedding based on Relevance Analysis with reduced computational complexity:Application to the detection of pathologies from biosignal recordings.Neurocomputing. 2014;132(Supplement C):148 – 158.6. Marwan N. How to Avoid Potential Pitfalls in Recurrence Plot Based DataAnalysis. International Journal of Bifurcation and Chaos. 2011;21:1003.doi:10.1142/S0218127411029008.7. Stergiou N, Decker LM. Human movement variability, nonlinear dynamics, andpathology: Is there a connection? Human Movement Science. 2011;30(5):869 –888.8. Bradley E, Kantz H. Nonlinear time-series analysis revisited. Chaos: AnInterdisciplinary Journal of Nonlinear Science. 2015;25(9):097610.doi:10.1063/1.4917289.9. Takens F. Detecting strange attractors in turbulence. Dynamical Systems andTurbulence, Warwick 1980: Proceedings of a Symposium Held at the Universityof Warwick 1979/80. 1981; p. 366–381.10. Iwanski JS, Bradley E. Recurrence plots of experimental data: To embed or notto embed? Chaos: An Interdisciplinary Journal of Nonlinear Science.1998;8(4):861–871. doi:10.1063/1.166372.11. Garland J, Bradley E, Meiss JD. Exploring the topology of dynamicalreconstructions. Physica D: Nonlinear Phenomena. 2016;334:49 – 59.12. Kantz H, Schreiber T. Nonlinear Time Series Analysis. New York, NY, USA:Cambridge University Press; 2003.October 23, 2018 32/403. Packard NH, Crutchfield JP, Farmer JD, Shaw RS. Geometry from a Time Series.Phys Rev Lett. 1980;45:712–716. doi:10.1103/PhysRevLett.45.712.14. Aguirre LA, Letellier C. Modeling Nonlinear Dynamics and Chaos: A Review.Mathematical Problems in Engineering. 2009;11(12):1–22.15. Guastello SJ, Gregson RAM. Nonlinear Dynamical Systems Analysis for theBehavioral Science Using Real Data. CRC Press; 2011.16. Quintana-Duque J, Saupe D. Evidence of Chaos in Indoor Pedaling Motion usingNon-linear Methods. July. Routledge; 2013.17. Krakovsk´a A, Mezeiov´a K, I HB. Use of False Nearest Neighbours for SelectingVariables and Embedding Parameters for State Space Reconstruction. Journal ofComplex Systems. 2015;2015.18. Casdagli M, Eubank S, Farmer JD, Gibson J. State space reconstruction in thepresence of noise. Physica D: Nonlinear Phenomena. 1991;51(1):52 – 98.19. Gibson JF, Farmer JD, Casdagli M, Eubank S. An analytic approach to practicalstate space reconstruction. Physica D: Nonlinear Phenomena. 1992;57(1):1 – 30.20. Uzal LC, Grinblat GL, Verdes PF. Optimal reconstruction of dynamical systems:A noise amplification approach. Phys Rev E. 2011;84:016223.doi:10.1103/PhysRevE.84.016223.21. Uzal LC, Verdes PF. Optimal Irregular Delay Embeddings. In: Dynamics DaysSouth America 2010. National Institute for Space Research; 2010.Available from: http://urlib.net/8JMKD3MGP7W/3824342 .22. Xochicale MP. How Well You Move Dataset; 2018. Available from: https://github.com/mxochicale/hwum-dataset .23. Kennel MB, Brown R, Abarbanel HDI. Determining embedding dimension forphase-space reconstruction using a geometrical construction. Phys Rev A.1992;45:3403–3411. doi:10.1103/PhysRevA.45.3403.24. Cao L. Practical method for determining the minimum embedding dimension of ascalar time series. Physica D: Nonlinear Phenomena. 1997;110:43–50.25. Emrani S, Gentimis T, Krim H. Persistent Homology of Delay Embeddings andits Application to Wheeze Detection. IEEE Signal Processing Letters.2014;21(4):459–463. doi:10.1109/LSP.2014.2305700.26. Garcia SP, Almeida JS. Nearest neighbor embedding with different time delays.Phys Rev E. 2005;71:037204. doi:10.1103/PhysRevE.71.037204.27. Fraser AM, Swinney HL. Independent coordinates for strange attractors frommutual information. Phys Rev A. 1986;33:1134–1140.28. Kabiraj L, Saurabh A, Wahi P, Sujith RI. Route to chaos for combustioninstability in ducted laminar premixed flames. Chaos: An InterdisciplinaryJournal of Nonlinear Science. 2012;22(2):023129. doi:10.1063/1.4718725.29. Garcia CA, Sawitzki G. nonlinearTseries: Nonlinear Time Series Analysis; 2016.Available from: https://github.com/constantino-garcia/nonlinearTseries .October 23, 2018 33/400. Pecora LM, Moniz L, Nichols J, Carroll TL. A unified approach to attractorreconstruction. Chaos: An Interdisciplinary Journal of Nonlinear Science.2007;17(1):013110. doi:10.1063/1.2430294.31. Small M, Tse CK. Optimal embedding parameters: a modelling paradigm.Physica D: Nonlinear Phenomena. 2004;194(3):283 – 296.doi:https://doi.org/10.1016/j.physd.2004.03.006.32. Marwan N, Romano MC, Thiel M, Kurths J. Recurrence plots for the analysis ofcomplex systems. Physics Reports. 2007;438(5):237 – 329.doi:https://doi.org/10.1016/j.physrep.2006.11.001.33. Eckmann JP, Kamphorst SO, Ruelle D. Recurrence Plots of Dynamical Systems.EPL (Europhysics Letters). 1987;4(9):973.34. Marwan N, Webber CL. Mathematical and Computational Foundations ofRecurrence Quantifications. In: Webber CL, Marwan N, editors. RecurrenceQuantification Analysis: Theory and Best Practices. 1st ed. Springer, Cham;2015. p. 3–43.35. Zbilut JP, Webber CL. Embeddings and delays as derived from quantification ofrecurrence plots. Physics Letters A. 1992;171(3):199 – 203.doi:https://doi.org/10.1016/0375-9601(92)90426-M.36. Trulla LL, Giuliani A, Zbilut JP, Webber CL. Recurrence quantification analysisof the logistic equation with transients. Physics Letters A. 1996;223(4):255 – 260.doi:https://doi.org/10.1016/S0375-9601(96)00741-4.37. Marwan N. A historical review of recurrence plots. The European PhysicalJournal Special Topics. 2008;164(1):3–12. doi:10.1140/epjst/e2008-00829-1.38. Apthorp D, Nagle F, Palmisano S. Chaos in Balance: Non-Linear Measures ofPostural Control Predict Individual Variations in Visual Illusions of Motion.PLOS ONE. 2014;9(12):1–22. doi:10.1371/journal.pone.0113897.39. Duran ND, Fusaroli R. Conversing with a devil’s advocate: Interpersonalcoordination in deception and disagreement. PLOS ONE. 2017;12(6):1–25.doi:10.1371/journal.pone.0178140.40. Rhea CK, Silver TA, Hong SL, Ryu JH, Studenka BE, Hughes CML, et al. Noiseand Complexity in Human Postural Control: Interpreting the DifferentEstimations of Entropy. PLOS ONE. 2011;6(3):1–9.doi:10.1371/journal.pone.0017696.41. Gouaillier D, Hugel V, Blazevic P, Kilner C, Monceaux J, Lafourcade P, et al.Mechatronic Design of NAO Humanoid. In: Proceedings of the 2009 IEEEInternational Conference on Robotics and Automation. ICRA’09. Piscataway, NJ,USA: IEEE Press; 2009. p. 2124–2129. Available from: http://dl.acm.org/citation.cfm?id=1703775.1703795 .42. Comotti D, Galizzi M, Vitali A. neMEMSi: One step forward in wireless attitudeand heading reference systems. In: 2014 International Symposium on InertialSensors and Systems (ISISS); 2014. p. 1–4.43. Shoaib M, Bosch S, Incel OD, Scholten H, Havinga PJM. Complex HumanActivity Recognition Using Smartphone and Wrist-Worn Motion Sensors. Sensors.2016;16(4). doi:10.3390/s16040426.October 23, 2018 34/404. Anguita D, Ghio A, Oneto L, Parra X, Reyes-Ortiz JL. A Public DomainDataset for Human Activity Recognition Using Smartphones. In: 21th EuropeanSymposium on Artificial Neural Networks, Computational Intelligence andMachine Learning, ESANN 2013; 2013.45. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C(2Nd Ed.): The Art of Scientific Computing. New York, NY, USA: CambridgeUniversity Press; 1992. Available from: .46. Savitzky A, Golay MJE. Smoothing and differentiation of data by simplified leastsquares procedures. Analytical Chemistry. 1964;36:1627–1639.47. Schafer RW. On the frequency-domain properties of Savitzky-Golay filters. In:2011 Digital Signal Processing and Signal Processing Education Meeting(DSP/SPE); 2011. p. 54–59.48. signal R developers. signal: Signal processing; 2014. Available from: http://r-forge.r-project.org/projects/signal/ .49. de Vries JIP, Visser GHA, Prechtl HFR. The emergence of fetal behaviour. I.Qualitative aspects. Early Human Development. 1982;7(4):301 – 322.doi:https://doi.org/10.1016/0378-3782(82)90033-0.50. Mori H, Kuniyoshi Y. Is the developmental order of fetal behaviors self-organizedin an uterine environment? In: 2012 IEEE International Conference onDevelopment and Learning and Epigenetic Robotics (ICDL); 2012. p. 1–2.51. Flash T, Hogan N. The coordination of arm movements: an experimentallyconfirmed mathematical model. Journal of Neuroscience. 1985;5(7):1688–1703.doi:10.1523/JNEUROSCI.05-07-01688.1985.52. Eager D, Pendrill AM, Reistad N. Beyond velocity and acceleration: jerk, snapand higher derivatives. European Journal of Physics. 2016;37(6):065008.53. Letellier C. Estimating the Shannon Entropy: Recurrence Plots versus SymbolicDynamics. Phys Rev Lett. 2006;96:254102. doi:10.1103/PhysRevLett.96.254102.October 23, 2018 35/40 N REC H F V N V F DETRATIOENTR H N H F V N V F RECDETRATIOENTR H S R S Fig 28. 3D surfaces of RQA metrics for sensors and activities.
3D surfaceswith increasing embedding parameters and recurrence thresholds are for HS01 and RS01sensors of HN, HF, VN and VF activities. RQA metrics values are for time series ofparticipant p01 for sensors (HS01 and RS01), activities (HN, HF, VN and VF) and forsg0zmuvGyroZ axis with 500 samples window length. R code to reproduce the figure isavailable from [22].October 23, 2018 36/40 R E C w250 w500 w750 D E T R A T I O E N T R Fig 29. 3D surfaces for RQAs metrics with four window lengths.
3D surfacesof RQAs metric values with increasing embedding parameters and recurrence thresholdsare for four window lengths (w100, w250, w500 and w750). RQA metrics values are fortime series of participant p01 using HS01 sensor, HN activity and sg0zmuvGyroZ axis.R code to reproduce the figure is available from [22].October 23, 2018 37/40 g REC s g s g DETRATIOENTR
Fig 30. 3D surfaces for RQA metrics with three levels of smoothness. REC p02p03
DETRATIOENTR
Fig 31. 3D surfaces for RQA metrics with three participants.
3D surfaces ofRQAs metric values for participants p01, p02 and p03 with increasing embeddingparameters and recurrence thresholds. RQA metrics values are for time series of HS01sensor, HN activity and 500 samples window length. R code to reproduce the figure isavailable from [22].October 23, 2018 39/40 ig 32. Time series for horizontal movement trajectory. (A). Triaxialaccelerometer (in red) is moved repetitively across a line of 251 mm from point a to b and then from b to a . The points a and b indicate when a click sound is produced. (B).Person’s hand holding and moving the sensor horizontally across the line. (C). Timeseries for the triaxial accelerometer ( A x ( n ), A y ( n ), A z ( n )) for ten repetitive horizontalmovements across a line. The top time series only shows A y axis which corresponds toone cycle of the horizontal movement and the black arrows represent the movement’sdirection of the accelerometer with respect to the produced time series. R code toreproduce the figure is available from [22]. A B C D