Multivariate Analysis of Surface Electromyography Signals by Joyce Hsien-yin Chiang B.A.Sc , The University of British Columbia, 2004 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF Master of Applied Science in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University of British Columbia August, 2007 © Joyce Hsien-yin Chiang 2007 Abstract As the primary method of measuring muscle activation, the surface electromyography (sEMG) is of great importance in the study of motor deficits seen in patients with brain injuries and neuromuscular disorders. While clinicians have long intuitively understood that deficits in motor control are related to inappropriate recruitment of muscle synergies across several muscles, sEMG recordings are still typically examined in a univariate fash-ion. However, most traditional univariate techniques are unable to quantitatively capture the complex interactions between muscles during natural movements. To address this is-sue, multivariate signal processing techniques are employed in this thesis to study muscle co-activation patterns in patient populations. A method for classification of multivariate sEMG recordings between stroke and healthy subjects is proposed. The proposed classification scheme utilizes the eigenspectra of time-varying covariance patterns between sEMG channels as feature vectors and the support vector machines (SVM) as classifiers. Despite the minimal differences in the RMS profiles of individual muscles, the proposed scheme is able to effectively differentiate between healthy and stroke subjects. Moreover, the classification rate is shown to be monotonically related to the severity of motor impairment. This simple, biologically-inspired approach is able to quantitatively capture the subtle differences in muscle recruitment patterns between two populations and appears to be a promising means to measure motor performance. The other approach to modeling multivariate sEMG utilizes the HMM-mAR framework, which combines hidden Markov models (HMMs) and multivariate autoregressive (mAR) models. Different forms of sEMG data are analyzed, including raw sEMG, amplitude sEMG and carrier sEMG. The classification between healthy and stroke subjects is performed ii using structural features derived from estimated model parameters. Both the raw and carrier data produce excellent classification performance. The proposed method represents a fundamental departure from most existing classification methods where only amplitude sEMG is analyzed or mAR coefficients are directly used as feature vectors. In contrast, our analysis shows that the structural features of the carrier sEMG can enhance the classification performance and provide additional insights into motor control. iii Table of Contents Abstract ii Table of Contents iv List of Tables vii List of Figures viii Acknowledgements xi 1 Introduction 1 1.1 Background and Motivation 1 1.2 Background on Surface Electromyography 3 1.3 Mathematical Representation of sEMG signals 5 1.4 Literature Survey 6 1.4.1 Traditional Techniques 6 1.4.2 Multivariate Time-series Analysis Techniques 8 1.5 Research Objectives and Methodology 10 1.6 Bibliography 13 2 Time-varying Eigenspectrum 15 2.1 Introduction 15 2.2 Methods 17 2.2.1 Feature Extraction 17 iv 2.2.2 Pattern Classification 18 2.3 Results and Discussion • • • 19 2.3.1 Experiment 1: Stationary Targets 20 2.3.2 Experiment 2: Looming Virtual Targets 21 2.4 Conclusion 22 2.5 Bibliography 27 3 HMM-mAR Model 28 3.1 .Introduction 28 3.2 Methods 31 3.2.1 HMM-mAR model 31 3.2.2 Parameter Estimation via E M Algorithm 32 3.2.3 Amplitude-Modulated sEMG Representation 34 3.2.4 Analysis of Residuals and Carrier Signals 35 3.2.5 Construction of Muscle Networks and 3-node Subnets 36 3.2.6 Feature Extraction and Classification 37 3.3 Results 39 3.3.1 Data Collection and Preprocessing 39 3.3.2 Fitting of HMM-mAR 40 3.3.3 Muscle Network 42 3.3.4 Classification Results 43 3.4 Further Discussion .„ 45 3.5 Conclusion 46 3.6 Bibliography 58 4 Supplemental Analysis for Windowed Eigenspectrum Method 60 4.1 Introduction 60 4.2 Methods 60 v 4.2.1 Whitening of Raw sEMG Data . . . 60 4.2.2 Analysis of Amplitude and Carrier Components 62 4.2.3 Cross Validation 62 4.3 Results 62 4.4 Conclusion 63 4.5 Bibliography , 67 5 Conclusion 68 5.1 Conclusions 68 5.2 Future Work 69 Appendices A Derivation of EM for HMM-mAR Model 70 vi List of Tables 2.1 Classification rates for different severities of impairment using various pre-processing techniques 24 3.1 Cross-subject Classification Error Rates Based on Structural Feature Sets . . 47 3.2 Cross-subject Classification Error Rates Using Combined Structural Feature Sets 47 3.3 Cross-subject Classification Error Rates Using A R Coefficients 47 3.4 Cross-subject Classification Error Rates between Different Handiness and Motor Impairment using Structural Feature Sets Extracted from A R Coeffi-cient Matrices 48 3.5 Cross-subject Classification Error Rates between Different Handiness and Motor Impairment using Structural Feature Sets Extracted from Covariance Matrices 48 4.1 Classification rates for different severities of impairment using different forms of whitened data 65 4.2 Classification rates for different severities of impairment using different forms of unwhitened data 65 vii List of Figures 1.1 Motor Neuron and Its Innervating Muscle Fiber (based on [7]) 12 1.2 Illustrations of sEMG Signal Models (based on [9]) 12 2.1 Covariance of sEMG recoridngs. A histogram of g(X) (Eqn. 2.1), where X is a number of channels by timepoint matrix of sEMG recordings. Note that in both healthy and stroke subjects, a significant cross-correlation (g{X) > 0) between muscles can be detected 24 2.2 (a): Time-varying eigenspectral patterns obtained from non-paretic side(solid line) and paretic side(dotted line) of a stroke subject over five trials. Note the consistency across reaching trials, (b): Moving-RMS profiles obtained from non-paretic side(solid line) and paretic side(dotted line) of the same subject. 25 2.3 Top panel (a): Histogram of decision values for Healthy (open circles), Se-vere (solid line), Moderate (dashed line), and Mild (dotted line). Bottom panel (b): The linear relationship between classification decision values and the severity of stroke 26 2.4 Mean eigenspectral vectors between control(solid line) and experimental con-ditions(dashed line with open circle) 26 3.1 Three forms of sEMG signal: raw sEMG data, amplitude data, and carrier data 49 3.2 Flowchart summarizing processing steps performed on various form of sEMG data • 50 viii 3.3 The group average muscle networks extracted from A R coefficient matrices of the State-1 HMM-mAR model fitted to raw data. Only top 20 edges were selected. Frequency of occurrences of an edge is reflected by its width, (a): Average muscle network of the healthy subjects, (b): Average muscle network of the stroke subjects 51 3.4 The group average muscle networks extracted from residual covariance matri-ces of the mAR model fitted to carrier data. Only top 12 edges were selected. Frequency of occurrences of an edge is reflected by its width, (a): Average muscle network of the healthy subjects, (b): Average muscle network of the stroke subjects 52 3.5 (a): Labels for directed triple subnetworks, (b): Labels for undirected triple subnetworks 53 3.6 Model order selection using the Schwarz's Bayesian Criterion versus model order P. The optimal order is the one which gives the lowest SBC, and it is around P = 4 in this case 53 3.7 Fitting a 4-th order HMM-mAR model to raw sEMG data obtained from stroke subject over five consecutive reaching trials. The top panel, (a), shows the estimated state sequence of a 2-state hidden Markov chain. The bottom panel, (b), shows the envelopes of the raw sEMG data (solid line) and the residual of fitted HMM-mAR model (dotted line) . 54 3.8 Fitting HMM-mAR model to the amplitude component of sEMG data. The top panel, (a), shows the estimated state sequence which is much more sta-ble than the raw data case as shown in Fig. 3.7. Moreover, the estimation performance has greatly improved as indicated by the significant reduction in the magnitude of the residual (dotted line in the bottom panel) 54 ix 3.9 The best classification trees with edges as classification features, (a): Classi-fication tree based on A R coefficient matrix of rriAR model fitted to carrier data. Top 20 elements of the matrix were selected to construct the muscle connectivity network, (b): Classification tree based on residual covariance matrix of mAR model fitted to carrier data. Top 12 elements of the matrix were selected 55 3.10 The best classification trees with triples as classification features, (a): Clas-sification tree based on A R coefficient matrix of State-1 HMM-mAR model fitted to raw data. Top 20 elements of the matrix were selected, (b): Clas-sification tree based on residual covariance matrix of mAR model fitted to carrier data. Top 12 elements of the matrix were selected 56 3.11 Fitting HMM-Gaussian to the detrended residual. The top panel, (a), shows the estimated state sequence. The bottom panel, (b), shows the power spec-trum of posterior likelihood of the estimated data. Notice the significant increase in power around 20 Hz (indicated by the arrow) 57 4.1 Illustration of leave-one-subject-out cross validation scheme 65 4.2 Average eigenspectral vectors extracted from whitened data of different sub-ject groups. The vertical error bars represent the standard deviation 66 x Acknowledgements I would like to take this opportunity to convey my great appreciation to my supervisors, Dr. Jane Wang and Dr. Martin McKeown, whose patient guidance, persistent encouragement and deep insight in the research area have helped me immeasurably throughout the course of my thesis research. This thesis would never have been written without their assistance. I would like to thank my thesis examination committee members for their time and efforts. I would also like to acknowledge Dr. Janice Eng for providing us with the data and helpful discussions regarding current analysis techniques. I am deeply indebted to my family for their endless love, constant support, and immense encouragement over the years. This work was supported by Vancouver Coastal Health Research Institute(VCHRI) In-terdisciplinary Grant. xi Chapter 1 Introduction 1.1 Backg round and M o t i v a t i o n Recent advances in electronics and signal processing have allowed the wide use of biomedical signals such as electroencephalograms (EEG), electrocardiograms (ECG), and electromyo-grams (EMG) in both clinical and research settings. E M G is the study of muscle functions through the electrical signals that muscle cells generate during contractions [1]. E M G has been widely used in rehabilitation, ergonomics, and prostheses as it provides clinicians and researchers with a means to assess motor performance and also an opportunity to gain deeper understanding of the motor control mechanisms. In the past, the acquisition of E M G signals involved the insertion of needle electrodes through the skin into the muscle tissues. This operation can only be performed by trained medical professionals [2]. As the E M G instrumentation continues to improve, the surface recording technique emerged. Surface E M G (sEMG) has become more widely used as it provides an easy, non-invasive way to measure muscle activity. In addition, sEMG allows simultaneous recording of a large number of muscles and thereby provides additional spatial information which is not easily accessible by needle E M G [2]. A longstanding question in the study of motor control is how the central nervous sys-tem (CNS) efficiently coordinates the large number of muscles to produce various complex movements [3]. Several studies on the organization of the motor system have hypothesized that the CNS simplifies this complicated control problem by organizing the motor system into a modular, hierarchical structure, which is comprised of a set of basis movements called 1 "muscle synergies". Muscle synergies are considered the smallest functional units in a motor task, and they are realized by the coordinated activations of a group of, possibly spatially distributed, muscles. Therefore, instead of controlling individual muscles separately, the CNS constructs complex movements by recruiting and coordinating appropriate muscle synergies. Based on this observation, recent research on sEMG has started to focus on the characterization of co-activation patterns in synergistic muscles. Yet traditionally, simul-taneously recorded sEMG signals are analyzed individually, in a univariate fashion. While these univariate techniques provide valuable information about activity of each individual muscle, they are unable to quantitatively capture the complex interaction between muscles during natural movements. To address this issue, recent research has explored multivari-ate time-series analysis techniques in an effort to devise new methods for studying muscle synergies. Pattern classification of sEMG signals has been an active research topic. It has been applied in many areas such as functional electrical stimulations [4], prosthetic device control [5], and the study of motor impairment after spinal cord injuries [6]. A typical classification scheme consists of two components: a feature extraction algorithm and a classifier. Over the last few decades, numerous classifiers have been proposed in the literature, such as self-organizing maps, support vector machines (SVM), and neural networks. On the other hand, extracting appropriate features from sEMG signals that can succinctly capture the unique characteristics of each class is a challenging task. This problem is further complicated by inter-subject variability when cross-subject classification is performed. The main objective of this research is to develop new classification schemes for discriminating subjects with different severities of motor impairment using sEMG signals. The extracted features should be invariant to inter-subject variability within a group, yet still sensitive to detect the subtle differences between different patient populations. The rest of this chapter is organized as follows: Section 1.2 and 1.3 present an intro-duction to the E M G signals and their mathematical representations. This is followed by a 2 literature survey of various sEMG analysis techniques, including univariate and multivariate techniques, in Section 1.4. Finally, the research objective and methodology are presented in Section 1.5. 1.2 Background on Surface Electromyography Movements are generated by a combination of the spatial and temporal patterns of muscle contractions coordinated by the CNS. A muscle is comprised of many small muscle fibers. These muscle fibers are controlled by a set of lower motor neurons located in the spinal cord and brainstem. A motor neuron and the group of fibers it innervates form a motor unit (MU) which is the smallest functional entity that can be activated to produce movement. The contraction of a motor unit is initiated by an electrical pulse originating from the cor-responding motor neuron. As the pulse propagates along the axon of the motor neuron and reaches the neuromuscular junction, action potentials (AP) are induced in the innervated muscle fibers and consequently, the M U contracts (see Fig. 1.1). The combination of APs from the muscle fibers in a contracting M U is called a motor unit action potential (MUAP). To sustain the force generated by muscle contractions, the M U is repeatedly stimulated by the CNS and a train of MUAPs is produced. The number of recruited MUs, together with the rate at which they are stimulated, determines the amount of force generated [7]. E M G measures the electrical activity generated by active MUs during muscle contrac-tions. There are two ways of acquiring E M G signals. One of them is needle E M G which involves insertion of needle electrodes through the skin into the muscle tissue. Needle E M G has a small pickup range which enables the detection of individual MUAPs generated by a single muscle fiber. The most common use of the needle E M G is to diagnose disorders of peripheral nervous system and muscle. However, due to its invasive nature, needle E M G is less commonly used in non-clinical settings. The other method of measuring muscle ac-tivity is surface E M G , which is the focus of this thesis. sEMG measures the activity of 3 superficial muscles with surface electrodes placed on the skin surface. sEMG signal is a summation of electrical activity of active MUs within the detection range of the surface electrode. sEMG is more commonly used than needle E M G because it provides a relatively simple, non-invasive way to assess muscle activity. Furthermore, sEMG allows simultaneous recording of the activity of a large group of muscles and thereby provides additional spatial information about muscle recruitment [2]. Despite the various benefits offered by sEMG, the major problem plaguing sEMG signals is the possibility of crosstalk. sEMG crosstalk refers to the electrical signal detected over the muscle of interest that is actually generated by nearby muscles. The most common countermeasure is to use differential or bipolar electrodes [8]. They work on the principle that any signal (such as crosstalk) originating from distant muscles will appears as a common signal to both adjacent electrodes. As the recorded signals are subtracted, the common signals including crosstalk and other ambient noise will be canceled out. The recorded sEMG signals are continuous waveforms with amplitude ranging from 0 to 6 millivolts (mV) peak to peak. For analysis purposes, the sEMG signals are first amplified and filtered by an analog bandpass filter. The low-frequency cutoff of the bandpass filter, typically ranging between 5 and 20Hz, removes any DC offset and baseline drift caused by movement and perspiration. The high-frequency cutoff removes high frequency noise, and its range can vary from 200Hz to 1 kHz. The analog sEMG signals are converted to digital signals with a sampling rate at least two times higher than the high-frequency cutoff of the bandpass filter. The resulting digital sEMG signals (subsequently referred to as "raw" sEMG signals) are further post-processed depending on the application. A more detailed description of the acquisition of sEMG signals can be found in [1]. 4 1.3 Mathematical Representation of sEMG signals During most voluntary muscle contractions, a sufficiently large number of MUs contract in an asynchronous pattern. The contraction of each M U produces a train of MUAPs whose inter-arrival times are randomly distributed. As the MUAPs propagate and reach the skin surface, the surface electrodes pick up a noise-like interference pattern which represents the superposition of individual MUAPs. In general, the resulting sEMG signal can be mathematically modeled as a stochastic process. Two different representations have been proposed [9]. A. Static Contractions In static (i.e. constant) contractions, the number of active MUs is constant. The recorded sEMG signal y(k) can be represented as an algebraic sum of the detected M U A P trains and it can be written as N y(k) = ^Si(k)*Pi(k), (1.1) where N denotes the number of active MUs, Si(t) is an impulse train representing the electrical stimulations originating from the motor neuron of i-th M U , Pi(k) is a single AP of i-th. M U , and (*) denotes the convolution operator. Each term Sj(fc) * Pi(k) represents a M U A P train. As the number of active MUs, N, gets large, by the central limit theorem, the sEMG signal can be modeled as a Gaussian process. B. Dynamic Contractions In dynamic muscle contractions, the number of active MUs and their firing rate vary with time. Thus, the recorded sEMG signal is non-stationary. To model the time-varying characteristics of the signal, the following mathematical representation is often used [10]: y(k) = x(k)m(k) + n(k), (1.2) where m(k) is the time-varying amplitude of sEMG, which is proportional to the 5 number of active MUs and their firing rate. x(k) is a wide-sense stationary stochastic process of unit variance, and n(k) denotes the measurement noise. The resulting sEMG signal, y(k) is a zero-mean Gaussian process with time-varying variance. In this thesis, we will analyze each term of Eqn. 1.2, namely the raw sEMG, y(k), the carrier sEMG, x(k), and the amplitude sEMG, m(k). Illustrations of these two sEMG signal models are shown in Fig. 1.2. 1.4 L i te ra tu re Survey In the literature, numerous techniques have been proposed for the analysis of sEMG signals. There are two main categories: one is traditional univariate techniques which focus on the estimation of amplitude and frequency information of single-channel sEMG signals, and the other one is multivariate time-series analysis techniques which exploit the spatial correlations in multi-channel sEMG signals. This section presents a brief survey of the commonly used techniques of each category. 1.4.1 Traditional Techniques Traditional techniques for analyzing sEMG signals typically focus on the amplitude and frequency information of a single-channel sEMG signal. An extensive amount of work has been done to improve the estimation of the amplitude and frequency features. The most commonly used estimators are described below. Estimation of amplitude information In its raw form, the sEMG signal is a noise-like waveform which oscillates randomly above and below the zero value. Its amplitude is typically associated with the amount of muscle energy expended. In order to quantify this information, some form of amplitude estimation is usually applied. The most widely used amplitude estimator are average rectification value 6 and root-mean-square. Average Rectification Value (ARV): A R V represents the arithmetic mean of the rec-tified (or absolute value). s E M G over a period of time. It can be mathematically expressed as follows: , i+N ARV(i) = -J2\y(k)l (1-3) k=i where y(k) is the signal sample at time k and N is the width of the time window over which the averaging is performed. The time window is shifted across the signal by a fixed step size and this procedure results in a time series of A R V s , one for each time position of the window. Typically, in order to increase the temporal resolution of the A R V s , the step size is chosen such that the time windows are mostly overlapping. Root-mean-square (RMS): R M S , as its name suggests, computes the square root of the mean of the squares of the s E M G signal over a period of time. It can be obtained as follows: RMS{i) = \ 1 i+N k=i R M S has become more widely used because it measures the power of a s E M G signal whereas A R V has no clear physical meaning [8]. Several important variables can be derived from the amplitude information, including onset time of muscle activations, the amount of force exerted by the muscle, etc. Estimation of frequency information The frequency information of a s E M G signal is typically determined using fast Fourier transform (FFT) which decomposes the signal into its various frequency components. The resulting frequency spectrum provides information about muscle activations, such as the muscle fiber conduction velocity and dominance of fast-twitch or slow-twitch muscle fibers [7]. Two spectral parameters are commonly used to characterize the frequency spectrum: 7 mean frequency (MNF) denoted by fmean and median frequency (MDF) denoted by fmed-They are denned as follows: fmean = * ^ ~~~ and Pi = Pi, (1-5) ^i=lFi i= l i=fmed where Pi is the i-th component of the power spectrum and N is the highest frequency considered. M N F and M D F are used to monitor changes in the frequency spectrum during sustained muscle contractions. A typical example is the detection of muscle fatigue [7]. In the fatigued muscle, due to the changes in volume conduction velocity and a shift in dominance from fast-twitch muscle fibers to slow-twitch fibers, the shape of the frequency spectrum is skewed to the left. This spectral change is reflected in a decrease in MNF and MDF. Although F F T is useful for determining frequency information of the sEMG signals, its applicability is limited due to its requirement of the stationarity of the signals. Since sEMG signals are non-stationary during dynamic muscle contractions, F F T may not be suitable in this case. To address this issue, various time-frequency analysis methods have been proposed to estimate the time-varying frequency spectrum. Examples include short-time Fourier transform, wavelet transform, and Wigner-Ville distribution. A comparative study of these techniques can be found in [11]. 1.4.2 Multivariate Time-series Analysis Techniques One of the most widely used time-series techniques for analyzing sEMG signals is the au-toregressive (AR) models. The A R models have been extensively applied in various areas such as economics and speech processing. The basic idea of the A R model is that it repre-sents each sample of a time series as a weighted sum of its previous samples. The simplest form of A R models is the scalar A R (sAR) which strictly handles univariate time-series. Multivariate autoregressive models (mAR) extend sAR models to multivariate signals and 8 they can be formulated as follows: p y(k) = Z~2ary(k ~p">+ w(k)> (L6) where y(k) is an M-dimensional column vector denoting the values of the M variables at time k. The term ap denotes the M x M A R coefficient matrix at time lag p. The coefficient a,ij in ap represents the weight of variable j in determining the value of variable i after p time points. w(k) is an M-dimensional Gaussian noise with the covariance matrix R. mAR models have been widely applied to biophysiological signals mainly due of its ability to capture temporal and spatial dependence relationships. In [12, 13], the authors explored the use of sAR and mAR models to extract features from electroencephalographic signals (EEG) for the purpose of mental task classification. In [14], Harrison et al. applied mAR models to functional MRI data. Based on the concept of Granger causality, Harrison et al. established the effective connectivity between brain regions through finding non-zero elements in the estimated A R coefficient matrices. The possibility of modeling sEMG signals with A R models was first demonstrated by Graupe et al. in [15] for prosthesis control purposes. They showed that sEMG signals recorded over major muscles (biceps or triceps) approximately follow a Gaussian distribu-tion and are piece-wise (locally) stationary and, thus, A R models are well suited for analysis of sEMG signals. Furthermore, Graupe et al. exploited the spectral variations reflected in the A R coefficients to discriminate between different limb functions of interest. With a sin-gle channel of sEMG signal, the results in [15] showed an 85% of accuracy in discriminating between five limb functions. In [5], Doerschuk et al. argued that during natural move-ments, several muscles are brought into play and thus, significant information might exist in the cross-correlation between sEMG signals measured at spatially distributed muscles. By using only one channel of sEMG signal, Graupe et al.'s method was unable to utilize this spatial information and therefore had limited operational ranges. Doerschuk et al. ex-9 tended the work of Graupe et al. by exploiting both the spatial and temporal correlations in multi-channel sEMG signals and used mAR models for limb function classification. Sev-eral subsequent works adopted a similar approach to characterize sEMG signals [16], [17]. However, the validity of A R models was questioned by Sherif because of the nonstationary nature of sEMG signals during dynamic muscle contractions [18]. To address to this issue, Sherif et al. proposed modeling the sEMG signals with autoregressive integrated moving-average (ARIMA) models. Moser and Graupe [4] considered a variant of regular A R models in which the A R coefficients are time-varying and their dynamics are governed by another A R process. The other commonly used approach to modeling sEMG signals is matrix factorization. The central idea behind all matrix factorization algorithms is that they assume the dy-namics of measured signals are governed by some unobserved (latent) variables and their common goal is to estimate the latent variables in the presence of noise. Several algorithms have been explored to identify muscle synergies in raw sEMG signals. They include fac-tor analysis (FA), principal component analysis (PCA), independent component analysis (ICA), and nonnegative matrix factorization (NMF). An excellent comparative survey of these algorithms can be found in [19]. 1.5 Research Objectives and Me thodo logy In this thesis, we examine the sEMG signals collected during reaching movements from healthy subjects and subjects recovering from differing severity of stroke. The objectives of this research are: 1. to develop new classification schemes for discriminating subjects with different severity of motor impairment, and 2. to compare the muscle recruitment patterns used by stroke and healthy subjects. 10 To achieve these objectives, we explore the multivariate signal processing techniques and propose the following two approaches to characterizing sEMG signals: 1. Time-varying Eigenspectrum We exploit the spatial correlations in the multi-channel sEMG signals and present a novel feature extraction technique which utilizes the time-varying covariance pattern between muscles. The feature extraction technique is applied to the sEMG data collected from stroke and healthy subjects. The spectral feature vectors are classified with S V M classifiers based on the severity of motor impairment of the subjects. 2. Hidden-Markov model with mAR To characterize the nonstationarity of sEMG signals, we extend the aforementioned mAR models and propose combining the hidden Markov model (HMM) and the mAR model into the HMM-mAR framework for modeling sEMG signals. The proposed HMM-mAR model is applied to different forms of the sEMG signals: raw sEMG, am-plitude sEMG and carrier sEMG. The estimated model parameters including mAR coefficient matrices and residual covariance matrices are utilized to construct mus-cle connectivity networks and classification features that represent muscle coactivity. These feature vectors are classified between stroke and healthy subjects using decision trees. The rest of the thesis is organized as follows: In Chapter 2, we present a feature extrac-tion technique which exploits the time-varying covariance patterns between muscles. The performance of the proposed feature extraction technique in discriminating different severi-ties of motor impairment is evaluated. In Chapter 3, we present the HMM-mAR framework for modelings sEMG signals. We propose a classification scheme which utilizes the muscle connectivity network extracted from the estimated model parameters. Its performance in classifying between healthy and stroke subjects is evaluated. Chapter 4 contains the addi-tional analysis performed to supplement the classification scheme described Chapter 2. In Chapter 5, we present conclusions of this work along with suggestions for future work. 11 Action Potential Axon of Motor Neuron Neuromuscular Junction Muscle Fiber' Action Potential Reaches Junction ma x Figure 1.1: Motor Neuron and Its Innervating Muscle Fiber (based on 1 ". 1 1 • ( - f - r t t t J I I l_L sN(k) f = II I ) I I' PN(k) (a) Static Contractions y(k) x(k) + n(k) = m(k) (b) Dynamic Contractions y(k) Figure 1.2: Illustrations of sEMG Signal Models (based on 1.6 B ib l iog raphy [1] J . V . Basmajian and C . J . D . Luca, Muscles alive. Baltimore, M D : Wil l iams & Wilkins , 1985. [2] M . J . Zwarts and D . F . Stegeman, "Multichannel surface E M G : Basic aspects and clinical utility," Muscle Nerve, vol. 28, no. 1, pp. 1-17, 2003. [3] A . d 'Avella , P. Saltiel, and E . B i z z i , "Combinations of muscle synergies in the construction of a natural motor behavior." Nat. Neurosci., vol. 6, no. 3, pp. 300-308, 2003. [4] A . T . Moser and D . Graupe, "Identification of nonstationary models wi th application to myoelectric signals for controlling electrical stimulation of paraplegics," IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 5, pp. 713-719, 1989. [5] P. C. Doerschuk et ai, "Upper extremity limb function discrimination using E M G signal analysis," IEEE Trans. Biomed. Eng., vol. 30, no. 1, pp. 18-29, 1983. [6] R. Maksimovic and M . Popovic, "Classification of tetraplegics through automatic movement evaluation," Med. Eng. Phys., vol. 21, no. 5, pp. 313-27, 1999. [7] J . R. Cram, G . S. Kasman, and J . Holtz, Introduction to surface electromyography. Aspen, Gaithersburg, 1998. [8] C. J . D . Luca, "The use of surface electromyography in biomechanics," J. Appl. Biomech., vol. 13, no. 2, pp. 135-163, 1997. [9] K . C . M c G i l l , "Surface electromyogram signal modelling," Med. Biol. Eng. Comput., vol. 42, no. 4, pp. 446-454, 2004. [10] E . A . Clancy, S. Bouchard, and D . Rancourt, "Estimation and application of E M G amplitude during dynamiccontractions," IEEE Eng. Med. Biol. Mag., vol. 20, no. 6, pp. 47-54, 2001. [11] S. Karlsson, J . Y u , and M . Akay, "Time-frequency analysis of myoelectric signals during dynamiccontractions: a comparative study," IEEE Trans. Biomed. Eng., vol. 47, no. 2, pp. 228-238, 2000. [12] Z. A . Ke i rn and J . I. Aunon, " A new mode of communication between man and his surround-ings," IEEE Trans. Biomed. Eng., vol. 37, no. 12, pp. 1209-1214, 1990. [13] C . W . Anderson, E . A . Stolz, and S. Shamsunder, "Multivariate autoregressive models for clas-sification of spontaneous electroencephalographic signals during mental tasks," IEEE Trans. Biomed. Eng., vol. 45, no. 3, pp. 277-286, 1998. [14] L . Harrison, W . D . Penny, and K . Friston, "Multivariate autoregressive modeling of f M R I time series," Neuroimage, vol. 19, no. 4, pp. 1477-1491, 2003. [15] D . Graupe, J . Magnussen, and A . Beex, " A microprocessor system for multifunctional control of upper-limb prostheses via myoelectric signal identification," IEEE Trans. Autom. Control, vol. 23, no. 4, pp. 538-544, 1978. [16] R. J . Triolo and G . D . Moskowitz, "The theoretical development of a multichannel time-series myoprocessor for simultaneous limb function detection and muscle force estimation," IEEE Trans. Biomed. Eng., vol. 36, no. 10, pp. 1004-1017, 1989. 13 [17] X . H u and V . Nenov, "Multivariate A R modeling of electromyography for the classification of upper arm movements," Clin. Neurophysiol, vol. 115, no. 6, pp. 1276-87, 2004. [18] M . H . Sherif, R . J . Gregor, and J . Lyman, "Effects of load on myoelectric signals: the A R I M A representation," IEEE Trans. Biomed. Eng., vol. 28, no. 5, pp. 411-416, 1981. [19] M . C. Tresch, V . C . K . Cheung, and A . d 'Avella , "Mat r ix factorization algorithms for the identification of muscle synergies: evaluation on simulated and experimental data sets," J. Neurophysiol, vol. 95, no. 4, pp. 2199-2212, 2006. 14 Chapter 2 Time-varying Eigenspectrum / S V M Method for sEMG Classification 1 2.1 In t roduc t ion Pattern classification of surface electromyographic (sEMG) signals from the simultaneous recordings of several muscles has applications for neural prostheses development and diag-nosis of diseases involving the motor system [20, 21]. Despite changes in arm dynamics after stroke obvious to trained clinicians, the abnormalities in the complex interference pattern of sEMG recordings during stroke can be subtle. In this chapter, we concentrate on the classification of reaching movements, rather than, e.g. unnatural sustained contraction, in healthy and stroke subjects using sEMG data. A typical classification application is comprised of two components: features selection and decision-making algorithm (classifier). Many different classifiers have been proposed in the literature, which, in general, can be divided into two categories: data-driven and model-driven. The former is more widely used due to its simplicity and generality. It includes methods such as self-organizing maps and machine learning based schemes including K-nearest neighbors, support vector machines (SVM) and neural network analysis [22]. In this study, we use the widely-applied linear S V M algorithm since it is a powerful tool in lA version of this chapter has been published/submitted for publication. J . Chiang, Z. J. Wang, and M . J . McKeown, "A time-varying eigenspectrum/SVM method for s E M G classification of reaching move-ments in healthy and stroke subjects," in Proc. 31th IEEE Int. Conf. Acoustics, Speech, and Signal Proc. (ICASSP31), vol. 2, 2006, pp.1188-1191. J. Chiang, Z. J. Wang, and M . J . McKeown, "A windowed eigen-spectrum method for multivariate s E M G classification during reaching movements," IEEE Signal Processing Letters, submitted for publication. 15 classification and pattern recognition and has been shown to provide excellent classification performance [23]. A clear challenge is to extract appropriate features from sEMG data for subsequent clas-sification. We therefore appeal to contemporary developments in theoretical neuroscience for guidance. Recent work has suggested that complex movements are implemented by low-dimensional basis movements encoded in the spinal cord [24]. These basis movements or "muscle synergies" are distributed across several muscles. Moreover, examination of muscle synergies may provide a fruitful avenue to succinctly summarize the often complex changes in muscle activation that are seen in disease states, such as motor impairments after stroke. Yet traditionally, sEMG recordings are examined individually, in a univariate fashion. Based on the following statistic testing where X denotes the sEMG recording matrix and | • | represents 1-norm, a simple inspection of basic covariance properties of sEMG recordings (see Fig. 2.1) suggests strong correlation between muscles during natural movements. This observation motivates us to use features based on the covariance patterns. A number of factors influence the amplitude of the sEMG: exact positioning of the electrodes, movement of the muscle with respect the electrodes, the amount of subcutaneous fat, and the impedance of the skin [25]. These factors can potentially affect the reliability of amplitude-based univariate classification schemes. On the other hand, examining the relationships between muscles may prove more robust and reliable in monitoring muscle functions. In this chapter, we propose a classification scheme which extracts features from the time-varying covariance between sEMG recordings in healthy subjects and subjects recovering from differing severity of stroke. We demonstrate that the proposed method is more reliable than examining muscles individually, and is monotonically related to the severity of stroke, 16 as assessed by traditional subjective clinical scales. 2.2 M e t h o d s In this section, we propose a feature extraction technique which exploits the covariance patterns between muscles. We then describe the S V M classifier used for the sEMG classifi-cation. 2.2.1 Feature Extraction Conventionally, the classification of sEMG signals has been performed by using the linear envelope of a single channel as the input feature. One widely used processing technique for extracting the linear envelope is the moving root-mean-square (RMS) method. This feature extraction technique examines each muscle in isolation, and does not take into account the correlation between muscle recordings. As previously mentioned, the studies in muscle synergies suggest that the muscles are activated in a coordinated fashion. One of our recent studies [26] and Fig. 2.1 suggest that sEMGs recorded from spatially distributed muscles may be correlated with each other. These observations lead us to a new feature extraction method which examines the covariance patterns between sEMG recordings. The method is outlined as follows: 1. Normalization: To minimize trial-to-trial variability, we first normalize the raw sEMG signals to zero mean and unit variance. In addition, to account for the slight variability in reaching times towards the target, the normalized sEMG are resampled such that all recordings are of equal length . 2. Extraction of feature vectors: Let X denote the sEMG recording matrix for one trial whose dimension is number of channel, M, by number of timepoints, N. For each trial, a moving window of width T points is slid across all M channels with a step size of 5 points. For each time position of the moving window, the covariance matrix 17 and its largest eigenvalue, A, are computed from the data points inside the moving window. As the window slides across, a vector of eigenvalues, v, is generated and it can be written as follows v = [A(0), X{6), X(25), \{36),\(nS)} G R n + 1 , (2.2) where X(iS) denotes the largest eigenvalue of the covariance matrix calculated using the data from t — i6 to t = id + T. While a range of features based on the covariance matrix could be extracted (e.g. the entire eigenspectrum), our preliminary investigations [27] revealed that the largest eigen-value most reliably distinguished between the non-paretic side and paretic side in stroke subjects. It therefore motivated us to employ it as an efficient feature to classify healthy and stroke subjects. 2.2.2 Pattern Classification The classifier used in this study is the linear S V M which is widely applied in many different areas due to its robustness and reliability [28]. The goal of the S V M is to find a hyperplane in the input space which can maximally separate members of two classes. Specifically, the hyperplane has the following form x : w -x + b = 0, (2.3) where re is a vector in the feature space, w is a weight vector normal to the hyperplane, and b is the bias. The term w • x represents the dot product of w and x. In most cases, there are infinitely many hyperplanes which can separates two classes. However, in order to avoid the potential problem of overfitting, the optimal hyperplane is defined as the one with maximal margin of separation between two classes, ft can be found by solving the 18 following quadratic optimization problem, minimize \w 2 (2.4) subject to Vi(xi -w + b) — 1 > 0, V i where £j is i-th training sample and y* G { — 1,1} is the corresponding class label. The resulting decision function for the test sample, x, has the following form where h(x) is the decision value for x and o^ is the Lagrange multiplier associated with i-th inequality constraint in Eqn. 2.4. In the solution, only the training samples lying on the margin have non-zero a, and theses samples are called support vectors. The support vectors carry all relevant information required to specify the classifier. The linear SVMs can be easily extended to nonlinear SVMs to handle non-separable data by using kernel functions. Interested readers are referred to [28] for further details. In this study, inputs to the SVM classifier are the eigenspectral vectors v generated by applying the feature extraction technique described in Section 2.2.1 to the sEMG data collected from stroke and healthy subjects. The classification performance is evaluated using the leave-one-out cross-validation technique. 2.3 Results and Discussion In this section, we study the performance of the proposed classification scheme by examining the sEMG data collected from stroke and healthy subjects. Two different experiments were conducted. In the first experiment, subjects were asked to reach a stationary target, and the objective is to classify between healthy and stroke subjects. In the second experiment, subjects were asked to catch "looming" targets (balls coming towards them) in a simulated (2.5) 19 virtual environment. The objective is to compare the condition where only target balls were presented, and the condition where target balls were intermixed with distracter balls. The experiment protocol and classification results for each experiement are described below. 2.3.1 Experiment 1: Stationary Targets In this experiment, subjects were first seated in a chair with their hands resting on the thigh and were asked to reach and touch a fixed target upon hearing an auditory cue. The target was located in the subject's mid-sagittal plane at shoulder height, and its distance with the subject was adjusted such that the target was just within the workspace of paretic arm of a stroke subject or the non-dominant arm of a healthy subject. The reaching task was repeated five times for each subject. The electrical activities of the following seven muscles were recorded using surface elec-trodes: the anterior and lateral deltoid, long head and lateral triceps, biceps brachium, latissimus dorsi, and brachioradialis. A bipolar montage was used to minimize the effect of crosstalk. The 7-channel sEMG signals were amplified, sampled at 600 Hz, and high-pass filtered at 20Hz to reduce movement-related artifact. Pleaser refer to [29] for further details on the experimental procedures. The proposed feature extraction technique described in Section 2.2.1 is applied to the sEMG data. The resulting time-varying eigenspectral vectors from a stroke subject are shown in Fig. 2.2(a). The eigenspectral vectors are relatively consistent across trials on the same side, and the difference between the paretic and non-paretic arms is evident. On the other hand, although the lateral deltoid was reported as revealing significant differences between healthy and stroke subjects using the same data set [29], the moving-RMS profiles from this single muscle of the same subject are almost indistinguishable between the paretic and non-paretic sides, as shown in Fig. 2.2(b). To measure the correlation between the severity of stroke and the classification perfor-mance, the stroke subjects are divided into three classes based on their Fugl-Mayer (FM) 20 scores: Severe with F M score below 25, Moderate with F M score between 25 and 50, and Mild with F M score above 50. Classification is then performed comparing Severe, Moderate and Mild separately to normal subjects. The results are summarized in Table 2.1 (without preprocessing). Note that the classification rates are monotonically related to the severity of motor impairment assessed by F M scores. When all stroke and healthy subjects are classified jointly, the histogram of the classi-fication decision values computed using Eqn. 2.5 confirms the correspondence between the classification performance and severity of impairment, as shown in Fig. 2.3(a). In general, the average decision value is linearly related to severity of stroke (see Fig. 2.3(b)). This observation suggests that the proposed method may provide a quantifiable metric of motor performance. In the literature, several nonlinear sEMG preprocessing techniques have been explored [30]. Examples include the widely-applied rectification | • |, cubic function (-)3, and hy-perbolic tangent tanh(-). To evaluate their effects on the classification performance, we compare the classification rates produced by these preprocessing techniques. As shown in Table 2.1, the proposed classification method is relatively insensitive to the preprocessing techniques applied to the raw data. The results suggest that there is no universal prepro-cessing technique which can consistently improve the classification performance for all three cases. 2.3.2 Experiment 2: Looming Virtual Targets In this experiment, subjects were asked to catch looming balls rapidly approaching them in a simulated virtual environment. The experiment compared two conditions: experimental and control conditions. The experimental condition consists of distracter balls and target balls. Subjects were asked to catch only the target balls and ignore the distracter balls. The control condition consists of only target balls. Surface E M G electrodes were placed on: the anterior deltoid, biceps, triceps, trapesius, pectoralis, wrist flexors and wrist extensors. 21 Motion sensors were placed on the subject's forehead, upper arm, forearm, and back of the hand to record the trajectories of the hand motion. Detailed experimental setup can be found in [31]. When the experimental and control conditions are compared, the trajectories of the hand motion are found similar. On the other hand, the eigenspectral vectors extracted from the sEMG data reveal a noticeable difference in the sEMG co-activation patterns between the experimental and control conditions, especially during the early part of the movement as illustrated in Fig. 2.4. To assess the significance of S V M classification, p-values are computed using the permutation test described in [31]. The eigenspectral patterns are found significantly different between the experimental and control conditions (p < 0.05) in 25/29 data sets of the stroke subjects and in 33/40 data sets of the healthy subjects. The difference between the control and experimental conditions in the early part of the movement is suggested to be related to the secretion of the brain neuromodulator nore-pinephrine immediately after the decision to move. The results suggest that the proposed classification scheme is sufficiently sensitive to detect the subtle difference in muscle co-activation patterns despite the minimal differences in the RMS profiles and hand trajecto-ries. 2.4 Conc lus ion Performing a reaching task is a complex interplay between brain regions subserving motor planning, vision, attention, and motor execution, and the impairment of some of these will not be captured by the conventional clinical scale. Finding features that are relatively invariant to inter-individual differences, yet still sensitive to detect severity of impairment is a challenge. The proposed method is a first step towards this goal, as the exact combination of muscles may vary across subjects, but this may not affect the eigenvalues (as opposed to the eigenvectors) of the sEMG recordings over a specific time window. Also, by looking 22 at the eigenvalues, the results would be expected to be relatively insensitive to the exact positioning of the electrodes, a problem plaguing amplitude-based classification schemes. Moreover, since the classification decision values are monotonically related to the severity of stroke (as estimated by subjective clinical scales), it suggests that this method could be extended into a quantifiable assay of motor performance. 23 Table 2.1: Classification rates for different severities of impairment using various prepro-cessing techniques Without Preprocssing With Preprocessing |x tanh(x) Severe vs. Healthy 96.72 91.80 85.25 98.36 Moderate vs. Healthy 73.77 80.33 84.97 73.77 Mild vs. Healthy 67.21 85.25 75.41 68.85 Histogram of g(X) Figure 2.1: Covariance of sEMG recoridngs. A histogram of g{X) (Eqn. 2.1), where X is a number of channels by timepoint matrix of sEMG recordings. Note that in both healthy and stroke subjects, a significant cross-correlation (g(X) > 0) between muscles can be detected. 24 % Movement (a) Eigenspectral Patterns Moving-RMS Profile of Lateral Deltoid 3 .5 , , , , , , , , r % Movement (b) RMS Profiles Figure 2.2: (a): Time-varying eigenspectral patterns obtained from non-paretic side(solid line) and paretic side(dotted line) of a stroke subject over five trials. Note the consistency across reaching trials, (b): Moving-RMS profiles obtained from non-paretic side(solid line) and paretic side(dotted line) of the same subject. 25 Histogram of Decision Values Decision Value Average Decision Value vs. Severity of Stroke (b) Healthy Mid Stroke Moderate Stroke Severe Stroke Figure 2.3: Top panel (a): Histogram of decision values for Healthy (open circles), Se-vere (solid line), Moderate (dashed line), and Mild (dotted line). Bottom panel (b): The linear relationship between classification decision values and the severity of stroke. (0 > C Lu to Control Condition e - Experimental Condition 1 10 15 20 % Movement 25 Figure 2.4: Mean eigenspectral vectors between control(solid line) and experimental condi-tions (dashed line with open circle) 26 2.5 B ib l i og raphy [20] X . H u and V . Nenov, "Multivariate A R modeling of electromyography for the classification of upper arm movements," Clin. Neurophysiol, vol. 115, no. 6, pp. 1276-87, 2004. [21] R. Maksimovic and M . Popovic, "Classification of tetraplegics through automatic movement evaluation," Med. Eng. Phys., vol. 21, no. 5, pp. 313-27, 1999. [22] D . Kumar , N . M a , and P. Burton, "Classification of dynamic multi-channel electromyography by neural network," Electromyogr. Clin. Neurophysiol., vol. 41, no. 7, pp. 401-08, 2001. [23] T . S. Furey et al, "Support vector machine classification and validation of cancer tissue samples using microarray expression data," Bioinformatics, vol . 16, no. 10, pp. 906-14, 2000. [24] A . d 'Avel la and E . B izz i , "Shared and specific muscle synergies in natural motor behaviors," Proc. Natl. Acad. Sci. U.S.A., vol. 102, no. 8, pp. 3076-81, 2005. [25] D . Farina, R . Merlet t i , and R. M . Enoka, "The extraction of neural strategies from the surface E M G , " J. Appl. Physiol., vol. 96, no. 4, pp. 1486-1495, 2004. [26] M . J . M c K e o w n and R. Radtke, "Phasic and tonic coupling between E E G and E M G demon-strated wi th independent component analysis," J. Clin. Neurophysiol, vol. 18, no. 1, pp. 45-57, 2001. [27] J . Chiang, Z. J . Wang, and M . J . McKeown, " A time-varying eigenspectrum/SVM method for s E M G classification of reaching movements in healthy and stroke subjects," in Proc. 31th IEEE Int. Conf. Acoustics, Speech, and Signal Proc. (ICASSP31), vol. 2, 2006, pp. 1188-1191. [28] C . J . Burges, " A tutorial on support vector machines for pattern recognition," Data Mining and Knowledge Discovery, vol. 2, no. 2, pp. 121-167, 1998. [29] P. H . McCrea , J . J . Eng, and A . J . Hodgson, "Saturated muscle activation contributes to compensatory reaching strategies after stroke," J. Neurophysiol, vol. 94, no. 5, pp. 2999-3008, 2005. [30] H . Mathis and S. C. Douglas, "On the existence of universal nonlinearities for blind source separation," IEEE Trans. Signal Process., vol. 50, no. 5, pp. 1007-16, 2002. [31] Y . Wei et al, "Stimulation of endogenous Norepinephrine: motor effects of looming visual target stimuli in normal and stroke subjects," Restor. Neurol. Neurosci., 2006, under review. 27 Chapter 3 Markov mAR Modeling Framework and Network Analysis for sEMG data 2 3.1 In t roduc t ion The surface electromyogram (sEMG) signal is a stochastic process which can be represented as a zero-mean, band-limited and wide-sense stationary stochastic process (referred to as "carrier data" here) modulated by the amplitude of the sEMG signals [32]. As it is assumed that the amplitude data represent overall muscle activity, a common practice in the sEMG literature is to focus on the amplitude data, which is obtained by rectifying and low-pass filtering the raw sEMG signal, while discarding the carrier data. Only a limited number of studies were reported working on raw, unrectified sEMG signals. For instance, in [33], linear ICA was applied to noisy raw sEMG data and revealed meaningful interactions between muscles. In [34], a multivariate autoregressive (mAR) model was utilized to model multi-channel raw sEMG signals. To the best of our knowledge, no previous studies have focused on the carrier data of sEMG signals, yet our very recent study shows that the carrier data 2 A version of this chapter has been published/accepted for publication/submitted for publication. J . Chiang, Z. J . Wang, and M . J . McKeown, "Hidden Markov multivariate autoregressive (HMM-mAR) mod-eling of dynamic muscle association patterns in reaching movements," in Proc. 29th Can. Medical and Biological Soc. Conf. (CMBES29), 2006. J. Chiang, Z. J . Wang, and M . J . McKeown, "Hidden Markov multivariate autoregressive (HMM-mAR) modeling framework for surface electromyography (sEMG) data," in Proc. 29th Annu. Int. Conf. IEEE EMBS(EMBS29), 2007. J. Chiang, Z. J . Wang, and M . J. McKeown, "Markov multivariate autoregressive modeling framework and network analysis for surface E M G data," IEEE Trans. Signal Proc, submitted for publication. 28 may also be informative [35]. In this study, to have a more thorough understanding of sEMG signals, we will use a fundamentally different approach and will model different forms of the sEMG signals: raw sEMG, amplitude sEMG and carrier sEMG. A topic of increasing interest in motor control is to determine how different muscles effi-ciently collaborate or interact together during natural movements. A number of studies have suggested that complex movements are implemented by low-dimensional basis movements called muscle synergies, which are coordinated activations of a group of, possibly spatially distributed, muscles [36]. Based on this observation, several different methods have been proposed to infer synergistic action between muscles, such as Principal Component Analysis (PCA), linear Independent Component Analysis (ICA) [33] [37], and Nonnegative Matrix Factorization (NMF) [38]. A l l of these methods involve the characterization of observed sig-nals by latent variables. For example, ICA assumes that the observations are combinations of statistically independent components, and its objective is to find the underlying sources. Thus, the multivariate sEMG data can be represented by subspace projections. A common assumption implied in the above signal models is that the muscle activity patterns are stationary during the recording period. However, it is widely accepted that the sEMG signals recorded during natural movements, as opposed to isometric contractions, are non-stationary. For example, reaching movements are typically considered to consist of an initial phase, under an "open loop" control, where feedback is less critical, and a "closed loop" control phase, where feedback is important [39]. In particular, in people with diseases of the cerebellar hemispheres, initial portions of reaching movements are smooth, but an intention tremor becomes prominent towards the end of movements. The time-varying and non-stationary properties implied by open-loop and closed-loop phases of reaching suggests that techniques assuming stationarity (e.g. coherence, linear component analysis) may result in misleading interpretations. To address these concerns, a common approach is to incorporate a sliding time window into the original signal models, such that the stationarity assumption holds for the segment of data in the window. However, selecting an appropriate 29 window length is non-trivial, and it can have a significant effect on the analysis results [40]. Based on the aforementioned observations, in this paper we propose combining hidden Markov models(HMMs) and multivariate autoregressive mAR models into a HMM-mAR framework for modeling nonstationary multivariate sEMG time series and determining dy-namic muscle activity patterns during reaching movements. HMM-mAR is advantageous for modeling sEMG data due to a number of reasons. First, mAR models allow full represen-tation of multivariate sEMG signals, so it naturally captures the dependency relationships between different channels. The results can be further explored to reveal muscle associations and synergies. Furthermore, given that the data are nonstationary and we have no prior knowledge on the timing of state transition of the mAR models, the incorporation of H M M provides a probabilistically tractable and robust way of modeling the dynamic changes of state (i.e. different mAR models) [41]. While the HMM-mAR has been widely used in econometrics, target tracking, and statistical signal processing [42], only a limited number of studies applied HMM-mAR to biophysiological data. In [43], Cassidy and Brown applied HMM-mAR to electrophysiological signals for the purpose of spectral estimation. In the present study, we extend prior approaches to investigate the suitability of HMM-mAR for sEMG signals and study dynamic interactions among muscles during reaching. The main contributions of this paper are as follows: • Presents a HMM-mAR framework for modeling the muscle activities during reaching movements; • Investigates the results of model-fitting when the proposed analysis is applied to raw sEMG, amplitude sEMG and carrier sEMG; • Shows that the structural features derived from estimated model parameters can en-hance the classification performance; and • Suggests that both the carrier data and the residuals after model fitting are informative and provide additional insights into the underlying muscle activity patterns. 30 The paper is organized as follows. In Section 3.2, we describe the proposed HMM-mAR framework and discuss different forms of sEMG signal. A classification scheme is proposed for classifying sEMG signals recorded from healthy and stroke subjects. The classification results are discussed in Section 3.3. Finally,.we conclude our paper. In this section, we first present the HMM-mAR framework and introduce the expecta-tion maximization (EM) algorithm for estimating model parameters. We then present the amplitude-modulated signal model for sEMG signals. Finally, we describe the process of constructing muscle networks and propose a classification scheme based on the learned HMM-mAR model. We propose modeling the multi-channel sEMG signals using a hidden Markov-model mul-tivariate autoregressive (HMM-mAR) process. HMM-mAR model can be considered as a variant of a regular mAR model. The key difference is that in HMM-mAR, the model parameters, including mAR coefficients and noise covariance, are no longer time-invariant, but are modulated by an unobserved Markov chain. In other words, the HMM-mAR model switches between different "sub-models", each of which has its own set of parameters, and thus the resulting time series is piecewise stationary. A good introduction to a first-order HMM-mAR model can be found in [44], where it is referred to as a switching AR model. In our study, we consider a more general HMM-mAR model of order P. For simplicity, we assume that the hidden Markov chain is discrete and is of first order. The mAR model is with order P. Specifically, the sEMG data are formulated as 3.2 M e t h o d s 3.2.1 HMM-mAR model p 1 < k < T (3.1) P = I 31 where T denotes the total number of time points, yk e R M denotes the observed M-channel sEMG signal at time k, and in this study, M = 7. Sk is the hidden Markov chain at time k, and it can be in one of the distinct states {1, 2 , N s } . The state transition of the Markov chain is governed by the following two parameters: (i) transition probability matrix A 6 R N s x N s , with elements Aij — Pr(sk+i = j\sk = i), and (ii) initial state probability 7r € RNs, where 7Tj = Pr(si = i). ap(sk = i) & ] R M x M i s the A R coefficient matrix at time lag p given that the sk is in state i. The residual term, Wk{sk = i) ~ N(0, Ej) € E M is assumed to be an M-dimensional Gaussian with zero mean, and its covariance matrix, Ej, is determined by the hidden state, sk- In summary, the complete parameter set of the HMM-mAR process can be specified as follows A = (A,ir,ap(i),Hi), l<i<Ns. (3.2) 3.2.2 Parameter Estimation via E M Algorithm Given that we only observe Yi-T = (2/1 ,2/2, — ,VT) and the state transition of s is unknown, we use the expectation maximization (EM) algorithm to obtain the maximum likelihood (ML) estimate of the model parameters. A good introduction to the E M algorithm can be found in [45]. The key idea is that the unobserved data are first estimated by the initial parameters in the expectation step (E-step) and these data are used,- as if they were observed, to update the estimate of parameters in the maximization step (M-step). The updated parameter estimate are fed back to the E-step and the process repeats. The iterative E M algorithm keeps alternating between E-step and M-step until convergence is reached in a M L sense. Dey et al. [42] has shown the derivation of the E M algorithm for univariate H M M - A R models. We extended their work and derived the E M algorithm for the HMM-mAR model. The derivation for the HMM-mAR model is similar to the classic H M M model and is shown in Appendix A. Here, we only list the major results as follows: E-step: This step evaluates the auxiliary likelihood function, Q(-). To facilitate the 32 description of the auxiliary likelihood function, we first introduce two new variables, ak(i) and /3k(i). The forward variable ak(i) is defined as ak(i) = Pr(Yi:k, sk — i|A), and the backward variable Bk(i) is Bk(i) = Pr(Yk+i:T\sk = i, A). They can be recursively computed via the weli-known forward-backward procedure as follows: Ns ak+1(j) = 'Y^ak(i)Aijbj{yk+l), k = 1, 2, ...,T - 1 and 1 < i, j < Ns, (3.3) i=i Ns Bk{i) = J2AMyk+i)Pk+i(j), k = 1,2,...,T-I and 1 < i, j <Ns, (3.4) where bj(yk) denotes the observation symbol probability distribution in state j, and is defined as bj(yk) = Pr(yk\sk = j, Yk_P.k_i, A). Given that the noise follows a Gaussian distribution, bj(yk) can be written as N(yk + 2~2p=iap{J)tfk-p\Q^j)- Next, we compute the Q(-) function by taking the expectation of the complete data log likelihood as follows Q(A) = £ { l o g P r ( Y 1 ; T , 5 l : T | A ) | F 1 : T , A o W } Ns T Ns Ns i=i *:=i i=i j=i ^ T Ns P P where r\ A D / ., > N ak(i)Bk(i) E J = S i « f c ( j ' ) A 0 ' ) 6 ^ , j ) = P r ( s f c = i,sk+1 = j\y1:T, A) _ afc(i)Aj&j(Vfc+i)/?fc+i(i) E i J i E j = i afc(*)Aj&j(»fc+i)/?fc+i(i) M-step: This step re-estimates the model parameter A by computing argA maxQ(A) and 33 it yields E L i &(*>.?) ELi7fc(*) ^ 1 f r S i = — Y — — y^7fc(oi»fc+s^(oyfc-p]ii/fc+y)op(*)yfc-p]r» Efc=i 7*W fc=i p=l p=l B(i) — ( fll(i) a 2 (z) o P ( i ) ) fc=i fe=i where 2 /4 -2 By maximizing the Q(-) function at each iteration, the likelihood function will increase and eventually converge to a local maximum. Once the maximum likelihood estimates of the model parameters have been obtained, the optimal state sequence {si, s%,ST} can be found via the Viterbi algorithm [46]. 3.2.3 Amplitude-Modulated sEMG Representation Clancy et al. [32] have suggested that a sEMG signal y(k) recorded during voluntary dy-namic contractions usually can be considered as a zero-mean, band-limited and wide-sense stationary stochastic process, x(k), modulated by the E M G amplitude, m(k), such that y(k) = x(k)m(k), (3.5) 34 where k indicates time. To distinguish x(k) from y(k), we refer to x(k) as the "carrier data". In the literature, there have been several techniques proposed for accurate amplitude es-timation [47]. One that is used in this paper is root-mean-square (RMS), with a moving window applied. The optimal window length may be estimated by the sampling frequency and tasks performed [47]. For this study, we used RMS with a window length of 50 mil-liseconds, or equivalently, 30 data points at the sampling frequency of 600Hz. Based on the estimated amplitude m(k), the carrier stochastic process x(k) can be calculated by dividing the raw sEMG signal y(k) by the estimated amplitude m(k). Fig. 3.1 shows a segment of raw sEMG data and the estimated amplitude and carrier components. The resulting carrier data, x(k), can be approximately represented as a zero-mean Gaussian process [32] and is wide-sense stationary. To our knowledge, almost all existing sEMG analysis techniques focus on amplitude data, and the carrier is simply ignored, but our very recent work suggests that the carrier can also be informative and may provide novel insights into the underlying reaching movements [48]. Therefore, in this study, we propose applying HMM-mAR framework to the following three forms of sEMG data: the raw data y(k), the amplitude data m(k), and the carrier data x(k). 3.2.4 Analysis of Residuals and Carrier Signals Based on our preliminary results, we observed that the raw sEMG and the carrier signals were not well predicted from the mAR model. This observation motivates us to further explore the statistical characteristics of the residuals. Let yk denote the observed sEMG signals and yk denote the signal estimated by the HMM-mAR model. yk can be computed as yk = — Ep=i a>p(h)yk-P, where ap denotes the A R coefficient matrices estimated by the E M algorithm, and sk denotes the state sequence of the hidden Markov chain found by the 35 Viterbi algorithm. The residual is defined as follows rk = Vk-Vk p = Vk + y^QP(sfc)yfc-p-P = i To minimize the effect of the amplitude on the residual,we divide the residual r*, by the amplitude m(k), and call the resulting signal "normalized residual" which is denoted by T V The next step is to fit a model to fk. Given that the normalized residual follows a Gaussian distribution with covariance between channels varying with time, we propose fitting it with an HMM-Gaussian model which is a special case of HMM-mAR with order P = 0. Specifically, HMM-Gaussian model is formulated as fk = Wk(sk), where tu^is an M-dimensional Gaussian with zero mean and covariance modulated by an unobserved hidden Markov chain Sk-As mentioned in Section 3.2.3, carrier signals are wide-sense stationary, enabling a sta-tionary mAR model to be used. To summarize the processing steps performed on the various forms of sEMG, a flowchart is shown in Fig. 3.2. 3.2.5 Construction of Muscle Networks and 3-node Subnets mAR coefficients and residual covariance matrix of HMM-mAR and HMM-Gaussian mod-els essentially represent temporal and spatial correlation between muscle signals. From a biological perspective, they also reveal the dynamic association patterns between muscles during natural movements. To explore the muscle association patterns, we employ ideas from the graph theory and use directed and undirected graphs to represent the connectivity between muscles observed in mAR coefficients and covariance matrices, respectively. In the graph, each muscle is represented by a node variable. A directed graph can be constructed from mAR coefficient matrices by selecting the top n largest off-diagonal elements of the matrix and drawing directed edges between the corresponding muscle pairs. For instance, 36 if element of the matrix is selected, a directed edge is drawn from node j to node i to indicate that muscle j is exerting influence on muscle i. Examples of directed graphs constructed from A R coefficient matrices with n = 20 are shown in Fig. 3.3. An undirected graph can be constructed from a standardized covariance matrix in a similar fashion by forming edges between muscle pairs corresponding to top n coefficients. Since a covariance matrix is symmetric, connections between selected muscle pairs will always be bidirectional. To simplify the graph, the edges between selected muscle nodes are made undirected (see Fig. 3.4). A natural movement such as a reaching movement is typically considered to consist of a sequence of submovements. Each individual submovement can involve the co-activation of multiple spatially-distributed muscles. This observation motivates us to look at sub-networks which comprise a smaller set of muscles nodes in the original muscle connectivity graph. Given a graph with M nodes, CkM distinct subnetworks can be extracted from it, where C denotes the binomial coefficient and k is the number of nodes in subnetworks. In this study, similar to [49], we will focus on sub-networks consisting of 3 muscle nodes, , which we call triple subnets. These subgraphs can be succinctly characterized by the con-nection patterns between muscle nodes. Specifically, a directed triple subgraph derived from a mAR coefficient matrix can have 16 different connection patterns (see Fig. 3.5(a)), whereas an undirected subgraph derived from a covariance matrix has 4 different patterns (see Fig. 3.5(b)). Based on the connection pattern between muscle nodes, each triple sub-graph is assigned a label, and these labels can be further used as classification features for differentiating muscle networks from two different groups. Details are explained in the next section. 3.2.6 Feature Extraction and Classification Our goal in this study is to classify between stroke and healthy subjects based on the es-timated parameters of the sEMG models. One possible approach is to extract structural 37 features from the muscle connectivity graphs derived from mAR coefficient and/or covari-ance matrices. Two types of structural features are considered: edges and triples. The edge feature is essentially a string of l's and O's, each of which denotes the presence of a connectivity between a certain muscle pair (1 for connected pair and 0 for disconnected pair). To obtain edge features, we first convert the connectivity graph of m nodes into an adjacency matrix, T, where element tij = 1 if there is an edge going from j to node i and Uj = 0 otherwise. Since the connectivity graphs are constructed to be loop-free, the diagonal elements of T are all zeros and are neglected for the classification purpose. The off-diagonal elements of the adjacency matrix are then concatenated to form a binary vector of M(M — 1) elements and this is referred to as edge feature. The other structural feature is based on the triple subgraphs extracted from the muscle connectivity graphs (see Section 3.2.5). For an M-node connectivity graph, C\ triple subgraphs can be extracted, and each subgraph is assigned a label according to its connection pattern, as described in the previous section. For classification purpose, the labels of all Cj^ subgraphs are put together to form a categorical feature vector or triple feature vector. To classify the structural features extracted from the muscle connectivity graphs, we choose classification trees as classifiers in this study. Classification trees are chosen over other classifiers, such as support vector machine (SVM), because they are more suitable for categorical data. Furthermore, tree-based classifiers are known for their ability in selecting features which produce largest separation of objects from different classes. Therefore, the classification results can help us identify the muscle recruitment patterns which are most distinctive between stroke and healthy subjects. On the other hand, since SVM implicitly map feature vectors into a high dimensional space, the results are harder to interpret. In addition to structural features, there are alternative ways of extracting features from the learned mAR model. In Hu and Nenov's work [34], they proposed constructing feature vectors directly from A R coefficient matrices by concatenating all the entries in the matrix into a vector of PM2 elements, where P denotes the order of the A R model and M denotes 38 the number of muscles being looked at. A maximal likelihood classifier was employed to classify these feature vectors. For comparison purposes, we adopted a similar approach. First, we concatenate all entries of an A R coefficient matrix to form a M2-element feature vector, assuming that the A R model is of order 1. Because of the high dimension of the feature space, we use principal component analysis (PCA) to perform dimension reduction. The new feature vectors composed by principal component coefficients are classified using SVM [50]. Later, we will compare this classification scheme with the one based on structural features. To evaluate the performance of cross-subject classification, leave-one-subject-out cross-validation scheme is used. During the experiment, each subject performed ./V reaching trials. In each round of classification, all N trials of a subject are left out to serve as test samples while the remaining trials of other subjects are used to train the classifier. Each test trial is classified to either stroke or healthy. Once all test trials are classified, the class label of the test subject is determined based on the majority rule: the subject is labeled to the class that the majority of test trials belong to. This procedure is repeated until all trials of each subject are used as test samples once. 3.3 Results In this section, we first outline the experimental procedures for collecting sEMG data and discuss the data preprocessing techniques used to denoise the raw data. Then we evaluate the performance of the proposed HMM-mAR modeling for the sEMG data and discuss the classification results. 3.3.1 Data Collection and Preprocessing Nine healthy and nine stroke subjects were recruited from the community to perform reach-ing tasks. During the experiment, each subject was first seated in a chair with their hands 39 on the thigh and was instructed to reach and touch a fixed target after hearing an auditory cue. The target was located in the subject's mid-sagittal plane at shoulder height, and its distance with the subject was adjusted such that it is just within the workspace of the paretic arm of the stroke subject or non-dominant arm of the healthy subject. For every subject, the reaching tasks were performed five times on paretic side of stroke subjects and non-dominant side of healthy subjects, resulting in 90 trials in total. The electrical activities of seven muscles (the anterior and lateral deltoid, long head and lateral head triceps, biceps brachium, latissimus dorsi, and the brachioradialis) were recorded using surface electrodes. A bipolar montage was used to minimize the effect of crosstalk. The 7-channel sEMG signals were amplified, sampled at 600 Hz, and high-pass filtered at 20 Hz to reduce movement-related artifact. More details on the sEMG experimental procedures can be found in [51]. After normalizing the sEMG signals to zero mean and unit variance, we denoised the sEMG data using the Empirical Mode Decomposition (EMD) technique. EMD decomposes a univariate time series into a sum of intrinsic mode functions (IMF). The IMFs are extracted from the original signal sequentially via the sifting algorithm [52]. For the purpose of this study, we used the sum of the first three IMF's of each channel in place of the original sEMG data. 3.3.2 Fitting of H M M - m A R First, we study the dynamic changes of state during reaching movements by fitting an H M M -mAR model to the raw multi-channel sEMG signals. The order of the model is determined by Schwarz's Bayesian Criterion(SBC), which is shown by Liitkepohl [53] that among several order selection criteria, SBC yields the smallest mean-squared prediction error for fitted A R models, fn our case, a fourth-order HMM-mAR model is chosen based on SBC (see Fig. 3.6) to fit to the raw sEMG data. An example of the estimated state sequence of the hidden Markov chain from a stroke subject is shown in Fig. 3.7(a). It can be seen that the proposed 40 approach segments each reaching movement into two sections: the initial acceleration phase (State 1), and the terminal deceleration phase (State 2), and the pattern is consistent across all five trials. This result suggests that the proposed technique is robust enough to segment data regardless of the different onset latencies in individual muscles. Also, it is encouraging to note that the states of the Markov chain remain relatively stable in each phase despite the fluctuation in the amplitudes of the sEMG signals. The state change may indicate that the underlying muscle association pattern changes stably from one pattern to another during the reaching movement. When we looked closely at the residuals after HMM-mAR fitting, we noticed that re-gardless of the A R model order P, HMM-mAR provides very poor prediction performance when applied to the raw sEMG signals, as shown in Fig. 3.7(b) where large residuals were observed. Based on this observation, we believe that the residuals may also be informative and should not be naively ignored. Later, we will show that the structural features extracted from covariance matrices of the residuals actually yield good classification performance. Further, we applied the proposed approach to study the amplitude data. One example is shown in Fig. 3.8, where a fourth-order model was chosen based on SBC. It is clear that the residuals are sufficiently small, indicating that sEMG amplitude data are well modeled by the HMM-mAR approach. We also examined the carrier data. Since the carrier data are shown to be wide-sense stationary, instead of using an HMM-mAR model, we fitted a regular mAR model to the carrier data. The results are similar to the ones observed when HMM-mAR was applied to raw sEMG data in the sense that both signals are both poorly predicted by the models and, therefore, large residuals are observed. Our later classification results from using the carrier data show that, similar to the residuals, the carrier can be informative too, despite the fact that it is ignored by traditional analysis approaches. In summary, we observe that the HMM-mAR model can provide very accurate prediction of the amplitude of the sEMG signals, but yield very large residuals when applied to the raw and the carrier sEMG signals. It suggests that HMM-mAR, and thus mAR, could be 4 1 reasonably good for modeling the amplitude sEMG signals, but may not be suitable for modeling the raw and the carrier signals. 3.3.3 Muscle Network As we described in Section 3.2.5, both mAR coefficient and covariance matrices of the fitted HMM-mAR or regular mAR can be used to construct muscle connectivity networks. For illustration purpose, here, we focus on using the first-order model. To visualize the difference in muscle connectivity patterns between healthy and stroke subjects, we construct an average muscle network for each group by averaging individual networks of the subjects within the group. Examples of average group networks are shown in Fig. 3.3 and Fig. 3.4. fn Fig. 3.3, only the top 20 weighted edges are shown, whereas in Fig. 3.4, top 12 edges are shown, ft can be seen that in both figures, the healthy subject have different muscle connection patterns from that of the stroke subject, fn particular, in Fig. 3.3, anterior deltoid has substantially more connections in the healthy subject than in the stroke one. Also, we notice that the connection between anterior deltoid and lateral head triceps is weaker while the connection between lateral deltoid and lateral head triceps is stronger in the stroke subject than in the healthy one. In both Fig. 3.3 and Fig. 3.4, the importance of the deltoid (anterior and lateral) is consistent with the previous study of these data [51] where Eng et al. used traditional univariate analysis technique and observed that the deltoids' activity was significantly altered after stroke. In Fig. 3.4, the connectivity between the deltoids and long head triceps was found more prominently in the stroke subjects (also indicated later by the classification tree results). This observation may suggest a more traditional stroke synergy, where there is breakdown in the normal independent activation of muscles involving the shoulder girdle and those involved in movement of the elbow [54]. 42 3.3.4 Classification Results Our central objective is to classify between healthy group and stroke group, and in partic-ular, we focus on the comparing the nondominant arm of healthy subjects with the paretic arm of stroke subjects. The dominant side of healthy subjects may also be used. However, as suggested in [51], because dominant arm in healthy individuals is usually over-trained for simple tasks such as reaching, the nondominant arm is more comparable with the paretic arm of stroke subjects in terms of their motor capability. For classification purpose, a first-order, two state HMM-mAR model is fitted to the raw sEMG data, and a first order mAR. model is fitted to the carrier data. To avoid over-fitting, we restrict both classification trees and SVMs from using more than three elements of the input features. In our proposed scheme, classification is performed on structural features, i.e. edges an triples, extracted from A R coefficient and residual covariance matrices of the fitted models as described in Section 3.2.6. We let the classification tree do an exhaustive search over all possible subsets of structural features and select the ones which give the lowest error rate. Based on the model parameters from which structural features are derived, six different sets of features are explored. Table 3.1 summarizes the best error rates produced by the classification tree classifiers, where each row represents one set of feature choices. The corresponding clas-sification trees are shown in Fig. 3.9 and 3.10. In these illustrations, Model 1 refers to the estimated model parameters associated with State 1 of the hidden Markov chain, and Model 2 corresponds to State 2. In the rest of the paper, Model 1/2 and State 1/2 are used interchangeably. In general, the classification error rates drop as we increase the number of features used in classification from 1 to 3. Three features are sufficient for providing excellent classification results (e.g. 10% error rate). When considering the raw data, we noticed that the classification performances from using residual covariance matrices are comparable, if not better, to using A R coefficient matrices of either Model 1 or Model 2. This observation suggests that the residuals are useful in differentiating different subject groups (i.e. the 43 control group and the stroke group in our study). It leads us to believe that the residual is also informative and thus should not be simply ignored as noise. When we compare the classification performance from using Model 1 with that of Model 2, we note that they both produce good, comparable performances. Since Model 1 and Model 2 represent different portions of the reaching movement, these 2 representations should provide different but complementary insights into the underlying reaching movements. Such intuition motivates us to check whether combining these structural features from Model 1 and Model 2 could improve the overall performance, fn our study, as we combined features from these two models either using the A R coefficient matrix or the covariance matrix, the performance (see Table 3.2) is noticeably better than any individual model, fn the future research, we will further explore different combinations of input features. When comparing the results from using the carrier and the raw sEMG signals, we note that the classification performance of carrier data is comparable to that of raw data, ft suggests that the carrier could be informative in characterizing reaching movements. To evaluate the performance of the proposed classification scheme where the structural features and classification trees are used, we compare it with other classifiers. The first one is similar to the one proposed by Hu et al. [34] and it was outlined in Section 3.2.6. fn this method, the input features are directly obtained from mAR coefficient matrices by concatenating all elements in the matrices into an input vector. These features are then projected to a lower dimension using P C A before feeding to the S V M classifier. This method gives very poor classification performance, and the average error rate is around 0.5. To have a more fair comparison with the proposed structure-based method, we modified Hu's approach such that the number of elements used for classification are also limited to at most three. The results are reported in Table 3.3. ft can be seen that the proposed method significantly outperforms the modified version of Hu's method. The key difference between the above two approaches, in terms of feature selection, is that Hu's approach directly uses mAR coefficients as features, whereas the proposed approach uses structural features which 44 are based on the connection patterns between muscle nodes. The results shown in Table 3.1 and Table 3.3 suggest that the structural features are more robust and less sensitive to inter-subject variability than the mAR coefficient values themselves. 3.4 Fur ther Discuss ion During our study, we fit an HMM-Gaussian to the detrended residuals from the raw sEMG following the procedure explained in Section 3.2.4. The resulting estimated state sequence from a typical healthy subject is shown in Fig.3.11(a), where frequent oscillation is observed. One thing that is worth mentioning is that when we plot the power spectrum of posterior likelihoods of the state sequence, significant components are observed around 20 Hz band. A similar pattern is observed when we fit the HMM-Gaussian model to the residual from the carrier data. We speculate that this might represent cortical influence over the muscle as cortico-muscular coherence typically finds activity at around this frequency, but further study is required to verify this. So far, as in [51], we only performed the comparison between the non-dominant hand in healthy subjects and paretic hand in stroke subjects. Because stroke condition and hand dominance can both have influence on sEMG patterns during reaching movements, we thoroughly examined the differences by performing different group comparisons. More specifically, sEMG recordings were grouped into four types of experimental groups: healthy dominant hand (HD), healthy non-dominant hand (HN), stroke dominant hand (SD) and stroke non-dominant hand (SN). To focus on the effect of one factor at a time, we fixed the state of one factor each time and thus performed four group comparisons: HD v.s. HN, SD v.s. SN, HN v.s. SN and HD v.s. SD. In particular, we are interested in whether the carrier data is useful in distinguishing different groups, as it was shown to provide the best classification performance in earlier results. Table 3.4 and 3.5 respectively summarize the results of cross-subject validation error rates, based on A R coefficient matrix and residual 45 covariance matrix of the mAR model fitted to the carrier. We can see that classification trees based on up to 3 structural features can classify well the four types of experimental groups. For instance, the error rate is as low as 5% for the comparison of HN vs. SN. 3 . 5 Conc lus ion We have proposed a HMM-mAR framework for studying the dynamic muscle association patterns during reaching movements using sEMG signals. It seems that the proposed frame-work is able to effectively segment multivariate sEMG series and identify the dynamic mus-cle association patterns over time. Further, our sEMG results suggest that the proposed HMM-mAR model fits the amplitude signals well. While applying it to the raw or carrier data, the resulting residuals are informative and should be further investigated. Our sEMG analysis presents a fundamental departure from most existing methods and suggests that the properties of the sEMG carrier data or the residuals after model fitting are essential in characterizing reaching movements. Our future work will focus on inter-subject variabil-ity in muscle association patterns during reaching movements using the carrier data and residuals. 46 Table 3.1: Cross-subject Classification Error Rates Based on Structural Feature Sets Structural Feature Set Edges Triples 1 2 3 1 2 3 Raw: A R Coef., Model 1 0.2222 0.1667 0.1111 0.2778 0.1667 0.1111 Raw: Residual Cov., Model 1 0.3333 0.2222 0.1111 0.2778 0.1111 0.1111 Raw: A R Coef., Model 2 0.2222 0.1667 0.1111 0.2778 0.2778 0.1667 Raw: Residual Cov., Model 2 0.2778 0.1667 0.1111 0.2778 0.1667 0.0556 Carrier: A R Coef. 0.2778 0.2222 0.0556 0.2778 0.1667 0.1111 Carrier: Residual Cov. 0.1667 0.1111 0.1111 0.1111 0.1111 0.0556 Raw and carrier represent the form of sEMG signals to which the HMM-mAR/mAR model is fitted. AR coef. and Residual Cov. denote the model parameters, A R coefficient matrix and residual covariance matrix, respectively. Note that for every A R coefficient matrix, we used its top 20 largest elements (n = 20) to construct the muscle network, whereas for the residual covariance matrix, top 12 elements were used. Table 3.2: Cross-subject Classification Error Rates Using Combined Structural Feature Sets Structural Feature Set Edges Triples 1 2 3 1 2 3 Raw: A R Coef., Model 1 + Model 2 0.2222 0.0556 0.0556 0.2778 0.1667 0.0556 Raw: Res. Cov., Model 1 + Model 2 0.2778 0.1667 0.1111 0.2778 0.1111 0.0556 The error rates of cross-subject classification between healthy and stroke groups. Classifi-cation features were derived from A R coefficient/residual covariance matrices by combining Model 1 and Model 2 extracted from raw sEMG signals. Table 3.3: Cross-subject Classification Error Rates Using A R Coefficients Feature Set Edges Triples 1 2 3 1 2 3 Raw: A R Coef., Model 1 0.5000 0.2222 0.1667 0.5000 0.3333 0.3333 Raw: A R Coef., Model 2 0.3889 0.2222 0.0556 0.3333 0.2222 0.2222 Carrier: A R Coef. 0.1667 0.2222 0.1111 0.2778 0.2222 0.1667 The error rates of cross-subject classification between healthy and stroke groups with SVM as the classifier. Here, the classification features were directly composed from elements of mAR coefficient matrices, in contrast with Table 3.1 and 3.2 where structural features (i.e. edges and triples) were used. 47 Table 3.4: Cross-subject Classification Error Rates between Different Handiness and Motor Impairment using Structural Feature Sets Extracted from A R Coefficient Matrices Groups Edges Triples 1 2 3 1 2 3 HD vs. HN 0.2222 0.2222 0.1111 0.2778 0.1667 0.1111 SD vs. SN 0.2222 0.2222 0.0556 0.2778 0.1667 0.1111 HD vs. SD 0.2778 0.1667 0.1111 0.2778 0.1667 0.1111 HN vs. SN 0.1111 0.0556 0.0000 0.2778 0.1111 0.0556 Classification features are composed of structural features derived from A R coefficient ma-trices of the mAR model fitted to carrier data. For every A R coefficient matrix, we used its top 20 largest elements to construct the muscle network. Table 3.5: Cross-subject Classification Error Rates between Different Handiness and Motor Impairment using Structural Feature Sets Extracted from Covariance Matrices Edges Triples Groups 1 2 3 1 2 3 HD vs. HN 0.2778 0.1667 0.1111 0.2222 0.1667 0.1111 SD vs. SN 0.2222 0.1667 0.1111 0.2222 0.1111 0.1111 HD vs. SD 0.2222 0.1667 0.1111 0.1111 0.1667 0.0556 HN vs. SN 0.2222 0.1667 0.1111 0.1667 0.1111 0.0556 Classification features are composed of structural features derived from residual covariance matrices of the mAR model fitted to carrier data. For every residual covariance matrix, we used its top 12 largest elements to construct the muscle network 48 Raw sEMG Data 0 0.5 1 1.5 Amplitude Data 2 2.5 3 0 0.5 1 1.5 Carrier Data 2 2.5 3 mmmmmmmm 0 0.5 1 '1.5 2 2.5 3 Time (sec) F i g u r e 3.1: Three forms of s E M G signal: raw s E M G d a t a , a m p l i t u d e data , and carrier data . 49 Raw Data, y(k) HMM-mAR RMS 1 Raw, State 1 > 2 R a w , State 1 1 Raw, State 2 > £ Raw, State 2 Residual, r(k) Amplitude Data, m(k) HMM-mAR Detrended Residual, 7(k) 1 Amp, State 1 > E Amp, State 1 ' Amp, State 2 • 2 Amp, State 2 Carrier Data, x(k) 3 X Carrier ' Carrier Figure 3.2: Flowchart summarizing processing steps performed on various form of sEMG data. Anterior Deltoid Lateral Head Triceps Brachioradialis (a) Healthy Anterior Deltoid Lateral Head Triceps Brachioradialis (b) Stroke Figure 3.3: The group average muscle networks extracted from A R coefficient matrices of the State-1 HMM-mAR model fitted to raw data. Only top 20 edges were selected. Frequency of occurrences of an edge is reflected by its width, (a): Average muscle network of the healthy subjects, (b): Average muscle network of the stroke subjects. 51 (a) Healthy-Anterior Deltoid Lateral Head Triceps Brachioradialis (b) Stroke Figure 3.4: The group average muscle networks extracted from residual covariance matrices of the mAR model fitted to carrier data. Only top 12 edges were selected. Frequency of occurrences of an edge is reflected by its width, (a): Average muscle network of the healthy subjects, (b): Average muscle network of the stroke subjects. 52 Figure 3.5: (a): Labels for directed triple subnetworks, (b): Labels for undirected triple subnetworks. 10 12 14 Model Order (P) 16 18 20 Figure 3.6: Model order selection using the Schwarz's Bayesian Criterion versus model order P. The optimal order is the one which gives the lowest SBC, and it is around P = 4 in this case. 53 (a) Fit H M M - m A R to Raw Data 5 « <" -3 a. ^ Q . C . o nj > »J LU O • I J l M i A If 1 j 1 r 8 10 12 Time (sec) Figure 3.7: Fitting a 4-th order HMM-mAR model to raw sEMG data obtained from stroke subject over five consecutive reaching trials. The top panel, (a), shows the estimated state sequence of a 2-state hidden Markov chain. The bottom panel, (b), shows the envelopes of the raw sEMG data (solid line) and the residual of fitted HMM-mAR model (dotted line). Fit H M M - m A R to Ampli tude Data Figure 3.8: Fitting HMM-mAR model to the amplitude component of sEMG data. The top panel, (a), shows the estimated state sequence which is much more stable than the raw data case as shown in Fig. 3.7. Moreover, the estimation performance has greatly improved as indicated by the significant reduction in the magnitude of the residual (dotted line in the bottom panel). 54 (a) Carrier: AR Coef. (b) Carrier: Residual Cov. Figure 3.9: The best classification trees with edges as classification features, (a): Classification tree based on AR coefficient matrix of mAR model fitted to carrier data. Top 20 elements of the matrix were selected to construct the muscle connectivity network, (b): Classification tree based on residual covariance matrix of mAR model fitted to carrier data. Top 12 elements of the matrix were selected. Cn Cn (2, 3, 5, 7, 10, 13, 14, 15; (1, 2,4, 5, 6, 7, 8, 10, 11, 14, 15, 16) Biceps ( Long Head Triceps \ A n t . Deltoid (10,11,13) (1,2,3,4,5,6,7,8, 9, 12, 14, 15, 16) SR: 6 HR: 2 (a) Raw: A R coef., Model 1 (b) Carrier: Residual Cov. Figure 3.10: The best classification trees with triples as classification features, (a): Classification tree based on A R coefficient matrix of State-1 HMM-mAR model fitted to raw data. Top 20 elements of the matrix were selected, (b): Classification tree based on residual covariance matrix of mAR model fitted to carrier data. Top f2 elements of the matrix were selected. 0 5 (a) (b) Fit HMM-Gaussian to Detrended Residual 6 8 10 12 Time (sec) PSD of Posterior Likelihood of Estimated Data 90 120 150 180 210 240 270 300 Frequency (Hz) Figure 3.11: Fitting HMM-Gaussian to the detrended residual. The top panel, (a), shows the estimated state sequence. The bottom panel, (b), shows the power spectrum of posterior likelihood of the estimated data. Notice the significant increase in power around 20 Hz (indicated by the arrow). 57 3.6 B ib l iog raphy E . A . Clancy, S. Bouchard, and D . Rancourt, "Estimation and application of E M G amplitude during dynamiccontractions," IEEE Eng. Med. Biol. Mag., vol. 20, no. 6, pp. 47-54, 2001. M . J . McKeown , "Cort ical activation related to arm-movement combinations." Muscle Nerve Suppl., vol. 9, pp. 19-25, 2000. X . H u and V . Nenov, "Multivariate A R modeling of electromyography for the classification of upper arm movements," Clin. Neurophysiol., vol. 115, no. 6, pp. 1276-87, 2004. J . L i , Z. J . Wang, and M . J . McKeown, "Bayesian network modeling for discovering depedeht synergies among muscles in reaching movements," IEEE Trans. Biomed. Eng., 2007. A . d 'Avel la , P. Saltiel, and E . B i z z i , "Combinations of muscle synergies in the construction of a natural motor behavior." Nat. Neurosci., vol. 6, no. 3, pp. 300-308, 2003. M . J . M c K e o w n and R. Radtke, "Phasic and tonic coupling between E E G and E M G demon-strated wi th independent component analysis," J. Clin. Neurophysiol., vol. 18, no. 1, pp. 45-57, 2001. A . d 'Avel la and E . B i z z i , "Shared and specific muscle synergies in natural motor behaviors," Proc. Natl. Acad. Sci. U.S.A., vol. 102, no. 8, pp. 3076-81, 2005. P. Davidson et al., "Simulating closed- and open-loop voluntary movement: a nonlinear control-systems approach," IEEE Trans. Biomed. Eng., vol. 49, no. 11, pp. 1242-52, 2002. D . Farina and R . Merlet t i , "Comparison of algorithms for estimation of E M G variables during voluntary isometric contractions." J. Electromyogr. Kinesioi, vol . 10, no. 5, pp. 337-349, 2000. Z. Ghahramani, " A n introduction to hidden Markov models and Bayesian networks," Intern. J. Pattern Recognit. Artif. Intell., vol. 15, no. 1, pp. 9-42, 2001. S. Dey, V . Krishnamurthy, and T. Salmon-Legagneur, "Estimation of Markov-modulated time-series v ia E M algorithm," IEEE Signal Process. Lett, vol. 1, no. 10, pp. 153-155, 1994. M . Cassidy and P. Brown, "Hidden Markov based autoregressive analysis of stationary and non-stationary electrophysiological signals for functional coupling studies," J. Neurosci. Meth-ods., vol. 116, no. 1, pp. 35-53, 2002. K . P. Murphy, "Switching Ka lman filters," Department of Computer Science, University of California, Berkeley, Tech. Rep., 1998. T . K . Moon , "The expectation-maximization algorithm," IEEE Signal Process. Mag., vol. 13, no. 6, pp. 47-60, 1996. L . R . Rabiner, " A tutorial on hidden Markov models and selected applications inspeech recognition," Proc. IEEE, vol. 77, no. 2, pp. 257-286, 1989. E . Clancy, E . Mor in , and R. Merlet t i , "Sampling, noise-reduction and amplitude estimation issues in surface electromyography," J. Electromyogr. Kinesioi, vol. 12, no. 1, pp. 1-16, 2002. 58 [48] J . Chiang, Z. J . Wang, and M . J . McKeown, "Hidden Markov multivariate autoregressive ( H M M - m A R ) m o d e l i n g of dynamic muscle association patterns in reaching movements," in Proc. 29th Can. Medical and Biological Soc. Conf. (CMBES29), 2006. [49] R. M i l o et ai, "Network motifs: simple building blocks of complex networks," Science, vol. 298, no. 5594, pp. 824-827, 2002. [50] C . J . Burges, " A tutorial on support vector machines for pattern recognition," Data Mining and Knowledge Discovery, vol. 2, no. 2, pp. 121-167, 1998. [51] P. H . McCrea , J . J . Eng, and A . J . Hodgson, "Saturated muscle activation contributes to compensatory reaching strategies after stroke," J. Neurophysiol., vol. 94, no. 5, pp. 2999-3008, 2005. [52] N . E . Huang et al., "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis," Proc. Math. Phys. Eng. Sci., vol. 454, no. 1971, pp. 903-995, 1998. [53] H . Li i tkepohl , "Comparison of criteria for estimating the order of a vector autoregressive process," Journal of Time Series Analysis, vol. 6, no. 1, pp. 35-52, 1985. [54] P. L . Gribble and D . J . Ostry, "Independent coactivation of shoulder and elbow muscles," Exp. Brain Res., vol. 123, no. 3, pp. 355-360, 1998. 59 Chapter 4 Supplemental Analysis for Windowed Eigenspectrum Method 4.1 In t roduc t ion The objective of this chapter is to extend and improve upon the windowed eigenspectrum classification scheme presented in Chapter 2. The following modifications/improvements are made to the original analysis technique: (1) spatial whitening of raw data, (2) analysis of amplitude sEMG and carrier sEMG, and (3) the use of leave-one-subject-out cross val-idation scheme combined with majority vote. These modifications are described in detail in Section 4.2 along with their justification. The analysis results based on the modified windowed eigenspectrum method are discussed in Section 4.3. 4.2 M e t h o d s This section discusses three modifications to be incorporated into the windowed eigenspec-trum classification scheme. 4.2.1 Whitening of Raw sEMG Data fn the original work presented in Chapter 2, a minimal amount of pre-processing was done in order to maximally preserve the information in raw sEMG data. However, further in-vestigation has shown that in some subjects, the sEMG signals recorded while the subject 60 was at rest are spatially correlated across different channels. The correlated noise can be caused by various factors, both physiological and non-physiological, and it can potentially contaminate the eigenspectral patterns extracted from the sEMG data. As a countermea-sure, a whitening filter is applied to remove any cross-correlation between channels observed during the non-movement part. Although the application of the whitening filter can effectively remove the spatial cor-relation in the noise, it can possibly also remove the correlation associated with muscle co-activation during reaching movements. To overcome this complication, we assume that the sEMG noise observed during the non-movement part remains stationary throughout the entire recording period and was superimposed on the desired sEMG signals. Thus, the whitening of a sEMG recording is done by first constructing the whitening filter using the non-movement part of the raw data. The filter is then applied to the entire sEMG recordings. Let Y denote the raw sEMG data and N denote the non-movement segment of the raw data. Assume that both Y and N are zero-meaned. The covariance matrix of N, denoted by C, is given by C = E{NNT}. C can be decomposed by singular value decomposition (SVD) into C = VLVT, where V is an orthogonal matrix containing columns of eigenvectors, and L is a diagonal matrix consisting of eigenvalues of C. The whiting filter, denoted by W, is given by W = VL-X'2VT = C-1'2. (4.1) The whitened sEMG data, Y, is obtained by Y = WY = C~1/2Y. (4.2) A l l subsequent analysis is performed on the whitened sEMG data, Y. 61 4.2.2 Analysis of Amplitude and Carrier Components In Chapter 3, the proposed HMM-mAR technique was applied to the amplitude sEMG and carrier sEMG of raw data and promising results were found. Therefore, in addition to the whitened sEMG data Y defined in Eqn. 4.2, we also apply the windowed eigenspectrum technique to amplitude and carrier sEMG. The procedure of extracting amplitude and carrier data is described in Section 3.2.3. 4.2.3 Cross Validation In the classification method described in Chapter 2, the performance of cross-subject clas-sification was evaluated using leave-one-trial-out cross validation (LOTOCV). That is, at each time, the recording from a single trial of a subject serves as the test set, and the remaining data were used to train the classifier. The other possible approach is to use leave-one-subject-out cross validation (LOSOCV). The main difference is that, instead of using a single trial each time, LOSOCV uses the recordings from all five trials of a subject as the test set. The test trials are classified one by one and the class label (i.e. healthy or stroke) of the test subject is determined by the majority vote of the classified test trials. Fig. 4.1 provides an illustration for the LOSOCV with majority voting. LOSOCV is more preferable for our study as it not only measures the classification accuracy in the presence of unseen data, but also provides a good indication of robustness of the proposed classification scheme to inter-subject variability. Therefore, LOSOCV is used to subsequent analysis. 4.3 Resul ts To evaluate the performance of the modified windowed eigenspectrum classification scheme, we apply it to the sEMG data collected in Experiment 1, which is described in Section 2.3.1. Based on the motor impairment assessed by Fugl-Meyer scale, participants in the experiment are divided into four subject groups: Healthy, Severe, Moderate, and Mild. Each subject 62 performs 5 reaching trials and a total of 105 trials of sEMG data are collected. For each trial, eigenspectral vectors are extracted from different forms of whitened sEMG data. Fig. 4.2 shows the average eigenspectral vectors of each subject group based on whitened data. The extracted eigenspectral pattern can be intuitively interpreted as the representation of changes in the strength of the dominant muscle synergy. As expected, Healthy and Mild groups have very similar eigenspectral patterns, suggesting that they employ similar muscle recruitment strategies during reaching movements. Furthermore, it can be seen that the magnitude of the eigenspectral patterns is inversely related to the severity of motor impairment as assessed by clinical scales. This observation suggests that in healthy subjects, the reaching movement was dominated by few strongly activated muscle synergies, whereas in severely impaired subjects, more muscle synergies were recruited to complete the task. This agrees with Eng et aVs observations, fn [55], Eng et al. suggest that in stroke subjects, additional muscles are typically recruited to compensate the weaker muscles in order to complete the reaching task. To quantitatively compare muscle activity of different subject groups, the eigenspectral vectors are classified between healthy and stroke subjects with a third order polynomial SVM. Based on LOSOCV, the classification rates given by whitened and unwhitened sEMG data are summarized in Table 4.1 and 4.2, respectively. In general, whitening increases the classification rates for raw and amplitude sEMG considerably. Furthermore, similar to the results reported in Chapter 2, classification rates increase with the severity of motor impairment assessed by Fugl-Meyer scale. 4.4 Conc lus ion fn Chapter 2, we introduced a novel feature extraction technique which exploits the time-varying covariance patterns between sEMG channels. The proposed scheme was able to effectively classify between stroke and healthy subjects. In this chapter, we applied a whiten-63 ing filter to the raw data in an attempt remove the spatial correlation detected during the non-movement part of the recordings. The classification scheme was applied to the whitened sEMG data, and it classification performance evaluated by LOSOCV was compared with the non-whitened data. The results show that whitening can enhance the classification of eigenspectral patterns. Furthermore, the proposed feature extraction technique is able to quantitatively characterize the subtle differences in muscle recruitment patterns between dif-ferent severities of motor impairment. This simple, biologically-inspired approach provides a new perspective to the understanding of motor control in various disease states. 64 Table 4.1: Classification rates for different severities of impairment using different forms of whitened data Whitened Amplitude Carrier sEMG sEMG sEMG Severe vs. Healthy 93.33% 93.33% 66.67% Moderate vs. Healthy 91.67% 100.00% 91.67% Mild vs. Healthy 71.43% 71.43% 28.57% Table 4.2: Classification rates for different severities of impairment using different forms of unwhitened data Unwhitened Amplitude Carrier sEMG sEMG sEMG Severe vs. Healthy 86.67% 80.00% 80.00% Moderate vs. Healthy 50.00% 75.00% 66.67% Mild vs. Healthy 71.43% 42.86% 57.14% i — Trial 1 Test Trial 1i Majority Vote V Subject 1 is Stroke / Healthy Figure 4.1: Illustration of leave-one-subject-out cross validation scheme 65 3.5 3 — Healthy - - Severe - Moderate Mild CD 1\ , , , , 1 0 20 40 60 80 100 % Movement Figure 4.2: Average eigenspectral vectors extracted from whitened data of different subject groups. The vertical error bars represent the standard deviation. 66 4.5 B ib l i og raphy [55] P. H . McCrea , J . J . Eng, and A . J . Hodgson, "Saturated muscle activation contributes to compensatory reaching strategies after stroke," J. Neurophysiol, vol. 94, no. 5, pp. 2999-3008, 2005. 67 Chapter 5 Conclusion 5 .1 Conclus ions In this research, we presented two multivariate approaches to analyzing multi-channel sEMG signals collected from healthy and stroke subjects during reaching movements. In Chapter 2, we exploited the time-varying covariance patterns between muscles, and proposed a feature extraction technique which utilizes the eigenspectrum of windowed covariances. The pro-posed technique was applied to raw sEMG signals, and the extracted eigenspectral feature vectors were classified based on the severity of motor impairment, which was assessed by clinical scales. The proposed scheme was able to effectively differentiate between healthy and stroke subjects, in spite of the minimal differences in the RMS profiles of individ-ual muscles between two groups. Furthermore, we showed that the classification rate is monotonically related to the severity of motor impairment. This result suggested that the proposed technique was able to quantitatively capture the subtle difference in muscle re-cruitment patterns between two populations and, thus, can be extended to a quantifiable metric of motor performance. In Chapter 3, we extended prior approaches of modeling multivariate sEMG signals with mAR models and proposed a dynamical mAR framework for characterizing the non-stationarity of sEMG signals. The proposed algorithm was applied to different forms of sEMG data including raw sEMG, amplitude sEMG, and carrier sEMG. Overall, the pro-posed model was able to consistently segment reaching movements into initial acceleration phase and the terminal deceleration phase. Furthermore, we showed that the structural 68 features extracted from the estimated model parameters yielded significantly better cross-subject classification performance than the conventional approaches where m A R coefficients are directly used as classification features. More interestingly, the carrier data produced clas-sification results comparable to the raw data, suggesting that carrier s E M G is informative and should not be discarded. 5.2 Future Work In the course of the investigations presented in this thesis, a few interesting questions have been raised which merit further research, fn Chapter 3, we proposed the H M M - m A R framework for modeling s E M G signals and extracting muscle connectivity networks. In the analysis, several parameters are chosen based on experimentation, such as the the number of states in the hidden Markov model and the thresholds used for determining the connec-tivity between two muscles. While the parameters suggested in Chapter 3 are suitable for the current study, more rigorous parameter selection methods must be derived before the proposed analysis technique can be widely used in other applications. 69 Appendix A Derivation of E M for H M M - m A R Model Considering the following signal model for HMM-mAR framework: p Vk + 'Y^ap(sk)yk^p = wk(sk), (A.l) where yk G E M denotes the observed M-channel signal and P is the order of the mAR model. sk is the hidden Markov chain at time k, and it can be in one of the distinct states {1 ,2 ,Ns} . The state transition of the Markov chain is governed by the following two parameters: (i) transition probability matrix A G ~RNsxNs, with elements Aij = Pr(sk+i = j\sk = i), and (ii) initial state probability 7T G H N S , where TTJ = Pr(si = i). The term Op(sfc = i ) G E M x M is the A R coefficient matrix at time lag p given that the sk is in state i. The noise term is denoted by wk(sk = i) ~ JV(0,£j). fn the state space notation, Eqn. A . l can be expressed as follows Mz~\ sk)yk = wk{sk), (A.2) where p A{z-\sk) = l + ^2ap{sk)z-\ (A.3) In summary, the complete parameter set of the HMM-mAR model can be specified as follows 6 = (i4,fl-,.a p(i) )E i) ) l.<i<Ns. (AA) To find the maximum likelihood estimate of the model parameters, the Expectation 70 Maximization (EM) algorithm is used. The complete data log likelihood can be found as follows T - l fc=i logPr(Yi : T )Si::r) = log Pr(sj) + ^log Pr(sk+1\sk) + E l o g P r(y f c | s f c , F f c - P : f c - i ) fe=i Ns = ^2S(s! = i)\ogPr(si = i) i= l T - l JVs Ns + EEE 5( S f c+ 1 =3,sk = i) log Pr(sk+1 = j\sk = i) k=l i=l j=l T Ns P + EE 5 ( S f c = *) l o g P r^(yfc + Eap(Sfc = i)yfc-P|0,Ei) fc=l i= l p=l JVs T - l iVs Ns = E 5(si = i) log 7Ti + E E E <5(sfc+i = j , sk = i) log Ay 1=1 fc=l t=l j = l T Ns P + E E 5 ( S f c = l o g ^ (W* + E A P ^ S F C = *)2/fc- P |0, E i ) fc=l i= l p=l Ns T - l Ns Ns = E 5(s i = *) l o§ n i + E E E s ( S k + i = s* = *) los i= l fc=l i= l j = l T Ns +E E 5 ^ = l o § i ^ i - « ^ z _ 1 > s*)y*]} *:=i i=i T - l Ns Ns - E 5 ( s i = o i ° g 7 r t + EEE(5(Sfc+i=j.** = o i o g ^ i j i= l fc=l i= l j=l ^ Ns T ~2/Zz~25(Sk = ^ ) log l^ij i= l fe=l j ATs T " o E E ^ » = i)[A (2 _ 1 ,s f c = i)y f c] tEr 1[A( z- 1 ) S f c = i)v*]} i=i fe=i Define 7fc(z) = Pr(s f c = i) &(*, j) - Pr{sk = i, sk+1 = j) (A.5) (A.6) 71 The auxiliary likelihood function, Q(-), is found by taking the expectation of the com-plete data log likelihood as follows Q{9) = E[\ogPr(Y1:T,s1:T)] Ns T-l Ns Ns i=l k=l i=l j=l j Ns T i=l fc=l ^ Ns T ~ 2 ^ E P r ^ = i)[A(z-\sk = i)yk\tT1rl[A{z-\sk = i)yk]} i=l k=l Ns T-l Ns Ns = 5^1og7Ti7i(i) + ^2l°S^k(i,J) i=l k=l i=l j=l ^ Ns T i=l fc=l ^ Ns T i=l k=l The model parameters can be computed by setting dQ/d9 — 0 and solving for the param-Ns eters 9. In particular, for the Markov chain parameters Atj and 7^, the constraints 2~2 ni = 1 i=i Ns and 2~Z A-ij = 1 must be satisfied to preserve the statistical properties of the parameters. 3 = 1 This can be accomplished by augmenting the Q function with Lagrange multipliers, A. 72 Therefore, 7Tj and A^ can be solved as follows: Ns 1 3 = 1 ^UA = 0 Ns Ns xzZ^ = - X ^ i i ? ) 3=1 3 = 1 Ns Ns A = ~ 7i(j) ( since ^ ^ = 1 ) 3=1 3=1 1X. = = 7i ( 0 A E ^ i 7 i ( j ) a i = l ,'=1 1=1 j = l E L I ^ , J ) + A I = 0 Aij Ns T Ns 3=1 fc=l 3=1 T k=l Ns Ns (since ^ A^ = 1 and ^^k{hj) = 7fc(0 ) j=i j=i - E L i ^ ( ^ i ) A, Efc=i&(»,j) E L i 7fc(») 73 U s i n g the identit ies = ( X " 1 ) ' a n d log \X\ = - log I * - 1 ! , we get dT l o g | E , | = - - ^ - j t o g l E r 1 ! d E r 1 0 1 " d E j = - E A l s o , using the i d e n t i t y ^ - v ^ = aa ' , we get d=i\A(z \sk = i)yk]%1[A(z 1, sk = i)yk] = [A(z 1, sk = i)yk][A(z \sk = i)yk}' Therefore dQ{6) 1 T 1 r d E r 1 2 , K V " / " T 2 1 k=l k=l E, ; .= 1 1 *;=i fc=i E L i 7 f c ( 0 [ ^ - 1 ' . sfc = « ) » f c ] [ ^ ( 2 _ 1 , « f c = i)Vk]' ELi7fc(«) Since the equat ion for the A R coefficient ap(i) also depends o n other A R coefficients, we need to estimates a l l A R coefficients jo int ly . W e can do th is by concatenat ing a l l A R coefficients ap(i) to create B(i) a n d concatenating yt-i,yt-i, --Vi-p to create zt-\. B(i) ± ( 0 l ( i ) a2(i) . . . o P ( i ) ) zt-i = y * - 2 \Vt-pJ G i v e n the i d e n t i t y d[(Xa + b)'C(Xa + b)} dX = (C + C')(Xa + b)a', 74 we can evaluate the derivat ive of Q function w i t h respect to B(i) as follows T dQ(9) _ d . { _ l ^ 7 t W [ ^ + j B ( , ) ; 2 t _ i r E - i [ y 4 + j B ( ^ _ i ] } dB(i) dB(i) 1 2 X = ~\£ ^ 7wr> + s(0*t-i]'Er1^ + £(0*-i] = -IzZ vox?:1 + sr1)^ * + * ( ^ - i K - i T t=i B y equat ing the above expression to 0, we get 0 = dQ{6) dB(i) T T t=i t=i B(i) = - E T ^ W - I I E ^ ^ M ^ , ] - 1 t=l t=l
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multivariate analysis of surface electromyography signals
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multivariate analysis of surface electromyography signals Chiang, Joyce Hsien-yin 2007
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Multivariate analysis of surface electromyography signals |
Creator |
Chiang, Joyce Hsien-yin |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | As the primary method of measuring muscle activation, the surface electromyography (sEMG) is of great importance in the study of motor deficits seen in patients with brain injuries and neuromuscular disorders. While clinicians have long intuitively understood that deficits in motor control are related to inappropriate recruitment of muscle synergies across several muscles, sEMG recordings are still typically examined in a univariate fashion. However, most traditional univariate techniques are unable to quantitatively capture the complex interactions between muscles during natural movements. To address this issue, multivariate signal processing techniques are employed in this thesis to study muscle co-activation patterns in patient populations. A method for classification of multivariate sEMG recordings between stroke and healthy subjects is proposed. The proposed classification scheme utilizes the eigenspectra of time-varying covariance patterns between sEMG channels as feature vectors and the support vector machines (SVM) as classifiers. Despite the minimal differences in the RMS profiles of individual muscles, the proposed scheme is able to effectively differentiate between healthy and stroke subjects. Moreover, the classification rate is shown to be monotonically related to the severity of motor impairment. This simple, biologically-inspired approach is able to quantitatively capture the subtle differences in muscle recruitment patterns between two populations and appears to be a promising means to measure motor performance. The other approach to modeling multivariate sEMG utilizes the HMM-mAR framework, which combines hidden Markov models (HMMs] and multivariate autoregressive (mAR) models. Different forms of sEMG data are analyzed, including raw sEMG, amplitude sEMG and carrier sEMG. The classification between healthy and stroke subjects is performed using structural features derived from estimated model parameters. Both the raw and carrier data produce excellent classification performance. The proposed method represents a fundamental departure from most existing classification methods where only amplitude sEMG is analyzed or mAR coefficients are directly used as feature vectors. In contrast, our analysis shows that the structural features of the carrier sEMG can enhance the classification performance and provide additional insights into motor control. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065582 |
URI | http://hdl.handle.net/2429/31587 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-0358.pdf [ 4.03MB ]
- Metadata
- JSON: 831-1.0065582.json
- JSON-LD: 831-1.0065582-ld.json
- RDF/XML (Pretty): 831-1.0065582-rdf.xml
- RDF/JSON: 831-1.0065582-rdf.json
- Turtle: 831-1.0065582-turtle.txt
- N-Triples: 831-1.0065582-rdf-ntriples.txt
- Original Record: 831-1.0065582-source.json
- Full Text
- 831-1.0065582-fulltext.txt
- Citation
- 831-1.0065582.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065582/manifest