Electromyography and mechanomyography signal recognition: Experimental analysis using multi-way array decomposition methods

biocybernetics and biomedical engineering 37 (2017) 103–113

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/bbe

Original Research Article

Electromyography and mechanomyography signal recognition: Experimental analysis using multi-way array decomposition methods Andrzej Wołczowski, Rafał Zdunek * Faculty of Electronics, Wroclaw University of Science and Technology, Wroclaw, Poland

article info

abstract

Article history:

In this study, we considered the problem of controlling a prosthetic hand with noisy

Received 29 February 2016

electromyography (EMG) and mechanomyography (MMG) signals. Several dimensionality

Received in revised form

reduction methods were analyzed to assess their efficiency at classifying these signals,

29 July 2016

which were registered during the performance of grasping movements with various objects.

Accepted 2 September 2016

Using the cross-validation technique, we compared various dimensionality reduction

Available online 20 January 2017

methods, such as principal components analysis, nonnegative matrix factorization, and some tensor decomposition models. The experimental results demonstrated that the high-

Keywords:

est classification accuracy (exceeding 95% for all subjects when classifying 11 grasping

Biosignal processing

movements) and lowest computational complexity were obtained when higher-order sin-

Electromyography

gular value decomposition was applied to a multi-way array of multi-channel spectrograms,

Mechanomyography

where the temporal EMG/MMG signals from all channels were concatenated.

Multi-array decomposition method

© 2017 Nałęcz Institute of Biocybernetics and Biomedical Engineering of the Polish

Supervised classification

1.

Introduction

The measurement and analysis of biosignals is used widely in medical diagnostics, sports, and rehabilitation. Computer recognition of biosignals provides significant support in these applications, thereby allowing the automation of human tasks. A new application of biosignal recognition is the control of machines, particularly processes in the human body that depend on the human will, such as the contraction of individual skeletal muscles. The accompanying myoelectric and myomechanic signals can be used to control a machine's

Academy of Sciences. Published by Elsevier B.V. All rights reserved.

action after recognition. In contrast to diagnostics, biosignal classification can be used to interpret human intentions. However, for both automatic diagnosis and automatic control, the efficiency of recognition depends on the methods used to analyze (representation and dimensionality reduction) and classify biosignals. A particular challenge is the control of prosthetic hands. The human hand is capable of performing a variety of movements, thereby allowing gesticulation and the playing of musical instruments, but mainly grasping and manipulating objects. The hand is also a source of complex tactile sensations. The loss of a hand dramatically reduces the

* Corresponding author at: Faculty of Electronics, Wroclaw University of Science and Technology, Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland. E-mail address: [email protected] (R. Zdunek). http://dx.doi.org/10.1016/j.bbe.2016.09.004 0208-5216/© 2017 Nałęcz Institute of Biocybernetics and Biomedical Engineering of the Polish Academy of Sciences. Published by Elsevier B.V. All rights reserved.

104

biocybernetics and biomedical engineering 37 (2017) 103–113

possibility of full human functioning in personal and professional life, and the loss of both hands practically excludes independent functioning. The aim of the prosthetic hand is to partially restore the function of the lost limb, especially its manipulating and prehensile functions. To achieve these needs, the prosthesis must have the ability to perform a variety of movements and to adopt different finger configurations. This requires the transmission of a wide range of control decisions to the prosthesis. Modern transradial prostheses are generally controlled via the recognition of hand stump muscle-derived signals. This type of control is possible after an upper limb is amputated below the elbow, where the forearm stump retains a significant part of the muscle that moves the fingers in the healthy hand. The movements of these muscles are still subject to the human will. During contraction, the muscles produce both electrical (myopotentials) and mechanical (myovibration) signals, which are referred to as electromyography (EMG) and mechanomyography (MMG) signals, respectively. These signals are associated with the movements of the healthy hand and fingers, and they can clearly represent the will of the user with respect to control of the prosthesis. It is essential that these signals are recorded non-invasively by suitable sensors placed directly above the muscles on the skin of a forearm. Therefore, prosthetic biocontrol involves measuring surface signals from the hand stump muscles and then recognizing the type of intended prosthesis action by classifying these observations [26,11,15,13]. For obvious reasons, the basic quality criterion for this type of control is the reliability of intent recognition (because the prosthesis cannot perform unintended movements). It is easy to show that this criterion is closely related to a fixed set of movements. Changing the repertoire of movements significantly affects the reliability of their recognition. Another criterion is the time required for recognition, which affects the prosthesis response time to control biosignals and it strongly influences agility. The response delay (also due to psychological reasons) results in a significant decrease in the precision of motion control. The third criterion is based on the requirements of the computer system resources (CPU and memory usage) used for running the algorithm. Clearly, this is not related to the synthesis stage of control (e.g., classifier training), but it is important for the action stage of control. Unfortunately, the number of recognition errors increases rapidly if more classes need to be recognized. This is partly because the information regarding human intentions contained in the registered signal is usually distorted by noise, and because the form in which it appears in the signal loads the classification process (disrupts the efficiency and reliability). The preceding step of classification involves the process employed for analyzing biosignals. In our experiments, we measured the muscular activity using both EMG and MMG signals. To extract useful features for classifying various muscular activities, the temporal waveforms of the signals were transformed into the time-frequency domain. This transformation is usually performed with the short-time Fourier transform (STFT) but other transformations, such as the discrete wavelet transform (DWT) can be also used [5,16]. An appropriate transformation of the waveform into feature vectors can significantly decrease interference (in the raw

signal). The set of features should retain information regarding human intention (human control decision). This part of the signal analysis process is referred to as feature extraction. Most feature extraction methods produce redundant sets of features, but we should apply the next part of the signal analysis process to increase the efficiency of classification, i.e., selection/reduction of features, which yields a new optimal set of features (optimal in the sense of the minimum classification errors). Many of the methods for feature dimensionality reduction and selection can be applied to improve the classification of EMG/MMG signals [20]. One of the best-known techniques is principal components analysis (PCA) [9], which captures the orthogonal directions in the input data space where the variance is maximized. PCA works very well for timeseries, but two important problems occur when it is applied to 2D data, such as spectrograms. First, the low-dimensional factors obtained with PCA contain negative entries, which substantially hinder their physical interpretation. Second, PCA loses information regarding the 2D interactions between the variables in a spectrogram or scalogram. Thus, the lowdimensional factors provided by PCA do not contain information regarding the local smoothness of spectrograms, which can considerably reduce the efficiency of classification. Several approaches have been proposed to address these problems. To avoid the first problem, nonnegativity constraints can be imposed on all of the estimated factors to obtain the nonnegative matrix factorization (NMF) model [14,3]. Yazama et al. [28] first applied the NMF to the recognition of EMG signals. Subsequently, Theis and Garcia [22] compared various NMF algorithms in terms of their effectiveness at decomposing EMG signals. Recent advances in the use of the NMF for processing these signals include: the identification of EMG finger movements [1], separating ECG from a single-channel EMG signal [17], analysis of motor modules of human locomotion [18], coherence analysis of elbow and shoulder muscles [21], extraction of neural control information from EMG signals [8], and the classification of hand grip configurations [15]. Despite the undoubted advantages of the NMF for processing EMG signals, the other problem cannot be alleviated if a 2D representation is vectorized, but multi-way array decomposition methods [3] can be used to avoid this issue, which provide multi-linear features that capture the nature of the high-dimensional data along each mode. After the decomposition of a multi-way array of ordered spectrograms, the separate features reflect the temporal and frequency structures, as well as the discriminant features corresponding to the mode of ordering in the spectrograms. If nonnegativity constraints are imposed on each estimated factor, these methods can be regarded as extensions of NMF. In particular, nonnegative tensor factorization is a simple extension of NMF to multi-way nonnegative data. The strategies used in its application to model dimensionality reduction and the classification of multi-way array data were presented by [3,19]. According to a strategy based on core tensor classification, Xie and Song [27] performed experiments with twochannel EMG signals, where they successfully classified two hand actions: full elbow flexion and full elbow extension. EMGbased experiments with multi-way array decomposition methods were also performed by Kim et al. [11], where EMG

biocybernetics and biomedical engineering 37 (2017) 103–113

observations were represented by multi-linear decomposition models to extract the synergies required to control a robot hand. The multi-factor models were trained using the fundamental alternating least squares (ALS) algorithm [3]. In this study, we analyzed several computational strategies for feature dimensionality reduction using both EMG and MMG signals. We compared PCA, NMF, and four multi-factor decomposition approaches in terms of their usefulness for classifying 11 grasping movements. Various numerical algorithms were used for training and various constraints, such as nonnegativity and orthogonality, were imposed on the estimated factors. In the case of the Tucker decomposition model [3], classification can be performed using multi-linear features represented by the core tensor or by using the features of a factor that corresponds to the mode of ordering in the spectrograms, and we studied both approaches. The biosignals were collected from an arm while performing various grasping movements, and the muscular activity was monitored using both EMG and MMG signals. The remainder of this paper is organized as follows. The first section introduces the problem of controlling a prosthetic hand with biosignals. The concepts of EMG and MMG signals, and the factors that affect the quality of their recognition are discussed in Section 2. The dimensionality reduction methods used for extracting the most meaningful features from biosignals are presented in Section 3. Section 4 explains the classification methods used to perform the experiments. The results obtained after classifying various grasping movements are described in Section 5. Finally, we give our conclusions in Section 6.

2. Biosignals and factors that affect the quality of their recognition 2.1.

Physiology of biosignals

When skeletal muscles are activated, their cells produce an electric potential, which is called the myopotential due to the movement of ions (Na+, K+-ions) flowing out and into the interior of the myocytes (muscle fibers) that are stimulated to contract. Due to progressive depolarization of the stimulated cells, Ca2+-ions are released and they penetrate deep into the contractile protein structures of the cell myofibrils, thereby changing the spatial structure of the protein and allowing the contraction process. During the muscle fiber contraction process, the chemical energy from myocytes is converted into mechanical energy of motion (chemo-mechanical transduction). The contraction of a single fiber is accomplished according to the ‘‘all-or-nothing’’ principle, where the strength and shrinkage rate of the whole muscle depend on the number of fibers activated. The individual muscle fibers are activated by the incoming axons of the motoneurons. Each motoneuron can have synapses with multiple muscle fibers. The muscle fibers activated simultaneously by the same motoneuron (along with it) are called a muscle motor unit. The activity of individual motor units in a contracting muscle changes randomly over time, where they are activated (recruited) and released in turn. However, the number

105

of recruited units remains constant and proportional to the muscle activation level, and it increases when the required muscle power is higher. Myocyte activation is accompanied by a change in the electrical potential between the interior of the cell and its exterior from
60 mV to +50 mV. This change in the potential can propagate through the surrounding tissue and reach up to the surface of the limb. The superposition of changes in potential from all the active muscle cells can be recorded on the skin over the active muscles as an EMG signal. The root mean squared (RMS) value of an EMG signal from human skeletal muscles clearly depends on the degree of contracting muscle activation, as well as depending on the locations of the measurement electrodes relative to the muscles, with values in the range of 0–1.5 mV. The effective energy of the signal falls into the band of 0–500 Hz and the highest values are in the range of 50–150 Hz. The potentials change in a similar manner for individual myocytes. However, their superposition on a body surface due to spatial filtration of the signals in body tissues (different degree of dampening frequency components) depends significantly on the spatial positions of the active fibers relative to the locations of the measurement electrodes. Hence, the form of an EMG signal contains information about the muscles that are active. Various movements of fingers involve different muscle groups in the forearm so ultimately the EMG signals carry information about the type of hand movement being performed. Mechanical vibrations comprise another type of signal associated with skeletal muscle activity, where they are formed during the activation of individual fibers in the muscle motor units due to changes in their outer geometry (shortening and swelling) during contraction, the friction of muscle membranes moving relative to each other, as well as the movements of tendons and joints. Similar to myopotentials, these vibrations spread along a limb tissue and reach the surface where their superposition can be recorded as a MMG signal using a microphone or a vibrometer placed on the skin. Like an EMG signal, the form of a MMG signal also depends on the type of active muscle and thus it carries information about the type of movement. Due to differences in the nature of these two types of signals, the information contained in them may complement each other, thereby yielding a synergistic effect.

2.2.

Factors that disrupt the recognition of biosignals

Biocontrol is based on the classification of measured biosignals, where a particular decision that controls the motion of a biomanipulator, bioprosthesis, or other technical device is assigned to each of the recognized classes. To ensure the reliability of recognition, the shape of the signals recorded during the same hand movement should always be the same (unchangeable). However, both of the physiological phenomena discussed above, which underlie the formation of biosignals as well as their (biosignals) measurement process itself, do not provide signals with unchanging forms. Physiologically, the changes in EMG and MMG signals can be explained by variations in the patterns of recruited muscle motor units due to muscle fatigue or weakness. Fatigue is affected by the concentration of electrochemical metabolites within the surrounding muscles (K+-ions accumulate outside

106

biocybernetics and biomedical engineering 37 (2017) 103–113

the muscle membrane), which can be detected by monitoring the mean or median frequency, or based on the decay of the signal envelope (RMS). The signal also changes due to a decrease in the muscle conduction velocity. In addition, during their registration, the biosignals are affected by electrical and mechanical phenomena that accompany the signal measurement process, which for EMG signals comprise:  external electromagnetic fields superimposed on the signal, and  imperfect contact between the electrode and the skin, which may be due to changes in skin conductance caused by perspiration, the composition of sweat, the adhesion of electrodes, etc.). The electromagnetic smog (usually at the frequency of 50 Hz in a power network, but also in mobile phone, radio, and television bands) may induce skin interference with a magnitude that is 103 times greater than the desired signal. The electrode-skin contact also has a great effect on the signalto-noise ratio. For MMG signals, the primary source of interference is outside noise as well as (in the case of prosthesis) sounds from the prosthesis drives and interactions with the grasped object. These phenomena increase the intraclass dispersion in the measured signals, thereby impairing the reliability of signal recognition and hindering the decision-making process required for control. Ambient electromagnetic noise can be eliminated effectively by using a differential measurement system [26], which comprises two active electrodes placed directly on the skin at the test point (above the examined muscle) and an additional reference electrode placed as far as possible from and above an electrically neutral tissue (e.g., directly above a bone or a joint). The signals obtained from the active electrodes are subtracted from each other and then amplified. Thus, the common components of the signals, including the surrounding noise, are subtracted and the differential signal is amplified. The other type of interference (due to skin-electrode contacts and movements of the cables connecting the electrodes with the measuring amplifier) can be reduced significantly by applying an amplifier for high input impedance (much higher than the impedance of the skin) and low output impedance, which is positioned directly on the electrodes. Thus, in the concept of an active electrode, the electrode is integrated with a pre-amplifier to give a signal that resists external disturbances. For MMG signals, the basic problem is isolating the microphone-sensor from external sound sources and ensuring the acquisition of the muscle sound propagating in the patient's tissue. External interference can be reduced effectively by:  using sound insulating materials to produce the probe housing, and the acquisition of a good signal can be ensured by:  forming the probe housing from an appropriate acoustic chamber, matching the acoustic impedance of the

microphone to the impedance of the sound source (the skin and tissue stump); and  reducing the noise generated by the artificial hand (using appropriate materials for the gears and transmission parts, as well as the covering prosthetic fingers with suitable materials). The last factor that needs to be considered is the selection of the site(s) for recording signals. In commercial prosthetic hands, measurement systems based on two electrodes are the most common [8,13,15,26]. An important part of the registration process is optimizing the bio-signals, which requires the determination of appropriate positions for the electrodes on the forearm. In this study, we analyzed the experimental selection process based on the number of redundant electrodes (channel measurement system) and we determined the best locations for the electrodes by considering the effects on the recognition quality of information obtained from the individual channels.

3.

Dimensionality reduction

Each grasping movement was measured with an array of EMG/ MMG sensors. Multi-way array data can be processed in various ways to extract the features for classification. The EMG and MMG modalities can be combined in the model dimensionality reduction stage to estimate the common features. In another approach, the modalities can be treated separately and the features can then be combined. In this study, we treated the two modalities separately. Let Y ¼ ½yi1 ;i2  2 RI1 I2 be a 2D discrete representation of any observed EMG/MMG signal y(t), defined at discrete time instants {t0, t1, . . ., tT
1}, where i1 = 1, . . ., I1 and i2 = 1, . . ., I2. By applying the STFT to y(t), Y can be expressed by the spectrogram, i.e., an estimate of the power spectral density of the signal in time segments. The indices i1 and i2 refer to the frequency bin and time segment, respectively. In the case of the DWT, Y can be regarded as the scalogram representing the wavelet energy density, which contributes to the signal energy on the scale i1 and the time shift i2. The 2D representations of the observed signals in each modality can be collected in a 4D multi-way array Y ¼ ½yi1 ;i2 ;...;i4  2 RI1 I2 I3 I4 , where I3 is the number of EMG or MMG sensors, and I4 is the number of grasping movements. According to the notation in [12], Y is a fourth-order tensor. Therefore, the measurement of the kth grasping movement is modeled by a 3D subtensor Y i4 ¼k 2 RI1 I2 I3 , which belongs to one of C classes. Assuming that each movement is measured within a time window of 2 s and the sampling frequency is 1 kHz, the number of time instances amounts to T = 2000. Using the STFT with a 64-sample window, 50% overlapping, and 128 frequency points, we have: I1 = 65 and I2 = 61. When the DWT with 32 scales is applied, the size of the scalogram is determined by I1 = 32 and I2 = 2000. Our measurement system had eight channels for EMG signals and the same number for the MMG signals, i.e., I3 = 8 for each modality. Thus, each Y i4 ¼k has about 31720 variables if the STFT is applied. Moreover, the variables come from various subspaces and they have different meanings, which motivated

107

biocybernetics and biomedical engineering 37 (2017) 103–113

us to employ dimensionality reduction methods to extract the most meaningful information from each Y i4 ¼k . We also considered another model for multi-channel signal representation, where the channels were concatenated. In this model, the temporal waveforms contained 16,000 samples (eight channels time 2000 samples) in each modality. With respect to the spectral characteristics of the observed signals, the number of samples can be reduced by using decimation. This approach is justified by observations, which showed that a large amount of the spectral energy of the registered signals f was located below 4s , where fs is the sampling frequency. Thus, 2nd order decimation was applied to each concatenated signal. After using the STFT with the same parameters mentioned above, we obtained the third-order tensor Y 2 RI1 I2 I3 , where I1 = 65, I2 = 249, and I3 is the number of grasping movements. In this case, Y i3 ¼k had only 16185 elements but dimensionality reduction was still needed.

3.1.

PCA

Several dimensionality reduction methods are available but PCA [9] is probably the most popular in many research areas. Using orthogonal transformations, PCA extracts a given number of principle components (PCs), which capture the largest possible variance and they are linearly uncorrelated. Hence, a few of these components may explain a considerable amount of the total variance in the data. To apply PCA to our tensorial data, the tensor had to be reshaped into a matrix. This can be achieved with an unfolding operation [12], which transforms an Nth order tensor Y ¼ Q In  I p 6¼ n p along mode n, ½yi1 ;...;iN  2 RI1   IN into a matrix YðnÞ 2 R where n = 1, . . ., N. In our case, N = 3 or N = 4. By unfolding our 4th order tensor Y ¼ ½yi1 ;i2 ;...;i4  2 R

I1 I2 I3 I4

along the 4th mode,

we obtained: Yð4Þ ¼ ½yi4 ;j  2 RI4 I1 I2 I3 , where j = i1 + (i2
1)I1 + (i3
1)I1I2. The third order tensor can be unfolded in a similar manner. In the following, we present model dimensionality reduction methods only for the 4D tensor and we assume that the 3D tensor can be treated in a similar manner. PCA of the data matrix Y(4) can be expressed by the following matrix factorization model: Yð4Þ ¼ ZVTJ ;

(1) I4 I1 I2 I3

is the centralized data mawhere Yð4Þ ¼ Yð4Þ
EðYð4Þ Þ 2 R trix, the column vectors of Z 2 RI4 J determine the PCs, and VJ ¼ ½v1 ; . . .; vJ  2 RI1 I2 I3 J is a submatrix created from the J largest eigenvectors of the covariance matrix obtained from Y ð4Þ . The eigenvectors are mutually orthogonal and they are referred to as the feature vectors. By using PCA for model dimensionality reduction, we have: J  I1I2I3. The number of PCs, i.e., the number J, can be roughly estimated by observing the behavior of the eigenvalues {li}. The ratio of the variance explained by J PCs relative to the total variance is given by: PJ li j ¼ PI1i¼1 100%: I2 I3 li i¼1

(2)

Hence, J should be as small as possible but it should also be selected to maximize j. The problem of determining the optimal number of PCs has been discussed widely, e.g., [10,25]. In our approach, we determined J by setting j to a constant value equal to 80%.

3.2.

Tucker decomposition

Extracting uncorrelated information from the tensor Y unfolded along only one mode may result in the loss of spatial information for neighboring pixels in a 2D representation, as well as regarding the information contribution of each sensor to a given class. To avoid this problem, model dimensionality reduction can be performed using tensor decomposition methods, such as the Tucker decomposition [24], CANDECOMP/PARAFAC (CP) [2,7], or various tensor decomposition methods with nonnegativity constraints [3]. Kim et al. [11] used tensor factorization models to extract the multi-linear features from EMG signals. The features were then used to analyze the synergies among grasping movements for a robot hand, which was remotely controlled by human movements. Xie and Song [27] attempt to extract similar features but they used the Tucker decomposition with nonnegativity constraints. In our approach, we used the Tucker model with orthogonality constraints to extract lowdimensional multi-linear features from the array of EMG and MMG signals for the supervised classification of grasping movements. The basic Tucker decomposition model has the following form:

Y ffi G1 Uð1Þ 2 Uð2Þ 3    N UðNÞ ¼

J1 X J2 X j1 ¼1j2 ¼1



JN X

ð1Þ

ð2Þ

ðNÞ

1

2

N

gj1 ;j2 ;...;jN uj uj    uj ;

jN ¼1

(3)

where G ¼ ½gj1 ;j2 ;...;jN  2 RJ1 J2   JN for Jn ≤ In and n = 1, . . ., N is the ðnÞ

ðnÞ

ðnÞ

core tensor, UðnÞ ¼ ½u1 ; . . .; uJn  ¼ ½uin ;j  2 RIn Jn is the factor n

matrix, which contains the feature vectors along the nth mode of the tensor Y. The operator n denotes the tensor-matrix product along the mode n and the operator is the outer product between the vectors. The tensor G has the (J1, J2, . . ., JN)-rank, i.e., 8n : rank(G(n)) = Jn, where G(n) is obtained by unfolding the tensor G along the nth mode. If no constraint is imposed on the factors {U(n)} or the core tensor G, then the decomposition in (3) is generally not unique. Under the assumption that the column vectors in {U(n)} are mutually orthogonal, the Tucker decomposition can be regarded as an extension of PCA to multi-way arrays. In T particular, when 8n : Jn = In and ðUðnÞ Þ UðnÞ ¼ I 2 RIn In , the Tucker decomposition is higher-order singular value decomposition (HO-SVD) [4]. In this case, the core tensor is not diagonal but it satisfies the following conditions: hGin ¼a ; Gin ¼b i ¼ 0 if a 6¼ b, and hGin ¼a ; Gin ¼b i > 0 for a = b, as well as: jjGin ¼1 jjF jjGin ¼2 jjF    jjGin ¼N jjF . If the condition that Jn  In holds for any mode n, then HO-SVD can be interpreted as a truncated version, which leads to strong model dimensionality reduction. In our approach, we assumed that: Jn ≤ In for n = 1, T . . ., 3, J4  I4, and 8 n : ðUðnÞ Þ UðnÞ ¼ I 2 RJn Jn .

108

biocybernetics and biomedical engineering 37 (2017) 103–113

By unfolding the multi-way model in (3) with respect to the mode n, we have: 
n T

ðnÞ

YðnÞ ¼ U GðnÞ ðU

Þ ;

(4)

where Q Q Ip  J p 6¼ n p . U
n ¼ UðNÞ     Uðnþ1Þ Uðn
1Þ     Uð1Þ 2 R p 6¼ n The symbol  denotes the Kronecker product. Based on the orthogonality conditions for the core subtensors, we know that G(n) has orthogonal rows. Thus,

3.3.

CP decomposition

The model in (3) can be regarded as a general form of tensor decomposition. Assuming that J = J1 =    = J4 and G ¼ I 2 RJJJJ is a super-diagonal identity tensor, then the model (3) is the following CP decomposition: Y ffi I 1 Uð1Þ 2 Uð2Þ 3    N UðNÞ ¼

J X ð1Þ ð2Þ ðNÞ uj uj    uj ;

(7)

j¼1 ðnÞ

ðnÞ

where 8 n : UðnÞ ¼ ½u1 ; . . .; uJ  2 RIn J . By unfolding the model (7) along the nth mode, we obtain:

T

T

YðnÞ ðYðnÞ ÞT ¼ UðnÞ GðnÞ ðU
n Þ U
n ðGðnÞ ÞT ðUðnÞ Þ T

¼ UðnÞ LðUðnÞ Þ 2 RIn In ;

(5)

Jn Jn

is a diagonal matrix. Assuming that rank(Y(n)) where L 2 R = In, then the matrix Y(n)(Y(n))T is symmetric and positive definite. It is easy to see that the column vectors of U(n) are eigenvectors of Y(n)(Y(n))T. Hence, similar to PCA, the factor U(n) can be computed by applying eigenvalue decomposition (EVD) to Y(n)(Y(n))T, and then selecting the eigenvectors that to the Jn largest eigenvalues. Since correspond T 8 n : ðUðnÞ Þ UðnÞ ¼ I, then the core tensor can easily be computed by:

T

T

T

G ¼ Y1 ðUð1Þ Þ 2 ðUð2Þ Þ 3    N ðUðNÞ Þ :

(6)

The complete algorithm for computing the truncated HOSVD is given by Algorithm 1. The function eigs(X, r) computes the r most significant eigenvectors of the matrix X. Algorithm 1. HO-SVD

T

Y ðnÞ ¼ UðnÞ ðU
n Þ ;

(8) Ip J

. The where U
n ¼ UðnÞ . . . Uðnþ1Þ Uðn
1Þ . . . Uð1Þ 2 R p 6¼ n symbol denotes the Khatri-Rao product. If nonnegativity In J can be obtained constraints are imposed, the factor UðnÞ 2 Rþ by solving the nonnegatively constrained least squares problem: 1 T UðnÞ ¼ arg min jjYðnÞ
UðU
n Þ jj2F ; U 2

s:t: U 0

(9)

Various numerical algorithms can be applied in this case, such as that described by Cichocki et al. [3]. In the experiments, we used two methods: the well-known projected ALS and hierarchical ALS (HALS), which is regarded as one of the most efficient approaches in many applications. The rank J can be estimated in a similar manner to the number of PCs in PCA, i.e., by setting the explained variance in (2) to a given level. This task can also be achieved by using the DIFFIT method (difference in fit) [23], which is particularly useful when the curve of the explained variance exhibits a breakdown with respect to the number of latent components.

4.

The column vectors in the factor U(n) represent the feature vectors with respect to the mode n. By decomposing our fourth order tensor obtained from EMG/MMG measurements, the ð1Þ vectors fuj g contain information on the frequency profiles. ð2Þ The temporal structures are represented by the vectors fuj g. The activity of the EMG and MMG sensors can be inferred from ð3Þ the sets of the vectors fuj g. The row vectors of U(4) correspond to the grasping movements. If J4  I4, then there is strong model dimensionality reduction along the fourth mode. Therefore, each grasping movement can be represented in the low-dimensional space RJ4 . Classification can be performed in the space of the vectors ð4Þ fuj g, or in the space of the core tensor G.

Q

Classification

Let Y r 2 RI1 I2 I3 IR be a tensor containing IR training samples. Similarly, let Y t 2 RI1 I2 I3 IT contain IT testing samples. Both tensors are constructed in the same way, as described in Section 3. The training samples are assigned to C classes and each testing sample is assumed to belong to one of these classes. When using PCA for mapping from the high-dimensional observation subspace to the low-dimensional latent component subspace, both tensors Y r and Y t must be unfolded along ðrÞ the fourth mode to the corresponding matrices Y ð4Þ 2 RIR I1 I2 I3 ðtÞ IT I1 I2 I3 and Yð4Þ 2 R . After centralization, we obtain the matrices ðrÞ ðtÞ Y ð4Þ 2 RIR I1 I2 I3 and Y ð4Þ 2 RIT I1 I2 I3 . The feature eigenvectors are ðrÞ ðrÞ then extracted only from Yð4Þ . Let VJ 2 RI1 I2 I3 J comprise J ðrÞ eigenvectors of the covariance matrix for Yð4Þ . By applying the ðrÞ orthogonal mapping with respect to the basis V J to both datasets, we obtain the matrices ZðrÞ 2 RIR J and ZðtÞ 2 RIT J , which contain the training and testing samples in the J-dimensional subspace, respectively. Each row vector from Z(t) can be readily classified using any classifier, e.g., k-nearest neighbor (k-NN), linear discriminant analysis (LDA) or support vector machine (SVM). In our experiments, we used 1-NN with the distance metrics determined by the cosine measure.

biocybernetics and biomedical engineering 37 (2017) 103–113

After applying the Tucker decomposition to the tensor Y r , ð4Þ we obtain the set of factor matrices fUð1Þ r ; . . .; Ur g and the core ðrÞ ðtÞ tensor Gr . Let Gð4Þ 2 RJ4 J1 J2 J3 and Yð4Þ 2 RIT I1 I2 I3 be the unfolded versions of Gr and Y t , respectively. The multi-way samples in Y t should have the same multi-linear features as the multiway samples in Y r (the samples come from the same observation subspace), thus, classification can be performed in two ways. In the first approach, we may assume that the observation samples belonging to the same class can be represented by very similar coefficients of a linear combination of multi-linear features. If this is the case, the classification process can be performed in the space RJ4 , which contains the row vectors of ð4Þ ð4Þ the factors Uð4Þ r and Ut . The factor Ut containing the vectors of a linear combination of the multi-linear features can be estimated from (4). By formulating the normal equations, the ð4Þ estimate for Ut can be calculated as: ð4Þ

Ut

ðtÞ

 
1 ðrÞ T ðrÞ ðrÞ T
4 T
4 Gð4Þ ðU Þ U ðGð4Þ Þ ; r r


4 ¼ Y ð4Þ U ðGð4Þ Þ r

Q

Ip 

Q

(10)

J

p 6¼ 4 p 6¼ 4 p
4 ¼ Uð3Þ     Uð1Þ 2 R where U r  r r 
1 . Note that the ðrÞ T ðrÞ ðrÞ T 
4 
4 T 
4 can be precommatrix Ur ðGð4Þ Þ Gð4Þ ðUr Þ Ur ðGð4Þ Þ ð4Þ puted, and the computation of each row vector of Ut only has computational complexity of O(I1I2I3). The row vectors of ð4Þ Ut can be classified with respect to the classes of the row vectors in Uð4Þ r by using similar classifiers to those mentioned above. In the following, this classification method is referred to as HO-SVD(1).

In the other approach, the core tensors can be regarded as objects containing relevant features. Previous studies by [3,19], as well as by [27], reported that the classification process can be performed directly in the space of the core tensors. Thus, we calculate partial cores that are reduced by all but one mode. These tensors are computed as follows: T

T

T

T

T

T

Z ðrÞ ¼ Y r 1 ðUð1Þ Þ 2 ðUð2Þ Þ 3 ðUð3Þ Þ 2 RJ1 J2 J3 IR ; Z ðtÞ ¼ Y t 1 ðUð1Þ Þ 2 ðUð2Þ Þ 3 ðUð3Þ Þ 2 RJ1 J2 J3 IT :

(11) (12)

ðrÞ Z in ¼k

Each subtensor for k = 1, . . ., IR contains the latent components related to the kth training grasping movement. Note that the column vectors in the factors {U(1), U(2), U(3)} have the same permutation order and scaling; hence, there is no problem with permutation and scaling ambiguity in the tensors Z ðrÞ and Z ðtÞ . After unfolding both tensors with respect to the fourth ðrÞ ðtÞ mode, we obtain the matrices Zð4Þ 2 RIR J1 J2 J3 and Zð4Þ 2 RIT J1 J2 J3 , which contain the training and testing latent components, respectively. If the ranks J1, J2, and J3 are sufficiently small, then all of the aforementioned classifiers can be used. This dimensionality reduction technique is called HO-SVD(2). If nonnegative CP decomposition is used for model dimensionality reduction, classification can be performed in ð4Þ a similar manner to HO-SVD(1). In this case, the factor Ut is determined from the nonnegatively constrained least square problem: ð4Þ

Ut

1 ðtÞ
4 T 2 ¼ arg min jjY ð4Þ
UðU Þ jjF ; r U 2

s:t: U 0;

(13)

109

Q I J p 6¼ 4 p .
4 ¼ Uð3Þ    Uð1Þ 2 R where U The factors r r r ð4Þ fUð1Þ r ; . . .; Ur g are obtained by factorizing Y r with the nonnegative CP decomposition. Classification can be performed in the J-dimensional subspace in a similar manner to that described above.

5.

Experiments

In this section, we present empirical comparisons of the dimensionality reduction methods discussed in Section 3 in terms of their efficiency at classifying grasping movements.

5.1.

Setup

In the experiments, we used the surface EMG/MMG signals measured using the original measurement system 1, which allowed us to simultaneously collect EMG and MMG signals from eight sensors placed on the forearm. Two subjects (healthy males aged from 20 to 30 years) participated in the experiments. Each subject was asked to perform 200 repetitions of each of 11 grasping movements with various objects, i.e., a pen, small screw, shaft of a potentiometer, cylindrical glass, cup with a handle, kettle, keycard, mobile phone, computer mouse, and a suitcase positioned horizontally. The mobile phone was moved in two different ways, with and without supination. Examples of these movements are illustrated in Fig. 1. Note that the images in Fig. 1(h) and (i) represent two phases of one grasping movement. Thus, 11 different grasping movements had to be classified correctly based on the EMG/MMG observations. The measurement electrodes were arranged in two rows (upper and lower side of the forearm) with four electrodes in a row, where they formed a straight line running from the thumb to the elbow and from the small finger to the elbow. The raw temporal signals were transformed directly into 2D representations using the STFT or DWT. Preprocessing was not applied to the raw data. Thus, we obtained one multi-way array according to the rules given in Section 3 for each modality (EMG and MMG). In the case of STFT, we formed the multi-way array based on the log-scale power spectral density of the signals observed in each segment. The logarithmic scale was necessary to emphasize the low-energy components, especially in the high frequency range. Next, the 2D representations were locally averaged using moving Gaussian filters. When DWT was used, the scalograms had the same temporal resolution as the original temporal signals, and thus the corresponding tensor had many more entries along the second mode. We experimentally verified that the accuracy of scalogram classification was comparable to that of the spectrograms, but there was a large difference in computational complexity, even after decimation of the second-order to temporal signals. Hence, due to space limitations, we only present the results obtained using the STFT. We analyzed two examples of multi-channel signals. First, each channel had its own 2D representation, which yielded a 1 The measurement system was designed and developed in the Department of Cybernetics and Robotics, Wroclaw University of Science and Technology, Poland.

110

biocybernetics and biomedical engineering 37 (2017) 103–113

Fig. 1 – Examples of grasping movements using the following objects and their assignments to classes: (a) pen – class 1, (b) screw – class 2, (c) shaft of a potentiometer – class 3, (d) cylindrical glass – class 4, (e) cup – class 5, (f) kettle – class 6, (g) keycard – class 7, (h) mobile phone with supination (rotation upward – lift phase) – class 8, (i) mobile phone with supination (rotation upward – rotation phase) – class 8, (j) computer mouse – class 9, (k) suitcase positioned horizontally – class 10, (l) mobile phone without supination (rotation) – class 11.

4D tensor. In this model, the first and second modes of the tensor referred to the frequency and temporal profiles, respectively, and the third mode represented the channels. The fourth mode corresponded to the grasping movements. In the other case, the channels were concatenated along the second mode, which yielded a 3D tensor with more entries in this mode. We compared the following model dimensionality reduction methods: PCA (Section 3.1), NMF with the HALS algorithm [3], CP decomposition with the projected ALS and HALS algorithms (Section 3.3), and the two versions of HO-SVD (see Sections 3.2 and 4). In PCA, setting J = 20 gave about 80% of the explained variance. Obviously, each model with reduced dimensionality can only give a rough fit to the observed data, so exact factorization was not possible. Thus, the rank of factorization or decomposition in each model could be selected only as a tradeoff. Many methods can be used for estimating the rank of each decomposition method. To allow comparisons with PCA, we set J = 20 for NMF and CP factorization, and J = [10, 10, 8, 20] for the HO-SVD methods. For iterative methods such as NMF or CP decomposition, the iterations were terminated when the relative residual error dropped below the threshold 10
5 or the maximum number of iterations exceeded 100. Classification was performed according to the rules described in Section 4 and using the fivefold cross-validation technique. Moreover, each fold was repeated 10 more times,

which was especially important for NMF due to its non-convex optimization. The classification results were evaluated quantitatively using various measures. Fig. 2 shows the average classification accuracy and standard deviation (marked by the whiskers) for subject 1. The two modalities (EMG and MMG) were analyzed separately using all of the aforementioned dimensionality reduction methods. Fig. 2(a) and (b) show the 3D and 4D data models, respectively. All of the algorithms were coded in Matlab 2012 and run on a computational server equipped with two CPUs [Intel Xeon(R) X5650, 2.66 GHz]. Their computational complexity was empirically validated based on the average runtime, as shown in Fig. 3, for the same measurement scenarios described above. Fig. 4 shows Hinton diagrams of the confusion matrices for the classification of grasping movements performed by both subjects. Only the EMG data were used, and dimensionality reduction was performed using only the HO-SVD(2) algorithm. In the titles, P denotes the average accuracy.

5.2.

Discussion of the results

The results illustrated in Fig. 2 show that the EMG signals could be classified with higher accuracy than the MMG signals. This difference might be explained by many reasons but the acoustic disturbance in the MMG signals (mentioned in Section 2.2) appeared to affect the accuracy considerably.

biocybernetics and biomedical engineering 37 (2017) 103–113

111

Fig. 3 – Runtime of the training process with various dimensionality reduction methods (EMG data, subject 1).

Fig. 2 – Average classification accuracy (standard deviations are denoted by the whiskers) for the grasping movements performed by subject 1 (11 classes). Various dimensionality reduction methods were applied to the following models: (a) 3D and (b) 4D.

Greater accuracy was also obtained when we used the 3D model, i.e., when the channels were concatenated. This observation is rather surprising and it suggests that the 2D representations of the signals received simultaneously from all of the channels for a given grasping movement were not strongly correlated. Hence, their concatenation along the second mode yields informative and discriminant temporal profiles. Indeed, when J3 = 1 (i.e., averaging along the third mode) in the HO-SVD, the classification results were considerably worse. This lower performance may also have been due to a numerical problem, which was particularly clear for CP decomposition using the ALS algorithm (see Fig. 2). A simple projection onto the nonnegative orthant may have resulted in the non-monotonic convergence behavior, which had a greater effect on the quality of the multi-linear features when a higher-order tensor was processed. After comparing the

algorithms, we found that HO-SVD(2) and PCA gave the most accurate classifications of the EMG signals. The 3D model was also preferable in terms of its computational complexity, as confirmed by Fig. 3. This observation is fully justified by theoretical considerations because one mode (corresponding to channels) in the 4D tensor is absorbed into the covariance matrix of the second mode in the 3D model. By applying decimation of the secondorder to temporal instances in the 3D model, the length of the samples in the second mode was reduced by about half. Hence, the covariance matrix (particularly in PCA) in the 3D model was half the size of that for the 4D case, which affected the runtime considerably. Similar differences were found for the covariance matrices obtained using HO-SVD. For the 4D model, the covariance matrix with PCA was of order I1I2I3, whereas HO-SVD yielded three matrices of orders I1, I2, and I3. Therefore, it is obvious why the HO-SVD algorithms operated much faster than PCA, as confirmed empirically in Fig. 3. After comparing all of the results, we may conclude that the fastest algorithms belong to the HO-SVD family. In terms of overall efficiency, the HO-SVD(2) algorithm provided the best results with respect to the accuracy and computational complexity. This observation also suggests that the reduction of the model's dimensionality does not need to be as great as usually assumed. The k-NN classifier had a low computational cost and for this application, it was better to reduce the computational complexity of the HO-SVD model than the classifier. Thus, further experiments were performed only using the HO-SVD(2) algorithm. Fig. 4 illustrates the classification accuracy for the grasping movements performed by both subjects, which were measured using the EMG modality. The result suggests that subject 2 was better trained or more concentrated during the experiments. The main difficulties were with the classification of grasping movements 2 and 3, i.e., grasping a screw and the shaft of a potentiometer. Indeed, both grasping movements appeared to be very similar and they involved similar muscles. However, subject 2 could still perform this exercise better than subject 1. Other errors that led to the misclassification of EMG signals occurred rarely.

112

biocybernetics and biomedical engineering 37 (2017) 103–113

measurements. The STFT was used for feature extraction because it had much lower computational complexity (in the batch mode) than DWT. Due to the computational cost and accuracy, we also suggest concatenating the temporal signals obtained from all channels, before subjecting them to decimation. This approach was more efficient for representing the channels than using a separate mode in a multi-way array of observations. We also compared the EMG and MMG modalities, which showed that the classification quality was much better for the EMG signals. We also attempted to combine both modalities. The EMG and MMG signals were registered simultaneously using double-function sensors, so they should have represented the same muscle activity. In one trial, both modalities were represented by separate modes in the multi-way array of spectrograms. When the channels were not concatenated, the combined modalities yielded a 5D tensor, which could be decomposed in a similar manner to that described in Section 3. In the other trial, higher-order orthogonal iterations [3] were used to combine the normalized features. In both cases, the results were not as good as those obtained using only the EMG modality. Further research will be performed to develop better tools for combining both modalities.

Acknowledgements The main research was partially supported by the grant 2015/ 17/B/ST6/01865, funded by National Science Center (NCN) in Poland. The other part, related to preparation of the measurements, was supported by the statutory research grant from the Ministry of Science and Higher Education.

references

Fig. 4 – Hinton diagrams of the confusion matrices obtained by the HO-SVD(2) algorithm applied to the EMG data: (a) subject 1; (b) subject 2.

6.

Conclusions

This study considered the control of a prosthetic hand by EMG/ MMG signals, where our experiments focused on various aspects of the classification of these signals, which were registered during the performance of various grasping movements. We compared several dimensionality reduction methods in terms of their efficiency with this application. We proposed a new classification model based on HO-SVD and experiments demonstrated that it was the most efficient method, with the highest classification accuracy, i.e., exceeding 95% for all the subjects. The experiments also showed that the HO-SVD algorithms had the lowest computational complexity. We also performed other comparisons of the modalities, feature extraction methods, and models representing the

[1] Alkan A, Gunay M. Identification of EMG signals using discriminant analysis and SVM classifier. Expert Syst Appl 2012;39(1):44–7. [2] Carroll JD, Chang JJ. Analysis of individual differences in multidimensional scaling via an n-way generalization of Eckart-Young decomposition. Psychometrika 1970;35:283– 319. [3] Cichocki A, Zdunek R, Phan AH, Amari S-I. Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. Wiley and Sons; 2009. [4] De Lathauwer L, de Moor B, Vandewalle J. A multilinear singular value decomposition. SIAM J Matrix Anal Appl 2000;21:1253–78. [5] Gokgoz E, Subasi A. Comparison of decision tree algorithms for EMG signal classification using DWT. Biomed Signal Process Control 2015;18:138–44. [7] Harshman RA. Foundations of the PARAFAC procedure: Models and conditions for an ‘‘explanatory’’multimodal factor analysis.UCLA working papers in phonetics. vol. 16. 1970;p. 1–84. [8] Jiang N, Englehart KB, Parker PA. Extracting simultaneous and proportional neural control information for multipleDOF prostheses from the surface electromyographic signal. IEEE Trans Biomed Eng 2009;56(4):1070–80.

biocybernetics and biomedical engineering 37 (2017) 103–113

[9] Jolliffe IT. Principal component analysis. Springer series in statistics. 2nd ed. Springer; 2002. [10] Josse J, Husson F. Selecting the number of components in principal component analysis using cross-validation approximations. Comput Stat Data Anal 2012;56(6):1869– 79. [11] Kim S, Kim M, Lee J, Park J. Robot hand synergy mapping using multi-factor model and EMG signal. Proc. 14-th International Symposium on Experimental Robotics (ISER 2014), volume 109 of Springer Tracts in Advanced Robotics; 2016. pp. 671–83. [12] Kolda TG, Bader BW. Tensor decompositions and applications. SIAM Rev 2009;51(3):455–500. [13] Kurzynski M, Krysmann M, Trajdos P, Wolczowski A. Multiclassifier system with hybrid learning applied to the control of bioprosthetic hand. Comput Biol Med 2016;69 (1):286–97. [14] Lee DD, Seung HS. Learning the parts of objects by nonnegative matrix factorization. Nature 1999;401:788–91. [15] Lee J, Kim J, Park J. Classification of grip configuration using surface EMG. Proc. 13th International Conference on Control, Automation and Systems (ICCAS 2013); 2013. pp. 671–83. [16] Lucas M-F, Gaufriau A, Pascual S, Doncarli C, Farina D. Multi-channel surface EMG classification using support vector machines and signal-based wavelet optimization. Biomed Signal Process Control 2008;3(2):169–74. [17] Niegowski M, Zivanovic M. ECG-EMG separation by using enhanced non-negative matrix factorization.. Proc. 36-th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). August 2014. pp. 4212–5. [18] Oliveira AS, Gizzi L, Farina D, Kersting UG. Motor modules of human locomotion: influence of EMG averaging,

[19]

[20]

[21]

[22]

[23]

[24] [25]

[26]

[27]

[28]

113

concatenation and number of gait cycles. Front Hum Neurosci 2014;8(335). Phan AH, Cichocki A. Tensor decompositions for feature extraction and classification of high dimensional datasets. IEICE Nonlinear Theory Appl 2010;1(1):37–68. Phinyomark A, Phukpattaranont P, Limsakul C. Feature reduction and selection for EMG signal classification. Expert Syst Appl 2012;39(8):7420–31. Stephen M. EMG-EMG coherence analysis on the elbow and shoulder muscles (Master's thesis). Saint Louis, Missouri: Washington University in St. Louis; December 2013. Theis FJ, Garcia GA. On the use of sparse signal decomposition in the analysis of multi-channel surface electromyograms. Signal Process 2006;86(3):603–23. Timmerman ME, Kiers HAL. Three mode principal components analysis: choosing the numbers of components and sensitivity to local optima. Br J Math Stat Psychol 2000;53(1):1–16. Tucker LR. Some mathematical notes on three-mode factor analysis. Psychometrika 1966;31:279–311. Ulfarsson MO, Solo V. Selecting the number of principal components with SURE. IEEE Signal Process Lett 2015;22 (2):239–43. Wolczowski A, Kurzynski M. Human-machine interface in bioprosthesis control using EMG signal classification. Expert Syst 2010;27(1):53–70. Xie P, Song Y. Multi-domain feature extraction from surface EMG signals using nonnegative tensor factorization. Proc. 2013 IEEE International Conference on Bioinformatics and Biomedicine; 2013. pp. 322–5. Yazama Y, Mitsukura Y, Fukumi M, Akamatsu N. Recognition system for EMG signals by using non-negative matrix factorization. Proc. International Joint Conference on Neural Networks, vol. 3; 2003. pp. 2130–3.