**3. Analysis of the surface EMG signals based on chaos theory**

The discovery of chaotic phenomena is the third major breakthrough in the 20th century physics scientific community following the creation of relativity and quantum mechanics. It organically combines the two major theoretical systems of determinism and probabilism that have long been debated to create a scientific model of a new paradigm, so that people can use some simple rules to explain seemingly stochastic information in the past[17-19]. The practical significance that finds chaotic phenomena is to recognize that a deterministic nonlinear system can have inherent uncertainty. Perhaps a system has only a few degrees of freedom, but it can produce complex, similar to the random output signal. In the past, one could only denote a random-looking data as a random process from the view of the traditional time series analysis. The statistical methods or random time series models were used to analyze the data. Since the chaotic phenomenon was discovered by Lorenz[17], people have begun to reunderstand and restudy these random-looking signals so as to reveal the inherent deterministic mechanisms of these signals. That is, it is to explore that the systems which generate these signals may contain essentially deterministic characteristics. Chaos phenomenon breaks the path that the regularity is found in a lot of completely different systems. This will lead to a revolution in the field of influence of various disciplines. It is chaos to lead people to explore the complexity in nature.

At present, the idea of Chaos has been introduced into the analysis of time series to create the field of chaotic time series analysis. Since the inception of chaotic time series analysis, it has quickly been penetrated into other disciplines and engineering fields. Thus it becomes the most active branch of the modern nonlinear dynamics. This section describes the chaos definition and the phase space reconstruction of chaotic time series, discusses some parameters that are used to analysis chaotic time series, such as the correlation dimension and Lyapunov exponent, study the principal component analysis methods based on SVD, and propose the symplectic principal component method based on symplectic geometry. Then we use these methods to investigate the surface EMG signals.

#### **3.1. Chaos and its definition**

Computational Intelligence in Electromyography Analysis – 130 A Perspective on Current Applications and Future Challenges

*C* (*r*) *nl orig*

signal, *C* (*r*) *nl*

nonlinear dynamic properties.

a. The action surface EMG signal b. The fatigue surface EMG signal

In order to detect the nonlinearity of the action surface EMG signal, 39 FT-based surrogate data are used according to the null hypothesis 3. The generated surrogate data contain the linear properties of the raw data. Figure 9 is the analysis of surface EMG signal based on VWK with surrogate data. We can see that no matter whether it is the action or fatigue EMG

significance. The results illuminate that the action and fatigue surface EMG signals contain

The discovery of chaotic phenomena is the third major breakthrough in the 20th century physics scientific community following the creation of relativity and quantum mechanics. It organically combines the two major theoretical systems of determinism and probabilism that have long been debated to create a scientific model of a new paradigm, so that people can use some simple rules to explain seemingly stochastic information in the past[17-19]. The practical significance that finds chaotic phenomena is to recognize that a deterministic nonlinear system can have inherent uncertainty. Perhaps a system has only a few degrees of freedom, but it can produce complex, similar to the random output signal. In the past, one could only denote a random-looking data as a random process from the view of the traditional time series analysis. The statistical methods or random time series models were used to analyze the data. Since the chaotic phenomenon was discovered by Lorenz[17], people have begun to reunderstand and restudy these random-looking signals so as to reveal the inherent deterministic mechanisms of these signals. That is, it is to explore that the systems which generate these signals may contain essentially deterministic characteristics. Chaos phenomenon breaks the path that the regularity is found in a lot of completely different systems. This will lead to a revolution in the field of influence of

*surr* , \* is

*surr* . The null hypothesis 3 can be rejected in 95%

**Figure 9.** VWK combined with surrogate analysis of surface EMG signal. solid is *C* (*r*) *nl*

*2.2.3. Analysis of SEMG based on VWK method with surrogate data[16]* 

**3. Analysis of the surface EMG signals based on chaos theory** 

various disciplines. It is chaos to lead people to explore the complexity in nature.

*orig* is always smaller than *C* (*r*) *nl*

Chaos is "order in disorder". The order means its deterministic nature. The disorder means that the final results can be unpredictable for a long time. As a scientific concept, chaos generally denotes that the long-term dynamical behavior of a deterministic nonlinear system manifests as a random-looking behavior. Mathematically speaking, "chaos" has not been a unified strict definition. For the definition of chaos, there are at least nine different definitions, where the three definitions given by Li-Yorke, Devaney, Marotto are more commonly used. Here describes the definition of chaos by Li-Yorke[18].

**Li-Yorke Theorem:** Let *f* (*x*) as a continuous self-map in[*a*, *b*] . If *f* (*x*) has a periodic point with period 3, then for any positive integer *n*=1, 2, 3, …, there is a periodic point with period *n*.

This is the famous period 3 theorem. It becomes a milestone in the development history of chaos theory and promotes the creation and development of chaos theory. From this theorem, the first formal mathematical definition of the chaos is given.

**Chaos definition:** Let *f* (*x*) as a continuous self-map in closed interval *I*, i.e.

$$f: I \to I \subset \mathbb{R}^m, \quad \mathbf{y} = f(\mathbf{x}) \tag{20}$$

where *x* , *y*∈*I* . If *m* =1, *f* is one-dimensional mapping. If *m* ≠1, *f* is multi-dimensional mapping. Denote the *n* times iteration of *f* as *f* (*x*) *<sup>n</sup>* . If Eq.20 satisfies the following conditions, then it has chaotic motion:


$$\lim\_{n \to \infty} \sup \left| f''(\mathbf{x}) - f''(\mathbf{y}) \right| > 0, \qquad \forall \mathbf{x}, \mathbf{y} \in S, \quad \mathbf{x} \neq \mathbf{y} \tag{21}$$

$$\lim\_{n \to \infty} \inf \left| f''(x) - f''(y) \right| = 0, \qquad \forall x, y \in S \tag{22}$$

$$\lim\_{n \to \infty} \sup \left| f''(\mathbf{x}) - f''(p) \right| > 0, \qquad \forall \mathbf{x} \in S, \quad \forall p \in P(f) \tag{23}$$

where *P*( *f* )Δ {*x* | *x* is a periodic point of *f*}.

This definition explains "existence" of chaos in mathematics. According to the above theorem and the definition, the description of chaotic motion is different from the general periodic and quasi-periodic motion. Its motion is not a single periodic orbit but an envelope for a bunch of tracks, where the infinite number of countable stable periodic orbits and uncountable stable aperiodic orbits are embedded densely. Meanwhile, there is at least one unstable aperiodic orbit. Overall, the chaos not only contains some inherent regularity, but also shows that the system has ergodicity. That is, the system has a long-term unpredictability. In other words, the long-term behavior of the system can not be predicted if the system displays the so-called "sensitive dependence on initial conditions". The meaning of this definition is that the aperiodicity of chaotic system is exhibited accurately. For a dynamical system, the observable behaviour was called stochastic in the past. In fact, it can be random-looking, i.e. "stochastic behaviour occurring in a deterministic system". Therefore, it is challenging to quantitatively describe the nature of chaotic dynamics and distinguish between the so-called random and chaotic motions from a time series, especially from an experimental time series. At present, chaotic time series analysis methods have been widely attention in fields of mathematics, physics, biology, biomedicine, robotics, geology, engineering, economics, finance, and so on.

Nonlinear Analysis of Surface EMG Signals 133

, …, ( ( )) <sup>2</sup>

*t m <sup>d</sup>* − τ=

<sup>2</sup> *x x y <sup>m</sup>*

ϕ.

is a smooth diffeomorphism. *<sup>t</sup> y* denotes the state of the

. If the embedding dimension is greater than 2

*x*(*t*) = *h*(*Y*) (24)

τ

*T*

)) (25)

, called as the lag time or delay

ϕ

is an embedding of → <sup>2</sup>*m*+<sup>1</sup> *M R* . The manifold *M* is

is the delay time. The signal observed in *M* at time *t* consist of

→ <sup>2</sup>*m*+<sup>1</sup> *M R* . For a practical time series { }*<sup>t</sup> x* , the state of the original system is equivalent to the *m*-dimensional manifold *M*. In fact, { }*<sup>t</sup> x* is a signal observed in the *m*-dimensional

times the dimension *m* of the attractor of the original system, the phase space with the base of the practical signal delay time coordinates is equivalent to the state space of the original system. That is, Takens' embedding theorem states that if the time series is indeed composed of scalar measurements of the state from a dynamical system, then under certain genericity assumptions, a one-to-one image of the original set {*x*} is given by the time-delay embedding, provided *d* is large enough. At present, the delay coordinate method has been

*<sup>d</sup> <sup>d</sup> <sup>d</sup> X <sup>t</sup>* = (*x*(*t*), *x*(*t* −

τ

time. *d* is and embedding dimension, *d*≥2*m*+1. *m* is an attractor dimension of the original

Takens' embedding theorem offers in the absence of noise, the possibility of reconstructing *n*-dimensional dynamics from one-dimensional infinite data of one observable-measurable system. This means that in the case of any delay time, a time series can always be embedded into the state space of the system, and when the embedding dimension is sufficiently large, reconstructed space and embedded space is almost one-to-one correspondence. Therefore, one can reconstruct a phase space from an experimental time series so as to estimate dynamical invariants of the time series, such as dimensions, Lyapunov exponents, entropies[21, 23, 24] and so on. However, the embedding theorem does not directly answer how to choose embedding dimension *d* and delay time *t*. In practical application, the experimental data is always limited and noisy so that the estimation of the above parameters presents some difficulties[25, 26]. Accuracy of the phase space reconstruction is critically important to the estimation of invariant measures characterizing system behavior.

),, *x*(*t* − (*d* − 1)

τ

and embedding dimension *d* always has a great impact on the

 , where ( ) *<sup>t</sup> <sup>t</sup> x* = *x y* , *x x*( *y* ) *x*( ( *y*)) *<sup>d</sup> <sup>t</sup>*<sup>−</sup>τ = *<sup>t</sup>*<sup>−</sup>τ=

τ

widely used to give the phase space reconstruction from the original signal.

τ), *x*(*t* − 2

is an integer multiple of the sampling interval

manifold *M*. If let *<sup>d</sup> <sup>t</sup> <sup>t</sup> y y* ϕ

( , , , ) ( , ) <sup>2</sup> *<sup>d</sup> <sup>m</sup> <sup>d</sup> <sup>x</sup> <sup>t</sup> <sup>t</sup> <sup>t</sup> x x x*

τ<sup>−</sup>

τ

system *M* at time *t*. *<sup>d</sup>*

{ , , , } <sup>2</sup> *<sup>d</sup> <sup>m</sup> <sup>d</sup> <sup>t</sup> <sup>t</sup> <sup>t</sup> x x x* <sup>−</sup>τ<sup>−</sup>

= <sup>−</sup>

where *<sup>d</sup>* τ

system.

φ ϕ → <sup>−</sup>

τ

τ<sup>−</sup>

For a time series *x*(*t*) observed by the measure function *h*, i.e.

τ

diffeomorphic with { , , , } <sup>2</sup> *<sup>d</sup> <sup>m</sup> <sup>d</sup> <sup>t</sup> <sup>t</sup> <sup>t</sup> x x x* <sup>−</sup>

the vector *X <sup>t</sup>* can be constructed as follows,

*3.2.2. Problems in phase space reconstruction* 

τ

The choice of delay time *<sup>d</sup>*

phase space reconstruction.

τ: ,

ϕ

### **3.2. Phase space reconstruction theory**

#### *3.2.1. Phase space reconstruction*

Phase space reconstruction is generally the first step of chaotic time series analysis from a time series data. The dynamic characteristic of the system can be explored through phase space reconstruction of the original time series so that the mechanism of the original system can be revealed from the original time series[20]. It has been proved by the so-called Takens' embedding theorem[21]. According to the theorem, the reconstructed phase space can maintain the invariance of geometry for the original dynamical system[22], such as the characteristic value of the fixed point, the fractal dimension of the attractor, the Lyapunov exponent in the phase space orbit, and so on.

**Definition 1:** Let (*N*, ρ) and ( , ) *N*<sup>1</sup> ρ<sup>1</sup> as two metric spaces. If the mapping <sup>1</sup> ϕ :*N* → *N* satisfies ①ϕ is a surjection; ② ( , ) ( , ) <sup>1</sup> ρ *x y* = ρ ϕ*x* ϕ*y* , then (*N*, ρ) and ( , ) *N*<sup>1</sup> ρ<sup>1</sup> are called as the isometric isomorphism.

**Definition 2:** If ( , ) *N*<sup>1</sup> ρ<sup>1</sup> is isometric isomorphism with a subspace ( , ) *N*<sup>0</sup> ρ <sup>2</sup> of another metric space ( , ) *N*<sup>2</sup> ρ <sup>2</sup> , then ( , ) *N*<sup>1</sup> ρ<sup>1</sup> can be embedded into ( , ) *N*<sup>2</sup> ρ<sup>2</sup> .

**Theorem 1:** Let *M* be a compact manifold of dimension *m*. For pairs (ϕ, *x*), ϕ : *M* → *M* a smooth diffeomorphism and *x* : *M* → *R* a smooth function, it is a generic property that the map <sup>2</sup> <sup>1</sup> ( , ) : → *<sup>m</sup>*<sup>+</sup> φ ϕ *<sup>x</sup> M R* is an embedding, where ( ) { ( ), ( ( )), , ( ( ))} <sup>2</sup> ( , ) *y x y x y x y <sup>m</sup>* φ ϕ *<sup>x</sup>* = ϕ ϕ.

In terms of the above definitions and theorem, ( ) { ( ), ( ( )), , ( ( ))} <sup>2</sup> ( , ) *y x y x y x y <sup>m</sup>* φ ϕ *<sup>x</sup>* = ϕ ϕ is a subspace of <sup>2</sup>*m*+<sup>1</sup> *R* . φ (ϕ,*x*) and the subspace are isometric isomorphic. Then the manifold *M* of dimension *m* can be embedded into <sup>2</sup>*m*+<sup>1</sup> *R* . In other words, φ (ϕ, *<sup>x</sup>* ) is an embedding of

→ <sup>2</sup>*m*+<sup>1</sup> *M R* . For a practical time series { }*<sup>t</sup> x* , the state of the original system is equivalent to the *m*-dimensional manifold *M*. In fact, { }*<sup>t</sup> x* is a signal observed in the *m*-dimensional manifold *M*. If let *<sup>d</sup> <sup>t</sup> <sup>t</sup> y y* ϕ → <sup>−</sup>τ : , ϕ is a smooth diffeomorphism. *<sup>t</sup> y* denotes the state of the system *M* at time *t*. *<sup>d</sup>* τ is the delay time. The signal observed in *M* at time *t* consist of { , , , } <sup>2</sup> *<sup>d</sup> <sup>m</sup> <sup>d</sup> <sup>t</sup> <sup>t</sup> <sup>t</sup> x x x* <sup>−</sup>τ <sup>−</sup> τ , where ( ) *<sup>t</sup> <sup>t</sup> x* = *x y* , *x x*( *y* ) *x*( ( *y*)) *<sup>d</sup> <sup>t</sup>*<sup>−</sup>τ = *<sup>t</sup>*<sup>−</sup>τ = ϕ , …, ( ( )) <sup>2</sup> <sup>2</sup> *x x y <sup>m</sup> t m <sup>d</sup>* − τ = ϕ . ( , , , ) ( , ) <sup>2</sup> *<sup>d</sup> <sup>m</sup> <sup>d</sup> <sup>x</sup> <sup>t</sup> <sup>t</sup> <sup>t</sup> x x x* φ ϕ = <sup>−</sup>τ <sup>−</sup> τ is an embedding of → <sup>2</sup>*m*+<sup>1</sup> *M R* . The manifold *M* is diffeomorphic with { , , , } <sup>2</sup> *<sup>d</sup> <sup>m</sup> <sup>d</sup> <sup>t</sup> <sup>t</sup> <sup>t</sup> x x x* <sup>−</sup>τ <sup>−</sup> τ . If the embedding dimension is greater than 2 times the dimension *m* of the attractor of the original system, the phase space with the base of the practical signal delay time coordinates is equivalent to the state space of the original system. That is, Takens' embedding theorem states that if the time series is indeed composed of scalar measurements of the state from a dynamical system, then under certain genericity assumptions, a one-to-one image of the original set {*x*} is given by the time-delay embedding, provided *d* is large enough. At present, the delay coordinate method has been widely used to give the phase space reconstruction from the original signal.

For a time series *x*(*t*) observed by the measure function *h*, i.e.

$$\mathbf{x}(t) = h(Y) \tag{24}$$

the vector *X <sup>t</sup>* can be constructed as follows,

Computational Intelligence in Electromyography Analysis – 132 A Perspective on Current Applications and Future Challenges

engineering, economics, finance, and so on.

**3.2. Phase space reconstruction theory** 

exponent in the phase space orbit, and so on.

ρ

φ (ϕ

<sup>2</sup> , then ( , ) *N*<sup>1</sup>

ρ

ρ

) and ( , ) *N*<sup>1</sup> ρ

 is a surjection; ② ( , ) ( , ) <sup>1</sup> ρ*x y* =

ρ

of dimension *m* can be embedded into <sup>2</sup>*m*+<sup>1</sup> *R* . In other words,

**Theorem 1:** Let *M* be a compact manifold of dimension *m*. For pairs (

ρ ϕ*x* ϕ

*3.2.1. Phase space reconstruction* 

**Definition 1:** Let (*N*,

ϕ

isometric isomorphism.

**Definition 2:** If ( , ) *N*<sup>1</sup>

metric space ( , ) *N*<sup>2</sup>

map <sup>2</sup> <sup>1</sup> ( , ) : → *<sup>m</sup>*<sup>+</sup>

subspace of <sup>2</sup>*m*+<sup>1</sup> *R* .

φ ϕ

satisfies ①

This definition explains "existence" of chaos in mathematics. According to the above theorem and the definition, the description of chaotic motion is different from the general periodic and quasi-periodic motion. Its motion is not a single periodic orbit but an envelope for a bunch of tracks, where the infinite number of countable stable periodic orbits and uncountable stable aperiodic orbits are embedded densely. Meanwhile, there is at least one unstable aperiodic orbit. Overall, the chaos not only contains some inherent regularity, but also shows that the system has ergodicity. That is, the system has a long-term unpredictability. In other words, the long-term behavior of the system can not be predicted if the system displays the so-called "sensitive dependence on initial conditions". The meaning of this definition is that the aperiodicity of chaotic system is exhibited accurately. For a dynamical system, the observable behaviour was called stochastic in the past. In fact, it can be random-looking, i.e. "stochastic behaviour occurring in a deterministic system". Therefore, it is challenging to quantitatively describe the nature of chaotic dynamics and distinguish between the so-called random and chaotic motions from a time series, especially from an experimental time series. At present, chaotic time series analysis methods have been widely attention in fields of mathematics, physics, biology, biomedicine, robotics, geology,

Phase space reconstruction is generally the first step of chaotic time series analysis from a time series data. The dynamic characteristic of the system can be explored through phase space reconstruction of the original time series so that the mechanism of the original system can be revealed from the original time series[20]. It has been proved by the so-called Takens' embedding theorem[21]. According to the theorem, the reconstructed phase space can maintain the invariance of geometry for the original dynamical system[22], such as the characteristic value of the fixed point, the fractal dimension of the attractor, the Lyapunov

<sup>1</sup> as two metric spaces. If the mapping <sup>1</sup>

ρ<sup>2</sup> .

( , ) *y x y x y x y <sup>m</sup>*

ϕ

φ (ϕ

( , ) *y x y x y x y <sup>m</sup>*

ϕ

) and ( , ) *N*<sup>1</sup> ρ

ρ

*y* , then (*N*,

<sup>1</sup> is isometric isomorphism with a subspace ( , ) *N*<sup>0</sup>

φ ϕ*<sup>x</sup>* =

,*x*) and the subspace are isometric isomorphic. Then the manifold *M*

<sup>1</sup> can be embedded into ( , ) *N*<sup>2</sup>

smooth diffeomorphism and *x* : *M* → *R* a smooth function, it is a generic property that the

*<sup>x</sup> M R* is an embedding, where ( ) { ( ), ( ( )), , ( ( ))} <sup>2</sup>

In terms of the above definitions and theorem, ( ) { ( ), ( ( )), , ( ( ))} <sup>2</sup>

φ ϕ*<sup>x</sup>* = ϕ:*N* → *N*

ρ

ϕ.

ϕ

, *<sup>x</sup>* ) is an embedding of

ϕ, *x*), ϕ

<sup>1</sup> are called as the

<sup>2</sup> of another

: *M* → *M* a

is a

$$\overline{X}\_{d} = \left(\mathbf{x}(t), \mathbf{x}(t-\tau\_{d}), \mathbf{x}(t-2\tau\_{d}), \dots, \mathbf{x}(t-(d-1)\tau\_{d})\right)^{T} \tag{25}$$

where *<sup>d</sup>* τ is an integer multiple of the sampling interval τ , called as the lag time or delay time. *d* is and embedding dimension, *d*≥2*m*+1. *m* is an attractor dimension of the original system.

#### *3.2.2. Problems in phase space reconstruction*

Takens' embedding theorem offers in the absence of noise, the possibility of reconstructing *n*-dimensional dynamics from one-dimensional infinite data of one observable-measurable system. This means that in the case of any delay time, a time series can always be embedded into the state space of the system, and when the embedding dimension is sufficiently large, reconstructed space and embedded space is almost one-to-one correspondence. Therefore, one can reconstruct a phase space from an experimental time series so as to estimate dynamical invariants of the time series, such as dimensions, Lyapunov exponents, entropies[21, 23, 24] and so on. However, the embedding theorem does not directly answer how to choose embedding dimension *d* and delay time *t*. In practical application, the experimental data is always limited and noisy so that the estimation of the above parameters presents some difficulties[25, 26]. Accuracy of the phase space reconstruction is critically important to the estimation of invariant measures characterizing system behavior. The choice of delay time *<sup>d</sup>* τ and embedding dimension *d* always has a great impact on the phase space reconstruction.

Some researchers have studied the choice of delay time *<sup>d</sup>* τ [26-30]. If the delay time *<sup>d</sup>* τ is too small, the reconstructed attractor will be crowded around the main diagonal, which is called as redundance. If *<sup>d</sup>* τ is too large, the dynamic shape of the attractor will be broken, which is called as irrelevance, and the phase space reconstruction is no longer representative of the true dynamics in the real system[28]. In normal circumstances, in order to make the elements of *Yt* independent, *<sup>d</sup>* τ is the same for all the embedding dimension *d*[27-29]. The autocorrelation function method and mutual information technique[30] have been most commonly used to give the delay time *<sup>d</sup>* τ although the issue of the delay time choice has still been completely resolved.

Nonlinear Analysis of Surface EMG Signals 135

(27)

μ

, the correlation

μ

, then for

*<sup>l</sup>*<sup>→</sup> <sup>=</sup> (28)

= − − = = *n i*

with *l*→0, there are a scaling relation <sup>2</sup> ( ) *<sup>D</sup> C l* ∝ *l* between the correlation integral *C*(*l*) and *l*. A

In practical computation, *D*2 is the slope of log *C* vs log *l* curve over a selected straight line

**Theorem 2:** Let a map *<sup>n</sup> <sup>n</sup> G* : *R* → *R* . *A* is an attractor of map *G* that has only a finite number

( ) [ ( ), ( ( )), , ( ( ))] <sup>1</sup> ( 1) *F X h X h G X h G X <sup>d</sup>*

The theorem says that with the embedding dimension increasing, the slope of corresponding correlation integral curve will converge to the correlation dimension *D*<sup>2</sup> of the original system attractor. Therefore, the optimal embedding dimension can be estimated by using the correlation dimension *D*<sup>2</sup> . That is, if the embedding dimension *D*<sup>2</sup> *d* ≥ , the slope of the correlation integral curve is equal to the correlation dimension. This also indicates that the dimension estimation actually does not have to meet the requirements of the embedding theorem on the embedding dimension 2 1 *d* ≥ *D*<sup>2</sup> + . When the embedding dimension *D*<sup>2</sup> *d* ≥ , the reconstructed attractor can contain the fractal structural feature of the original system attractor to reflect the chaotic characteristics of the original system.

Correlation dimension has been widely used in the analysis of chaotic time series.

Chaos has a fractal structure so that the corresponding correlation dimension *D*2 is a fractional value. The estimation of correlation dimension *D*2 from a time series can be used to determine whether the time series is chaotic. If *D*<sup>2</sup> is fractional, the original time series can have chaotic features, otherwise, it cannot be chaotic. According to the correlation dimension theorem, when the embedding dimension *d* of the reconstructed phase space is increased to a certain value, the correlation dimension *D*2 will be saturated. Then, the

is a natural probability measure in *Rd* of ( ) *hF A* . If 2 *d D*≥ ( )

. For a measure function *h* of *A*, *h R R* : *<sup>n</sup>* → , define a delay

<sup>−</sup> <sup>−</sup> <sup>−</sup> = (29)

of periodic points of period *P*. Under a natural probability measure

ln( ) ln ( ) lim <sup>0</sup> <sup>2</sup> *l <sup>C</sup> <sup>l</sup> <sup>D</sup>*

*n j <sup>X</sup> <sup>i</sup> <sup>X</sup> <sup>j</sup> <sup>l</sup> <sup>n</sup> <sup>C</sup> <sup>l</sup>* 1 1 <sup>2</sup> ( ) <sup>1</sup> ( ) θ

correlation dimension *D*2 is defined.

*3.3.2. Correlation dimension theorem[41]* 

μ

*h*

μ.

μ=

*3.3.3. Chaotic test based on correlation dimension* 

dimension of *A* is ( ) *D*<sup>2</sup>

where *d P* ≥ . ( ) *hF*

coordinate map *<sup>n</sup> <sup>d</sup> Fh* : *R* → *R* as

μ

almost every *h*, 2 2 ( ( )) ( ) *DF D <sup>h</sup>*

range.

For the embedding dimension *d*, there are three methods that are usually used to choose the appropriate embedding dimension, including the correlation dimension, singular value decomposition(SVD), the false neighbors[21, 31, 32]. The correlation dimension method is to estimate appropriate dimension *d* in terms of the correlation theorem[8, 21, 33]. By increasing the embedding dimension, one notes an appropriate dimension d when the value of the correlation dimension stops changing. Broomhead and King[31] used the singular value decomposition (SVD) technique to determine an appropriate embedding dimension *d* directly from the raw time series. The false neighbor method is based on the fact that choosing a too low embedding dimension results in points, which are far apart in the original phase space, being moved closer together in the reconstruction space[32]. Besides, there are also some other methods and modified extensions developed based on the above methods. However, there are still problems on how to determine the appropriate embedding dimension from a scalar time series[34-38].

#### **3.3. Correlation dimension**

If a system is chaotic, the strange attractor in a region of the phase space constitutes an infinite hierarchy of self-similar structure, i.e. a fractal structure. One can use quantitative measures to define the fractal nature. The correlation dimension is a useful measurement. Grassberger and Proccacia give a kind of computation method, called GP algorithm[33, 39, 40].

#### *3.3.1. GP algorithm of correlation dimension*

Let *X*1, *X2*, ..., *Xn* be a point of the attractor in phase space. *Cl*(*Xj*) is denoted as a hypersurface sphere with the radius *l* at the reference point *Xj*. μ[*Cl*(*Xi*)] is the probability that *Xi* (*i*=1, ..., *n*) falls into *Cl*(*Xj*), as follows.

$$\mu[C\_I(X\_{\\_j})] = \frac{1}{n} \sum\_{i=1}^n \theta(l - \left\|X\_i - X\_{\\_j}\right\|) \tag{26}$$

where • is Euclidean norm. θ(*r*) is Heaviside function whose value is 1 if *r*≥0, otherwise, zero.

Then a correlation integral function is defined as

Nonlinear Analysis of Surface EMG Signals 135

$$C(I) = \frac{1}{n^2} \sum\_{i=1}^{n} \sum\_{j=1}^{n} \Theta(I - \left\|{X\_i - X\_j}\_j\right\|)\tag{27}$$

with *l*→0, there are a scaling relation <sup>2</sup> ( ) *<sup>D</sup> C l* ∝ *l* between the correlation integral *C*(*l*) and *l*. A correlation dimension *D*2 is defined.

$$D\_2 = \lim\_{l \to 0} \frac{\ln C(l)}{\ln(l)} \tag{28}$$

In practical computation, *D*2 is the slope of log *C* vs log *l* curve over a selected straight line range.

#### *3.3.2. Correlation dimension theorem[41]*

Computational Intelligence in Electromyography Analysis – 134 A Perspective on Current Applications and Future Challenges

as redundance. If *<sup>d</sup>*

τ

commonly used to give the delay time *<sup>d</sup>*

elements of *Yt* independent, *<sup>d</sup>*

still been completely resolved.

**3.3. Correlation dimension** 

falls into *Cl*(*Xj*), as follows.

where • is Euclidean norm.

40].

zero.

Some researchers have studied the choice of delay time *<sup>d</sup>*

τ

embedding dimension from a scalar time series[34-38].

*3.3.1. GP algorithm of correlation dimension* 

sphere with the radius *l* at the reference point *Xj*.

Then a correlation integral function is defined as

μ

θ

[ ( )]

τ

is too large, the dynamic shape of the attractor will be broken, which is

is the same for all the embedding dimension *d*[27-29]. The

although the issue of the delay time choice has

[*Cl*(*Xi*)] is the probability that *Xi* (*i*=1, ..., *n*)

(26)

small, the reconstructed attractor will be crowded around the main diagonal, which is called

called as irrelevance, and the phase space reconstruction is no longer representative of the true dynamics in the real system[28]. In normal circumstances, in order to make the

autocorrelation function method and mutual information technique[30] have been most

For the embedding dimension *d*, there are three methods that are usually used to choose the appropriate embedding dimension, including the correlation dimension, singular value decomposition(SVD), the false neighbors[21, 31, 32]. The correlation dimension method is to estimate appropriate dimension *d* in terms of the correlation theorem[8, 21, 33]. By increasing the embedding dimension, one notes an appropriate dimension d when the value of the correlation dimension stops changing. Broomhead and King[31] used the singular value decomposition (SVD) technique to determine an appropriate embedding dimension *d* directly from the raw time series. The false neighbor method is based on the fact that choosing a too low embedding dimension results in points, which are far apart in the original phase space, being moved closer together in the reconstruction space[32]. Besides, there are also some other methods and modified extensions developed based on the above methods. However, there are still problems on how to determine the appropriate

If a system is chaotic, the strange attractor in a region of the phase space constitutes an infinite hierarchy of self-similar structure, i.e. a fractal structure. One can use quantitative measures to define the fractal nature. The correlation dimension is a useful measurement. Grassberger and Proccacia give a kind of computation method, called GP algorithm[33, 39,

Let *X*1, *X2*, ..., *Xn* be a point of the attractor in phase space. *Cl*(*Xj*) is denoted as a hypersurface

μ

= − − = *n <sup>i</sup> <sup>l</sup> <sup>j</sup> <sup>X</sup> <sup>i</sup> <sup>X</sup> <sup>j</sup> <sup>l</sup> <sup>n</sup> <sup>C</sup> <sup>X</sup>* <sup>1</sup>

θ

( ) <sup>1</sup>

(*r*) is Heaviside function whose value is 1 if *r*≥0, otherwise,

τ

[26-30]. If the delay time *<sup>d</sup>*

τ

is too

**Theorem 2:** Let a map *<sup>n</sup> <sup>n</sup> G* : *R* → *R* . *A* is an attractor of map *G* that has only a finite number of periodic points of period *P*. Under a natural probability measure μ , the correlation dimension of *A* is ( ) *D*<sup>2</sup> μ . For a measure function *h* of *A*, *h R R* : *<sup>n</sup>* → , define a delay coordinate map *<sup>n</sup> <sup>d</sup> Fh* : *R* → *R* as

$$F\_h(X) = [h(X), h(G^{-1}(X)), \dots, h(G^{-(d-l)}(X))]\tag{29}$$

where *d P* ≥ . ( ) *hF* μ is a natural probability measure in *Rd* of ( ) *hF A* . If 2 *d D*≥ ( ) μ , then for almost every *h*, 2 2 ( ( )) ( ) *DF D <sup>h</sup>* μ = μ.

The theorem says that with the embedding dimension increasing, the slope of corresponding correlation integral curve will converge to the correlation dimension *D*<sup>2</sup> of the original system attractor. Therefore, the optimal embedding dimension can be estimated by using the correlation dimension *D*<sup>2</sup> . That is, if the embedding dimension *D*<sup>2</sup> *d* ≥ , the slope of the correlation integral curve is equal to the correlation dimension. This also indicates that the dimension estimation actually does not have to meet the requirements of the embedding theorem on the embedding dimension 2 1 *d* ≥ *D*<sup>2</sup> + . When the embedding dimension *D*<sup>2</sup> *d* ≥ , the reconstructed attractor can contain the fractal structural feature of the original system attractor to reflect the chaotic characteristics of the original system. Correlation dimension has been widely used in the analysis of chaotic time series.

#### *3.3.3. Chaotic test based on correlation dimension*

Chaos has a fractal structure so that the corresponding correlation dimension *D*2 is a fractional value. The estimation of correlation dimension *D*2 from a time series can be used to determine whether the time series is chaotic. If *D*<sup>2</sup> is fractional, the original time series can have chaotic features, otherwise, it cannot be chaotic. According to the correlation dimension theorem, when the embedding dimension *d* of the reconstructed phase space is increased to a certain value, the correlation dimension *D*2 will be saturated. Then, the optimal embedding dimension *d* will be given from a time series. The corresponding correlation dimension *D*2 is called as the correlation dimension of this time series.

Lorenz chaotic time series is given by the state variable *x* of Lorenz system as follows.

$$\begin{cases} \frac{d\mathbf{x}}{dt} = \sigma(\mathbf{y} - \mathbf{x})\\ \frac{d\mathbf{y}}{dt} = r\mathbf{x} - \mathbf{y} - \mathbf{x}\mathbf{z} \\ \frac{d\mathbf{z}}{dt} = -bz + \mathbf{x}\mathbf{y} \end{cases} \tag{30}$$

Nonlinear Analysis of Surface EMG Signals 137

is chosen as the sampling interval. With the

curves of data in Fig. a

the recontructed phase space, the delay time *<sup>d</sup>*

τ

integral curves will tend to be parallel and keep unchange in the range.

**Figure 10.** The correlation dimension analysis of surface EMG signal

dimension is over 6.

data and all data of *F* ∈ *F*

When the sampling interval

φ

increase of embedding dimension, the straight line segments of the computed correlation

a. The original data b. The correlation integral

The corresponding slope value is the correlation dimension of the surface EMG signal, about 3.8050± 0.0560(see Table 3). The result indicates that the surface EMG signal during movement has fractal feature to reveal the implied chaotic motion behavior. Table 3 shows that the reconstructed attractor contains the nature of the raw system when the embedding

*d* 2 3 4 5 6 7 8 9 10 11 12 *D*<sup>2</sup> 1.9315 2.7037 3.2808 3.4837 3.8047 3.7597 3.8511 3.8864 3.8129 3.7160 3.8039

The correlation dimension is a quantitative index that describes the fractal structure of chaotic attractor. It measures the freedom degree and complexity of the system. For the raw

test results will be obviously significant[42]. Here, the correlation dimension is used as a test statistic to analyze the surface EMG signal. According to the null hypothesis 3, 39 sets of surrogate data are produced in confidence level 95%. In order to quickly obtain the correlation dimension in a variety of circumstances, the linear parts of all correlation integral curves are taken as the same. Though these values are not accurate, this test algorithm is

analysis for Lorenz chaotic time series based on correlation dimension. There are significant differences between the original data and its surrogate data in *m* = 5 (see Fig. 11a). This result explains that the null hypothesis can be rejected in confidence level 95%. The

, in the case of the same embedding dimension, the corresponding

=1, Figure 11a and b show the results of the surrogate data test

**Table 3.** The analysis of correlation dimension of surface EMG signal during movement

*3.3.5. Study of surrogate data test method based on correlation dimension* 

great effective to test the chaotic fractal feature of the experimental data.

τ

where σ=10, *b*=8/3, γ=28, initial conditions: *x*(0)=5, *y*(0)=5, *z*(0)=15. The sampling intervalτ=0.1. The sampling points *N*=1000. For delay time τ*d*=τ, the corresponding correlation dimension values are given in Table 1 when the embedding dimension *d* is increased from 2 to 12. From this table, we can see that the correlation dimension of the time series is about 2.07. The result shows that the reconstructed attractor has a fractal structure to reflect the chaotic feature of the system. The time series can reconstruct the state space of the original system when the embedding dimension *d*=6.

Logistic chaotic time series *<sup>N</sup> <sup>n</sup> <sup>n</sup> x* <sup>1</sup> { } <sup>=</sup> is given by Logistic system in Eq.16 with α=3.9 and *e* =0. The length *N* is 1000 points. Here,τ*<sup>d</sup>*=1 (i.e. discrete time series interval). With increasing the embedding dimension *d*, the corresponding correlation dimension is 0.97 (see Table 2). The optimal embedding dimension is 2.

For finite sampling number (e.g. *N*=1000), the reconstructed attractor will be broken when the embedding dimension *d* is increased continuously to a higher value. The estimation of correlation dimension will fail during computation. Therefore, embedding dimension *d* should not be unlimitedly increased.


**Table 1.** The analysis of correlation dimensions of Lorenz chaos time series


**Table 2.** The analysis of correlation dimensions of Logistic chaos time series

#### *3.3.4. The analysis of surface EMG signal based on correlation dimension*

From the above analysis, we can see that the surface EMG signal has deterministic nonlinear component. Here, the correlation dimension is further used to study whether its nonlinear component are chaotic. Figure 10a shows a raw data for forearm pronation. Figure 10b gives the correlation integral curve of the data under the embedding dimension from 2 to 12. In the recontructed phase space, the delay time *<sup>d</sup>* τ is chosen as the sampling interval. With the increase of embedding dimension, the straight line segments of the computed correlation integral curves will tend to be parallel and keep unchange in the range.

**Figure 10.** The correlation dimension analysis of surface EMG signal

Computational Intelligence in Electromyography Analysis – 136 A Perspective on Current Applications and Future Challenges

where σ

=10, *b*=8/3,

γ

when the embedding dimension *d*=6.

Logistic chaotic time series *<sup>N</sup>*

The length *N* is 1000 points. Here,

optimal embedding dimension is 2.

should not be unlimitedly increased.

The sampling points *N*=1000. For delay time

optimal embedding dimension *d* will be given from a time series. The corresponding

= − +

*bz xy dt*

τ*d*=τ

values are given in Table 1 when the embedding dimension *d* is increased from 2 to 12. From this table, we can see that the correlation dimension of the time series is about 2.07. The result shows that the reconstructed attractor has a fractal structure to reflect the chaotic feature of the system. The time series can reconstruct the state space of the original system

embedding dimension *d*, the corresponding correlation dimension is 0.97 (see Table 2). The

For finite sampling number (e.g. *N*=1000), the reconstructed attractor will be broken when the embedding dimension *d* is increased continuously to a higher value. The estimation of correlation dimension will fail during computation. Therefore, embedding dimension *d*

*d* 2 3 4 5 6 7 8 9 10 11 12 *D2* 1.8009 1.9284 1.9718 2.0389 2.0737 2.0966 2.086 2.0788 2.0705 2.0760 2.0753

*d* 1 2 3 4 5 6 7 8 9 10 11 *D2* 0.9598 0.9718 0.9689 0.9713 0.9718 0.9591 0.9774 0.9621 0.9827 0.9826 0.9839

From the above analysis, we can see that the surface EMG signal has deterministic nonlinear component. Here, the correlation dimension is further used to study whether its nonlinear component are chaotic. Figure 10a shows a raw data for forearm pronation. Figure 10b gives the correlation integral curve of the data under the embedding dimension from 2 to 12. In

=28, initial conditions: *x*(0)=5, *y*(0)=5, *z*(0)=15. The sampling interval

*<sup>n</sup> <sup>n</sup> x* <sup>1</sup> { } <sup>=</sup> is given by Logistic system in Eq.16 with

= − −

*rx <sup>y</sup> xz dt*

(30)

, the corresponding correlation dimension

*<sup>d</sup>*=1 (i.e. discrete time series interval). With increasing the

α

τ=0.1.

=3.9 and *e* =0.

= −

*<sup>y</sup> <sup>x</sup> dt*

correlation dimension *D*2 is called as the correlation dimension of this time series.

 

*dz*

 

*dy*

*dx* σ( )

τ

**Table 1.** The analysis of correlation dimensions of Lorenz chaos time series

**Table 2.** The analysis of correlation dimensions of Logistic chaos time series

*3.3.4. The analysis of surface EMG signal based on correlation dimension* 

Lorenz chaotic time series is given by the state variable *x* of Lorenz system as follows.

The corresponding slope value is the correlation dimension of the surface EMG signal, about 3.8050± 0.0560(see Table 3). The result indicates that the surface EMG signal during movement has fractal feature to reveal the implied chaotic motion behavior. Table 3 shows that the reconstructed attractor contains the nature of the raw system when the embedding dimension is over 6.


**Table 3.** The analysis of correlation dimension of surface EMG signal during movement

#### *3.3.5. Study of surrogate data test method based on correlation dimension*

The correlation dimension is a quantitative index that describes the fractal structure of chaotic attractor. It measures the freedom degree and complexity of the system. For the raw data and all data of *F* ∈ *F*φ , in the case of the same embedding dimension, the corresponding test results will be obviously significant[42]. Here, the correlation dimension is used as a test statistic to analyze the surface EMG signal. According to the null hypothesis 3, 39 sets of surrogate data are produced in confidence level 95%. In order to quickly obtain the correlation dimension in a variety of circumstances, the linear parts of all correlation integral curves are taken as the same. Though these values are not accurate, this test algorithm is great effective to test the chaotic fractal feature of the experimental data.

When the sampling interval τ=1, Figure 11a and b show the results of the surrogate data test analysis for Lorenz chaotic time series based on correlation dimension. There are significant differences between the original data and its surrogate data in *m* = 5 (see Fig. 11a). This result explains that the null hypothesis can be rejected in confidence level 95%. The

a. \* is *D*2 value of raw data in *m* = 5, τ = 1 histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

c. \* is *D*2 value of raw data in *m* = 4, τ=0.1 histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

e. \* is *D*2 value of raw data in *m* = 2, τ=0.005 histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

Nonlinear Analysis of Surface EMG Signals 139

=0.1. The differences between the raw data and its

=0.005. Figure 11e shows the surrogate

b. The curves of surface EMG and surrogate data with m, abscissa is m=2~8, ordinate is *D*<sup>2</sup>

differences disappear between the raw data and its surrogate data in *m* = 10 (see Fig. 11b). This illustrates that the reconstructed attractor appears broken. The reconstructed phase space is similar to that of the surrogate data with linear stochastic noise characteristics.

surrogate data can be seen in Fig. 11c (*m* = 4). When *m*>2, these differences become larger as

data test histogram in *m* = 2. The correlation dimension curves of the raw data and its surrogate data are given in *m* = 2 ~ 10 (see Fig. 11f). Even in the case of oversampling, the correlation dimension as test statistic can also make the surrogate data method very

Figure 12 shows the surrogate data analysis for the surface EMG signal in Fig. 12a based on correlation dimension. When *m* = 6, the correlation dimension value of the raw data is different from those of its surrogate data generated by the null hypothesis 3. The correlation dimension curves of the raw data and its surrogate data are given when *m* = 2 ~ 8 in Fig. 12b. We can see the differences between the original data and its surrogate data. The null hypothesis 3 can be rejected in confidence level 95%. The result indicates that the surface

**Figure 12.** The surrogate data test analysis based on correlation dimension for surface EMG signal.

The Lyapunov exponent method is to directly identify whether a system is chaotic. If the system is chaotic, the Lyapunov exponent is positive. Otherwise, the Lyapunov exponent is negative. For this, the Lyapunov exponent can be used to test the chaotic feature of a signal under study. The first algorithms developed computed the whole Lyapunov spectrum by

τ

τ

*3.3.6. The surface EMG signal analysis based on the surrogate data and correlation* 

EMG signal has deterministic nonlinear components, even chaotic.

a. \* is *D*2 of surface EMG in *m* = 6, histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

**3.4. Largest Lyapunov exponent** 

Figure 11c and d show the results with

effective.

*dimension* 

*m* increases. Figure 11e and f show the results with

b. The *D*2 curves of raw data and surrogate data with m, abscissa is m=2~10, ordinate is *D*<sup>2</sup>

d. The *D*2 curves of raw data and surrogate data with m, abscissa is m=2~10, ordinate is *D*<sup>2</sup>

f. The *D*2 curves of raw data and surrogate data with m, abscissa is m=2~10, ordinate is *D*<sup>2</sup>

**Figure 11.** The surrogate data test analysis based on correlation dimension for Lorenz chaos time series by sampling intervals

differences disappear between the raw data and its surrogate data in *m* = 10 (see Fig. 11b). This illustrates that the reconstructed attractor appears broken. The reconstructed phase space is similar to that of the surrogate data with linear stochastic noise characteristics.

Figure 11c and d show the results with τ=0.1. The differences between the raw data and its surrogate data can be seen in Fig. 11c (*m* = 4). When *m*>2, these differences become larger as *m* increases. Figure 11e and f show the results with τ=0.005. Figure 11e shows the surrogate data test histogram in *m* = 2. The correlation dimension curves of the raw data and its surrogate data are given in *m* = 2 ~ 10 (see Fig. 11f). Even in the case of oversampling, the correlation dimension as test statistic can also make the surrogate data method very effective.

#### *3.3.6. The surface EMG signal analysis based on the surrogate data and correlation dimension*

Figure 12 shows the surrogate data analysis for the surface EMG signal in Fig. 12a based on correlation dimension. When *m* = 6, the correlation dimension value of the raw data is different from those of its surrogate data generated by the null hypothesis 3. The correlation dimension curves of the raw data and its surrogate data are given when *m* = 2 ~ 8 in Fig. 12b. We can see the differences between the original data and its surrogate data. The null hypothesis 3 can be rejected in confidence level 95%. The result indicates that the surface EMG signal has deterministic nonlinear components, even chaotic.

**Figure 12.** The surrogate data test analysis based on correlation dimension for surface EMG signal.

#### **3.4. Largest Lyapunov exponent**

Computational Intelligence in Electromyography Analysis – 138 A Perspective on Current Applications and Future Challenges

a. \* is *D*2 value of raw data in *m* = 5,

c. \* is *D*2 value of raw data in *m* = 4,

e. \* is *D*2 value of raw data in *m* = 2,

by sampling intervals

histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

histogram is *D*2 distribution of surrogate data, abscissa is *D*2, ordinate is histogram

τ= 1

τ=0.1

τ=0.005

**Figure 11.** The surrogate data test analysis based on correlation dimension for Lorenz chaos time series

b. The *D*2 curves of raw data and surrogate data with m, abscissa is m=2~10, ordinate is *D*<sup>2</sup>

d. The *D*2 curves of raw data and surrogate data with m, abscissa is m=2~10, ordinate is *D*<sup>2</sup>

f. The *D*2 curves of raw data and surrogate data with m, abscissa is m=2~10, ordinate is *D*<sup>2</sup>

The Lyapunov exponent method is to directly identify whether a system is chaotic. If the system is chaotic, the Lyapunov exponent is positive. Otherwise, the Lyapunov exponent is negative. For this, the Lyapunov exponent can be used to test the chaotic feature of a signal under study. The first algorithms developed computed the whole Lyapunov spectrum by Wolf et al. [43] and Sano et al. [44]. Meanwhile, the largest Lyapunov exponent is sufficient for assessing the presence of chaos. At present, there are many algorithms to estimate the largest Lyapunov exponent from a time series, such as an algorithm given by Rosenstein et al.[45]. This algorithm is aimed specifically at estimating the largest Lyapunov exponent from short data.

#### *3.4.1. Algorithm of largest Lyapunov exponent estimation[45]*

For a short time series, Rosenstein et al. present a robust estimation algorithm of the largest Lyapunov exponent. First, the attractor is reconstructed, refer to Eq. 25. Next, the algorithm locates the closest neighbor of each point *Xi* on the trajectory, with respect to the Euclidian distance. Then, one defines the distance between two neighboring points at instant *n*=0 by:

$$d\_i(0) = \min\_{X\_{\nearrow}} \left\| X\_{\nearrow} - X\_{\nearrow} \right\| \tag{31}$$

Nonlinear Analysis of Surface EMG Signals 141

b. *y*(*n*) curves of time series x of Logistic

α

dimension m=2~8

d. *y*(*n*) curves of gaussian noise embedding dimension m=2~8

f. *y*(*n*) curves of Lorenz chaos series in

= 1 , embedding dimension m=3~8

τ

=3.9, embedding

chaos system in

a. *y*(*n*) curves of time series x of Lorenz

embedding dimension m=2~8

c. *y*(*n*) curves of Henon chaos series embedding dimension m=2~8

e. *y*(*n*) curves of linear stochastic process embedding dimension m=2~8

**Figure 13.** *y*(*n*) curves of signals, abscissa is *n*, ordinate is *y*(*n*)

τ=0.01,

chaos system sampled by

where • is the Euclidian norm. Here, the temporal separation of the nearest neighbors should be greater than the mean period of the time series.

$$|i - j| > mean \quad period \tag{32}$$

According to time, the average distance between two neighboring vectors can be simply

$$d\_i(n) = \left\| X\_{\,\, :\_{i \leftrightarrow n}} - X\_{\,\, :\_{i \leftrightarrow n}} \right\|\tag{33}$$

Assume that the system is controlled by the largest Lyapunov exponent only. Then, the distance between two neighbor points obey the following relationship:

$$d(t) = Ce^{\lambda t} \tag{34}$$

For *t* = *n*Δ*t* , there is:

$$d\_{\cdot}(n) \simeq C\_{\cdot} e^{\lambda n \Delta t} \tag{35}$$

$$\mathcal{A} = \frac{1}{n \cdot \Delta t} \ln \frac{d\_i(n)}{C\_i} \tag{36}$$

$$
\ln \cdot \mathcal{X} = \frac{1}{\Delta t} (\ln d\_i(n) - \ln C\_i) \tag{37}
$$

$$\begin{split} n \{ \dot{\mathcal{A}} \} &= \frac{1}{\Delta t} \left\langle \ln d\_{\cdot}(n) \right\rangle - \frac{1}{\Delta t} \left\langle \ln C\_{\cdot} \right\rangle \\ &= \frac{1}{\Delta t} \left\langle \ln d\_{\cdot}(n) \right\rangle - b \end{split} \tag{38}$$

where *Ci <sup>t</sup> <sup>b</sup>* ln <sup>1</sup> <sup>Δ</sup> <sup>=</sup> .

**Figure 13.** *y*(*n*) curves of signals, abscissa is *n*, ordinate is *y*(*n*)

Computational Intelligence in Electromyography Analysis – 140 A Perspective on Current Applications and Future Challenges

*3.4.1. Algorithm of largest Lyapunov exponent estimation[45]* 

should be greater than the mean period of the time series.

from short data.

For *t* = *n*Δ*t* , there is:

where *Ci <sup>t</sup> <sup>b</sup>* ln <sup>1</sup> <sup>Δ</sup> <sup>=</sup> .

Wolf et al. [43] and Sano et al. [44]. Meanwhile, the largest Lyapunov exponent is sufficient for assessing the presence of chaos. At present, there are many algorithms to estimate the largest Lyapunov exponent from a time series, such as an algorithm given by Rosenstein et al.[45]. This algorithm is aimed specifically at estimating the largest Lyapunov exponent

For a short time series, Rosenstein et al. present a robust estimation algorithm of the largest Lyapunov exponent. First, the attractor is reconstructed, refer to Eq. 25. Next, the algorithm locates the closest neighbor of each point *Xi* on the trajectory, with respect to the Euclidian distance. Then, one defines the distance between two neighboring points at instant *n*=0 by:

*<sup>j</sup> <sup>i</sup> <sup>X</sup> di <sup>X</sup> <sup>X</sup> <sup>j</sup>*

where • is the Euclidian norm. Here, the temporal separation of the nearest neighbors

According to time, the average distance between two neighboring vectors can be simply

Assume that the system is controlled by the largest Lyapunov exponent only. Then, the

*<sup>t</sup> d t Ce* λ

*<sup>i</sup> <sup>i</sup> d n C e* <sup>Δ</sup> ≈

*n t*

λ

λ

λ

*t*

*<sup>d</sup> <sup>n</sup> <sup>t</sup> <sup>n</sup>*

⋅ Δ

(ln ( ) ln ) <sup>1</sup> *<sup>i</sup> Ci <sup>d</sup> <sup>n</sup> <sup>t</sup> <sup>n</sup>* <sup>−</sup> <sup>Δ</sup> <sup>⋅</sup>

*d n b*

ln <sup>1</sup> ln ( ) <sup>1</sup>

*i*

<sup>Δ</sup> <sup>−</sup> <sup>Δ</sup> <sup>=</sup> ln ( ) <sup>1</sup>

<sup>−</sup> <sup>Δ</sup> <sup>=</sup>

*n t*

*i i C d n*

*t*

*i i*

λ

( ) ln <sup>1</sup>

distance between two neighbor points obey the following relationship:

(0) = min − (31)

*i* − *j* > *mean period* (32)

*<sup>i</sup> X <sup>j</sup> <sup>n</sup> Xi <sup>n</sup> d n* = <sup>+</sup> − <sup>+</sup> ( ) (33)

( ) = (34)

( ) (35)

= (36)

= (37)

(38)

*C*

Then, the Lyapunov exponent can be given by using a least-squares fit to the "average" line:

$$\text{hyp}(n) = \frac{1}{\Delta t} \left\langle \ln d\_{\cdot}(n) \right\rangle \tag{39}$$

Nonlinear Analysis of Surface EMG Signals 143

(40)

**3.5. Principal component analysis** 

Broomhead and King[31] proposed the idea of singular system analysis that determines an appropriate embedding dimension *d* directly from the raw time series. It provides its convenience for the further analysis of the given system. Numerical experience, however, led several authors to express some doubts about reliability of singular system analysis in the attractor reconstruction[46-48]. Palus and Dvorak[37] explain why singular-value decomposition(SVD), the heart of the singular system analysis and by nature a linear method, may become misleading technique when it is used in nonlinear dynamics studies that reconstruction parameters are time-delay, embedding dimension (or embedding windows). For this, we propose a novel nonlinear analysis method based symplectic

Let a time series *x*1, *x*2,..., *xn* be the measured signal by sampling interval *ts*, *n* is the number of samples. According to Takens' embedding theorem, a trajectory matrix *X* can be given by

> τ*<sup>d</sup>*=1):

1 2

where *d* is embedding dimension. *m*=*n*-*d*+1 is the number of points in *d*-dimension

there are a *m*×*d* orthogonal matrix *V* and a *d*×*d* orthogonal matrix. The matrix *X* can be

*S m i j d ij ij <sup>i</sup>* =

*<sup>T</sup>* (*V* ⋅*V* ) =

In order to facilitate the calculation, Broomhead et al. applies the covariance matrix *C* of the

*ij <sup>T</sup>* (*U* ⋅*U* ) = (*U* ⋅*U* ) =

*ij ij*

δ

*ij ij T*

δ

+

*m m n*

 

*x x x*

1

*x x x x x x*

2 3 1

 

=

 

*T m*

*T T*

2 1

 

*X*

δ σ

*X X*

=

*X*

where *S* is *d*×*d* diagonal matrix, whose elements are defined

matrix *X* to replace the matrix *X*. The details are as follows:

 

+

*<sup>i</sup>* , = 1,, , denotes a point in the attractor. For the matrix *X*,

*<sup>T</sup> X* =*VSU* (41)

, , =1,2,, (42)

(43)

(44)

*d d*

geometry, called symplectic principal component analysis(SPCA)[49].

*3.5.1. Principle and algorithm of principal component analysis* 

time delay coordinates method, refer to Eq. 25(

reconstruction attractor, *X i m <sup>T</sup>*

Since the matrix *V* is orthogonal, then

Meanwhile, for the matrix *U*, there are

decomposed as follows.

where *y*(*n*) = *n* λ + *b* . The largest Lyapunov exponent λ is the slope of *y*(*n*) in the above equation.

This method is deduced directly from the largest Lyapunov exponent definition. The accurate evaluation of λ depends on the full use of the data. In practice, the curve *y*(*n*) will tend to saturation. The largest Lyapunov exponent λ is given by computing the slope of the linear part in the curve *y*(*n*).

#### *3.4.2. Chaos test based on largest Lyapunov exponent*

In general, if the signal is chaotic, the slope of the curve *y*(*n*) will be independent of the embedding dimension. Otherwise, if the signal is not chaotic, the slope of the curve *y*(*n*) will depend on the embedding dimension. When the embedding dimension m is chosen from 2 to 8, the Lyapunov exponent of the curve *y*(*n*) of the signal is shown in Figure 13. For a chaotic signal, a good illustration is given (see Figure 13a, b and c). The *y*(*n*) curves are different from those of a non-chaotic signal (compare with Figure 13d and e). However, even for chaotic signals, the *y*(*n*) curves are not always parallel. For example, in the case of undersampling (τ = 1 ), the *y*(*n*) curves of Lorenz chaotic time series are similar to those of linear stochastic process(compare with Figure 13e and f). In the literature[23], the *y*(*n*) curves of the Ikeda chaotic time series are also not parallel.

#### *3.4.3. The analysis of surface EMG signal based on the largest Lyapunov exponent*

Figure 14 gives the curves of Lyapunove exponent *y*(*n*) for the surface EMG signal. The *y*(*n*) curves are not very parallel for the surface EMG signal. It is difficult to distinguish the curves of *y*(*n*) for the surface EMG signal from those of Figure 13d, f and g. The surface EMG signal can not be determined as chaotic, or as stochastic. But it can be a highdimensional system.

**Figure 14.** *y*(*n*) curves of surface EMG signal, embedding dimension *m*=2~8, abscissa is *n*, ordinate is *y*(*n*)

#### **3.5. Principal component analysis**

Computational Intelligence in Electromyography Analysis – 142 A Perspective on Current Applications and Future Challenges

λ

tend to saturation. The largest Lyapunov exponent

*3.4.2. Chaos test based on largest Lyapunov exponent* 

of the Ikeda chaotic time series are also not parallel.

where *y*(*n*) = *n*

undersampling (

dimensional system.

*y*(*n*)

τ

accurate evaluation of

linear part in the curve *y*(*n*).

equation.

λ

Then, the Lyapunov exponent can be given by using a least-squares fit to the "average" line:

ln ( ) <sup>1</sup> ( ) *<sup>d</sup> <sup>n</sup> <sup>t</sup>*

This method is deduced directly from the largest Lyapunov exponent definition. The

In general, if the signal is chaotic, the slope of the curve *y*(*n*) will be independent of the embedding dimension. Otherwise, if the signal is not chaotic, the slope of the curve *y*(*n*) will depend on the embedding dimension. When the embedding dimension m is chosen from 2 to 8, the Lyapunov exponent of the curve *y*(*n*) of the signal is shown in Figure 13. For a chaotic signal, a good illustration is given (see Figure 13a, b and c). The *y*(*n*) curves are different from those of a non-chaotic signal (compare with Figure 13d and e). However, even for chaotic signals, the *y*(*n*) curves are not always parallel. For example, in the case of

linear stochastic process(compare with Figure 13e and f). In the literature[23], the *y*(*n*) curves

Figure 14 gives the curves of Lyapunove exponent *y*(*n*) for the surface EMG signal. The *y*(*n*) curves are not very parallel for the surface EMG signal. It is difficult to distinguish the curves of *y*(*n*) for the surface EMG signal from those of Figure 13d, f and g. The surface EMG signal can not be determined as chaotic, or as stochastic. But it can be a high-

**Figure 14.** *y*(*n*) curves of surface EMG signal, embedding dimension *m*=2~8, abscissa is *n*, ordinate is

*3.4.3. The analysis of surface EMG signal based on the largest Lyapunov exponent* 

+ *b* . The largest Lyapunov exponent

*<sup>y</sup> <sup>n</sup> <sup>i</sup>* <sup>Δ</sup> <sup>=</sup> (39)

is the slope of *y*(*n*) in the above

is given by computing the slope of the

λ

depends on the full use of the data. In practice, the curve *y*(*n*) will

λ

= 1 ), the *y*(*n*) curves of Lorenz chaotic time series are similar to those of

Broomhead and King[31] proposed the idea of singular system analysis that determines an appropriate embedding dimension *d* directly from the raw time series. It provides its convenience for the further analysis of the given system. Numerical experience, however, led several authors to express some doubts about reliability of singular system analysis in the attractor reconstruction[46-48]. Palus and Dvorak[37] explain why singular-value decomposition(SVD), the heart of the singular system analysis and by nature a linear method, may become misleading technique when it is used in nonlinear dynamics studies that reconstruction parameters are time-delay, embedding dimension (or embedding windows). For this, we propose a novel nonlinear analysis method based symplectic geometry, called symplectic principal component analysis(SPCA)[49].

#### *3.5.1. Principle and algorithm of principal component analysis*

Let a time series *x*1, *x*2,..., *xn* be the measured signal by sampling interval *ts*, *n* is the number of samples. According to Takens' embedding theorem, a trajectory matrix *X* can be given by time delay coordinates method, refer to Eq. 25(τ*<sup>d</sup>*=1):

$$X = \begin{bmatrix} X\_1^T \\ X\_2^T \\ \vdots \\ X\_n^T \end{bmatrix} = \begin{bmatrix} \mathbf{x}\_1 & \mathbf{x}\_2 & \cdots & \mathbf{x}\_d \\ \mathbf{x}\_2 & \mathbf{x}\_3 & \cdots & \mathbf{x}\_{d+1} \\ \vdots & \vdots & \cdots & \vdots \\ \mathbf{x}\_n & \mathbf{x}\_{n+1} & \cdots & \mathbf{x}\_n \end{bmatrix} \tag{40}$$

where *d* is embedding dimension. *m*=*n*-*d*+1 is the number of points in *d*-dimension reconstruction attractor, *X i m <sup>T</sup> <sup>i</sup>* , = 1,, , denotes a point in the attractor. For the matrix *X*, there are a *m*×*d* orthogonal matrix *V* and a *d*×*d* orthogonal matrix. The matrix *X* can be decomposed as follows.

$$X = VSU^{\top} \tag{41}$$

where *S* is *d*×*d* diagonal matrix, whose elements are defined

$$S\_{\circ} = \mathcal{S}\_{\circ} \sqrt{\sigma\_i m} \,, \qquad i, j = 1, 2, \dots, d \tag{42}$$

Since the matrix *V* is orthogonal, then

$$(\boldsymbol{V}^{\boldsymbol{T}} \cdot \boldsymbol{V})\_{\boldsymbol{\psi}} = \boldsymbol{\mathcal{S}}\_{\boldsymbol{\psi}} \tag{43}$$

Meanwhile, for the matrix *U*, there are

$$(\boldsymbol{U}^{\boldsymbol{T}} \cdot \boldsymbol{U})\_{\boldsymbol{\bar{\boldsymbol{\eta}}}} = (\boldsymbol{U} \cdot \boldsymbol{U}^{\boldsymbol{T}})\_{\boldsymbol{\bar{\boldsymbol{\eta}}}} = \boldsymbol{\delta}\_{\boldsymbol{\bar{\boldsymbol{\eta}}}} \tag{44}$$

In order to facilitate the calculation, Broomhead et al. applies the covariance matrix *C* of the matrix *X* to replace the matrix *X*. The details are as follows:

Computational Intelligence in Electromyography Analysis – 144 A Perspective on Current Applications and Future Challenges

$$C = \frac{1}{m} \cdot X^T \cdot X \tag{45}$$

Nonlinear Analysis of Surface EMG Signals 145

b. The principal components of Logistic

d. The principal components of Logistic

σ2=0.82

with exterior noise

σ2=0.0012

with interior noise

σ2=0.0012

*3.5.2. Influence of noise on the principal component spectrum of chaotic time series* 

principal component spectrum of Logistic attractor with the measurement noise

a. The principal components of Logistic with no noise

c. The principal components of Logistic with

SVD, d=3 : 2 : 23, abscissa is *d*, ordinate is log( / ( )) *<sup>i</sup> <sup>i</sup>*

σ2=0.0012

**Figure 15.** The principal component analysis of Logistic chaos series with different noises based on

σ *tr* σ

exterior noise

and σ

Figure 15a shows the principal component spectrum of Logistic attractor from a Logistic chaotic time series without noise. The principal component spectrum has not a significant noise level. When the interior noise is Gaussian noise with zero mean and 0.0012 variance, the principal component spectrum is given for Logistic attractor. Figure 15c and d give the

2=0.82, respectively. It can be seen that the Logistic attractor with the internal noise has the same principal component spectrum as the attractor without noise. The curves of the principal component spectrum are also slanting. The total energy is significantly distributed into each principal axes. The principal components are declining with the index *i* so that there is no noise floor. It is difficult to choose an appropriate embedding dimension *d*. For the larger measurement noise, the corresponding principal component spectrum of Logistic attractor slant into a floor area with increasing the embedding dimension. In the floor area, the principal components keep unchanged and do not decline with the index *i*, called noise floor. Broomhead and King[31] have suggested that this noise floor can be used to determine the embedding dimension and filter out noise from the data. The signal energy will be focused on the truncated principal components and the corresponding principal axes

$$C\_{\dot{y}} = \frac{1}{m} \sum\_{k=1}^{m} \mathbf{x}\_{i+k-1} \mathbf{x}\_{j+k-1} \,, \qquad i, j = 1, 2, \cdots, d \tag{46}$$

Its values reflect the degree of correlation between the time delay coordinate variable *i* and *j*.

$$C = \frac{1}{m} \cdot X^{\top} \cdot X = \frac{1}{m} \cdot U \cdot S \cdot S \cdot U^{\top} = \frac{1}{m} \cdot U \cdot S^2 \cdot U^{\top} \tag{47}$$

$$\mathbf{U}^{\top} \cdot \mathbf{C} \cdot \mathbf{U} = \frac{1}{m} \cdot (XU)^{\top} \cdot XU = \frac{1}{m} \cdot S^{\top} \tag{48}$$

Let *Y*=*XU*, then:

$$U^{\top} \cdot C \cdot U = \frac{1}{m} \cdot Y^{\top} \cdot Y = \frac{1}{m} \cdot S^{\bot} \tag{49}$$

where *UTCU* is the covariance matrix of the matrix *Y*. Its elements are zero, except that the diagonal elements are equal to σ*<sup>i</sup>*. This means that the variables *i* and *j* of the matrix *Y* are independent. The coordinate system is orthogonal, which is constituted from the variables of the matrix *Y* after the above transformed. The σ*<sup>i</sup>* is called the principal component or singular value in accordance with the order of the largest to the smallest. The orthogonal vector *Ui* corresponding to the principal component σ*<sup>i</sup>* is called the principal axis. The principal component describes the distribution of the signal energy. That is, the value of the principal component reflects the projection of the signal energy in the corresponding principal axis. In the different principal axes, a distribution value is given as <sup>=</sup> *d i i* 1 *i* σ / σ , where <sup>=</sup> *d i* 1σ *<sup>i</sup>* is the total energy of the signal. If σ *<sup>i</sup>*+<sup>1</sup> ≈ ≈ σ *<sup>d</sup>* , the distribution values are called a noise floor. The distribution can be used to estimate the dimension of the dynamical system that generates a time series or to filter out the noise. Let *U <sup>i</sup>* be the principal axes corresponding to the principal components over the noise floor. Zero vector describes the principal axes corresponding to the noise floor. Thus, for σ*<sup>i</sup>* > noise floor, a new coordinate transform matrix is made up of:

$$\overline{U} = [\overline{U}\_1, \overline{U}\_2, \dots, \overline{U}\_i, 0, \dots, 0] \tag{50}$$

In order to filter out noise, the trajectory matrix *X* is first projected into the coordinate system *U*.

$$Y = XU\tag{51}$$

The variables in the matrix *Y* are independent. Then, the original coordinate system is updated by using the matrix *Y*:

$$
\overline{X} = Y\overline{U}^{\overline{r}} = (XU) \cdot \overline{U}^{\overline{r}} \tag{52}
$$

That is, *X* is a new time delay coordinate system.

#### *3.5.2. Influence of noise on the principal component spectrum of chaotic time series*

Figure 15a shows the principal component spectrum of Logistic attractor from a Logistic chaotic time series without noise. The principal component spectrum has not a significant noise level. When the interior noise is Gaussian noise with zero mean and 0.0012 variance, the principal component spectrum is given for Logistic attractor. Figure 15c and d give the principal component spectrum of Logistic attractor with the measurement noise σ2=0.0012 and σ2=0.82, respectively. It can be seen that the Logistic attractor with the internal noise has the same principal component spectrum as the attractor without noise. The curves of the principal component spectrum are also slanting. The total energy is significantly distributed into each principal axes. The principal components are declining with the index *i* so that there is no noise floor. It is difficult to choose an appropriate embedding dimension *d*. For the larger measurement noise, the corresponding principal component spectrum of Logistic attractor slant into a floor area with increasing the embedding dimension. In the floor area, the principal components keep unchanged and do not decline with the index *i*, called noise floor. Broomhead and King[31] have suggested that this noise floor can be used to determine the embedding dimension and filter out noise from the data. The signal energy will be focused on the truncated principal components and the corresponding principal axes

Computational Intelligence in Electromyography Analysis – 144 A Perspective on Current Applications and Future Challenges

Let *Y*=*XU*, then:

where <sup>=</sup> *d i* 1σ

system *U*.

diagonal elements are equal to

transform matrix is made up of:

updated by using the matrix *Y*:

That is, *X* is a new time delay coordinate system.

= = <sup>=</sup> <sup>+</sup> <sup>−</sup> <sup>+</sup> <sup>−</sup>

Its values reflect the degree of correlation between the time delay coordinate variable *i* and *j*.

*<sup>T</sup> <sup>T</sup> <sup>T</sup> <sup>U</sup> <sup>S</sup> <sup>U</sup> <sup>m</sup> <sup>U</sup> <sup>S</sup> <sup>S</sup> <sup>U</sup> <sup>m</sup> <sup>X</sup> <sup>X</sup> <sup>m</sup> <sup>C</sup>* <sup>=</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>=</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>=</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>1</sup> <sup>1</sup> <sup>1</sup> <sup>2</sup>

where *UTCU* is the covariance matrix of the matrix *Y*. Its elements are zero, except that the

independent. The coordinate system is orthogonal, which is constituted from the variables

singular value in accordance with the order of the largest to the smallest. The orthogonal

principal component describes the distribution of the signal energy. That is, the value of the principal component reflects the projection of the signal energy in the corresponding principal axis. In the different principal axes, a distribution value is given as <sup>=</sup>

called a noise floor. The distribution can be used to estimate the dimension of the dynamical system that generates a time series or to filter out the noise. Let *U <sup>i</sup>* be the principal axes corresponding to the principal components over the noise floor. Zero vector describes the

In order to filter out noise, the trajectory matrix *X* is first projected into the coordinate

The variables in the matrix *Y* are independent. Then, the original coordinate system is

σ

σ*<sup>i</sup>*+<sup>1</sup> ≈ ≈

σ

*m*

σ

of the matrix *Y* after the above transformed. The

vector *Ui* corresponding to the principal component

*<sup>i</sup>* is the total energy of the signal. If

principal axes corresponding to the noise floor. Thus, for

*<sup>X</sup> <sup>X</sup> <sup>m</sup> <sup>C</sup> <sup>T</sup>* <sup>=</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>1</sup> (45)

(47)

*<sup>k</sup> ij <sup>i</sup> <sup>k</sup> <sup>j</sup> <sup>k</sup> <sup>x</sup> <sup>x</sup> <sup>i</sup> <sup>j</sup> <sup>d</sup> <sup>m</sup> <sup>C</sup>* <sup>1</sup> <sup>1</sup> <sup>1</sup> , , 1, 2, , <sup>1</sup> (46)

<sup>1</sup> <sup>2</sup> ( ) <sup>1</sup> *<sup>S</sup> <sup>m</sup> XU XU <sup>m</sup> <sup>U</sup> <sup>C</sup> <sup>U</sup> <sup>T</sup> <sup>T</sup>* <sup>⋅</sup> <sup>⋅</sup> <sup>=</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>=</sup> <sup>⋅</sup> (48)

<sup>1</sup> <sup>1</sup> <sup>2</sup> *<sup>S</sup> <sup>m</sup> <sup>Y</sup> <sup>Y</sup> <sup>m</sup> <sup>U</sup> <sup>C</sup> <sup>U</sup> <sup>T</sup> <sup>T</sup>* <sup>⋅</sup> <sup>⋅</sup> <sup>=</sup> <sup>⋅</sup> <sup>⋅</sup> <sup>=</sup> <sup>⋅</sup> (49)

*<sup>i</sup>*. This means that the variables *i* and *j* of the matrix *Y* are

σ

σ

[ , , , ,0, ,0] *U* = *U*<sup>1</sup> *U* <sup>2</sup> *Ui* (50)

*Y* = *XU* (51)

*<sup>T</sup> <sup>T</sup> X* = *YU* = (*XU* ) ⋅*U* (52)

*<sup>i</sup>* is called the principal component or

*<sup>i</sup>* is called the principal axis. The

*<sup>d</sup>* , the distribution values are

*<sup>i</sup>* > noise floor, a new coordinate

*d i i* 1 *i* σ/

σ,

a. The principal components of Logistic with no noise

exterior noise σ2=0.0012

b. The principal components of Logistic with interior noise σ2=0.0012

d. The principal components of Logistic with exterior noise σ2=0.82

**Figure 15.** The principal component analysis of Logistic chaos series with different noises based on SVD, d=3 : 2 : 23, abscissa is *d*, ordinate is log( / ( )) *<sup>i</sup> <sup>i</sup>* σ *tr* σ

when the principal components above noise floor are only held. The number of the principal components above the noise floor is the optimal embedding dimension.

Nonlinear Analysis of Surface EMG Signals 147

=1, each line is separated from each other and tends to

to those in the Figure 16a. When

noise in chaotic data[59, 60].

product:

**3.6. Symplectic principal component analysis** 

**3.7. Symplectic principal component method** 

τ

horizontal line in the case of different embedding dimensions(see Fig. 16c). It shows that the distribution of the total energy has little difference in each principal axis, like the Gaussian noise (see Fig. 16d). For the Gaussian noise, its principal component spectrum curves are horizontal lines, where *N*=10000. It shows that every principal component is equal to each other. The energy distributes into every principal axis averagely. Therefore, it can be seen that sampling interval affects the determination of embedding dimension. When the sampling interval is not undersampling, the determination of embedding dimension depends the amount of signal-noise-ratio. In the case of undersampling, the chaotic time

series is similar to noise so that the embedding dimension seems to be estimated as 1.

The symplectic geometry is a kind of phase space geometry. Its nature is nonlinear. It can describe the system structure, especially nonlinear structure, very well. It has been used to study various nonlinear dynamical systems[50-52] since Feng Kang[53] has proposed a symplectic algorithm for solving symplectic differential. However, from the view of data analysis, few literatures have employed symplectic geometry theory to explore the dynamics of the system. Our previous works have proposed the estimation of the embedding dimension based on symplectic geometry from a time series[49, 54-56]. Subsequently, Niu et al. have used our method to evaluate sprinter's surface EMG signals[57]. Xie et al[58] have proposed a kind of symplectic geometry spectra based on our work. Subsequently, we show that SPCA can well represent chaotic time series and reduce

In SPCA, a fundamental step is to build the multidimensional structure (attractor) in symplectic geometry space. Here, in terms of Taken's embedding theorem, we first construct an attractor in phase space, i.e. the trajectory matrix *X* from a time series. That is, for a measured data (the observable of the system under study) *x*1, *x*2, ..., *xn* recorded with sampling interval *ts*, the corresponding *d*-dimension reconstruction attractor, *Xm*×*<sup>d</sup>* can be given (refer to Eq.40).Then we describe the symplectic principal component analysis (SPCA)

SPCA is a kind of PCA approaches based on symplectic geometry. Its idea is to map the investigated complex system in symplectic space and elucidate the dominant features underlying the measured data. The first few larger components capture the main relationship between the variables in symplectic space. The remaining components are composed of the less important components or noise in the measured data. In symplectic space, the used geometry is called symplectic geometry. Different from Eulid geometry, symplectic geometry is the even dimensional geometry with a special symplectic structure. It is dependent on a bilinear antisymmetric nonsingular cross product——symplectic cross

based on symplectic geometry theory and give its corresponding algorithm.

Besides, the new coordinate system corresponding to the principal axes can eliminate the noise floor to reduce the noise from the data. However, the truncated position of the principal components depends on the signal-noise-ratio, especially for the measurement noise. The principal components of the chaotic time series based on SVD spectrum more easily subject to the measurement noise so that the embedding dimension estimation is directly affected. For the smaller noise, there is the more number of principal components above the noise floor. For the larger noise, the number of the corresponding principal components will be reduced. Here, the above calculation accuracy is 2.2204e-016, which does not consider the numerical calculation error.
