Detecting nonlinearity of action surface EMG signal

19 November 2001

Physics Letters A 290 (2001) 297–303 www.elsevier.com/locate/pla

Detecting nonlinearity of action surface EMG signal Min Lei a,∗ , Zhizhong Wang b , Zhengjin Feng a a Institute of Mechatronic Control System, Shanghai Jiao Tong University, Shanghai 200030, PR China b Department of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200030, PR China

Received 12 December 2000; received in revised form 4 October 2001; accepted 10 October 2001 Communicated by A.R. Bishop

Abstract The action surface EMG (ASEMG) signal contains the electrical properties of limb muscle contraction that undergo many complex transitions in limb different movement states; however, because of the small data nature of this kind signal, it is not clear whether its essence is stochastic or deterministic nonlinear (even chaotic). In this Letter, we show for the first time that ASEMG has deterministic nonlinear character by using the surrogate data method. Furthermore, we study the nonlinear dynamic features of ASEMG by computing its correlation dimension and applying correlation dimension as test statistics. These results indicate that ASEMG is a high-dimension nonlinear signal (even chaotic). In addition, this Letter improves the surrogate data method based on Fourier transform (FT) algorithm to avoid limitations of the previous FT algorithm.  2001 Elsevier Science B.V. All rights reserved. Keywords: Deterministic nonlinearity; Action surface EMG signal; Surrogate data method; Small time series

1. Introduction The surface EMG (SEMG) signal is currently widely used in several fields such as physiological muscle assessment [1], muscle metabolic activity [2], functional electrical stimulation [3] (FES), neurology to evaluate muscle fiber conduction velocity [4,5], rehabilitation, sport and geriatric medicine for its convenience and noninvasiveness [6]. The nature of SEMG plays an important role in neuromuscular disorder diagnosis, muscle fatigue monitoring, prosthesis control, etc. At present, although many researchers have still been studying the linear properties of SEMG, some authors have begun studying the nonlinearity of

* Corresponding author.

E-mail address: [email protected] (M. Lei).

SEMG during isotonic contraction [7,8], during isokinetic contraction [9], and during isometric contraction [10–12] by computing fractal dimension (FD). Yang et al. [7] studied the complexity and dynamic characteristic of SEMG during isotonic contraction by presenting the phase space reconstruction analysis of SEMG during muscle fatigue and by computing the correlation dimension and maximum Lyapunov exponent, the data length is 1024 points. Yang et al. [8] analyzed the nonlinear characteristics of SEMG in the muscle contracted state and relaxed state by computing complexity, correlation time, correlation dimension, Lyapunov exponent and entropy, the data length is 20000 points. During isokinetic contraction, Ref. [9] applied 10000 point data to analyze fractal dimension of SEMG. To sum up, one often first assumes whether SEMG is stochastic signal or nonlinear signal, then applies relevant analysis methods (linear or nonlinear ways) to

0375-9601/01/$ – see front matter  2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 5 - 9 6 0 1 ( 0 1 ) 0 0 6 6 8 - 5

298

M. Lei et al. / Physics Letters A 290 (2001) 297–303

deal with it. Especially for ASEMG, people only use various methods to extract movement information to decide the intended command of limb movements for prosthesis control [13–16]. It is known that the essential problem using which kind of the disposals is what the dynamics of ASEMG is. However, because of the small data number feature of ASEMG, few researchers have studied and detected whether ASEMG has deterministic nonlinear properties or not. How and by which method can one test whether ASEMG has linear or nonlinear properties (even chaotic)? This is an untouched problem. Here, we first propose to apply the method of testing data nonlinearity: the surrogate data method, to detect the nonlinear dynamic features of ASEMG, because this method is fit for the small data set. In this Letter, first we give an introduction to ASEMG to be analyzed. Then, the surrogate data methods based on the theories of chaotic time series analysis are reviewed. The two testing ways are applied to analyze ASEMG obtained from normal people in Section 4. Conclusions and discussions are given in the last section.

2. Method Surrogate analysis is currently important empirical way of testing nonlinearity of time series, that enables us to test whether the dynamic are consistent with linearly filtered noise or a nonlinear dynamical system. It has been widely and rapidly employed in the study of related chaotic time series [17–19]. Surrogate analysis includes two parts: the null hypotheses and test statistics. Different types of surrogate data are generated to test membership of specific dynamical system classes, referred to as hypotheses. The three types of surrogates described by Theiler et al. [20], referred to as algorithm 0, 1 and 2, address the three hypotheses: (0) linearly filtered noise; (1) linear autocorrelated Gaussian noise; (2) static nonlinear transform of linear Gaussian noise. The second null hypothesis is often applied to test nonlinearity of the raw data. This Letter employs it to detect nonlinearity of ASEMG. In this section, we first briefly review this method, then give an improved algorithm. First, the Fourier transform is computed for positive and negative frequencies f = 0, 1/N, 2/N , . . . , 1/2,

and without the benefit of windowing. The Fourier transform has a complex amplitude at each frequency; to randomize the phases, we multiply each complex amplitude by eiφ , where φ is independently chosen for each frequency from the interval [0, 2π]. In order to the inverse Fourier transform to be real (no imaginary components), one must symmetrize the phases, so that φ(f ) = −φ(−f ). Finally, the inverse Fourier transform is the surrogate data. However, there are two limitations in it [20]. One limitation of this algorithm is that it does not reproduce “pure” frequencies very well. The author thought that this may not be too surprising since it is difficult to make a linear stochastic process with a long coherence time. The second problem, which is most evident for highly sampled continuous data, is that spurious high frequencies can be introduced. This could be understood as an artifact of the Fourier transform which assumes the time series is periodic with period N . For the first problem, though several researchers have studied and improved it in some extent [20,26], those results are not still very well. This Letter thinks the reason of these limitations is the nonreasonableness of the used FT algorithm itself. Because, in randomizing the phase program, one not only must assure φ(f ) = −φ(−f ), but also φ(0) = φ(f/2) = 0. This point can be easily proved [29]. However, it did not seem to be noticed in the previous literature because if so, the above problems will not exist in the generated surrogate data. (We repeated to produce surrogates by the previous FT algorithm many times but “pure” frequencies are not always reproduced very well.) Here, we give an improved FT algorithm: random phase ϕ should be uniform random among the interval [−π, π] and antisymmetric except for the first and (or) middle phases (according to the number of the original data). If the data length N is odd, ϕ(f1 ) = 0, ϕ(fi ) = −ϕ(fk ), i = 2–(N + 1)/2, k = N–(N + 1)/2 + 1; if N is even, ϕ(f1 ) = 0, ϕ(fN/2+1 ) = 0, ϕ(fi ) = −ϕ(fk ), i = 2–N/2, k = N –N/2 + 2. Thus, there are no imaginary components (see Fig. 1(a); the result in it is very minute and can be regarded as computer errors) and accordingly the reproduced “pure” frequencies are very well. Fig. 1(b) shows that the previous FT algorithm cannot make the imaginary components of Fourier inverse transform to be zero (but from FT theory, its imaginary components should be zero). So, if one only uses the real part of Fourier inverse transform

M. Lei et al. / Physics Letters A 290 (2001) 297–303

299

where · is mean operator, m is time delay [20], C(l) =

n n  1   θ l − Xi − Xj , 2 n i=1 j =1

T = D2 = lim

l→0

(a)

(b) Fig. 1. The imaginary components of surrogate data by improved (a) and previous (b) FT algorithm.

ln C(l) , ln l

(3)

where θ (X) is Heaviside function whose value is 1 if X
0; otherwise, and · is the usual distance function in R d , l is the radius of hypersphere in R d , n is the number of X, n = N − (d − 1)τ , N is the length of time series, D2 is called correlation dimension [22]. The T value in Eq. (1) is invariant to change in translation or scale; therefore it has the same distribution for all Gaussians, regardless of µ and σ . The T value in Eq. (2) is a simple skewed difference statistic [20] that is both rapidly computable and often quite powerful. (Informally, this statistic indicates the asymmetry between rise and fall times in the time series.) The T value in Eqs. (3) is a very complicated statistic, correlation dimension defined by GP algorithm. Since correlation dimension is a measure of the complexity and the number of active degree of freedom of a system this is a relatively simple matter [23]. When using correlation dimension to test composite hypothesis, one can get some more convincing results.

3. Materials as surrogate data and omit its imaginary components, the obtained surrogates will have the above two limitations. To compare the data to its surrogates, a suitable test statistic must be selected. In a recent paper, Theiler et al. [21] suggest that a useful statistic should be pivotal and independent of the way surrogates are generated. That is, for every data set z and every realization zi of any Fi ∈ Fφ , T (z) = T (zi ), where Fφ represents the hypothesis. Here we give three discriminating statistics as follows:
2 ¯ 4 (x − x) ¯ 2 , T = (x − x)

(1)

where we write f (x) to denote the sample mean (1/n) ni=1 f (x(i)) [21],  
  T = (xt +m − xt )3 (xt +m − xt )2 ,

(2)

The data to be analyzed are got from physiological instruments. Humid surface electrode and AD1216LG collecting card of physiology signal made by Contec. are used in the whole experiment. The 100 sets of ASEMG are sampled at 1 kHz from each of four different postures: finger flexion, finger tension, forearm supination and forearm pronation. Each data set consists of two channels. Channel 1 is from forearm flexor of normal people. Channel 2 is from forearm extensor. All data are directly stored in the computer. The whole experiment process was done at Hua Shan Hospital in Shanghai. A typical waveform of ASEMG during action is shown in Fig. 2. This data is arbitrarily chosen from sampled signals. From Fig. 2, one can see that ASEMG belong to small data sets (other data also have the similar results). To reveal the mechanism underlying ASEMG, the frequency spectrum of ASEMG of Fig. 2 is given in Fig. 3. From Fig. 3, we can

300

M. Lei et al. / Physics Letters A 290 (2001) 297–303

Fig. 2. The typical surface EMG wave during forearm supination (the upper is channel 1 wave, the lower is channel 2 wave), abscissa is data number, ordinate is amplitude value of EMG (mv).

Fig. 4. Surrogate data constructed by the null hypothesis 1 (abscissa is number, ordinate is value mv).

Fig. 3. The frequency spectrum of EMG in Fig. 2 (ordinate is frequency spectrum value, abscissa is frequency).

see ASEMG is not period or quasi-period (other data also have the similar results). This shows ASEMG can be stochastic or certain deterministic nonlinear signal. This Letter proposes ASEMG could be deterministic nonlinear signal. Here, we regard each set data as an approximate stationary process, then time series analysis method such as the surrogate technique is applied to test ASEMG data. To clearly show the process of analyzing ASEMG signal, in next sections, we take channel 1 of ASEMG of forearm supination in Fig. 2 for an example to analyze ASEMG and abbreviate it to ASEMG data. The length of data for analysis selected is 1000 points during the beginning of action because this span contains the information of action data [13].

4. Analysis and results For ASEMG data recorded in our experiments, we first assume ASEMG data is consistent with the null hypothesis 1 in Section 2. Then, according to the null hypothesis 1, the 39 surrogates [21] are produced by the improved FT algorithm. Next, we analyze the non-

Fig. 5. T distribution of surrogate data sets generated by the null hypothesis 1. ∗ is T value of EMG signal, where T is calculated by Eq. (1).

linearity of ASEMG. One set of surrogate data is given in Fig. 4. The produced surrogates have not only the same mean value and variance as those of raw data, but also the same power spectrum property as that of the raw. T values of ASEMG and the surrogates calculated by Eq. (1) are given in Fig. 5 by histogram. From this figure, one can see that there is the significant difference between ASEMG and the surrogates, then the null hypothesis 1 is rejected in 95% degree of confidence. This result shows that ASEMG is not stochastic signal produced by linear stochastic process, but contains nonlinear components. However, this result could not be sure this nonlinearity must be nonlinearity of dynamic system. For this, we further assume that ASEMG is stochastic signal consistent with the null hypothesis 2 to test that its nonlinear components are intrinsic deterministic. Fig. 6 gives the T values of ASEMG and surrogates calculated by Eq. (2). This statistic indicates the asymmetry between rise and fall times in the time series. The most direct example we know is due to Brock et al. [24], who used technical

M. Lei et al. / Physics Letters A 290 (2001) 297–303

Fig. 6. T value of surrogate data by the null hypothesis 2 (∗) and EMG signal (+), where T is calculated by Eq. (2) (abscissa is m, ordinate is T value).

trading rules as discriminating statistics for financial data. From Fig. 6, we can see that there is the difference between data and surrogates, and the null hypothesis 2 is rejected in 95% credibility. This result shows that the nonlinearity of ASEMG is intrinsic and deterministic. In addition, ASEMG does not contain periodic or quasi-periodic characteristic (see Fig. 3), so other form of nonlinear dynamics (such as limit cycles or invariant loops) also does not exist. Then the above results indicate that the deterministic nonlinearity of ASEMG may be chaotic. In the following, we use correlation dimension to further test its deterministic nonlinearity. It is known that the fractal dimension of a strange attractor is very sensitive to noise, so it not only can detect whether ASEMG is linear or nonlinear, but also can test whether it is deterministic or not. Here, we apply correlation dimension to study whether the nonlinearity of ASEMG is chaotic or not. Fig. 7 gives correlation integral curves ln C(r)–ln r of the ASEMG data with respect to different embedding dimension d. One can see that they approximately parallel with the increase of embedding dimension and keep it invariant, the corresponding local slope is about 3.8050 ± 0.056 (see Table 1). The results show ASEMG have fractal feature and chaotic dynamic. From Table 1, one can see that an embedding dimension of reconstruction phase space should be at least 6. In other words, the high-dimension reconstruction attractor could reflect the main characteristics of the origin system. This result shows that it is difficult to classify ASEMG of different movement for prosthesis control by only using a single nonlinear parameter. Fig. 8(a) gives the correlation dimension curves of the ASEMG data and surrogates (39 sets of surrogate

301

Fig. 7. The correlation integral curves of ASEMG data.

data generated according to the null hypothesis 1), m = 2–8. Fig. 8(b) gives the correlation dimension values of the data and surrogates, m = 6. From the two figures, we can see that there is the obvious difference of correlation dimensions between the data and surrogates. It shows the null hypothesis 1 is rejected in 95%, and indicates the deterministic nonlinearity of ASEMG is chaotic. Here it is necessary to say that for the calculated dimension D2 by using correlation integral curves, the choice of linear parts of ln C(r)– ln r and their length have a great effect on calculated results. In order to get rapidly the D2 values on the various embedding dimension conditions, this Letter chooses the same linear part for all correlation integral curves. Though these values are not very exact in comparison with the real values, this test algorithm is great effective in nonlinear test methods. It is necessary to say that we can get similar results not only from ASEMG of the above subject but also from ASEMG of other two subjects. This proves the methods used in this Letter are suit to analysis ASEMG.

5. Conclusions In summary, we have shown that nonlinear dynamic characteristic can be detected from small data sets of different action surface EMG signals. To our knowledge, we are the first to report that surface EMG signals, which reflect certain limb movement process, contain deterministic nonlinear components (such as chaos) and high-dimension dynamics. So using single dimension is difficult to describe the natural feature of ASEMG, this is consistent with what Wang et al. reported [25]. In comparison with those previous

302

M. Lei et al. / Physics Letters A 290 (2001) 297–303

Table 1 The analysis of correlation dimension of surface EMG signal during movement d

2

3

4

5

6

7

8

9

10

11

D2

1.9315

2.7037

3.2808

3.4837

3.8047

3.7597

3.8511

3.8864

3.8129

3.7160

(a)

(b) Fig. 8. The test analysis based on correlation dimension for ASEMG data, where the surrogate data are generated by the null hypothesis 1, D2 is calculated by Eqs. (3). (a) The curves of surface EMG and surrogate data with m (abscissa is m = 2–8, ordinate id D2 ). (b) D2 distribution of surrogate data. ∗ is D2 of ASEMG in m = 6 (abscissa is D2 , ordinate is histogram).

works on analysis of movement information from real surface EMG signals [7–12,25], where much room of controversy was left, because nonlinear dynamics has not yet been conformed and the existence of chaos still remains in the phase of speculation, here, nonlinear detection analysis of the experimental surface EMG data first clearly shows that ASEMG signal contains intrinsic deterministic nonlinearity, and even chaotic firings in real muscle contraction experiments by the correlation dimension analysis and spectral analysis.

This can serve as another convincing evidence that chaos does exist in irregular firing activities of neurons [27]. The basic idea of the null hypothesis 1 is that generated surrogate data have the same linear correlation as the raw data by reconstructing spectral of the original data. If the used algorithm cannot generate the above surrogates very well, the generated surrogates will be difficult to reflect the linear properties of the raw, then the result will be inconvincible. So, this Letter gives the improved FT algorithm that can satisfy the requirement of the null hypothesis 1. And it proves effective by using it to detect the nonlinearity of ASEMG. It is also worthwhile noting to choose various kinds of the test statistics in nonlinear test analysis. For the T values of Eqs. (1) and (2), they can be fast calculated on the analysis of 39 surrogates and real EMG data. Therefore one can conveniently utilize them to detect the nonlinearity of time series. However, for value T of Eqs. (3) (correlation dimension), the speed of its calculation is rather slow although this method is very effective. So, if it is effective to use the T values of Eqs. (1) and/or (2) to test nonlinearity of time series, it is not necessary to apply the T value of Eqs. (3). Finally, to further determine ASEMG signal as a chaotic signal, its largest Lyapunov exponent L should also be computed, because, if L > 0, there is chaotic dynamic character in the system. We apply a Lyapunov exponent algorithm, which was designed by Rosenstein et al. [28] for small data set, to calculate the most Lyapunov exponent of ASEMG at d = 6. The calculated L values are different for different ASEMG signal but always positive. These results show that ASEMG has high-dimension chaotic dynamics, like correlation dimension analysis. With the high-dimension nonlinear dynamics, it is difficult to classify ASEMG by movement kinds for prosthesis control by only using a single nonlinear parameter. Some multiparameter methods maybe more suit to describe the movement kinds of ASEMG, such as multi-

M. Lei et al. / Physics Letters A 290 (2001) 297–303

fractal dimension. In Ref. [29], we apply multifractal dimension to classify ASEMG into four kinds of action and obtain better recognition rate than the previous.

Acknowledgements The present Letter is revised from part of the Ph.D. thesis of the first author [29]. This work is supported by National Natural Science Foundation Grant No. 69675002 of PR China.

References [1] C.J. De Luca, Crit. Rev. Biomed. Eng. 11 (1984) 251. [2] L. Brody, M. Pollock, S. Roy, C. De Luca, B. Celli, J. Appl. Physiol. 71 (1991) 1878. [3] D. Graupe, IEEE Trans. Biomed. Eng. 36 (1989) 711. [4] V. Dimanico, P. Filippi, M. Barra, P. Bertone, F. Laterza, R. Merletti, in: 11th Congress of ISEK, Enschede, The Netherlands, 1996, pp. 146–147. [5] T. Sadoyama, T. Masuda, H. Miyata, S. Katsuta, Eur. J. Appl. Physiol. 57 (1988) 767. [6] R. Merletti, M. Knaflitz, C.J. De Luca, J. Appl. Physiol. 69 (1990) 1810. [7] Z. Yang, G. Zhao, Acta Biophys. Sinica 14 (2) (1998) 257. [8] J. Yang, B. Liu, J. Peng, Z. Ma, Space Med. Med. Engrg. 12 (3) (1999) 185.

303

[9] V. Gupta, S. Suryanarayanan, N.P. Reddy, in: IEEE 17th Annual Conference Engineering in Medicine and Biology Society, Vol. 2, 1995, p. 1331. [10] Anmuth, G. Goldberg, N.H. Mayer, Muscle Nerve 17 (1994) 953. [11] Gitter, J.M. Czerniecki, D. DeGroot, Muscle Nerve 14 (1991) 884. [12] Gitter, J.M. Czerniecki, J. Neurosci. Methods 58 (1995) 103. [13] D. Graupe, W.K. Cline, IEEE Trans. Syst. Man Cyber. 5 (1975) 252. [14] G.N. Saridis, T.P. Gootee, IEEE Trans. Biomed. Eng. 29 (1982) 403. [15] W.-J. Kang, J.-R. Shiu, C.-K. Cheng, J.-S. Lai, H.-W. Tsao, T.-S. Kuo, IEEE Trans. Biomed. Eng. 42 (8) (1995) 777. [16] B. Hudgins, P. Parker, R.N. Scott, IEEE Trans. Biomed. Eng. 40 (1993) 82. [17] M. Barahona et al., Nature 381 (5) (1996) 215. [18] U. Parlitz et al., Int. J. Bifurc. Chaos 7 (2) (1997) 407. [19] J. Theiler et al., Phys. Lett. A 196 (1995) 335. [20] J. Theiler et al., Physica D 58 (1992) 77. [21] J. Theiler et al., Physica D 94 (1996) 221. [22] P. Grassberger, I. Procaccia, Physica D 9 (1983) 189. [23] M. Small, K. Judd, Int. J. Bifurc. Chaos 8 (6) (1998) 1231. [24] W.A. Brock, J. Lakonishok, B. LeBaron, J. Finance, to appear. [25] R. Wang, C. Huang, Y. Chang, N. Yang, Chinese J. Med. Instrum. 23 (3) (1999) 125, 138. [26] T. Schreiber, A. Schmitz, Phys. Rev. Lett. 77 (4) (1996) 635. [27] Y.F. Gong, Ph.D. thesis, Xi’an Jiaotong University (1998), in Chinese. [28] M. Rosenstein, J. Collins, C.D. Luca, Physica D 65 (1993) 117. [29] L. Min, Ph.D. thesis, Shanghai Jiaotong University (2000), in Chinese.