View source for Projects:2016s2-240 Electromyographic Signal Processing for Controlling an Exoskeleton

contents
== Introduction ==
=== Background ===
The exoskeleton is a wearable device which can reduce physiotherapist intervention but ensure the efficacy of rehabilitation in the medical and health area, which can also help soldiers fight better because it can help them get better protection and carry more weapons and equipment in the military area. Recent years, researchers from around the world are working hard so that the exoskeleton can be widely used in practice. 

EMG signal is a bio-medical signal which generated from muscle activities. EMG signal originates from motor neurons in the spinal cord, which is the part of the central nervous system. Motor neurons are in the cell body; their axons extend to the muscle fiber and couples fiber at the end plate area. Muscle movement is controlled by awareness, when the brain is excited and downward conduction, the somatosensory and dendrites of the motor neurons in the central nervous system generate electrical impulses under the stimulation from synapses. These electrical impulses are along the neurons of the axon conduction to the peripheral nerve and muscle contact. When the motor nerve contacts the muscle, the axon branches to many fibers, and each branch ends in the muscle fibers to form synapses, which are called motor end plate. The generated electrical impulses are broadcasting along muscle fibers, which cause a serious change inside muscle fibers, and then produce a contraction of muscle fibers. Therefore, a large number of muscle fibers produce muscle contraction. It can be seen that the spread of action potential of muscle fibers leads to the muscle contraction, while the spread of electrical signals in human soft tissue results in the current field, and the potential difference is represented between the detection electrodes, which is EMG signal [1]. Usually, the amplitude of EMG signal is between 50μV and 5mV. The bandwidth is of EMG signal is from 0—500Hz.

SEMG signal is a bioelectric signal released from the neuromuscular activity recorded by the electrode from the surface of the human skeletal muscle. Different muscle movement patterns are caused by different muscle contractions, and their accompanying SEMG signals are different. It can be used to research or detect the muscle bioelectric activity to determine muscle system skills as well as morphological changes and contribute to neuromuscular system research or provide diagnosis [2]. The frequency range of the SEMG signal is from 0 to 2000Hz, but from the perspective of intensity, SEMG signals are mainly distributed between 20 and 500Hz.

=== Aim ===
This project was based on signal processing and focused on developing algorithms for feature extraction and classification to control lower limbs of an exoskeleton. The feature extracted from the recorded EMG signal will be compared to the previous training recordings to identify what kind of movement the muscles are trying to implement.

The project built up a system that has following contents: 1) Surface EMG signals recording. 2) Surface EMG signal denoising. 3) Surface EMG feature extraction. 4) Feature comparing and recognizing with the help of previous data. 

The outcome of the project is a control system that extracts the feature of the input EMG signal and then delivers the muscle state as an output to the motor control of the exoskeleton. The system is accurate and fast in real-time processing.

 [[File:21.png|bottom|alt text]] 
       System breakdown

== Approach ==
=== Classifier ===
The neural network is a type of the computational model which is used in computer science and statistics widely. The neural network consists of many simple elements connected to each other called artificial neurons [3]. The neural network can be trained through examples, which are also known as training data, rather than explicitly programming, and perform will in those fields where traditional computer programs are difficult to express solutions or feature detection [4]. The neural network has been used to solve tasks in different areas, which are difficult to solve by rule-based programmings, such as computer version and image recognition. The neural network is a system that can be learned and summarized, which means that they can be trained from the existing data. Typically, the neural network is the connection layer, which is a media that the signal propagating from the input layer to the output layer. The input layer has many neurons, which have a large amount of non-linear input information. The hidden layer consists of layers that contain links between the input and output layers of neurons. The output layer generates the resulting output by analyzing the information that is converted between neurons. 
In this project, neural network method is used to classify the movement pattern, especially the BP (Back Propagation) neural network, which is widely used in training the neural network. It is usually combined with some optimization methods, such as gradient descent. The BP algorithm repeats two-stage loops, propagation and weight updates. When the vector is placed on the network, it propagates forward to the output layer. The output of the network is then compared to the desired output (reverse propagation requires the desired output for each input value) to find the error value for each neuron. The error value will then propagate backward until each neuron has an associated error value that toughly represents its contribution to the original output. Then, the neurons of the network can be organized since the neural network has been trained. In addition, via different neurons, the characteristics of the input vector can be learned [26]. when an incomplete or noisy vector is placed in the network, the neurons in the hiding layer will respond to the active output if the new input contains any pattern that is similar to the feature contained in the neuron training.

=== Algorithm ===
In this project, three algorithms are used to analyze the EMG signal.

==== Power spectrum one-threshold algorithm ====

==== Process ====
Usually, different actions generate different powers, therefore, the muscle action ratio can be used. Muscle action ratio is measured by the contraction of muscle rigidity, and the degree of this contraction depends on the frequency of nervous fiber transfer impulse. Therefore, power ratio can be used to represent different EMG signal which generated from the variety of actions by the way of intensity estimation [3]. 
The value K can be used to represent this ratio. The steps of calculating this value K are the following:
	By frequency spectrum, find central frequency P(f);
	Calculating the power of a segment which is [P(f)-20, P(f)+20];
	Dividing this calculated value by the whole segment value, the formula is 
       [[File:20.png|bottom|alt text]]
The reason for using the segment of 5Hz to 400Hz is the main frequency of SEMG signal is between 5Hz and 400Hz.
	Comparing this value K with reference signal recorded in the database, identify the threshold value.
As shown in the following two charts, the shadow area is the power of the center frequency±20Hz; furthermore, many signals with smaller amplitudes existed when the muscle is contracted, so there is much more power speared out when the muscle is contracted. Therefore, the value K of the relaxed muscle is bigger that the value K of the contracted muscle. This is verified through the off-line test.

During the off-line analysis, as for the 64-channels SEMG signal motioned in the previous part, extracting one channel signal (assuming channel 1), the center frequency got is 76.29Hz. Figure 9 is the power spectrum of channel 1 signal, and its center frequency has been marked. Then using proportional value K to identify muscle action. Figure 10 is the plot of the calculated value K. Comparing with the corresponding time domain signal, it is suggested that the value K changes significantly when muscle switches from tight situation to relax situation and vice versa.

In real-time analysis, considering the condition of the muscle activity, the value K will be calculated every half second. Because the sampling frequency is 4000Hz, the data segment is 2000. Usually, the first step is to identify the threshold value K, which is operated through visual observation, because this value depends on who the subject is, which channel is used and the electrode placed on which muscle. The chart below is the value K captured in the real-time test, the threshold value is 0.91 this time. Then typing this value into the MATLAB code, the system will do the classification work automatically in the real-time experiment and the result will be displayed on the screen. If the calculated value K is bigger than the threshold value, then muscle is relaxed; if the calculated value K is smaller than the threshold value, then muscle is contracted.

[[file:19.jpg|bottom|alt text]]

                                        The value K with corresponding denoised EMG signal

==== Result ====
There were 3 female subjects participated in the real-time test, and their ages were 24, 25, 26 respectively. Each participant was tested 10 times for 50 seconds each time, which means 100 segments of data each time was processed. The following table illustrates the error rate for each time.
[[File:18.png|bottom|alt text]]

=== Autoregressive (AR) model coefficient estimation algorithm ===

==== Process ====
AR model coefficient estimation method is based on analyzing SEMG signals in the time domain. In detailed, as for the time sequence X_t, where t = 1, 2, 3, …, N, the AR_((n)) model of which can be represented as the linear fit of X_(t-1), X_(t-2), …, X_(t-n). It can be written as X_t=φ_1 X_(t-1) + φ_2 X_(t-2) + … + φ_n X_(t-n) + a_t, where a_t  is the stationary white noise sequence with zero-mean and variance Var=σ_a^2. The coefficients of φ_1, φ_2 , …, φ_n and σ_a^2 are the parameters that need to be estimated. In this report, the least square method is used to estimate parameters, because it is simple, unbiased and high precision [29]. The matrix form of the signal can be expressed as the following equation. The order chosen for this AR model in this project is 15.

[[file:17.png|bottom|alt text]]

==== Result ====
As mentioned above, the SEMG signal is still stored and processed every 0.5s. Because the sampling frequency is 4000Hz, each set of data contains 2000 sampling points. Then, building the AR model for each set of data through MATLAB, the coefficients matrix can be obtained.
The following graph compares one set of the EMG signal in contracted pattern and one set of the EMG signal in the relaxed pattern.

[[file:16.png|bottom|alt text]]


According to the above tables and graph, what can be found is that as for the same pattern, the coefficients matrixes of the AR model are similar; as for different pattern, the coefficients matrixes of the AR model are significantly different. Utilizing different AR model matrixes to classify, the different pattern results can be obtained, and that can be done easily through the BP neural network.
During the real-time test, each participant was tested 10 times for 50 seconds each time, which means 100 segments of data each time was processed. Subjects are the same as mentioned above. The following table illustrates the error rate for each time.

[[file:15.png|bottom|alt text]]

=== Singular value decomposition algorithm ===
==== Process ====
Like stated above, real-time EMG signal is stored every 500ms. In this paragraph, another algorithm will be introduced to extract the feature of the EMG signal.
Singular value is one of the characters of a matrix. The following theorems state the benefits of extracting features by singular value decomposition.
Theorem 1: Stability of singular value
A^(m*n),B^(m*n)  ∈R^(m*n)
The singular values of two matrix A and B are σ_1≥σ_2≥⋯≥σ_n and τ_1≥τ_2≥⋯≥τ_n respectively. Then it can be proved that：
|σ_i-τ_i |≥|A-B|_2    i=1,2,…,n
Therefore, when a matrix A has weak vibration, the change of its singular values is not larger than the Euclid norm. Which means singular values of a matrix is fairly stable. For signals that contain noise or trend, the feature matrix after SVD does not vary a lot. This character reduces the requirement of pre-treatment to sEMG signal and guarantees the accuracy of feature matching.

Theorem 2: Proportion invariance
The singular values of matrix A^(m*n) areσ_1≥σ_2≥⋯≥σ_k. The singular values of matrix 〖α*A〗^(m*n) are〖σ_1〗^*≥〖σ_2〗^*≥⋯≥〖σ_k〗^*. Then 

[[File:11.png|bottom|alt text]] 

Therefore, if the two matrixes are proportional, then their singular values are also proportional, and the proportion remains unchanged. This means normalization which is used to get rid of the influence of amplitude does not change the relative amplitude of singular values.
It can be seen from the theorems above that singular value is an inherent character of a matrix, and it satisfies the requirement of SEMG feature extraction (stability and proportion invariance)[4]. It can be confirmed that singular value is an effective way to represent the feature of a matrix.
The matrix that is used to represent the segment is the sEMG signal after STFT. By calculating the singular values of the complex matrix containing EMG segment after STFT and put them into the neural network, the project team was able to extract the feature of muscle state as the output.
In this project, the signal after STFT was stored in a matrix, MATLAB code then calculated the singular value of the complex matrix. A plot of the singular values of an example EMG signal is shown below: 

[[File:14.jpg|bottom|alt text]]

This figure shows that there is a significant difference between the singular values of contracted segments and relaxed segments. Contracted segments have larger singular values and relaxed ones have smaller singular values. Moreover, the plot of two different states of muscle follows their own trend, which makes it possible to separate two features using a BP neural network.
It is clear that the singular value difference after the 5th order is too small compared to the first few orders. Therefore, all the small singular values of order 5 are omitted.
In this project, the first five orders of the singular values were selected to be the feature singular values. These singular values will be put into the training data matrix. The training matrix has m rows and n columns, where m equals twice the length of training time in second (segment length = half a second) and n is the number of order taken into consideration, which is 5 in this project. All the rows for contracted EMG signal in the training matrix correspond to a [0 1] output in the neural network. The rest rows contain singular values from relaxed EMG segment and lead to a [0 0] output. In this way, the neural network is trained with the EMG signal generated by the same person’s same muscle at the same position, which will minimize the error caused by individual difference. All the training session is completed automatically without any manually set value.

 
The raw EMG signals used in this algorithm were recorded from a 24 years old healthy male’s left forearm. Through a 64-channel electrode and amplified with REFA amplifier. The singular values below are from contracted and relaxed segments of channel 1. The program asks the tester to contract the muscle for 8 seconds to record sEMG for contracted muscle. The segments are then analyzed by STFT and a new complex matrix will be created for the signal after STFT. After that, singular values of the complex matrix are put into the first 16 rows of the training feature matrix. After contraction training finishes, the tested will be asked to relax the muscle for 8 seconds to record sEMG for relaxed muscle. The Same procedure runs again and fills the training feature matrix with the singular values of relaxed segments.
Four groups of singular values of contracted and relaxed segments are shown below respectively.

[[file:13.png|bottom|alt text]]

It can be seen that there is significant difference between the singular values of contracted and relaxed segments. Which means the separation of two states of muscle can be done well with the BP neural network. After training, the program starts real-time processing. The same procedure runs to calculate the singular values of the segments that are stored in the buffer. The singular values are then put into the neural network that was trained by the training feature in the preparation step. The output of the neural network is the state of the muscle in the period of 500ms. The results of muscle state are displayed on the screen.

==== Result ====
During the real-time test, each participant was tested 10 times for 50 seconds each time, which means 100 segments of data each time was processed. Subjects are the same as mentioned above. The following table illustrates the error rate for each time.

[[file:12.png|bottom|alt text]]

== Conclusion ==
=== Discussion ===
In general, this project has been successful in most measures, three algorithms were brought to extract the key features of the SEMG signal and classify different muscle activity. Through the above experimental results, it is not difficult to find that all the three algorithms perform well in analysis and classifying the SEMG signal. By comparing the average error rate of the three methods, obviously, the SVD algorithm has the lowest error rate, and the power spectrum one-threshold algorithm has the highest error rate, but it is still acceptable.
The power spectrum one-threshold algorithm is a simple algorithm, which directly analysis the SEMG signal from the frequency domain and calculating the power to capture the key features. It is convenient to classify obvious simple action, such as contracted pattern and relaxed pattern, while for some complex actions, the one-threshold is not enough. Furthermore, the selection of the one-threshold value is manual and the process is not operated in real time situation, which means it will bring some inconvenience in the practice operation, especially for controlling the exoskeleton. In addition, the existence of high error rate is also due to this situation. 
The AR model coefficient estimation algorithm and the SVD algorithm are similar, they both analyze the SEMG signal from the time domain, and the classification work is operated by the BP neural network, so the error rate of them is close. Unlike the power spectrum one-threshold method, both of them can be totally performed automatically. The reason for the different error rate is the coefficients of the key feature matrix, because some coefficients of the AR model matrix for distinguishing different action patterns are very close, while the distinction of coefficients of different patterns in SVD matrix is obvious.
However, the current three algorithms still have some shortcomings, especially for the detection of the soft muscle activity which plays a vital role in controlling the exoskeleton.   In practice, patients are hoped to operate the exoskeleton with soft muscle contraction. The reasons are the following:

1) The exoskeleton is an intelligent mechanical equipment which can be driven by soft muscle activities;

2) The exoskeleton is a wearable device that covers the limbs that control the exoskeleton. Before the exoskeleton is activated, the patient who wears it does not have enough space to contract the limbs violently and begin to imitate the patient’s movement;

3) Some patients have less muscle strength, so they cannot control the exoskeleton with a very intense contraction.

unfortunately, the results of current three algorithms for detecting mild muscle activity are not very ideal, because the error rate increases drastically when testing slight muscle activity. The error rate of the power spectrum one-threshold algorithm to detect soft muscle activity is 70%-80%; the error rate of the AR model coefficient estimation algorithm is around 50%; the error rate of the SVD algorithm is nearly 30%. As for the power spectrum one-threshold algorithm, the values K of the soft contraction and relaxation are close, and sometimes they are even the same, so manually choosing the threshold value is not a perfect method. The same situation happens in the AR model coefficients, coefficients for distinguishing different action patterns getting much closer which increases the error rate. At last, for the SVD algorithm, the coefficients for slight contraction is smaller than intensely contraction but they are still bigger than that for relaxation. However, this disparity cannot be recognized very well by the backpropogation neural network.

=== Future work ===
Based on the current status of the project, there will be lots of options for the future work. Still in order to control the exoskeleton with the EMG signal eventually, some suggestions are listed as follows:

1) Decreasing the error rate
As mentioned above, the result of the current algorithms for detecting slight muscle activity is unsatisfactory. Therefore, a more detailed study of the key features of the EMG signal is imperative, which means focusing on the difference between slight contraction and relaxation is necessary. Furthermore, another reason for existing errors is that all three algorithms cannot deal with the signal generated very well when the muscle switched from the contraction to the relaxation and vice versa. Therefore, suitable processing the signal of the switching state is important to decrease the error rate as well.

2) Achieving complex exoskeleton activities
The current status of the project is mainly focusing on one degree of freedom of the muscle activities – contraction and relaxation because the used equipment can only detect one muscle activity at one time, while analysis complex muscle activities must obtain the EMG signals from different muscles at the same time. Therefore, in the future, research multi degrees of freedom to achieve complex exoskeleton activities will be better for the project, but the premise is the better equipment which can extract EMG signals from multi muscles at the same time.

3) Developing better neural network algorithm
The BP neural network MATLAB code used in this project is provided by the MathWorks Ltd. That is a simple classifier and the accuracy is not high, which is a reason for increasing the error rate. Therefore, it would be a good choice for developers to write their own neural network algorithm in the future.

== Team members ==
Hanzhi Wang & Junyao Wang

== Supervisiors ==
A. Prof. Mathias Baumert

Dr. Tien-fu Lu

== References ==
[1] Y. H. Chiou, et al, “The comparison of electromyographic patter classifications with active and passive electrode,” Medical Engineering and Physics, Vol.26, pp.605-610, 2004.

[2] D. Rodolfino, et al, “Effects of viewing affective pictures on sEMG activity of masticatory and postural muscles,” Neuroscience Letters, Vol.554, pp.10-14, Jun. 2013.

[3] M. Saduf and M. A. Wani, “Comparative study of back propagation learning algorithms for neural networks,” Dept. of Computer Science, Univercity of Kashmir, Srinagar, J&k, India, International journal of advanced research in computer science and software engineering, Vol3, Issue 12, Decemeber 2013.

[4] G. Deng and M. C. Ferris, “Neuro-dynamic programming for fractionated radiotherapy planning,” Springer Optimization and Its Applications, Vol.12, pp.47-70, 2008.

[5] C. Fleischer and G. Hommel, “A Human–Exoskeleton Interface Utilizing Electromyography,” Comp. Eng. Elec.,Univ.Berlin, Berlin, Rep.TRO.2008.926860, Feb 2008.

[6] G.Golub and W.Kahan, “Calculating the singular values and pseudo-inverse of a matrix,” in Journal of the Society for Industrial and Applied Mathematics: Series B, Numerical Analysis, Vol.2, No.2, 1965, pp.205-224.