**2. Detecting nonlinearity of the surface EMG signals**

In many areas of science and engineering, it is a critical issue to determine whether an observed time series is purely stochastic, or deterministic nonlinear, even chaotic. One may know about the intrinsic properties of the observed phenomenon by distinguishing between nonlinear deterministic dynamics and noisy dynamics from a time series. In this section, we review and discuss the surrogate data test method[1] and Volterra-Wiener-Korenberg (VWK) model test method[2] for identifying the nonlinearity of a time series. These methods have been successfully used to detect and characterize nonlinear dynamics from recordings in biology and medicine[2-5].

Surrogate analysis is currently an important empirical technique of testing nonlinearity for a time series. The aim is to test whether the dynamics are consistent with linearly filtered noise or a nonlinear dynamical system[1, 6]. The basic idea of the surrogate data method is to first specify some kind of linear stochastic process as a null hypothesis that mimics "linear properties" of the original data. According to the null hypothesis, surrogate data sets are generated. Then, a discriminating statistic is calculated for the original and for each of the surrogate data sets. If the statistic of the original data is significantly different from those of surrogate data sets, the null hypothesis can be rejected within a certain confidence level. It shows that the original data is from a nonlinear dynamical system. The method is demonstrated for numerical time series generated by known chaotic systems and applied to the analysis of SEMG.

VWK test method is a kind of nonlinear detection of time series based on linear and nonlinear Volterra-Wiener-Korenberg model [2, 5]. That is, it first produces the linear and nonlinear predicted data from the original time series and then compares their information criterions to detect the nonlinearity of the raw data. VWK test technique is capable of robust and highly sensitive statistical detection of deterministic dynamics, including chaotic dynamics, in experimental time series. This method is superior to other techniques when applied to short time series, either continuous or discrete, even when heavily contaminated with noise or in the presence of strong periodicity. Later, an extension of the Volterra algorithm (called the numerical titration algorithm) was given to detect and quantify chaos in noisy time series[7]. Here, the surrogate data method and VWK test approach are used to analyze the nonlinearity of surface EMG signals.

#### **2.1. Surrogate data test method**

Surrogate data method includes two parts: a null hypothesis and a test statistic. The null hypothesis is a specific process which may or may not adequately explain an origin of the data. The test statistic provides a quantitative description to demonstrate the data sources.

### *2.1.1. Null hypotheses and algorithms[1]*

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

**2. Detecting nonlinearity of the surface EMG signals** 

dynamical theory.

in biology and medicine[2-5].

the analysis of SEMG.

analyze the nonlinearity of surface EMG signals.

**2.1. Surrogate data test method** 

sources.

build on a broad and strong foundation of nonlinear time series analysis and chaotic

In many areas of science and engineering, it is a critical issue to determine whether an observed time series is purely stochastic, or deterministic nonlinear, even chaotic. One may know about the intrinsic properties of the observed phenomenon by distinguishing between nonlinear deterministic dynamics and noisy dynamics from a time series. In this section, we review and discuss the surrogate data test method[1] and Volterra-Wiener-Korenberg (VWK) model test method[2] for identifying the nonlinearity of a time series. These methods have been successfully used to detect and characterize nonlinear dynamics from recordings

Surrogate analysis is currently an important empirical technique of testing nonlinearity for a time series. The aim is to test whether the dynamics are consistent with linearly filtered noise or a nonlinear dynamical system[1, 6]. The basic idea of the surrogate data method is to first specify some kind of linear stochastic process as a null hypothesis that mimics "linear properties" of the original data. According to the null hypothesis, surrogate data sets are generated. Then, a discriminating statistic is calculated for the original and for each of the surrogate data sets. If the statistic of the original data is significantly different from those of surrogate data sets, the null hypothesis can be rejected within a certain confidence level. It shows that the original data is from a nonlinear dynamical system. The method is demonstrated for numerical time series generated by known chaotic systems and applied to

VWK test method is a kind of nonlinear detection of time series based on linear and nonlinear Volterra-Wiener-Korenberg model [2, 5]. That is, it first produces the linear and nonlinear predicted data from the original time series and then compares their information criterions to detect the nonlinearity of the raw data. VWK test technique is capable of robust and highly sensitive statistical detection of deterministic dynamics, including chaotic dynamics, in experimental time series. This method is superior to other techniques when applied to short time series, either continuous or discrete, even when heavily contaminated with noise or in the presence of strong periodicity. Later, an extension of the Volterra algorithm (called the numerical titration algorithm) was given to detect and quantify chaos in noisy time series[7]. Here, the surrogate data method and VWK test approach are used to

Surrogate data method includes two parts: a null hypothesis and a test statistic. The null hypothesis is a specific process which may or may not adequately explain an origin of the data. The test statistic provides a quantitative description to demonstrate the data The null hypotheses usually specify some certain properties of the original data that reflect some structure characteristics of the dynamical system, such as mean and variance, and possibly also the Fourier power spectrum. Different null hypotheses describe different specific dynamical systems. In terms of the corresponding null hypothesis, the surrogate data can be generated so as to test the corresponding specific dynamical system class.

**Null hypothesis 1** The observed data is produced by independent and identically distributed (IID) random variables.

For this hypothesis, the corresponding surrogate data can be generated by shuffling the time-order of the original time series so that it has the same mean, variance and amplitude distribution as the original data. But any temporal correlations of the original data are destroyed in the surrogate data. Schienkman and LeBaron[8] applied this hypothesis to analyze stock market returns. Breeden and Packard also used this hypothesis to demonstrate that a time series of quasar data which were sampled nonuniformly in time has some dynamics structure[9].

The algorithm of the null hypothesis is that one first create gaussian random numbers from 1 to *N*, where *N* is the length of the original data *x*. Then, the original data *x* is permuted by the random numbers to generate the surrogate data.

**Null hypothesis 2** The observed data is produced by the Ornstein-Uhlenbeck process.

The surrogate data generated by the Ornstein-Uhlenbeck process is a sequence that has the simplest time correlation. The Ornstein-Uhlenbeck process can be given as follows.

$$\mathbf{x}\_{t} = a\_{0} + a\_{1}\mathbf{x}\_{t-1} + \mathbf{oc}\_{t} \tag{1}$$

where *et* is a Gaussian random with zero mean and unit variance. The coefficients *a*0, *a*1, and σ work together to determine the mean, variance, and autocorrelation time of the time series *xt*. In this case, its autocorrelation function is exponential form. Let <sup>1</sup> λ = − log *a* , ⋅ denotes an average over time *t*. That is,

$$A(\mathbf{r}) \equiv \frac{\left\langle \mathbf{x}\_t \cdot \mathbf{x}\_{t-\tau} \right\rangle - \left\langle \mathbf{x}\_t \right\rangle^2}{\left\langle \mathbf{x}\_t^2 \right\rangle - \left\langle \mathbf{x}\_t \right\rangle^2} = e^{-\lambda \|\mathbf{r}\|}\tag{2}$$

In order to generate the surrogate data that is consistent with this hypothesis, the algorithm is that one first calculates the mean μ , variance γ and autocorrelation *A*(1) (in Eq. (2)) from the original data *x*. Then, the coefficients in Eq. (1) can be estimated: (1) *a*<sup>1</sup> = *A* , (1 ) <sup>0</sup> <sup>1</sup> *a* = μ − *a* , and (1 ) 2 1 2 σ = γ − *a* . The Gaussian *et* can be generated by a pseudorandom number generator. Finally, the surrogate data can be produced by iterating Eq. (1).

**Null hypothesis 3** The observed data is produced by the linear autocorrelated gaussian process with the mean and variance of the original time series.

The hypothesis has been usually used to test whether the original time series contains nonlinear components. It can be described by using a linear autoregressive (AR) model.

Computational Intelligence in Electromyography Analysis –

122 A Perspective on Current Applications and Future Challenges

$$\mathbf{x}\_{t} = a\_{0} + \sum\_{k=1}^{q} a\_{k} \mathbf{x}\_{t-k} + \sigma \mathbf{e}\_{t} \tag{3}$$

Nonlinear Analysis of Surface EMG Signals 123

(*fk*), *i*=2~(*N*+1)/2, *k*=*N*~(*N*+1)/2+1; If *N*

(10)

(*fk*), *i*=2~*N*/2, *k*=*N*~*N*/2+2. Thus, the surrogate data *x'*(*n*)

φ<sup>=</sup> 0)0( ,

given by the inverse Fourier transform is a sequence of real numbers.

*<sup>N</sup> nx*

In practical, if the data length *N* is odd,

(*fN/*2+1)=0,

ϕ(*fi*)= ϕ

is even,

ϕ(*f*1)=0, ϕ

two limitations[1, 12].

gaussian process.

φ

ϕ(*f*1)=0, ϕ(*fi*) =ϕ

′ <sup>=</sup> ′ <sup>⋅</sup> <sup>−</sup> = 1 0 /2 )( <sup>1</sup> )( *<sup>N</sup> k*

a b

Thus, there are no imaginary components (see Fig. (1a)). The values of the imaginary parts are very little (magnitude 10-14) so that they can be regarded as computing precision errors. The surrogate data has the same Fourier transform spectrum as the original data by using this algorithm. The reproduced "pure" frequencies are very well. Fig.(1b)shows that the previous FT algorithm[1] cannot make the imaginary components of Fourier inverse transform to be zero. So, if one only uses the real part of Fourier inverse transform as surrogate data and omit its imaginary components, the obtained surrogates would have the

**Null hypothesis 4** The observed data is produced by the static nonlinear transform of linear

The static nonlinear transform is that the observation or measure function is nonlinear. The static means that the measure data *xt* only depends on the state *yt* of the dynamic process at

The generated surrogate data not only contain the linear correlated characteristic, but also can reflect the static, monotonic nonlinearity of the original data. Strictly speaking, time series in this class are nonlinear. But this nonlinearity is not from the dynamics. This hypothesis can be used to indicate whether the nonlinearity is from the dynamical system or

( )*tt* = *yhx* (11)

the time *t* , not on derivatives or values in the past. Let *h* be a measure function, then

the amplitude distribution (i.e. the measure process).

**Figure 1.** The imaginary components of surrogate data by our (a) and previous (b) FT algorithm

*<sup>N</sup>* <sup>=</sup> 0)2/( .

*Nink ekX*

π

There are the two algorithms to produce the surrogate data in accord with this hypothesis. One algorithm is to directly use Eq.3. That is, the coefficients are firstly identified by using the original data. Then, the surrogate data is generated by repeatedly iterating Eq.3. However, the performance of this algorithm is very unstable. If the values of the coefficients are mis-estimated slightly, this algorithm may lead to the iterative results which easily diverge to infinity. The alternative algorithm is that a surrogate data is generated by randomizing the phases of a Fourier transform. According to the Weiner-Khintchine theorem, the two algorithms are equivalent in essence[1, 10]. The surrogate data has the same Fourier spectrum as the original data. Meanwhile, the algorithm based on the Fourier transform is stabler in the numerical calculation than the first algorithm. The following is the steps of this algorithm.

Let an observed data as *x*(*n*) . The Fourier transform of *x*(*n*) is computed as follows[3]:

$$X(k) = \sum\_{n=0}^{N-1} x(n)e^{-2\pi ink/N} \tag{4}$$

The Fourier transform has a complex amplitude at each frequency. One can randomize the phases of the Fourier transform by multiplying *<sup>i</sup>*ϕ*e* ,

$$X'(k) = X(k) \cdot e^{\prime \varphi} \tag{5}$$

where ϕ is independently chosen for each frequency from the [0, 2 π ]. In order to the inverse Fourier transform to be real (no imaginary components), the phases must satisfy the antisymmetric condition ϕ() ( ) *k* = −ϕ *N* − *k* . Meanwhile, ϕ( ) 0 = 0 , ϕ( ) *N* / 2 = 0 (when *N* is even), so that

$$X'(k) = \overline{X'(N-k)}\tag{6}$$

This point can be easily proved[11].

Proof:

According to the nature of DFT of a real time series *x*(*n*), if *x*(*n*)∈*R* , then

$$
\phi(k) = -\phi(-k)\tag{7}
$$

where φ(*k*) is the phase angle of *X* (*k*) . *k* = 0, 1, … *N*-1. *N* is the period of Fourier Transform. Then, for *k* = 0, *k* = *N* /2 (*N* is even), there are

$$
\phi(0) = -\phi(0),
\tag{8}
$$

$$
\phi(N/2) = -\phi(-N/2) = -\phi(N/2) \tag{9}
$$

In order to ensure that the inverse Fourier transform results are real values, there must be

$$
\phi(0) = 0 \quad\_{\cdot} \phi(N/2) = 0 \quad\_{\cdot}
$$

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

phases of the Fourier transform by multiplying *<sup>i</sup>*

ϕ

steps of this algorithm.

where

Proof:

where φ

ϕ

even), so that

antisymmetric condition

This point can be easily proved[11].

Then, for *k* = 0, *k* = *N* /2 (*N* is even), there are

= + + <sup>=</sup> <sup>−</sup> *q <sup>k</sup> <sup>t</sup> <sup>k</sup> <sup>t</sup> <sup>k</sup> <sup>t</sup> <sup>x</sup> <sup>a</sup> <sup>a</sup> <sup>x</sup> <sup>e</sup>* <sup>1</sup> <sup>0</sup>

There are the two algorithms to produce the surrogate data in accord with this hypothesis. One algorithm is to directly use Eq.3. That is, the coefficients are firstly identified by using the original data. Then, the surrogate data is generated by repeatedly iterating Eq.3. However, the performance of this algorithm is very unstable. If the values of the coefficients are mis-estimated slightly, this algorithm may lead to the iterative results which easily diverge to infinity. The alternative algorithm is that a surrogate data is generated by randomizing the phases of a Fourier transform. According to the Weiner-Khintchine theorem, the two algorithms are equivalent in essence[1, 10]. The surrogate data has the same Fourier spectrum as the original data. Meanwhile, the algorithm based on the Fourier transform is stabler in the numerical calculation than the first algorithm. The following is the

Let an observed data as *x*(*n*) . The Fourier transform of *x*(*n*) is computed as follows[3]:

<sup>−</sup> <sup>1</sup> 0 <sup>2</sup> / ( ) ( ) *<sup>N</sup> n ink <sup>N</sup> X k x n e*

The Fourier transform has a complex amplitude at each frequency. One can randomize the

inverse Fourier transform to be real (no imaginary components), the phases must satisfy the

*N* − *k* . Meanwhile,

ϕ*e* ,

(*k*) is the phase angle of *X* (*k*) . *k* = 0, 1, … *N*-1. *N* is the period of Fourier Transform.

(−*N* / 2) = −

φ

π

*i*ϕ

*X* ′(*k* ) = *X* (*k* ) ⋅ *e* (5)

*X* ′(*k*) = *X* ′(*N* − *k*) (6)

ϕ( ) 0 = 0 ,

= − =

is independently chosen for each frequency from the [0, 2

() ( ) *k* = −ϕ

According to the nature of DFT of a real time series *x*(*n*), if *x*(*n*)∈*R* , then

φ

φ(*k*) = −φ

φ(0) = −φ

φ

In order to ensure that the inverse Fourier transform results are real values, there must be

(*N* / 2) = −

σ

(3)

(4)

π

ϕ

(−*k*) (7)

(0), (8)

(*N* / 2) (9)

]. In order to the

( ) *N* / 2 = 0 (when *N* is

In practical, if the data length *N* is odd, ϕ(*f*1)=0, ϕ(*fi*) =ϕ(*fk*), *i*=2~(*N*+1)/2, *k*=*N*~(*N*+1)/2+1; If *N* is even, ϕ(*f*1)=0, ϕ(*fN/*2+1)=0, ϕ(*fi*)= ϕ(*fk*), *i*=2~*N*/2, *k*=*N*~*N*/2+2. Thus, the surrogate data *x'*(*n*) given by the inverse Fourier transform is a sequence of real numbers.

$$\mathbf{x}'(n) = \frac{1}{N} \sum\_{k=0}^{N-1} X'(k) \cdot e^{2\pi i nk/N} \tag{10}$$

**Figure 1.** The imaginary components of surrogate data by our (a) and previous (b) FT algorithm

Thus, there are no imaginary components (see Fig. (1a)). The values of the imaginary parts are very little (magnitude 10-14) so that they can be regarded as computing precision errors. The surrogate data has the same Fourier transform spectrum as the original data by using this algorithm. The reproduced "pure" frequencies are very well. Fig.(1b)shows that the previous FT algorithm[1] cannot make the imaginary components of Fourier inverse transform to be zero. So, if one only uses the real part of Fourier inverse transform as surrogate data and omit its imaginary components, the obtained surrogates would have the two limitations[1, 12].

**Null hypothesis 4** The observed data is produced by the static nonlinear transform of linear gaussian process.

The static nonlinear transform is that the observation or measure function is nonlinear. The static means that the measure data *xt* only depends on the state *yt* of the dynamic process at the time *t* , not on derivatives or values in the past. Let *h* be a measure function, then

$$\mathbf{x}\_t = h(\mathbf{y}\_t) \tag{11}$$

The generated surrogate data not only contain the linear correlated characteristic, but also can reflect the static, monotonic nonlinearity of the original data. Strictly speaking, time series in this class are nonlinear. But this nonlinearity is not from the dynamics. This hypothesis can be used to indicate whether the nonlinearity is from the dynamical system or the amplitude distribution (i.e. the measure process).

For generating surrogate data corresponding to this null hypothesis, an algorithm is described. The aim is to shuffle the time-order of the data *xt* and to preserve the linear correlations of the underlying time series *yt* = *h*-1(*xt*). The first step is to make a Gaussian time series *yt*, where each element is generated independently from a Gaussian pseudorandom number generator. Next, we rescale *yt* in accordance with the time-order of the original data *xt*. The re-ordered *yt* has a time series which "follows" the static, monotonic nonlinearity of the original data. Then, the data *<sup>t</sup> y* ′ is created by using the above FT algorithm to deal with the re-ordered *yt*. Finally, the raw data *xt* is rescaled in terms of the time-order of the data *<sup>t</sup> y* ′ to generate the surrogate data *<sup>t</sup> x*′ . The "underlying" time series (*yt* and *<sup>t</sup> y*′ ) are Gaussian and have the same Fourier power spectrum. The produced *<sup>t</sup> x*′ matches the amplitude distribution of the raw data *xt*.

#### *2.1.2. Test statistics[6]*

The test statistic is a value which estimates some certain aspects of the time series. To compare the raw data to its surrogate data sets, a suitable test statistic must be selected. A useful statistic should be pivotal and independent of the way that surrogate data sets are generated. In other words, for every data set *z* and every realization *zi* of any *Fi*∈*F*φ, their test statistics should be different, i.e.

$$T(z) \not\cong T(z\_i) \tag{12}$$

Nonlinear Analysis of Surface EMG Signals 125

*2.1.3. Performance of surrogate data method based on our FT algorithm[3, 13, 14]* 

*B*min *p* −−= 1)1(2 . For 95% confidence level, there should be 39 sets of surrogate data.

discrepancy between the test statistic *T* of the raw data *x* and those of surrogate data.

reflects the null hypothesis 3.

series as follows.

α

otherwise it is called as measure noise.

where

The surrogate data method is suitable to detect the nonlinearity of a short, noisy time series. Here, a Gaussian data and a Logistic chaotic time series are used to study the performance of surrogate data method. For a two-sided test, the probability of rejecting the null hypothesis is given by the confidence level *p*, the surrogate data sets B must be at least as large as

A Gaussian data *x* is a random time series with zero mean and unit variance produced by the pseudorandom generator. The data length is 1000 points. According to the null hypothesis 3, 39 sets of surrogate data are generated by using our above FT algorithm. The *T* value is calculated by Eq. (13) and Eq. (14), respectively. In the Figure (2a and b), there are no statistical

The statistic *T* values of the raw data are on the range of the empirical distribution of *T* given by the surrogate data. The results show that the generated surrogate data has the same Fourier transform spectrum as the raw data besides the same mean and variance as the raw data because the *T* value in Eq. (14) is a measure of the time irreversibility of the data. The null hypothesis 3 is accepted at the confidence level 95%. The raw data is consistent with the stochastic process of the null hypothesis 3. The surrogate data produced by the above FT algorithm is equivalent to the raw data. The generated surrogate data

a. *T* calculated by Eq. (13) b. *T* calculated by Eq. (14) **Figure 2.** The histogram is *T* distribution of surrogate data given by FT algorithm, \* is *T* value of the

To further test that the surrogate data based on the above FT algorithm can be used to detect nonlinearity of a time series, we apply the Logistic chaotic system to produce a chaotic time

> *exxxt*<sup>+</sup><sup>1</sup> α

evolution process of the above equation, it is called as interior noise (or dynamic noise);

= 9.3 , ]1,0[ *x*<sup>0</sup> ∈ , *e* is white noise with mean 0, variance 0.0012. If *e* takes part in the

*tt* )1( +−= (16)

original data, where abscissa is *T*, ordinate is the number of surrogate data sets

where *F*φ represents the null hypothesis process. Meanwhile, the distribution of *T* under the null hypothesis does not depend on μ or σ. Here we give two discriminating statistics as follows:

$$T = \overline{(\mathbf{x} - \overline{\mathbf{x}})^4} \sqrt{\left(\mathbf{x} - \overline{\mathbf{x}}\right)^2} \tag{13}$$

$$T = \frac{1}{n-1} \sum\_{i=1}^{n-1} \left| (\mathbf{x}(i) - \overline{\mathbf{x}}) \cdot (\mathbf{x}(i+1) - \overline{\mathbf{x}}) \right| \sqrt{\left(\mathbf{x} - \overline{\mathbf{x}}\right)^2} \tag{14}$$

where "¯" denotes the average of the data. The mean μ and varianceσ have no effect on the *T* value in Eq. (13). Therefore, some linear structure characteristics can be determined except for the mean and variance. The *T* value in Eq.(14) can judge if the surrogate data are consist with the raw data in the view of the correlation with the mean and variance. The *T* value in Eq. (15) is a simple skewed difference statistic that is both rapidly computable and often quite powerful.

$$T = \left\langle \left(\mathbf{x}\_{t+m} - \mathbf{x}\_t\right)^3 \right\rangle \left\langle \left(\mathbf{x}\_{t+m} - \mathbf{x}\_t\right)^2 \right\rangle \tag{15}$$

where ⋅ is mean operator, *m* is time delay. This statistic *T* provides a more significant rejection of the hypothesis of the static nonlinear filter of an underlying linear process. Informally, this statistic indicates the asymmetry between rise and fall times in the time series.

#### *2.1.3. Performance of surrogate data method based on our FT algorithm[3, 13, 14]*

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

distribution of the raw data *xt*.

statistics should be different, i.e.

null hypothesis does not depend on

*2.1.2. Test statistics[6]* 

where *F*φ

follows:

quite powerful.

series.

For generating surrogate data corresponding to this null hypothesis, an algorithm is described. The aim is to shuffle the time-order of the data *xt* and to preserve the linear correlations of the underlying time series *yt* = *h*-1(*xt*). The first step is to make a Gaussian time series *yt*, where each element is generated independently from a Gaussian pseudorandom number generator. Next, we rescale *yt* in accordance with the time-order of the original data *xt*. The re-ordered *yt* has a time series which "follows" the static, monotonic nonlinearity of the original data. Then, the data *<sup>t</sup> y* ′ is created by using the above FT algorithm to deal with the re-ordered *yt*. Finally, the raw data *xt* is rescaled in terms of the time-order of the data *<sup>t</sup> y* ′ to generate the surrogate data *<sup>t</sup> x*′ . The "underlying" time series (*yt* and *<sup>t</sup> y*′ ) are Gaussian and have the same Fourier power spectrum. The produced *<sup>t</sup> x*′ matches the amplitude

The test statistic is a value which estimates some certain aspects of the time series. To compare the raw data to its surrogate data sets, a suitable test statistic must be selected. A useful statistic should be pivotal and independent of the way that surrogate data sets are

represents the null hypothesis process. Meanwhile, the distribution of *T* under the

2

μ

*<sup>x</sup> <sup>i</sup> <sup>x</sup> <sup>x</sup> <sup>i</sup> <sup>x</sup> <sup>x</sup> <sup>x</sup> <sup>n</sup> <sup>T</sup>* (14)

φ

have no effect on the

( ) ( )*<sup>i</sup> T z* ≠ *T z* (12)

<sup>4</sup> <sup>2</sup> *T* = (*x* − *x*) (*x* − *x*) (13)

and variance

<sup>3</sup> <sup>2</sup> ( ) ( ) *<sup>t</sup> <sup>m</sup> <sup>t</sup> <sup>t</sup> <sup>m</sup> <sup>t</sup> T* = *x* − *x x* − *x* <sup>+</sup> <sup>+</sup> (15)

. Here we give two discriminating statistics as

σ

, their test

generated. In other words, for every data set *z* and every realization *zi* of any *Fi*∈*F*

μ or σ

= 1 1

1 *<sup>n</sup> i*

where "¯" denotes the average of the data. The mean

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

*T* value in Eq. (13). Therefore, some linear structure characteristics can be determined except for the mean and variance. The *T* value in Eq.(14) can judge if the surrogate data are consist with the raw data in the view of the correlation with the mean and variance. The *T* value in Eq. (15) is a simple skewed difference statistic that is both rapidly computable and often

where ⋅ is mean operator, *m* is time delay. This statistic *T* provides a more significant rejection of the hypothesis of the static nonlinear filter of an underlying linear process. Informally, this statistic indicates the asymmetry between rise and fall times in the time

<sup>2</sup> ( ( ) ) ( ( 1) ) ( ) <sup>1</sup>

The surrogate data method is suitable to detect the nonlinearity of a short, noisy time series. Here, a Gaussian data and a Logistic chaotic time series are used to study the performance of surrogate data method. For a two-sided test, the probability of rejecting the null hypothesis is given by the confidence level *p*, the surrogate data sets B must be at least as large as *B*min *p* −−= 1)1(2 . For 95% confidence level, there should be 39 sets of surrogate data.

A Gaussian data *x* is a random time series with zero mean and unit variance produced by the pseudorandom generator. The data length is 1000 points. According to the null hypothesis 3, 39 sets of surrogate data are generated by using our above FT algorithm. The *T* value is calculated by Eq. (13) and Eq. (14), respectively. In the Figure (2a and b), there are no statistical discrepancy between the test statistic *T* of the raw data *x* and those of surrogate data.

The statistic *T* values of the raw data are on the range of the empirical distribution of *T* given by the surrogate data. The results show that the generated surrogate data has the same Fourier transform spectrum as the raw data besides the same mean and variance as the raw data because the *T* value in Eq. (14) is a measure of the time irreversibility of the data. The null hypothesis 3 is accepted at the confidence level 95%. The raw data is consistent with the stochastic process of the null hypothesis 3. The surrogate data produced by the above FT algorithm is equivalent to the raw data. The generated surrogate data reflects the null hypothesis 3.

**Figure 2.** The histogram is *T* distribution of surrogate data given by FT algorithm, \* is *T* value of the original data, where abscissa is *T*, ordinate is the number of surrogate data sets

To further test that the surrogate data based on the above FT algorithm can be used to detect nonlinearity of a time series, we apply the Logistic chaotic system to produce a chaotic time series as follows.

$$\mathbf{x}\_{t+1} = \alpha \mathbf{x}\_t (1 - \mathbf{x}\_t) + \mathbf{e} \tag{16}$$

where α = 9.3 , ]1,0[ *x*<sup>0</sup> ∈ , *e* is white noise with mean 0, variance 0.0012. If *e* takes part in the evolution process of the above equation, it is called as interior noise (or dynamic noise); otherwise it is called as measure noise.

Nonlinear Analysis of Surface EMG Signals 127

*2.1.4. Nonlinear test of surface EMG signal based on surrogate data[3, 13, 14]* 

thing, the length of FSEMG data is also 1000 points when the arm has been fatigue.

ensure that this nonlinearity must be from the dynamic system.

**Figure 5.** The surface EMG signals

surrogate data

The nature of SEMG plays an important role in neuromuscular disorder diagnosis, muscle fatigue monitoring, prosthesis control, etc. Here the analyzed data are collected from physiological instruments. Humid surface electrode and AD12-16LG collecting card of physiology signal are used in the whole experiment that was done at Hua Shan Hospital in Shanghai. The data are sampled at 1kHz for the action surface EMG (ASEMG) [3]and the fatigue surface EMG (FSEMG) when one hand carries a 1kg heavy thing [15](see Fig. 5). The length of data for ASEMG is 1000 points during the beginning of action because this time span contains the information of the forearm movement. In the case of carrying a 1kg heavy

For these surface EMG signals, 39 surrogate data are produced by the null hypothesis 2. The surrogate data analysis is given for the action surface EMG signal and the fatigue surface EMG signal, respectively(see Fig. 6). The results show that for action surface EMG signal and fatigue surface EMG signal, their *T* values are obviously different from those of surrogate data in terms of Eq.13. The null hypothesis 2 can be rejected in 95% degree of confidence. The action surface EMG signal and fatigue surface EMG signal is not stochastic signal produced by a linear stochastic process, but contains nonlinear components. However, this result could not

a. a typical action surface EMG wave b. a typical fatigue surface EMG wave

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

**Figure 6.** Surrogate data analysis of surface EMG signal. \* is *T*orig, histogram is *T*surr distribution of

**Figure 3.** *T* calculated by Eq. (13), histogram is *T* distribution of surrogate data, \* is *T* value of the original data, where abscissa is *T*, ordinate is the number of surrogate data sets

For the case without noise *e* = 0, we use Eq. (16) to compute the two Logistic chaotic time series with the length of 5000 points and 500 points. 39 surrogate data generated by the above FT algorithm contain the linear properties of the original data in terms of the null hypothesis 3. In Fig.3, we can see the obvious difference between the original data and its surrogate data, regardless of the length of 5000 points or 500 points. The null hypothesis can be rejected in 95% confidence level. The original data has nonlinear components. The results show that the data length has little effect on the surrogate data method based on the above FT algorithm. In Figure 4, we study the nonlinear test of Logistic chaotic time series with measure noise and interior noise, respectively. The data length is 1000 points. According to the null hypothesis, 39 sets of surrogate data are generated. The statistic *T* for the original data is significantly different from the values obtained for the surrogate data sets. The null hypothesis 3 can also be rejected in 95% significance. The nonlinearity of the original data can be detected. To sum up, the length and noise has no impact on the surrogate data method based on our FT algorithm.

**Figure 4.** *T* calculated by Eq.(13), histogram is *T* distribution of surrogate data, \* is *T* value of the original data, where abscissa is *T*, ordinate is the number of surrogate data sets

#### *2.1.4. Nonlinear test of surface EMG signal based on surrogate data[3, 13, 14]*

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

a. The length N=5000 of the raw time series b. The length N=500 of the raw time series

For the case without noise *e* = 0, we use Eq. (16) to compute the two Logistic chaotic time series with the length of 5000 points and 500 points. 39 surrogate data generated by the above FT algorithm contain the linear properties of the original data in terms of the null hypothesis 3. In Fig.3, we can see the obvious difference between the original data and its surrogate data, regardless of the length of 5000 points or 500 points. The null hypothesis can be rejected in 95% confidence level. The original data has nonlinear components. The results show that the data length has little effect on the surrogate data method based on the above FT algorithm. In Figure 4, we study the nonlinear test of Logistic chaotic time series with measure noise and interior noise, respectively. The data length is 1000 points. According to the null hypothesis, 39 sets of surrogate data are generated. The statistic *T* for the original data is significantly different from the values obtained for the surrogate data sets. The null hypothesis 3 can also be rejected in 95% significance. The nonlinearity of the original data can be detected. To sum up, the length and noise has no impact on the surrogate data

a. Chaos time series with measure noise b. Chaos time series with interior noise

**Figure 4.** *T* calculated by Eq.(13), histogram is *T* distribution of surrogate data, \* is *T* value of the

original data, where abscissa is *T*, ordinate is the number of surrogate data sets

**Figure 3.** *T* calculated by Eq. (13), histogram is *T* distribution of surrogate data, \* is *T* value of the

original data, where abscissa is *T*, ordinate is the number of surrogate data sets

method based on our FT algorithm.

The nature of SEMG plays an important role in neuromuscular disorder diagnosis, muscle fatigue monitoring, prosthesis control, etc. Here the analyzed data are collected from physiological instruments. Humid surface electrode and AD12-16LG collecting card of physiology signal are used in the whole experiment that was done at Hua Shan Hospital in Shanghai. The data are sampled at 1kHz for the action surface EMG (ASEMG) [3]and the fatigue surface EMG (FSEMG) when one hand carries a 1kg heavy thing [15](see Fig. 5). The length of data for ASEMG is 1000 points during the beginning of action because this time span contains the information of the forearm movement. In the case of carrying a 1kg heavy thing, the length of FSEMG data is also 1000 points when the arm has been fatigue.

For these surface EMG signals, 39 surrogate data are produced by the null hypothesis 2. The surrogate data analysis is given for the action surface EMG signal and the fatigue surface EMG signal, respectively(see Fig. 6). The results show that for action surface EMG signal and fatigue surface EMG signal, their *T* values are obviously different from those of surrogate data in terms of Eq.13. The null hypothesis 2 can be rejected in 95% degree of confidence. The action surface EMG signal and fatigue surface EMG signal is not stochastic signal produced by a linear stochastic process, but contains nonlinear components. However, this result could not ensure that this nonlinearity must be from the dynamic system.

**Figure 6.** Surrogate data analysis of surface EMG signal. \* is *T*orig, histogram is *T*surr distribution of surrogate data

Nonlinear Analysis of Surface EMG Signals 129

where *r* ∈[1, *M* ] is the number of polynomial terms of the truncated Volterra expansions

For each data series, there is the following numerical procedure to search for the optimal

4. Compare *C*lin(*r*) and *C*nl(*r*), if *C*nl(*r*) is obviously smaller than *C*lin(*r*), then the original system dynamics is nonlinear, the obtained time series is nonlinear, even chaos;

Note that when *k*opt is rather large, *M* is quite large, too, then the computational time will rapidly go up. In this case, *k* and *d* should be adjusted synchronously to search for *k*opt and *d*opt so as to make *C*nl(*r*) < *C*lin(*r*). Furthermore, one can obtain the corresponding linear and nonlinear models for surrogate data generated by the FT algorithm according to the null

*orig*

Here, the VWK method is used to deal with the surface EMG signals in Fig.5. For the action surface EMG signal, *Clin*(*r*) is almost equal to *Cnl*(*r*), i.e. *C* (*r*) *C* (*r*) *lin nl* ≈ , this technique can hardly detect its nonlinearity (see Fig.8a). For the fatigue surface EMG signal, *Cnl*(*r*) is distinctly smaller than *Clin*(*r*) so that its nonlinear component can be detected (see Fig.8b). So VWK technique can effectively detect the nonlinear dynamic speciality of fatigue surface EMG signal but fails to test the nonlinearity of the action surface EMG signal. In other words, VWK technique can not be used directly to deal with

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

*orig* > in the statistical sense.

otherwise, the original system dynamics is linear, the raw data is linear.

*nl surr*

*<sup>N</sup> <sup>n</sup> <sup>n</sup> <sup>y</sup> <sup>y</sup>* <sup>1</sup> <sup>1</sup> . For *d*=1, VWK model is linear, whereas the model is nonlinear for *d*>1.

(*k* , *d* ) is a normalized variance of the error residuals,

from the given pair {*k*, *d*}, <sup>2</sup>

= <sup>=</sup> *N*

pair {*k*opt, *d*opt}:

ε

2. with *k*=*k*opt, increasing *d*>1, search for *d*opt which minimizes *C*(*r*). 3. calculate *C*lin(*r*) with *d*=1 and *k*=*M*-1, and *C*nl(*r*) with *d*=*d*opt and *k*=*k*opt.

1. when *d*=1, search for *k*opt which minimizes *C*(*r*).

hypothesis 3 so that *C* (*r*), *C* (*r*), *C* (*r*) *C* (*r*) *nl*

*2.2.2. Analysis of SEMG based on VWK method[15, 16]* 

*lin*

the action surface EMG signal.

**Figure 8.** VWK test analysis of surface EMG signal

*lin surr*

**Figure 7.** Surrogate data test of surface EMG signal during movement, where surrogate data sets are 39 sets; \* is *T* value of surrogate data by the null hypothesis 4, +is *T* value of EMG signal, where *T* is calculated by Eq. 15

In order to test that the nonlinear components are intrinsic deterministic, we further assume that ASEMG is stochastic signal consistent with the null hypothesis 4. Fig.7 gives the *T* values of ASEMG and surrogate data calculated by Eq. 15. This statistic indicates the asymmetry between rise and fall times in the time series. From this figure, we can see that there is the difference between data and surrogates, and the null hypothesis 4 is rejected in 95% credibility. This result shows that the nonlinearity of ASEMG is intrinsic and deterministic.

#### **2.2. Volterra-Wiener-Korenberg test method**

#### *2.2.1. Volterra-Wiener-Korenberg test model[2]*

For a dynamic system, an observed time series *<sup>N</sup> <sup>n</sup> <sup>n</sup> y* <sup>1</sup> { } <sup>=</sup> can be treated as a closed loop Volterra series by utilizing the feedback of *yn*. Suppose the time series is univariate. A discrete Volterra-Winner-Korenberg series can be calculated as follows:

$$\begin{aligned} \mathbf{y}\_n^{\text{cal}} &= a\_0 + a\_1 \mathbf{y}\_{n-1} + a\_2 \mathbf{y}\_{n-2} + \dots + a\_k \mathbf{y}\_{n-k} + a\_{k+1} \mathbf{y}\_{n-1}^2 + a\_{k+2} \mathbf{y}\_{n-1} \mathbf{y}\_{n-2} + \dots + a\_{M-1} \mathbf{y}\_{n-k}^T\\ &= \sum\_{m=0}^M a\_m z\_m(n) \end{aligned} \tag{17}$$

where the memory *k* and combination degree *d* correspond to the embedding dimension and the degree of nonlinearity of the model, respectively. The coefficients *am* are recursively estimated through a Gram-Schmidt procedure from linear and nonlinear autocorrelations of the data itself with a total dimension *M*=(*k*+*d*)!/(*d*!*k*!).

There is the following information criterion in accordance with the parsimony principle:

$$C(r) = \log \varepsilon(r) + r/N \tag{18}$$

$$\varepsilon(k,d)^2 \equiv \frac{\sum\_{n=1}^{N} \left(\chi\_a^{calc}(k,d) - \chi\_n\right)^2}{\sum\_{n=1}^{N} \left(\chi\_n - \overline{\chi}\right)^2} \tag{19}$$

where *r* ∈[1, *M* ] is the number of polynomial terms of the truncated Volterra expansions from the given pair {*k*, *d*}, <sup>2</sup> ε (*k* , *d* ) is a normalized variance of the error residuals, = <sup>=</sup> *N <sup>N</sup> <sup>n</sup> <sup>n</sup> <sup>y</sup> <sup>y</sup>* <sup>1</sup> <sup>1</sup> . For *d*=1, VWK model is linear, whereas the model is nonlinear for *d*>1.

For each data series, there is the following numerical procedure to search for the optimal pair {*k*opt, *d*opt}:

1. when *d*=1, search for *k*opt which minimizes *C*(*r*).

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

**2.2. Volterra-Wiener-Korenberg test method** 

For a dynamic system, an observed time series *<sup>N</sup>*

discrete Volterra-Winner-Korenberg series can be calculated as follows:

0 1 1 2 2 1 1

*2.2.1. Volterra-Wiener-Korenberg test model[2]* 

*a z n*

( )

the data itself with a total dimension *M*=(*k*+*d*)!/(*d*!*k*!).

=

*cala n*

=

0

*M <sup>m</sup> <sup>m</sup> <sup>m</sup>*

calculated by Eq. 15

deterministic.

**Figure 7.** Surrogate data test of surface EMG signal during movement, where surrogate data sets are 39 sets; \* is *T* value of surrogate data by the null hypothesis 4, +is *T* value of EMG signal, where *T* is

In order to test that the nonlinear components are intrinsic deterministic, we further assume that ASEMG is stochastic signal consistent with the null hypothesis 4. Fig.7 gives the *T* values of ASEMG and surrogate data calculated by Eq. 15. This statistic indicates the asymmetry between rise and fall times in the time series. From this figure, we can see that there is the difference between data and surrogates, and the null hypothesis 4 is rejected in 95% credibility. This result shows that the nonlinearity of ASEMG is intrinsic and

Volterra series by utilizing the feedback of *yn*. Suppose the time series is univariate. A

− − − + − + − − − −

2

*n n k n k k n k n n M n k*

= + + + + + + + +

*y a a y a y a y a y a y y a y*

where the memory *k* and combination degree *d* correspond to the embedding dimension and the degree of nonlinearity of the model, respectively. The coefficients *am* are recursively estimated through a Gram-Schmidt procedure from linear and nonlinear autocorrelations of

There is the following information criterion in accordance with the parsimony principle:

ε

−

*y y*

( )

2

2

( , ) (19)

=

1

−

*<sup>n</sup> <sup>n</sup> calc a*

*y k d y*

( ( , ) )

*C*(*r*) = log

≡

2

*k d*

ε

= *N <sup>n</sup> <sup>n</sup>*

1

*N*

*<sup>n</sup> <sup>n</sup> y* <sup>1</sup> { } <sup>=</sup> can be treated as a closed loop

2 1 2 1

(*r*) + *r N* (18)

(17)

*d*


Note that when *k*opt is rather large, *M* is quite large, too, then the computational time will rapidly go up. In this case, *k* and *d* should be adjusted synchronously to search for *k*opt and *d*opt so as to make *C*nl(*r*) < *C*lin(*r*). Furthermore, one can obtain the corresponding linear and nonlinear models for surrogate data generated by the FT algorithm according to the null hypothesis 3 so that *C* (*r*), *C* (*r*), *C* (*r*) *C* (*r*) *nl orig nl surr lin surr lin orig* > in the statistical sense.

#### *2.2.2. Analysis of SEMG based on VWK method[15, 16]*

Here, the VWK method is used to deal with the surface EMG signals in Fig.5. For the action surface EMG signal, *Clin*(*r*) is almost equal to *Cnl*(*r*), i.e. *C* (*r*) *C* (*r*) *lin nl* ≈ , this technique can hardly detect its nonlinearity (see Fig.8a). For the fatigue surface EMG signal, *Cnl*(*r*) is distinctly smaller than *Clin*(*r*) so that its nonlinear component can be detected (see Fig.8b). So VWK technique can effectively detect the nonlinear dynamic speciality of fatigue surface EMG signal but fails to test the nonlinearity of the action surface EMG signal. In other words, VWK technique can not be used directly to deal with the action surface EMG signal.

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

**Figure 8.** VWK test analysis of surface EMG signal

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

Nonlinear Analysis of Surface EMG Signals 131

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.

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

**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

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

*f x f y x y S x y <sup>n</sup> <sup>n</sup>*

*f x f y x y S <sup>n</sup> <sup>n</sup>*

lim sup *f* (*x*) *f* ( *p*) 0, *x S* , *p P*( *f* ) *<sup>n</sup> <sup>n</sup> <sup>n</sup>* <sup>−</sup> <sup>&</sup>gt; <sup>∀</sup> <sup>∈</sup> <sup>∀</sup> <sup>∈</sup> →∞

*<sup>n</sup>* <sup>−</sup> <sup>&</sup>gt; <sup>∀</sup> <sup>∈</sup> <sup>≠</sup> →∞ lim sup ( ) ( ) 0, , , (21)

*<sup>n</sup>* <sup>−</sup> <sup>=</sup> <sup>∀</sup> <sup>∈</sup> →∞ lim inf ( ) ( ) <sup>0</sup> , , (22)

(23)

*f* : *I I R* , *y f* (*x*) *<sup>m</sup>* → ⊂ = (20)

Then we use these methods to investigate the surface EMG signals.

commonly used. Here describes the definition of chaos by Li-Yorke[18].

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.

2. There is an uncountable set *S* ⊂ *I* , which satisfies the following conditions:

**3.1. Chaos and its definition** 

conditions, then it has chaotic motion:

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

1. The period of periodic point of *f* has no upper bound.

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

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

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 signal, *C* (*r*) *nl orig* is always smaller than *C* (*r*) *nl surr* . The null hypothesis 3 can be rejected in 95% significance. The results illuminate that the action and fatigue surface EMG signals contain nonlinear dynamic properties.
