Bispectrum-based features classification for myoelectric control

Biomedical Signal Processing and Control 8 (2013) 153–168 Contents lists available at SciVerse ScienceDirect Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc Bispectrum-based features classification for myoelectric control Eugenio C. Orosco a,∗ , Natalia M. Lopez b , Fernando di Sciascio a a b Instituto de Automática, Facultad de Ingeniería, Universidad Nacional de San Juan, Argentina Gabinete de Tecnología Médica, Facultad de Ingeniería, Universidad Nacional de San Juan, Argentina a r t i c l e i n f o Article history: Received 20 June 2012 Received in revised form 13 August 2012 Accepted 21 August 2012 Available online 7 September 2012 Keywords: EMG Robust bispectrum Continuous classification Myoelectric control a b s t r a c t Surface electromyographic signals provide useful information about motion intentionality. Therefore, they are a suitable reference signal for control purposes. A continuous classification scheme of five upper limb movements applied to a myoelectric control of a robotic arm is presented. This classification is based on features extracted from the bispectrum of four EMG signal channels. Among several bispectrum estimators, this paper is focused on arithmetic mean, median, and trimmed mean estimators, and their ensemble average versions. All bispectrum estimators have been evaluated in terms of accuracy, robustness against outliers, and computational time. The median bispectrum estimator shows low variance and high robustness properties. Two feature reduction methods for the complex bispectrum matrix are proposed. The first one estimates the three classic means (arithmetic, harmonic, and geometric means) from the module of the bispectrum matrix, and the second one estimates the same three means from the square of the real part of the bispectrum matrix. A two-layer feedforward network for movement’s classification and a dedicated system to achieve the myoelectric control of a robotic arm were used. It was found that the classification performance in real-time is similar to those obtained off-line by other authors, and that all volunteers in the practical application successfully completed the control task. © 2012 Elsevier Ltd. All rights reserved. 1. Introduction Surface electromyographic signal (EMG) is an electric manifestation of muscular activity associated with muscle contraction and force level, which can be detected on the skin surface [1,2]. This complex signal is influenced by several physiological and anatomical factors, as well as for recording situations. The EMG signal contains rich information about muscular activity, neuromuscular disorders, and motion intentionality, and has the potential to become a reference input signal to myoelectric control systems. For this purpose, signal processing is required for conditioning the EMG in two main steps, namely, feature extraction and classification. The aim of feature extraction techniques is to reduce the dimensionality of the digitalized EMG signals by mapping these large-dimension objects into a smaller-dimensional space. This is done by selecting from the original EMG signal a small set of features (or characteristics) that are arranged in a feature vector. Ideally, these features should be simple to extract, invariant to irrelevant transformations, insensitive to noise, and useful for characterizing and discriminating motion patterns. The information in these features represents the motion intention of the user. On this ∗ Corresponding author. Tel.: +54 02644286993; fax: +54 02644213672. E-mail addresses: [email protected] (E.C. Orosco), [email protected] (N.M. Lopez), [email protected] (F. di Sciascio). 1746-8094/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.bspc.2012.08.008 accounts, the correct selection and design of the feature vector is more critical than the choice of the classification method [3]. The abstraction provided by the feature-vector representation of the EMG signals allows using of a well-developed domain independent classification theory. A wide variety of different features for EMG classification purposes have been extensively reported in the literature [4]. Features can be classified according to several criteria, e.g., their specific feature domain (time [5], frequency [6], or time-frequency domain [7]); the linear or nonlinear nature of the feature extraction mapping [8,9]; or whether the feature extraction method is based on second or high-order statistics. High-order statistics (HOS), and the probabilistic models based on them, allow modelling non-Gaussian and non-linear signals, as pointed out by Nikias, Mendel and Swami [10–12] and in more extensive theoretical works [13–17]. In past years, several authors have proposed to use of HOS concepts to analyze and classify the EMG signals, e.g., high-order cumulant sequences, high-order spectra or polyspectra, skewness, kurtosis, and other concepts [18–24]. HOS are more extensively used in EMG analysis, to explore the relationships between the bicoherence (normalized bispectrum) and the force level [20], to estimate the amplitude and the number of newly motor unit action potential using the second and the fourth-order moments [21] and using the third-order spectrum [22–24]. However, few researchers have addressed the problem of HOS-based feature classification. For example, Nazarpour et al. [18] employ HOS of EMG signal to classify four motions, based on the 154 E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 fact that useful information can be extracted from the second, third and fourth-order cumulants using a K-nearest neighbor classifier. Another bispectrum feature extraction method can be found in [19] that performs an off-line classification of EMG signals, using segmented matrix integration (axial and radial direction) and Fisher’s linear discriminant analysis. They obtain the best classification rate off-line, without considering the transition between movements; although they suggest the possibility of online implementation. In this paper, a bispectrum (the double Fourier transform of the third-order cumulant sequences) from multiple EMG channels is proposed as a HOS-based feature classification method applied to the myoelectric control of a robotic device, as a hard real-time application. The bispectrum is estimated from four EMG channels, corresponding to biceps brachii, triceps brachii, pronator and brachioradialis. The bispectrum is a bidimensional complex function represented by a complex matrix. Therefore, in order to perform the classification, it is necessary to reduce this 2D complex array into a real feature vector. To this aim, two feature vectors were proposed and analyzed under criterions of complexity, computational cost and feasibility of real-time implementation. The most adequate feature vector is selected to be used as an input to an artificial neural network to classify five movements of the upper limb: elbow’s flexion–extension, forearm’s pronation–supination and arm’s rest. With these results, a robotic arm is controlled in real time. The paper is organized as follows: In Section 2, we explain the experimental procedure for EMG signal acquisition, processing, classification, and the myoelectric control of a robot arm. We also provide a brief overview of data segmentation, surface electromyogram models, bispectrum definition and properties, and bispectrum processing (estimation methods, feature extraction and reduction). In Section 3, we present and discuss our results in relation to: (i) the statistical analysis of the proposed robust bispectrum estimators, (ii) the classification performance of the bispectrumbased feature extraction, and (iii) the real-time implementation of the myoelectric control of a robot arm. Finally, Section 4 states the conclusions. interest in this work are flexion, extension, pronation, supination, and rest position. The volunteers are encouraged to perform the sequence of movements prompted by beepers. The database consists of the records of six volunteers, three male and three females, volunteers aged 28.6 ± 5.4, their height 1.77 ± 0.13 m and weight 66.2 ± 12.7 kg. All volunteers are healthy, with no history of muscle weakness, neurological diseases or drug therapy. One of them has a congenital malformation, i.e., unilateral phocomelia below his elbow (Vol. 1). The volunteers have approved and signed an informed consent form, according to the experiments to be performed. The database was collected in four experimental trials, completed on four different days, while avoiding the effects of muscular fatigue [25]. Each trial consisted of five series of sequential movements, holding each motion for three seconds and interspersing the rest position. Each trial has 100 s maximum duration. The first trial was used for data training for artificial neural networks. The remaining sessions were used for data validation. The volunteers performed an additional session for real-time myoelectric control of a robotic arm in. 2.2. Front-end processing and EMG acquisition The myoelectric control system based on features extracted from the EMG signal bispectrum is described in this section. The block diagram of the overall myoelectric control scheme is shown in Fig. 1. Following the sequence outlined in the block diagram, this section is organized as follows: First, the experimental protocol, the front-end processing and EMG acquisition, and data segmentation are presented. Then, the conceptual subsections explain the relevant concerns about the models of the surface EMG, bispectrum, bispectrum estimation and, finally, feature extraction, feature reduction, the classification and real time myoelectric control methods. Two data acquisition systems were used, both of which featuring a front end processing stage and an EMG acquisition stage. The front end processing stage includes instrumentation amplifiers (Gain of 1000 (V/V) and CMRR > 100 dB), band pass filters (10–500 Hz) and opto-isolated subsystems [1,2]. The EMG acquisition stage has an analog-to-digital converter and the digital processing software. The EMG signals from four channels were registered by the two different data acquisition systems with a sampling frequency fixed at 1 kHz. The reasons for using these two systems were; first, the database needs a precise and accurate high-resolution data acquisition system, because the signals are used for features analysis and neural network training. On the other hand, the real time control requires a fast hardware and a fast communication link in order to effectively control these devices. The commercial data acquisition system, used for EMG database, is a 15LT – Grass Technologies® . Signals are sampled with a 16-bit A/D converter (National Instruments DAQPad-6015® ), and Matlab® software. The real time control was developed with a data acquisition system designed in our institute. The design of the front end processing stage is explained with details in [26,27]. The EMG acquisition stage has a PIC16F876 microcontroller which samples the signals through a 10-bit A/D converter. The QNX® real time operating system records the signals and it executes the processing, classification and control algorithm tasks. Also, in this stage, the EMG data from the database is statistically normalized in order to have zero mean and unit variance EMG signals. 2.1. Experimental protocol 2.3. Data segmentation A movement-based protocol is designed for recording the upper limb EMG signals for biceps brachii, triceps brachii, pronator and brachioradialis muscles, as is shown in Fig. 2. The movements of A time segment (time-slot or time-window) is a time interval used to acquire myoelectric data needed for some kind of signal processing, such as detection, estimation, filtering, feature 2. Materials and methods Fig. 1. Block diagram of the overall myoelectric control system. E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 155 Fig. 2. Electrodes placement used with the experiments. The considered muscles are biceps and triceps brachii, and pronator and brachioradialis. The movements are flexion, extension, pronation, supination, and rest position. extraction and classification (as in our case), and the like. The idea behind data segmentation is that the signal can be considered quasi-stationary during these time intervals; then, the data can be batch-wise processed from segment to segment. This signal processing mode is also known as real-time data processing over sliding time windows (with or without overlapping). Two important parameters must be properly selected in order to design a good data segmentation strategy; namely, the segment length, and the amount of overlapping between two consecutives segments. The segment length (measured in the number of time samples), determines the trade-off balance between the accuracy of the classification and the response time. This is so because, on the one hand, the segment length should be large enough to avoid a degradation of the classification performance, because the bias and variance of feature estimators increase as segment length decreases and, on the other hand, the segment length should be small enough to satisfy the hard real-time constraints of the myoelectric control application [4,28]. The myoelectric control must supply the control commands in less than 300 ms, so that the user does not feel any delay in the expected response from the system. In this work a segment length of 256 samples was adopted. The amount of overlapping is measured as a percentage of the segment length. The degree of overlapping between two consecutives segments normally ranges from 0%, for non-overlapped or adjacent windowing, and almost 100%. This parameter affects the rate with which class decisions will be made and can be adjusted to optimize the use of the computing resources. Here, a 50% overlapping (128 samples) is adopted to obtain a continuous classification scheme (see Fig. 3). The post-processing technique is explained in Section 2.7. 2.4. Models of the surface electromyogram The standard assumption about the surface EMG signal recorded during voluntary constant-force contractions is that it can be modeled as a band-limited, correlation-ergodic, centered (zeromean) Gaussian process, modulated by muscle activity and corrupted by additive Gaussian white noise [29]. The sum of these two independent and centered Gaussian processes is a centered Gaussian process. This assumption provides an adequate mathematical tool for modelling and feature extraction, with acceptable results. Nevertheless, some authors have proposed other representations of the EMG signal, such as Laplacian process [30] or Gaussian Mixture Models [31]. Based on the implicit assumption of ergodicity, several authors have investigated the sample or empirical marginal probability density function (pdf) of EMG signal amplitudes [32–36]. It is widely accepted that the Gaussian amplitude distribution is more adequate to represent voluntary contractions above 50% of maximum voluntary contraction (MVC), but the amplitude distribution tends towards a Laplacian or a double exponential distribution as the force level decreases (specially below 30% of maximum voluntary contraction) [33,35,36]. 2.5. Bispectrum In the following three subsections, the definitions of cumulant sequences and bispectrum functions are presented and also the influence of the symmetries of the joint and marginal density functions of a random signal on their bispectrum are discussed. Finally, E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 156 Fig. 3. The EMG channel is segmented into sliding windows of 256 samples with overlapping of 128 samples. Tree segments are shown with 50% overlapped. This configuration is implemented to do a continuous feature extraction. in Section 2.5.3, we explain several bispectrum estimation methods. 2.5.1. Bispectrum definition and basic properties The bispectrum is a third-order frequency-domain measurement containing information that conventional spectral analysis techniques cannot provide [10–12,37]. This is especially useful for non-Gaussian signals. Its definition and some basic properties are recalled in the following. Let {X(k)} , k = 0, 1, 2, 3, . . . a real stationary discrete time random process whose first n-moment exists. Then the n-moment function [37] can be defined as mxn (1 , 2 , . . . , n−1 ) = E{X(k)X(k + 1 )· · ·X(k + n−1 )} (1) where E {·} is the statistical expectation and  i is time shift. The most widespread moment functions used are: the second moment or autocorrelation function mx2 (), the third-order moment functions mx3 (1 , 2 ) and the fourth-order moment function mx4 (1 , 2 , 3 ). For a zero mean stationary random process and only for n = 3, 4, the definition of cumulant functions is written as g cnx (1 , 2 , . . . , n−1 ) = mxn (1 , 2 , . . . , n−1 ) − mn (1 , 2 , . . . , n−1 ) cnx (1 , 2 , . . . , n−1 ) = E   X(k)X(k + 1 )· · ·X(k + n−1 ) − E(g(k)g(k + 1 )· · ·g(k + n−1 )) (2) where {g(k)} is a Gaussian random process with the same autocorrelation function as {X(k)}. For a given random process, the cumulant functions provide a measure of the distance to Gaussianity. If {X(k)} is a Gaussian process, cumulants higher than second order are zero. n = 1, 2 and 3 For between moments and given by the general relationship cumulant functions are c1x = mx1 = E{X(k)} (3a) 2 2 c2x (1 ) = mx2 (1 ) − (mx1 ) = mx2 (−1 ) − (mx1 ) = c2x (−1 ) (3b) c3x (1 , 2 ) = mx3 (1 , 2 ) − mx1 [mx2 (1 ) + mx2 (2 ) + mx2 (1 − 2 )] − 2(mx1 ) 3 (3c) where, c1x is the mean value of {X(k)} and c2x is the autocovariance sequence. The third-order cumulant function is c3x (1 , 2 ) and, when the random sequence {X(k)} is zero mean, the cumulant function is equivalent to the moment function, that is, mx3 (1 , 2 ) = c3x (1 , 2 ) = E{X(k)X(k + 1 )X(k + 2 )} (3d) Important properties justify the use of cumulants instead of moments [10,11]. First, high-order cumulant functions extract information deriving from deviations from Gaussianity. Second, the polyspectra of the white noise is a multidimensional flat function. Third, the cumulant functions of two statistically independent random processes are equal to the sum of the individual cumulant functions of each process. Therefore, the cumulant functions can be treated as an operator and, consequently, the discrete time Fourier transform is applicable. From this point, this work will focus on the third-order cumulant sequence and assume that c3x (1 , 2 ) is absolutely summable [37], the third-order spectra or the bispectrum can be defined as the E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 two-dimensional discrete time Fourier transform of the third-order cumulant c3x (1 , 2 ) as follows. B3x (ω1 , ω2 ) = ∞ ∞   c3x (1 , 2 ) exp{−j(1 ω1 + 2 ω2 )} (4) 1 =−∞2 =−∞ where |ω1 | ≤  ; |ω2 | ≤  ; |ω1 + ω2 | ≤  and c3x (1 , 2 ) is the thirdorder cumulant function of the random sequence {X(k)}. The properties of symmetries in third-order cumulants function and the corresponding bispectrum function can be found in [10,37], in example; c3x (1 , 2 ) = c3x (2 , 1 ) = c3x (−2 , 1 − 2 ) = c3x (2 − 1 , −1 ) = c3x (1 − 2 , −2 ) = c3x (1 , 2 − 1 ) (5) third-order cumulant estimates. It is also well known that for a fixed number of data samples, the variance of these estimates is greater than the second-order cumulant estimates (autocorrelation function). In the following subsections, three third-order cumulant estimators are presented, based on different approximations of the expectation operator E(·) in Eq. (3d): the arithmetic mean, the median or second quartile, and the truncated or trimmed mean. 2.5.3.1. Third-order cumulant estimators. Arithmetic mean thirdorder cumulant estimator: The standard unbiased estimator of the third-order cumulant sequence is equivalent to using the sample mean estimator, and; ĉ3x−MEAN (m, n) = B3x (ω1 , ω2 ) = B3x (ω2 , ω1 ) = B3x (−ω2 , −ω1 ) = B3x (−ω1 − ω2 , ω2 ) = B3x (ω1 , −ω1 − ω2 ) = B3x (−ω1 − ω2 , ω1 ) = B3x (ω2 , −ω1 − ω2 ) (6) 2.5.2. Bispectrum and symmetry of probability density functions The bispectrum of a random process is identically zero when their joint probability density functions are symmetric. Only for the special case where the process is white (i.e., an independent and identically distributed sequence), the bispectrum vanishes when their marginal probability density functions are symmetric [38,39]. The following equalities between joint PDF’s must be satisfied in order to verify that a zero-mean, stationary random process X(k) has identically zero third-order statistics [39]: 1 N − max(m, n) N−max(m,n)  x(l)x(l + m)x(l + n) (8) l=1 where, in our case XN = {x(k), k = 1, 2, . . ., N} is the EMG sequences (time segment) of N = 28 = 256 samples, l = 1, 2, . . . is the index of the sample segment, and m = 0, 1, 2, . . . , n = 0, 1, 2, . . . are the time lags. Evidently, between l, m, n and N it must be satisfied the relation 1 ≤ max(l + m, l + n) ≤ N. Therefore index l varies from 1 to N − max(m, n). A slight modification of Eq. (8) leads to a biased (but asymptotically unbiased) estimator of the third-order cumulant sequence. ĉ3x−MEAN-B (m, n) = 1 N N−max(m,n)  x(l)x(l + m)x(l + n) (9) l=1 Eqs. (8) and (9) can be written in a compact form as: N−m fX(k),X(l),X(m) (xk , xl , xm ) = fX(k),X(l),X(m) (xk , xl , −xm ) ĉ3x−MEAN (m, n) = = fX(k),X(l),X(m) (xk , −xl , xm ) = fX(k),X(l),X(m) (−xk , xl , xm ) 157 1 zm,n (l) K (10) l=1 (7) where k, l, and m are distinct integers. The above equalities are much more restrictive symmetry constraints than the one expressed in terms of just the marginal PDF. Furthermore, it would be more difficult to verify the symmetry conditions in terms of the joint PDF’s, than to directly investigate when the bispectrum is zero (or nonzero). Because it is well-known that the EMG signal is not white, the preceding paragraphs mean that the symmetry of the marginal (amplitude) probability distribution of the EMG signals by itself is not sufficient to conclude that the third-order statistics are zero. This clearly justifies the use of EMG bispectrum features for movement classification task. 2.5.3. Bispectrum estimation methods The problem of estimating the bispectrum from a single finite time realization of a stochastic process has been largely addressed in the literature. There are basically two approaches that have been used to estimate the bispectrum, namely, the parametric approach (which is based on linear parametric models, e.g., autoregressive, moving average, and ARMA models), and the conventional or ‘Fourier type’ approach [37]. A review of the commonly used techniques can be found in [40–44], and references therein. Considering the real-time application of myoelectric control, this paper is focused in three estimators, members of the family of indirect class of conventional bispectrum estimators, that is, estimators based on the approximation of bispectrum definition given by (3d) and (4). It follows from this definition that the accuracy of bispectrum estimation depends mainly on the accuracy of the where zm,n (l) = x(l)x(l + m)x(l + n) is defined for fixed values of m and n over the triangular region m = 0, 1, 2, . . . and 0 ≤ n ≤ m ≤ N − 1, the index l ranges from 1 to N − m, and K takes the values N − m or N for unbiased or asymptotically unbiased estimator respectively. In what follows, only unbiased estimators are considered, because our estimation task is actually a small sample estimation problem, and, as G.W. Brown pointed out as early as 1947 [45] “To ensure, in small-sample estimation, that an estimate bears some relation to the parameter it is estimating, it has become the custom to require that an estimate be unbiased”. Median third-order cumulant estimator: The sample median is a classical method for estimating the middle of a data set. This estimator provides a robust alternative to sample mean when the underlying distribution of zm,n (l) is approximately Laplacian or double exponential, which is indeed our case. Fig. 4 shows the plots of nearly Laplacian empirical density functions for several zm,n (l). In this case the cumulant estimates is bias robust in the presence of outliers, and with less variance than the sample mean estimator. The sample median third-order cumulant estimator is represented compactly by: ĉ3x−MED (m, n) = Med(zm,n (l)) (11) where Med(·) is the sample median operator and, likewise above, zm,n (l) = x(l)x(l + m)x(l + n) is defined for fixed values of m and n over the triangular region m = 0, 1, 2, . . ., 0 ≤ n ≤ m ≤ N − 1, and 1 ≤ l ≤ N − m. If the data zm,n (1), . . ., zm,n (N − m) are arranged in ascending order from smallest to largest and written as zm,n [1] ≤ · · · ≤ zm,n [N − m], then zm,n [i] is the ith order statistic and E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 158 Fig. 4. The estimated probabilistic density functions of the several EMG windows of 256 ms are near a Laplacian or double exponential density distribution. The estimated probabilistic density functions are based on a normal kernel function. the values zm,n [i] are called the order statistics. Using this notation, the sample median is defined as: ĉ3x−MED (m, n) = Med(zm,n (l)) = ⎧ N −m+1 ⎨ zm,n 2 ĉ3x−MED (m, n). if N − m is odd ⎩ zm,n [(N − m)/2] + zm,n [[(N − m)/2] + 1] 2 if N − m is even (12) Symmetric trimmed mean third-order cumulant estimator: The truncated or trimmed mean estimator given by (13) is another approach to robust estimation of third-order cumulant. The sensitivity to outliers is reduced by removing the smallest ks and largest kl values from the data set. ĉ3x−TRIMM (m, n, ks , kl ) = 1 N − m − ks − kl N−m−kl  zm,n [l] (13) l=ks +1 where zm,n [1] ≤ · · · ≤ zm,n [N − m] are the previously defined order statistics. When ks = kl = k, a symmetric k-trimmed mean estimator is obtained, given by (14). ĉ3x−TRIMM (m, n, k) = 1 N − m − 2k N−m−k  zm,n [l] (14) = ⎧ N−m−1 ⎪ ⎨ if N − m is odd 2 ⎪ ⎩ N − m − 1 if N − m is even 2 2.5.3.2. Ensemble average third-order cumulant estimators. Ensemble (or time segments) averaging is a form of data stratification. This is a standard statistical technique used in the analysis of stationary time series in order to obtain an acceptable trade-off between the bias of the estimator and its variance. In averaging over time segments, the finite sequence XN = {x(k), k = 1, 2, . . ., N} is divided into D non-overlapping subsequences, each one having a length of M samples, such that N = D × M. If x(i) (j), j = 1, 2, . . ., M, i = 1, 2, . . ., D denote the k = (i − 1)M + jth observation in the original sequence, then the subsequences X(i) are defined as, X(i) = {x(i) (j) = x((i − 1)M + j), j = 1, 2, . . ., M} , i = 1, 2, . . ., D. Thus, the following ensemble average versions of the three estimators defined previously are calculated as the ensemble average of the cumulant subsequences estimates on the segments. Ensemble average arithmetic mean third-order cumulant estimator: l=k+1 where k can take the values k = 0, 1, . . . , (N − m − 1)/2 , and the symbol ⌊ ⌋ denotes the floor or the integer part function that is given by (16). N−m−1 2 It should be noted that the unbiased sample mean estimator (8) is a special case of the symmetric k-trimmed mean estimator (14) when k = 0, i.e., ĉ3x−TRIMM (m, n, k = 0) = ĉ3x−MEAN (m, n); and that the sample median estimator (12) is another  special case of (14), when  k = (N − m − 1)/2 , i.e., ĉ3x−TRIMM m, n, k = (N − m − 1)/2 = ĉ3x−MEAN-D (m, n) = (i) (15) D M−m i=1 l=1 1  1  (i) zm,n (l) D M−m (16) zm,n (l) = x(i) (l)x(i) (l + m)x(i) (l + n), i = 1, 2, . . . , D are where, defined in each subsequence X(i) for fixed values of m and n over the triangular region m = 0, 1, 2, . . ., 0 ≤ n ≤ m ≤ M − 1, the index l varies from 1 to M − m. E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 Ensemble averaging median third-order cumulant estimator: ĉ3x−MED-D (m, n) = 1 D D   intention of the user. In the following sections, the procedures used for EMG bispectrum feature extraction and feature reduction will be explained.  Med zm,n (l) i=1 = (i) ⎧ D  M −m+1 ⎪ 1 ⎪ (i) ⎪ zm,n ⎪ 2 ⎨D if M − m is odd i=1 D (17) (i) (i)  zm,n ⎪ [(M − m)/2] + zm,n [[(M − m)/2] + 1] 1 ⎪ ⎪ ⎪ ⎩D 2 if M − m is even i=1 (i) where the values zm,n [·] are the order statistics defined in the subsequence X(i) , i = 1, 2, . . ., D. Ensemble averaging symmetric trimmed mean third-order cumulant estimator: D ĉ3x−TRIMM-D (m, n, k) = 1 1 D M − m − 2k i=1 M−m−k  (i) zm,n [l] (18) l=k+1 2.5.3.3. Bispectrum estimates. As it was pointed out previously, the bispectrum estimation is an approximation of the definition (4), calculated as: B̂3x−XX (ω1 , ω2 ) = 159 L L   ĉ3x−XX (m, n)w(m, n) m=−Ln=−L × exp{−j(ω1 m + ω2 n)} (19) where the super index XX in the third-order cumulant and bispectrum estimates can be: MEAN, MED, TRIMM, or the ensemble average versions MEAN-D, MED-D, TRIMM-D; and L ≤ N − 1 or L ≤ M − 1, as appropriate. The cumulant symmetries and the properties of its two-dimensional fast Fourier transform should be considered in order to reduce the computational cost [10,11,37]. As in the case of conventional power spectrum estimation, a bidimensional window function w(m,n) is used to find the smooth bispectrum estimates. Furthermore, it is well known that the selection of the window function w(m,n) affect the bias-variance trade-off in the bispectrum estimation. The bidimensional window function w(m,n) must satisfy the following conditions [37]: • Symmetry properties: w(m, n) = w(n, m) = w(−m, n − m) = w(m − n, −n) • Zero outside the region high-order statistics estimates, i.e., w(m, n) = 0; |m| > L • Normalizing condition at the origin, i.e., w(0, 0) = 1 • Nonnegative Fourier transform, i.e., W(ω1 , ω2 ) ≥ 0, for all (ω1 , ω2 ) These properties can be satisfied by generating the bivariate function w(m,n) using some standard univariate window functions, such as Daniell, Hamming, Parzen, Priestley or Sasaki windows. The window functions can be evaluated in terms of two quantities that are proportional to the bias and the variance of the estimated bispectrum. This issue is fully explained and analyzed in Nikias and Petropulu [37]. Among other things, the authors concluded that the bispectrum estimated with Sasaki window has the lowest bias, and that of Parzen window has the lowest variance. According to this, the Parzen window function was chosen. 2.6. Feature extraction The goal of feature extraction is to find a small set of characteristics, grouped into a feature vector of the digitalized EMG signals that is useful for characterizing and discriminating the motion patterns. The information contained in these features represents the motion 2.6.1. Bispectrum feature extraction The performance of the bispectrum estimations in real time situations determines the most adequate estimator. For analysis the estimators explained in Section 2.5.3 were considered. In this feature extraction stage, the robust bispectrum median estimator using (17) is implemented. As is explained later, the median estimator produces the lowest variance estimation and it is robust against outliers as well, which justifies the decision for choosing this estimator. The EMG data segmentation allows the continuous classification of movements in sliding time windows. Those windows have 256 samples with overlapping on 128 samples (50%). The segments are assumed as a stationary stochastic process and it also provides the suitable dimension, thus making faster the two-dimensional fast Fourier transform. The bispectrum function was calculated in each of the four EMG channels over the sliding time windows. In Eq. (17), the EMG segment x(k) was divided into four groups of 64 samples (D = 4 and M = 64), and then, the estimation of the third-order cumulant sequences was calculated as the average of estimates on the segments. Then, the discrete time double fast Fourier transform was calculated using the cumulant sequences and their symmetries. The bispectrum estimation results in a complex square matrix of order 128 which contains several regions of symmetry. The fist triangular region was used to compute the features vector. A schematic explanation is shown in Fig. 5. 2.6.2. Bispectrum feature reduction The bispectrum is a bidimensional complex function that is represented by a complex-valued square matrix of size 2M × 2M; where M is the size of the fist quadrant redundancy region. In this work, M = 64 for each of the four channels; then, the matrix size is 128 × 128. In order to perform the classification, it is not suitable to use the whole bispectrum matrix due to its high dimensionality. Therefore, it is necessary to reduce this 2D complex-valued array into a real-valued feature vector without losing much information along the process. This is possible because of the redundant information present in the bispectrum. To this purpose, various authors have proposed different reduction techniques, by taking into consideration several factors, such as complexity, computational cost and feasibility of real-time implementation. Kaplanis et al. [20] have examined seven parameters of the bispectrum: peak amplitude in the x-direction; slope peak slope in the x-direction; peak amplitude in the y-direction; slope peak slope in the y-direction; Maximum amplitude; Test for Gaussianity; Test for Linearity. Similarly, Hussain et al. [46] have investigated the Test for Gaussianity based on the mean bicoherence power (a squared power normalized version of the bispectrum). Chen et al. [19] used segmented matrix integration (axial and radial direction) and a non-linear logarithm transform of the bispectrum. For each EMG channel, two feature reduction methods for the complex bispectrum matrix are proposed. The first one estimates, as shown in Table 1, the three classic Pythagorean means from the module of the bispectrum matrix, i.e., arithmetic mean (AM), harmonic mean (HM), and geometric mean (GM). With these, the vector [AMs , HMs , GMs ]T is formed, where the subscript s = 1, 2, 3, 4 indicates the EMG channel number (see Fig. 6(a)). The final feature vector is obtained by applying a quarter-root transformation of each component of the vector [AMs , HMs , GMs ]T . The purpose of this last transformation is to obtain a more peaked feature’s probability distribution which has better numerical properties. 160 E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 Table 1 Feature vectors. Module feature vector [AMs , HMs , GMs ]T , s = 1, 2, . . ., 4 Arithmetic mean (AM) AMs = 1 M M  1 M−i+1 i=1 Harmonic mean (HM) HMs = M M   M  (M − i + 1) Real arithmetic mean (RM) RAMs = 1 M 1 M−i+1 RHMs = M  M  x , EMG subsegment. (i, j)| −1 −1 M  (s) MED-D ℜeal(B̂3x (i, j)) 2 j=i (M − i + 1)  M  j=i (s) MED-D ℜeal(B̂3x −1 −1 (i, j)) ⎛  1/(M−i+1) ⎞1/M M M   2 (s) ⎠ RGMs = ⎝ ℜeal(B̂3x MED-D (i, j)) i=1 (s) |B̂3 −1 j=i M  i=1 Real geometric mean (RGM) x(s) MED-D j=i i=1 Real harmonic mean (RHM)  M  ⎛  1/(M−i+1) ⎞1/M M M   (s) ⎠ GMs = ⎝ |B̂3x MED-D (i, j)| i=1 Square real part feature vector [RAMs , RHMs , RGMs ]T , s = 1, 2, . . ., 4 (i, j)| j=i i=1 Geometric mean (GM) (s) MED-D |B̂3x j=i Fig. 5. A 256 samples segment is used to estimate the third-order cumulant and, then, the bispectrum is obtained by the two-dimensional fast Fourier transform. The fist triangular redundancy region is chosen for the analysis of movement intention. Fig. 6. The feature vectors are defined by the quarter-root of the three Pythagorean means obtained from (a) the module, and (b) the square of the real part of the complex bispectrum matrix. E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 161 Fig. 7. The neural network is a two-layer feed-forward back propagation network. The artificial neural network input vector is based on the EMG feature vector, and the five elements output vector (M1–5) characterizes the movements (flexion, extension, pronation, supination and rest). A post processing is applied to the outputs. The second feature reduction method is summarized in Fig. 6(b) and Table 1. This method differs from the former one in that it uses the square of the real part of the complex bispectrum matrix instead of the module. In this case, the vector obtained before applying the quarter-root transformation is [RAMs , RHMs , RGMs ]T . This feature reduction method turns to be more effective, and faster than the first one. In both reduction methods, the size of the final feature vector is 12 (three features per channel). 2.7. Classification method It is widely recognized that the correct selection of components of the feature vector is more decisive than the choice of a specific classification method [3]. Thus, the selection of the classifier is not our primary goal in this paper. Several classifiers have been used for EMG classification [8,9,47]. For example, different Artificial Neural Networks, Support Vector Machine, multiple classifiers with competence function, linear Bayesian classifier, and so on. In this paper, a two-layer feedforward neural network (Fig. 7) was used, because of, among other things, its great simplicity, straightforward numerical implementation, acceptable classification rate, and good performance in real time. This neural network has one hidden layer and one output layer. The hidden layer contains 20 neurons with hyperbolic tangent sigmoid activation function and the output layer has 5 neurons with a linear activation function constrained to the interval [0,1]. The selection of the number of hidden neurons has been determined directly by minimizing the combination of squared errors and weights, and by maximizing the classification rate on the training stage [48]. For both feature reduction methods discussed in Section 2.6.2, the size of the final augmented feature vector was 12 (three features for each of four EMG channels). This augmented feature vector is used as the input to a neural network classifier and a 5-dimension vector is returned as an output that represents the classified movements, i.e., flexion, extension, pronation, supination and rest. Continuous classification of movements using overlapped sliding windows was applied for training and validation of the neural network classifier. As described in Section 2.3, a segment length of 256 samples and an overlapping of 128 samples were adopted. For training the classifier, the first trial of each volunteer was extracted from the (our) database of experiments, comprising the sequence of five movements (see Section 2.1). In the training stage, only the steady-state signal was considered. Consequently, the transition periods were removed from the training data [28]. A Bayesian regularization algorithm [49,50] was used for training the neural network classifier. This method uses a modified performance function designed to minimize over-training, and obtains a classifier with improved generalization properties. The remaining trials of our database of experiments without removing the transition periods were used as validation data. Fig. 8. The figure shows the post processing of the motion classification system. Four classified motions are averaged and then, the final motion is the maximum of this average. Average M1–5 characterizes the movements. The decision increment is 128 ms to the outputs. Finally, in order to improve the motion classification accuracy, a post-processing step is applied. The analysis window size is 512 ms (two segment), and the decision increment is 128 ms (half segment). Four classified motions are averaged and the intended motion is assigned to the maximum of this average (see Eq. (21), and Fig. 8). This post processing produces a decision once every 128 ms. MMAX = max(M̄j ) = max j j  i+3 1 Mji 4 i  , j = 1, . . . , 5 (21) Fig. 9. The basic software consist on all those C-codes necessary to run the ADC, buffering data, filtering data, windowing, saving data to disk, graphics user interface, and to allow the option to run the process and control algorithms. Those tasks run in four parallel threads, in different priorities. This configuration avoids the competition between the fundamental thread and the GUI. Therefore, the hard real time is guaranteed. 162 E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 Fig. 10. The figure shows the joint coordinates calculation scheme. Selectors M1–4 are the active motions and M5 is rest. When the Selector M5 is active, an exponential decreasing function (with K < 1) is set. Constants K1 to K4 transform parameters to joint coordinates for each volunteer. Saturations from −90◦ to 90◦ are applied to avoid any damage to the robot. where i is the sliding time window index, and j = 1, . . ., 5 indicates the classified class motion. 2.8. Real-time myoelectric control Myoelectric control is widely used in assistive technologies, including multifunction prosthesis, wheelchairs, grasping control, virtual keyboards, diagnoses and clinical applications, such as functional neuromuscular stimulation. The main goal of this work is to control a robotic arm (Cyton ARM7) using the EMG signals from the upper limb and user’s visual feedback. The myoelectric control was done with a dedicated hardware and software system, and the control loop was closed using the visual perception of the user. The second acquisition system was used in this stage (see Section 2.2). The real time processing software was executed under QNX® real time operating system. QNX allowed multi-threads programming with different priority levels, from 1 to 255 (low to high priority). The applications were developed with an integrated development environment (IDE) in C-language. The IDE was performed under Windows® and the applications were transferred over Ethernet to a QNX host computer. This configuration decreases the development time and increases the versatility. The software was executed in four parallel threads with different priorities. A block diagram of the basic module’s structure is shown in Fig. 9. Technical details about this procedure can be found in [26,27]. The robotic arm, a Cyton ARM7 manipulator, was controlled by mimicking a human arm configuration with seven degrees of freedom. The joints used in this experiment were the elbow and the wrist rotation, and their coordinates were calculated using the Root-Mean Square (RMS) of the rectified signal, as the maximum likelihood estimator of muscular force [30]. When a movement is detected and correctly classified, the RMS of the data segment is integrated over time and the final joint coordinates are obtained by subtracting the integrated RMS of the agonist–antagonist muscles pairs. If movement switches to rest, the joint coordinate is modeled by an exponential decreasing function. This procedure is equivalent to low-pass filtering or smoothing the robot arm dynamics, as shown in Fig. 10. 3. Results and discussion In the following subsections, we present and discuss our results, divided into three parts. The first one shows the evaluation of the bispectrum estimators proposed in Section 2.5.3.3, while the second presents the classification performance of the bispectrum-based feature extraction. In the third subsection, the real time implementation of the myoelectric control scheme herein proposed is presented, applied to a robotic arm. 3.1. Evaluation of the bispectrum estimation The bispectrum estimators proposed in the Section 2.5.3.3 were evaluated in terms of accuracy, robustness against outliers and computational time. In addition, the underlying probability distributions of the estimates have been considered in the analysis of the results. It is reasonably assumed that EMG signals are sufficiently fast mixing sequences; that is, its values at widely-separated times are asymptotically independent [51]. For signals that exhibit this weakly time-dependent behavior, the moving block bootstrap resampling technique [52] could be used to generate a large number of independent pseudo-random sequences from the observed signal. With these new segments, it is possible to estimate the value of statistical quantities from the original signal via Monte Carlo simulations, e.g., mean, variance, distribution functions, and so on. 3.1.1. Accuracy of bispectrum estimation The accuracy of the bispectrum estimation depends on three factors, namely, the amount of data (segment length L), the bidimensional smoothing window function W(· , ·) and, finally, the accuracy of the third-order cumulant estimates ĉ3x−XX (m, n). By recalling the equations of Section 2.5, we can see that these factors are clearly shown in definition (4) and estimation (19) of the bispectrum. B3x (ω1 , ω2 ) = ∞ ∞   c3x (1 , 2 ) exp{−j(1 ω1 + 2 ω2 )} (4) 1 =−∞2 =−∞ B̂3x−XX (ω1 , ω2 ) = L L   ĉ3x−XX (m, n)w(m, n) exp{−j(ω1 m + ω2 n)} m=−Ln=−L (19) As was already pointed out in Section 2.3, the segment length and overlapping have been fixed to 256 samples and 128 samples, respectively. These two values were determined by the trade-off balance between the accuracy of the classification and the response times that satisfy the real-time constraints of the myoelectric control application. As it is noted in Section 2.5.3.3, two-dimensional smoothing window functions have been extensively analyzed by E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 Table 2 Mean and variance of bispectrum estimators from 1000 Monte Carlo simulation experiments. Estimator Arithmetic mean Median Trimmed mean D Mean (23) Variance (24) 1 2 4 1 2 4 1 2 4 0.031 0.0066 0.0051 0.00039 0.00014 0.00016 0.012 0.0027 0.0022 0.64 0.17 0.19 0.027 0.0038 0.0048 0.36 0.075 0.12 B̂3x−XX (0, 0) = L L   ĉ3x−XX (m, n)W (m, n) Table 3 Mean and variance of bispectrum estimators from 1000 Monte Carlo simulation experiments with the data sets that include outliers. Estimator Arithmetic mean Median Trimmed mean other authors. Based on these studies, Parzen window was selected because the estimated bispectrum showed the lowest variance. Finally, it is clear from the bispectrum definition (as a double discrete Fourier transform of the third-order cumulants) that the accuracy of the third-order cumulant estimates is an important factor that directly affects the accuracy of the bispectrum estimation. The bispectrum accuracy is analyzed through the mean (23) and variance (24) of the bispectrum estimated at the origin B̂3x−XX (0, 0). This term is a real random variable that is calculated by the matrixweighted average of the third-order cumulant estimates (see Eq. (22)) and, therefore, it provides an indirect evaluation of the thirdorder cumulant estimates for a fixed weight matrix W. (22) 163 D Mean (23) Variance (24) 1 2 4 1 2 4 1 2 4 0.035 0.0078 0.0065 0.00035 0.00015 0.00018 0.012 0.0027 0.0023 1.10 0.22 0.21 0.015 0.0064 0.011 0.50 0.056 0.075 although the variance shows a slight increase when the data are divided into four segments (D = 4). 3.1.2. Robustness property of bispectrum estimators: The robustness properties of bispectrum estimators were evaluated by analyzing the accuracy of B̂3x−XX-D (0, 0), when a single normally distributed outlier is introduced at a random (uniformly distributed) time into the original EMG signal. In data sets that include outliers, the variance of mean estimator grows further than the others estimators, while the robustness property of the median estimator results in nearly the same variance (as without outliers), as is shown in Table 3. It is recognized that the trimmed mean estimator has similar robustness properties. Nevertheless the variance is relatively greater than the median estimator. m=−Ln=−L The sample mean of B̂3x−XX (0, 0) is given by: NB 1  x(b) −XX B̂3 (0, 0) B̄3 (0, 0) = NB (23) b=1 and the sample variance as: NB B̃3 (0, 0) = 2 1  x(b) −XX (B̂3 (0, 0) − B̄3 (0, 0)) NB (24) b=1 where {x(b) , b = 1, 2, . . ., NB } are a set of NB independent pseudorandom sequences (segments) generated from the EMG signal in the moving block bootstrap technique. 3.1.1.1. Results of Monte Carlo simulation experiments. The previously outlined moving block bootstrap procedure is used to generate 1000 independent pseudo-random segments from a segment of EMG signal during a muscle contraction of 256 samples. Then, we calculate the sample mean (23), and sample variance (24) of the bispectrum estimators proposed in Section 2.5.3.3, namely, the arithmetic mean estimator, the median estimator, the symmetric trimmed mean estimator, and their ensembles average versions for three subsegments D = 1, 2, and 4 (see Section 2.5.3.2). Results are summarized in Table 2 and Fig. 11 from which the following observations were extracted: (i) The arithmetic mean estimator has a notable higher variance than the median and the symmetric trimmed mean estimators. (ii) The median estimator has a substantial lower variance than the other estimators. This result is consistent with the fact that the underlying distribution of z(l) was approximately Laplacian (see Fig. 3) and, in this case, the median is the maximum likelihood estimator. (iii) As it would be expected, both the mean and the variance decrease when the number of sub segments D increase, 3.1.3. Average computational time The real time constraints of myoelectric control require that the computational time for the bispectrum estimation, the features vector construction and the movements’ classification be smaller than the sliding time window of 128 ms. The average computational time was calculated using 1000 segments randomly selected from the four channels EMG signal, in different trials and volunteers. The time that the computer takes to perform the features extraction, features vector construction and movements’ classification were obtained for each segment. This last classification task is straightforward to implement, so the computational time is spent mainly by the following two tasks: • Bispectrum estimation: The bispectrum is calculated using the estimator evaluated in the previous sections, i.e., the arithmetic mean estimator, the median estimator, the symmetric trimmed mean estimator, and their ensembles average versions for three sub segments D = 1, 2, and 4. • Features vector construction: The Pythagorean means from the real part of the complex bispectrum method was used to obtain the features vector. This feature reduction method results more effective and faster (Section 2.6.2). These calculations were performed in Matlab, providing an estimate of the processing time. As it can be expected, the real time implementation is faster on the QNX operating system (Section 2.8). The average computational times for the estimators are shown in Table 4. The time difference between the ensemble estimators with two (D = 2) and four (D = 4) sub-segments are small. In the case Table 4 Averaged computational time. Estimator Arithmetic mean Median Trimmed mean D=1 D=2 D=4 61 ms 126 ms 311 ms 14 ms 32 ms 121 ms 14 ms 30 ms 109 ms 164 E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 Fig. 11. (a) The top of the figure shows a bar-plot of the estimated sample means (23) of the arithmetic mean, median, and trimmed mean bispectrum estimators for three ensembles versions (D = 1, 2 and 4). The top of the figure also shows an enlargement of the median mean estimator. (b) The bottom of the figure shows a bar-plot of the estimated sample variances (24). All means and variances are estimated with 1000 independent pseudo-random segments obtained with the moving block bootstrap method. of one segment (D = 1), the average times increases significantly, and this time is even greater than the sliding time windows for the trimmed mean estimator. • The percentage of correct classification rates were computed as the ratio of correct decisions over the total number of decisions multiplied by 100. 3.2. Classification performance of the bispectrum-based features Some authors, such as in [19,53], give values of correct classification rates for steady state signals, i.e., the transition periods between successive movements have been removed. Clearly, this approach is not suitable for hard real-time applications such as myoelectric control. Nevertheless, for comparison purposes, the classification rates were calculated with and without taking into account transitions between movements. But, for both cases, the transition periods were removed from the training data set. For continuous classification (with transitions), the average of the mean and standard deviation of the correct classification rates are summarized in Table 5 and Fig. 12. For a 12-dimension feature vector obtained from the module of the complex bispectrum matrix, the average of the mean and standard deviation classification rates are 92.22% and 2.97, respectively; and for features obtained from the square of the real part of the bispectrum matrix are 92.12% and 3.11, respectively. By comparing the performance of the module and real part feature vectors, it was observed that the two averages of the mean rates were very similar. There was only a small decrease (less than 2% in the worst case) for real part feature vector. In addition, on account his extensive experience, the amputee Vol. 1 shows an acceptable performance in both feature vectors. In this section, the performance evaluation of the neural network classifier described in Section 2.7 is presented. The correct classification rate is used to measure the classification accuracy for the two feature vectors proposed in Section 2.6.2. Recall that these feature vectors are defined by the quarter-root of the three classic means obtained from the module and the square of the real part of the complex bispectrum matrix (see Table 1). In order to properly interpret the results, the following issues are relevant: • The data set of movements was extracted from our own database of experiments (see Section 2.7). • The ensemble median bispectrum estimator (for D = 4) was used on the basis of the advantages outlined in Section 3.1. • The stages of features extraction and reduction, as well as the classification of movements using a two-layer feedforward neural network were performed in continuous mode (see Sections 2.6 and 2.7). • The feature vectors were calculated for each volunteer, and then, the neural networks were trained (see Section 2.7). E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 165 Table 5 The classification rates using our database with transitions. Volunteer Feature vector obtained from the module of the complex bispectrum matrix Feature vector obtained from the square of the real part of the complex bispectrum matrix Mean % SD Mean % SD 1 2 3 4 5 6 96.97 93.04 93.97 89.35 90.04 89.97 0.71 1.96 2.40 5.25 2.83 4.42 97.03 92.93 93.85 90.40 90.07 88.43 1.05 3.37 2.26 2.86 1.59 3.26 Averages 92.22 2.97 92.12 3.11 Mean %, mean percentage of classification rate; SD, standard deviation; Vol. 1, amputee. Table 6 The classification rates using our database without transitions. Volunteer Feature vector obtained from the module of the complex bispectrum matrix Feature vector obtained from the square of the real part of the complex bispectrum matrix Mean % SD Mean % SD 1 2 3 4 5 6 99.60 96.59 96.96 93.47 93.99 93.86 0.35 1.95 2.09 4.95 2.74 3.96 99.53 96.41 96.83 94.95 94.11 92.60 0.39 3.09 1.91 2.45 1.73 2.78 Averages 95.75 2.40 95.74 2.41 Vol. 1, amputee. The results without taking into account transitions between movements are summarized in Table 6. For the feature vector obtained from the module of the complex bispectrum matrix, the average of the mean and standard deviation classification rates are 95.75% and 2.40, respectively, and for features obtained from the square of the real part of the bispectrum matrix are 95.74% and 2.41, respectively. These classification rates are consistent with results obtained in [19,53]. When these results were compared with those obtained by continuous classification (Table 5), a slight improvement was noted in the classification rates. This improvement was expected, because most of the classification errors occur at the transitions between movements. feature vectors with an external database provided by Chan [53]. In this database, EMG signals were collected from 30 volunteers from seven sites on the forearm and one site on the biceps brachii. Within each trial, the volunteer repeats each limb motion four times, holding each motion for three seconds each time. A total of six trials were completed in a session. In this paper, we only use the data from one session of the fourth volunteer. Data from the first trial were used as training data and data from the remaining five trials were used as testing data. EMG data motions are: open hand, closed hand, supination, pronation, wrist flexion, wrist extension, and rest. 3.2.1. Cross-validation with an alternative database The purpose of this subsection is to validate the effectiveness of the neural network classifier for the same previous two Fig. 12. Bar-plot of the classification rates summarized in Table 5. The mean and standard deviation of the classification rate from the module feature vector and the real part feature vector are shown. Fig. 13. The user had to control the robot through muscle contraction and visual feedback to confirm the correct movements. The task involved a sequence of pronation, supination, flexion, extension movements interspersing with short rest intervals. E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 166 Fig. 14. The top four figures (a)–(d) show four EMG signals generated by the subject with a congenital amputation below the elbow. The movements are pronation, supination, flexion and extension. EMG signals are (a) pronator muscle; (b) brachioradialis muscle; (c) biceps muscle and (d) triceps muscle. The bottom two figures (e) and (f) show the joint coordinates of the Cyton robot arm. Joint coordinates are obtained by the subtraction of the integrated RMS of the agonist–antagonist muscles pairs. The wrist joint coordinate is from (e) the pronator–supinator functional muscles pairs, and the elbow joint coordinate is from (f) flexion–extension functional muscles pairs. The transition to rest state is made through an exponential decreasing function. The correct classification rates over six trials are summarized in Table 7. For the feature vector obtained from the module of the complex bispectrum matrix, the mean and standard deviation are 96.23% and 1.81, respectively, and for features obtained from the square of the real part of the bispectrum matrix are 95.73% and 1.91, respectively. Table 7 Correct classification rates calculated using the external database (hand protocol). Trials Feature vector obtained from the module of the complex bispectrum matrix Feature vector obtained from the square of the real part of the complex bispectrum matrix 1 2 3 4 5 6 98.29 96.20 93.18 95.26 97.09 97.37 97.78 96.41 92.93 93.98 97.35 95.95 Mean % SD 96.23 1.81 95.73 1.91 3.3. Myoelectric Control of a robotic arm The myoelectric control system based on features extracted from the EMG signal bispectrum has been shown in the block diagram of Fig. 1. The algorithms were executed in the real time operating system QNX and the design specifications were discussed in Section 2.8. These specifications explain, among other things, how the joint coordinates of the robot arm are calculated (Fig. 9, Section 2.8). In an additional experimental session, three volunteers (two able-bodied and one trained amputee volunteer) perform a task in real-time. The experimental task is a sequence of pronation, supination, flexion, and extension movements interspersed with short rest intervals. As we have already pointed out in Section 2.8, the robot joints used in this experiment were the elbow and the wrist. Two issues were evaluated for their influence on the system suitability: (i) the correct (movement) classification rate, and (ii) the system response delay. All volunteers successfully completed the control task, especially the amputee volunteer, due to his experience with myoelectric control systems (Fig. 13). A sequence of movements of E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 four EMG signals and the joints coordinates are shown in Fig. 14. The movement classification rate of 92% was calculated from the volunteers’ collected data. A classification error occurs more frequently during transitions between movements. These classification errors are naturally low-pass filtered by the robot dynamics. Therefore, it is unlikely that a robot could be sensitive to transitory misclassifications [28]. From a practical point of view, the users of the myoelectric control equipment have instantaneous visual feedback. This is so because the system delay does not exceed 128 ms, thereby avoiding any human–machine compatibility problems which may deteriorate the safety and controllability of the system. 4. Conclusions In previous work, several authors have already proposed the use of higher order statistics for off-line classification of EMG signals. Our approach differs from these earlier works in that we introduce a bispectrum-based continuous classification of EMG signals in a hard real-time application, as the case of myoelectric control of a robot arm. The proposed myoelectric control scheme was structured in several stages and we use the most advantageous techniques for EMG acquisition, data segmentation, bispectrum estimation, feature extraction and reduction, classification and control. The data segmentation was implemented as done in several other works in the literature, allowing for continuous EMG classification and for improving the performance of the system. Two initial assumptions were empirically verified, (i) that EMG signals are not a purely Gaussian process, thus the third order statistics are non-zero; and (ii) that HOS techniques cannot be discarded when the marginal probability distribution of the EMG signals is symmetric. These now verified facts allowed us to extract discriminative features from the EMG bispectrum, and thus achieving a good classifier. This work is focused on three bispectrum estimators members of the family of indirect class of conventional bispectrum estimators: MEAN, MED, and TRIMM estimators; and their ensemble averaging versions: MEAN-D, MED-D, and TRIMM-D estimators. On the feature extraction stage, the bispectrum was estimated by using the MED-D cumulants estimator, offering robustness against outliers and low variance. This can be explained by the fact the underlying distribution of zm,n (l) was approximately Laplacian and in this case, the median is the maximum likelihood estimator. The two feature vectors, defined by the quarter-root of the three Pythagorean means obtained from the module, and the square of the real part of the complex bispectrum matrix, were analyzed and compared. The classification performance was validated with our own database of experiments, and with an external one. It was found that: (i) our results in real-time classification are quite similar to those obtained off-line by other authors; (ii) when comparing the performance of the module and real part feature vectors, it was observed that the two averages of the mean rates were very similar (there was only a small decrease, less than 2% in the worst case for real part feature vector). Also, all volunteers successfully completed the control task, especially the amputee volunteer due to his experience with myoelectric control systems. The computational time was low enough to allow the real-time implementation; and the delay does not exceed the segment time of 128 ms. Finally, we conclude that the bispectrum-based method provides a good balance between complexity and performance for real time applications, even when the bispectrum is a complex matrix function. References [1] C. De Luca, Physiology and Mathematics of Myoelectric Signals, IEEE Transactions on Biomedical Engineering BME-26 (6) (1979) 313–325. 167 [2] C.J. De Luca, Surface Electromyography: Detection and Delsys Inc., Boston, MA, 2002, Available in Recording, http://www.delsys.com/Attachments pdf/WP SEMGintro.pdf. [3] R.O. Duda, P.E. Hart, D.G. Stork, Pattern Classification, 2nd ed., WileyInterscience Publication, New York, 2001. [4] M. Oskoei, H. Hu, Myoelectric control systems: a survey, Biomedical Signal Processing and Control 2 (2007) 275–294. [5] M. Zardoshti-Kermani, B. Wheeler, K. Badie, R. Hashemi, EMG feature evaluation for movement control of uer extremity prostheses, IEEE Transactions on Rehabilitation Engineering 3 (1995) 324–333. [6] D. Roman-Liu, M. Konarska, Characteristics of power spectrum density function of EMG during muscle contraction below 30% MVC, Journal of Electromyography and Kinesiology 19 (2009) 864–874. [7] X. Zhang, Y. Wang, P.S. Ray, Han, wavelet transform theory and its alication in EMG signal processing, in: Proc. of IEEE Seventh International Conference on Fuzzy Systems and Knowledge Discovery, vol. 5, no. 8, 2010, pp. 2234–2238. [8] M. Kurzynski, T. Woloszynski, A. Wolczowski, Multiclassifiers with competence function alied to the recognition of EMG Signals for the control of bio-prosthetic hand, in: Proc. of the 9th International Conference on Information Technology and Alications in Biomedicine ITAB 2009, 2009, pp. 1–4. [9] P.K. Artemiadis, K.J. Kyriakopoulos, An EMG-based robot control scheme robust to time-varying EMG signal features, IEEE Transactions on Information Technology in Biomedicine 14 (3) (2010) 582–588. [10] C.L. Nikias, J.M. Mendel, Signal processing with higher-order spectra, IEEE Signal Processing Magazine 10 (3) (1993) 10–37. [11] C.L. Nikias, M.R. Raghveer, Bispectrum estimation, a digital signal processing framework, Proceedings of the IEEE 75 (7) (1987) 869–891. [12] A. Swami, J.M. Mendel, C.L. Nikias, Higher-Order Spectral Analysis Toolbox for use with Matlab, User’s Guide, MathWorks Inc., Natick, MA, 1998. [13] M.B. Priestley, Spectral Analysis and Time Series, 11th ed., Academic Press Inc., Great Britain, 1981, pp. I871–I874. [14] D.R. Brillinger, Time Series, Data Analysis and Theory, Siam ed., Holden Day Inc., San Francisco, 2001 (Chapter 2). [15] D.R. Brillinger, Some history of the study of higher-order moments and spectra, Statistica Sinica 1 (1991) 465–476. [16] D.R. Brillinger, An introduction to polyspectra, The Annals of Mathematical Statistics 36 (5) (1965) 1351–1374. [17] M.G. Kendal, The Advanced Theory of Statistics, I, 2nd ed., Charles Griffin & Company Limited, London, 1945 (Chapter 3). [18] K. Nazarpour, A.R. Sharafat, S.M. Firoozabadi, Alication of higher order statistics to surface electromyogram signal classification, IEEE Transactions on Biomedical Engineering 54 (10) (2007) 1762–1769. [19] X. Chen, X. Zhu, D. Zhang, A discriminant bispectrum feature for surface electromyogram signal classification, Medical Engineering & Physics 32 (2) (2010) 126–135. [20] P.A. Kaplanis, C.S. Pattichis, L.J. Hadjileontiadis, V.C. Roberts, Surface EMG analysis on normal subjects based on isometric volunteer contraction, Journal of Electromyography and Kinesiology 19 (1) (2009) 157–171. [21] K. Kanosue, M. Yoshida, K. Akazawa, K. Fujii, The number of active motor units and their firing rates in volunteer contraction of human brachialiis muscle, Japanese Journal of Physiology 29 (4) (1974) 427–443. [22] S. Shahid, J. Walker, G.M. Lyons, C.A. Byrne, A.V. Nene, Higher order statistics techniques alied to EMG signal, IEEE Transactions on Biomedical Engineering 52 (7) (2005) 1195–1209. [23] K. Yana, H. Mizuta, R. Kajiyama, Surface electromyogram recruitment analysis using higher order spectrum, in: Proc. of the IEEE 17th Annual Conference Engineering in Medicine and Biology Society, vol. 2, 1995, pp. 1345–1346. [24] E. Plévin, D. Zazula, Decomposition of surface EMG signals using non-linear LMS optimisation of higher-order cumulants, in: Proc. of the 15th IEEE Symposium on Computer-Based Medical Systems, 2002, pp. 5912–5915. [25] F. Molinari, M. Knaflitz, P. Bonato, M.V. Actis, Electrical manifestations of muscle fatigue during concentric and eccentric isokinetic knee flexion–extension movements, IEEE Transactions on Biomedical Engineering 53 (2006) 1309–1316. [26] E. Orosco, N. López, C. Soria, F. di Sciascio, Surface electromyogram signals classification based on bispectrum, in: Proc. of the 2010 Annual International Conference of the IEEE of Engineering in Medicine and Biology Society (EMBC), 2010, pp. 4610–4613. [27] E. Orosco, N. López, C. Soria, F. di Sciascio, Hard real time modular system for acquisition, processing and classification of EMG signals, in: Proc. of the RPIC Buenos Aires, 2011. [28] K. Englehart, B. Hudgins, A robust, real-time control scheme for multifunction myoelectric control, IEEE Transactions on Biomedical Engineering 50 (7) (2003) 848–854. [29] E. Clancy, E. Morin, R. Merletti, Sampling, noise-reduction and amplitude estimation issues in surface electromyography, Journal of Electromyography and Kinesiology 12 (2002) 1–16. [30] E.A. Clancy, N. Hogan, Multiple site electromyograph amplitude estimation, IEEE Transactions on Biomedical Engineering 42 (2) (1995) 203–211. [31] Y. Huang, K. Englehart, B. Hudgins, A. Chan, Optimized Gaussian mixture models for uer limb motion classification, IEEE Transactions on Biomedical Engineering 52 (11) (2005) 1801–1811. [32] E.A. Clancy, N. Hogan, Probability density of the surface electromyogram and its relation to amplitude detectors, IEEE Transactions on Biomedical Engineering 46 (6) (1999) 730–739. 168 E.C. Orosco et al. / Biomedical Signal Processing and Control 8 (2013) 153–168 [33] H. Al-Timemy, G. Bugmann, A. Jackson, K. Nazarpour, A note on the probability distribution function of the surface electromyogram signal, Medical & Biological Engineering & Computing, submitted for publication. Available in http://www.staff.ncl.ac.uk/k.nazarpour/papers/EMG pdf.pdf. [34] A.S. Cherniz, C.E. Bonell, C.B. Tabernig, Study of the SEMG probability distribution of the paretic tibialis anterior muscle, Journal of Physics: Conference Series 90 (1) (2007) 1–7. [35] G.R. Naik, D.K. Kumar, S.P. Arjunan, Kurtosis and negentropy investigation of myo electric signals during different MVCs, in: ISSNIP Biosignals and Biorobotics Conference, 2011, pp. 1–4. [36] K. Nazarpour, A.R. Sharafat, S.M. Firoozabadi, Negentropy analysis of surface electromyogram signal, in: IEEE/SP 13th Workshop on Statistical Signal Processing, 2005, pp. 974–977. [37] C. Nikias, A. Petropulu, Higher-order Spectral Analysis: a Nonlinear Signal Processing Framework, Prentice-Hall, Inc., NJ, 1993. [38] M.R. Raghuveer, Higher-order statistics: laying a myth to rest, in: Proceedings of the 28th Asilomar Conference on Signals, Systems and Computers, vol. 1, 1994, pp. 5–8. [39] M.R. Raghuveer, Third-order statistics: issue of PDF symmetry, IEEE Transactions on Signal Processing 43 (7) (1995) 1736–1738. [40] P.J. Huber, Robust Statistics, John Wiley and Sons, New York, 1981. [41] A.K. Nandi, Robust estimation of third-order cumulants in applications of higher-order statistics, IEEE Proceedings of Radar and Signal Processing 140 (6) (1993) 380–389. [42] A.K. Nandi, D. Mampel, Development of an adaptive generalized trimmed mean estimator to compute third-order cumulants, Signal Processing 57 (3) (1997) 271–282. [43] P.J. Rousseeuw, A.M. LeRoy, Robust Regression and Outlier Detection, Wiley, New York, 1987. [44] Y. Zhang, D. Hatzinakos, A.N. Venetsanopoulos, Bootstraing techniques in the estimation of higher-order cumulants from short data records, in: 1993 IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP93, vol. 4, 1993, pp. 200–203. [45] G.W. Brown, On small-sample estimation, The Annals of Mathematical Statistics 18 (4) (1947) 582–585. [46] M.S. Hussain, M.B.I. Reaz, F. Mohd-Yasin, M.I. Ibrahimy, Electromyography signal analysis using wavelet transform and higher order statistics to determine muscle contraction, Expert Systems 26 (1) (2009) 35–48. [47] G. Li, A.E. Schultz, T.A. Kuiken, Quantifying pattern recognition-based myoelectric control of multifunctional transradial prostheses, IEEE Transactions on Neural Systems and Rehabilitation Engineering 18 (2) (2010) 185–192. [48] B. Curry, P.H. Morgan, Model selection in neural networks: some difficulties, European Journal of Operational Research 170 (2006) 567–577. [49] D.J.C. MacKay, Bayesian interpolation, Neural Computation 4 (3) (1992) 415–447. [50] F.D. Foresee, M.T. Hagan, Gauss–Newton aroximation to Bayesian regularization, in: Proceedings of the International Joint Conference on Neural Networks, vol. 3, 1997, pp. 1930–1935. [51] P. Doukhan, Mixing: Properties and Examples, Lecture Notes in Stat. 85, Springer-Verlag, New York, 1994. [52] A.M. Zoubir, D.R. Iskander, Bootstrap Techniques for Signal Processing, Cambridge University Press, New York, 2004. [53] A.R. Goge, A.D.C. Chan, Investigating classification parameters for continuous myoelectrically controlled prostheses, in: 28th Conference of the Canadian Medical & Biological Engineering Society, 2004, pp. 141–144.