EMG Decomposition and Artefact Removal

Traditionally, in clinical electromyography (EMG), neurophysiologists assess the state of the muscle by studying basic units of an EMG signal, which are referred to as motor unit action potentials (MUAPs). Information regarding themorphology and rate of occurrence ofMUAPs is often used for diagnosis of neuromuscular disorders. In addition, recent studies have shown that the analysis of the energy content of MUAPs is a possible way for discriminating among normal, neurogenic, and myopathic MUAPs [38], illustrating, thus, the clinical value of the interpretation of MUAP information.


Introduction
Traditionally, in clinical electromyography (EMG), neurophysiologists assess the state of the muscle by studying basic units of an EMG signal, which are referred to as motor unit action potentials (MUAPs). Information regarding the morphology and rate of occurrence of MUAPs is often used for diagnosis of neuromuscular disorders. In addition, recent studies have shown that the analysis of the energy content of MUAPs is a possible way for discriminating among normal, neurogenic, and myopathic MUAPs [38], illustrating, thus, the clinical value of the interpretation of MUAP information.
A common way of obtaining such information is by observing MUAP activities on an oscilloscope and listening to their audio characteristics over the speakers. When doing this, the researcher is implicitly performing a time and frequency analysis of MUAPs. However, the results of this analysis are dependent on the experience of the investigator and on his ability to extract relevant information from the visual and auditory analysis. Furthermore, this procedure is time-consuming and prone to error.
The drawbacks related to the procedure described above have motivated the use of computer-based techniques for extraction of MUAPs from EMG signals [3,19,22,27,29,32,40]. Such methods, also known as EMG decomposition techniques, aim at classifying MUAPs generated by a common source into the same group. The results of this classification may provide information regarding the orchestration of the neuromuscular system, and therefore of the state of the muscle. A similar problem, often referred to as spike sorting, is found in the study of neuronal activities [25]. In this case, neuronal action potentials from the same source are classified into a common group. detection of MUAPs from superficial muscles [8]. This advancement has received widespread support among researchers and clinicians because of the ease of use, reduced risk of infection, and the greater number of motor unit action potential trains obtained compared to needle sensor techniques [47].
Currently, computer-based EMG has become an indispensable tool for investigations seeking to explain the state of the muscle. Different methodologies, ranging from simple quantitative measures to automatic systems that enable the assessment of neuromuscular disorders, have been developed [30]. Such tools are important for standardization of results and also they may reveal important features in the signals, which might be barely perceived from a manual analysis [35].
A typical system for extraction of MUAPs from EMG signals may require several stages of signal processing, for instance, signal detection and filtering (i.e., artefact removal) [45], feature extraction or selection [36], data clustering or classification [23]. Specific research may be carried out in each of these steps.

Strategies for EMG decomposition
One of the most rudimentary strategies for isolation of MUAPs belonging to a common group is by means of a thresholding scheme. The central idea of this technique is to use a voltage threshold trigger for detection of MUAPs that have similar height, which is represented by the amplitude of the highest peak of a MUAP. In this method the experimenter positions the recording electrode so that MUAPs of interest are maximally separated from the background activity (noise). MUAP activities are then measured with a hardware threshold trigger, which generates a pulse whenever the measured voltage crosses the threshold. Such pulses may be used for triggering data collection and further storage of MUAPs or their occurrence time in a computer. The main advantages of this method are that: (i) It is easy to apply since it requires minimal hardware and software; (ii) It is a good starting point for detection of the strongest MUAPs. The main drawbacks of this technique are that: (i) The threshold level, mainly dependent on the signal-to-noise ratio, determines the trade-off between missed spikes (false negatives) and the number of background events that cross the threshold (false positives); (ii) Only a single feature (the MUAP highest peak) is used for data classification. As a consequence of this strategy two MUAPs with different shapes might be grouped together because they have similar peak amplitude.
The drawbacks related to the technique discussed above highlight the importance of a pre-processing stage prior to grouping MUAPs. The presence of high levels of background activity in the signal suggests that the use of a filter for reduction of the background activity is necessary. Different digital filters may be used for this purpose. Typically high band-pass filters are used for attenuation of very low frequency components in the signal related to noise, which can be either inherent from the hardware used for data acquisition or the contribution of distant MUAPs from the detection point.
Another common approach is the use of differential filters. Such filters are low-pass filters and have been used in many investigations. The main drawback of such filters is that they may generate artificial spikes and modify the shape of MUAPs. This has motivated 262 Computational Intelligence in Electromyography Analysis -A Perspective on Current Applications and Future Challenges the use of alternative filtering procedures, for instance, wavelets and a spatial filter, known as a Laplacian filter. The implementation of these filters as well as the introduction of an alternative method for EMG signal filtering is further discussed in the chapter.
Another relevant pre-processing step is the one that segments the EMG signal into windows containing active and inactive segments. This is a signal detection stage which aims to identify the activities of single MUAPs or their combinations, which is known as MUAP overlaps. Furthermore, this step separates noise in inactive segments from useful information in active segments. The detection of active and inactive segments may be performed visually or manually, i.e. the researcher may classify regions of an electromyographic signal into one of those categories, however such a method is time-consuming and requires concentration, which may introduce inconsistency in the signal analysis.
Different approaches may be employed for automation of this pre-processing stage: (i) The use of the root-mean square of the EMG signal together with a pre-defined threshold; (ii) A threshold proportional to the maximum peak in the signal; (iii) A threshold which is manually adjusted; (iv) Wavelets. The main assumption of this method is that there might be similarity between the mother wavelet and action potentials, and when the correlation between an active window and the mother wavelet is high an active segment is detected.
As one may note there exists a variety of strategies that could be considered for automation of the detection of active and inactive segments. Techniques that make an assumption about a pre-defined height and width of MUAPs may be more susceptible to failures, i.e. active regions that do not fit the predefined window will not be detected. Note also that the level of the background activity in the signal will also influence the determination of the beginning and end of active segments. The detection of active regions is further discussed in this chapter.
After detection of active regions it is possible to group them into logical units (clusters), however, it is very common to obtain features from those regions prior to data clustering. For example, morphological features of MUAPs, i.e. duration, amplitude, area, number of phases (number of baseline crossings) and number of turns (number of positive and negative peaks) have been employed. Other successful approaches include the use of coefficients of the Fourier transform, the coefficients of the Wavelet transform, the use of time samples of the band-pass filtered signal and low-pass differentiated signal, and the use of autoregressive and cepstral coefficients. Some other approaches, known as feature extraction and selection procedures, try to obtain features that maximize cluster separability. Examples of the use of such techniques include the application of Principal Component Analysis (PCA) and Independent Component Analysis (ICA).
The grouping of MUAPs is commonly performed by means of at least three distinct strategies: (i) Template matching: raw MUAPs, referred to as MUAP templates, are first classified or identified, and then used for classification of new MUAPs. The initial MUAP templates may be manually selected from the EMG signal or be chosen automatically from a clustering procedure. During data classification, usually MUAP templates are modified by an update rule, which takes into account the variability of MUAP shapes in the data set; (ii) Clustering: a clustering technique is used for grouping patterns represented by features selected from active regions; (iii) Hybrid: in this approach, first a clustering technique is used for grouping part of the data set (normally the first 3 to 5 seconds), and then the non-classified data set is grouped into one of the classes defined in the first step.
Both the template matching and hybrid techniques require a priori identification of patterns in the data set before classification of the entire data set. The main disadvantage of these methods is that if new MUAP classes appear they will not be identified. The main advantage is that extra information regarding MUAP activities, e.g. the study of the firing time of MUAPs, may be taken into account in the final classification. The clustering approach has the advantage that it makes no assumptions about the data set to be grouped.
The main processing steps discussed so far form the basis of a complete EMG decomposition system. When the final application of such system is to study the firing behaviour of sources that generate MUAPs, it may be necessary to include an extra stage that deals with a problem known as MUAP overlaps. Overlapping spikes occur when two or more spikes fire simultaneously. When using the clustering technique for grouping active segments it may be possible to detect such overlaps as outliers. There are at least three strategies for dealing with overlaps: (i) Once a spike is classified it is subtracted from the active segment, in the hope that this will improve the classification of subsequent spikes. This approach requires a template of the spike. It yields reasonable results when two spikes are separated well enough so that the first can be accurately classified, but fails when the spikes are close together. Another problem with this approach is that the subtraction can introduce more noise into the waveform if the spike model (template) is not accurate. Also subtraction-based approaches may introduce spurious spike-like shapes if the spike occurrence is not accurately estimated; (ii) Another approach is to compare all possible combinations of two or more spike models. However, for some applications the computation time for performing this comparison may be prohibitive; (iii) The use of multiple electrodes or an array of electrodes may reduce the problem of overlapping spikes, because what appears as an overlap on one channel might be an isolated unit on another. Since the main aim of solving the overlapping problem is to increase the accuracy of estimators (e.g. mean) obtained from the firing of motor units, an alternative option might be to work directly with a precise estimate of the estimator considering missing data points (i.e that some MUAPs are missing).
Finally, once the system is designed and implemented it is important to test its accuracy. At least three methods are well accepted for this purpose: (i) Synthetic signals: artificial EMG signals are generated and employed for testing the stages of the system. The main advantage of this approach is that the characteristics of the analyzed signal are totally known; (ii) Manual classification of MUAPs: MUAPs are visually classified by the researcher and the results of this classification are used as reference for evaluation of the automatic classification; (iii) Comparison between MUAP activities from different channels: the consistency of the decomposition data of the same units from two different electrodes provides an indirect measure of the accuracy in real data decomposition.

Artefact removal from EMG signals
The detection of electromyographic signals is a very complex process, which is affected not only by the muscle anatomy and the physiological process responsible for the signal generation but also by external factors, for instance, the inherent noise of the hardware employed in the signal amplification and digitalization. As a result EMG signals are often corrupted by noise.
It may be very difficult, if possible at all, to extract useful information from very poor signal-to-noise ratio EMG signals. In some applications, for example, the decomposition of electromyographic signals, a high level of background activity could impede the accurate segmentation of the signal into regions of activity that may represent the activity of single motor unit action potentials, influencing thus the final results of the EMG decomposition.

Low-pass differential filter
Since its introduction, the low-pass differential filter (LPD) [44] has been widely employed in EMG signal processing [15,19,32,[40][41][42]. This filter is implemented in the time-domain as: where x k is the discrete input time-series and y k is the filtered output. N is the window width to adjust the cut-off frequency. Increasing N will reduce the cut-off frequency of the filter. This may be easily perceived if Equation 1 is studied in the frequency domain.
For this, consider the following difference equation obtained from Equation 1: Its representation in the frequency domain may be obtained via its Z transform as follows, is the filter transfer function H(z), f is the frequency in Hz and f sr is the sampling frequency in Hz.  The main advantages of using the LPD filter are that it is very easy to implement, and it is considerably fast for real-time applications, but some drawbacks regarding this filter are discussed in [46]: • Under conditions of low signal-to-noise ratio the relatively strong high frequency noise (background activity) may be accentuated; • Since the LPD filter is not an ideal low-pass filter, there will exist severe Gibbs phenomenon, that is, the leakage of energy frequency out of the filter pass-band as illustrated in the Bode diagrams shown in Figure 1. As a result, many high frequency noise components will pass through the filter.

Weighted low-pass differential filter
The weighted low-pass differential filter (WLPD) was proposed in [46] as an alternative to the LPD filter. The main difference is that an appropriately weighted window is included in Equation 2 for reduction of the Gibbs effect. The WLPD filter is implemented in the time-domain as: where w(n) is an N point windowing function. Several windows such as Barlett, Hamming and Hanning may be employed. If a rectangular window is used then Equation 3 is an LPD filter. Similarly to the LPD filter the WLPD is very easy to implement and considerably fast for real-time applications, but results presented in [17,46] show that phase distortion may  be present in the filtered signal and depending on the level of noise false spikes may be generated.

Adaptive 50/60 Hz noise cancellation
When the source of noise is well known it is possible to employ methods for adaptively reduce its influence over the desired signal. This is, for instance, the aim of adaptive 50/60 Hz noise cancellation based on the least mean square (LMS) algorithm [2].
The block diagram shown in Figure 2 depicts the main idea of the adaptive noise cancelling algorithm. The signal s[n] is corrupted by the noise y[n] yielding the corrupted signal x[n].
The adaptive algorithm will estimate parameters for a filter capable of attenuating the desired components and maintaining relevant information from the signal. The piece of code below is a Matlab function created for attenuating 50 Hz noise from signals. It was used to illustrate the application of the technique on an EMG signal which was strongly corrupted by 50 Hz noise. The plots given in Figure 3 show the input signal, the estimated noise and the filtered signal. The power spectra of these signals are given in the graph at the bottom. Their analysis allows us to conclude that the application of the adaptive filter attenuated the undesired 50 Hz component. Note that the 50 Hz component is also part of the EMG signal, therefore its complete elimination is not desired.

Signal filtering based on wavelets
One of the main applications of wavelets is noise removal from corrupted signals. A general procedure for signal de-noising involves the following three steps that have been successfully applied [14, 16,  3. Signal reconstruction.
In the first step, both a wavelet prototype and the decomposition level (N) are chosen, and the wavelet decomposition of the signal at level N is performed. In the next step, for each level from 1 to N, a threshold, t N , is selected and soft-thresholding is applied to detail coefficients as shown in Equation 4.
where Y N is the de-noised version of the Nth detail coefficients and the function (x) + is defined as Finally, in the third step the denoised signal is estimated by using the original approximation coefficients of level N and the modified detail coefficients of levels from 1 to N. The application of this procedure is illustrated in the following example.

Example of application
The procedure for noise removal based on wavelets is applied to an experimental surface EMG signal in this section. The main aim is to attenuate the background activity present in the signal and at the same time to preserve its shape.
The Matlab Wavelet Toolbox 1 is used for data processing (see Figure 4). The EMG signal is decomposed into three levels (N = 5) and the db5 wavelet prototype is employed. The choice for the decomposition level may be justified by the fact that this signal is mostly contaminated by high frequency noise that will be highlighted by short time-scale components (e.g. detail coefficients of the first and second level) and therefore further decomposition would not contribute significantly to the final form of the filtered signal. The criteria for selection of the wavelet family (Daubechies) is mainly because of its reported success in removing/attenuating noise from EMG signals, and in its use for feature extraction from MUAPs [18,37,48].
The results of the signal decomposition are depicted in Figure 4. A window of the input EMG signal together with its filtered version are presented in Figure 5.

EMG signal filtering based on Empirical Mode Decomposition
In this section a novel algorithm for signal de-noising based on the Empirical Mode Decomposition (EMD) method is presented [7]. The developed filter may be employed as an alternative to the approach based on wavelets described above.

The Empirical Mode Decomposition method
Huang et al. (1998) [21] described a new technique for analyzing nonlinear and non-stationary data. The key part of the method is the Empirical Mode Decomposition (EMD) method, in which any complicated data set can be adaptively decomposed into a finite, and often small, number of Intrinsic Mode Functions (IMFs). The name Intrinsic Mode Function is adopted because those components represent the oscillation modes embedded in the data. In the case of Fourier analysis, oscillation modes (i.e. components) in a signal are defined in terms of sine and cosine waves. The EMD thus defines oscillation modes in terms of IMFs, which are functions that satisfy two conditions [21]: 1. In the whole time-series, the number of extrema and the number of zero crossings must be either equal or differ at most by one. Note that extrema are either local minima or local maxima. Furthermore, a sample g i in a time-series is a local maximum if g i > g i−1 and g i > g i+1 , and a sample q i is a local minimum if q i < q i−1 and q i < q i+1 , where i is a discrete time. The definition above is empirical and currently there is no explicit equation for estimating IMFs, thus any arbitrary time-series that satisfies conditions 1 and 2 is an IMF. Furthermore, by means of the analysis of the power spectra of IMFs it is possible to verify that these functions represent the original signal decomposed into different time-scales (or frequency bandwidths). This is illustrated in [4,6]. Thus, both wavelets and the Empirical Mode Decomposition provide the decomposition of a signal into different time-scales. The main difference is that one method performs the signal decomposition adaptively and based solely on the available data, whereas the other normally uses a set of pre-fixed filters.

The algorithm for signal filtering based on the EMD
The success of the general procedure for noise removal using wavelets is based on the fact that it is possible to filter signal components individually instead of filtering the original signal. This is desirable because some components may highlight the noise and thus it may be easier to attenuate its presence.
Similarly, the Empirical Mode Decomposition also provides the decomposition of a signal into different time-scales or IMFs. This means that it is also possible to filter signal components individually instead of the original signal. This suggests that the strategy for signal de-noising based on wavelets may also be applied to intrinsic mode functions. Thus the following procedure was proposed for EMG signal filtering [3,7]: 1. Decompose the signal into IMFs; 2. Threshold the estimated IMFs;

Reconstruct the signal.
This procedure is practical, mainly due to the empirical nature of the EMD method, and it may be applied to any signal as the EMD does not make any assumption about the input time-series. A block diagram describing the steps for its application is shown in Figure 6. First, the EMD method is used for decomposing the input signal into intrinsic mode functions, IMF 1 , . . . , IMF N , where N is the number of IMFs. These IMFs are then soft-thresholded, yielding tI MF 1 , . . . , tI MF N , which are thresholded versions of the original components. The filtered signal is obtained as a linear summation of thresholded IMFs.
A very common strategy used in the filtering procedure based on wavelets is to use the soft-thresholding technique described in [14]. The same idea is used for thresholding IMFs. For each IMF from 1 to N a threshold, t n , n = 1, . . . , N, is selected and soft-thresholding is applied to individual IMFs.
The threshold t n is estimated by using the following strategy: a window of noise is selected from the original signal and then the boundaries of this window are used to extract regions of noise from IMFs. The standard deviation of each of those regions is then estimated, and  these are regarded as the required thresholds (t 1 , ..., t N ). Figure 7 illustrates the procedure for estimating t n .

An example of application of the filter based on the EMD method
In this section, an example illustrating the application of the procedure for filtering EMG signals based on the EMD is provided. For this, the experimental surface electromyographic signal shown in Figure 8 is used. Figure 8. EMG signal and its filtered version. The residue, which is the difference between those signals is also shown.
The first step of the algorithm is to use the EMD (or sifting process) to decompose the experimental signal into intrinsic mode functions. Eight IMFs (IMF1, . . . , IMF8) are obtained and they are presented in Figure 9(a).
The subsequent step is to threshold the components IMF1, . . . , IMF8. Equation 4 was employed for denoising individual IMFs. The results of this procedure are shown in Figure  9(b). Note that in order to estimate thresholds for IMFs the boundaries of the window of noise (0.03 s and 0.07 s) indicated in Figure 8 were selected.
In the last step, the resulting (de-noised) components were combined to generate a filtered version of the original signal as shown in Figure 8. In the same figure, the residue, which is the difference between the original and filtered signals, is also presented. The random nature of this component and the attenuation of the noise in the EMG signal is apparent.

A system for extraction and visualization of MUAPs
The visual inspection of data is an important stage in signal analysis. It may help the investigator to identify relevant features, outliers and noise in the signal. Although visualization is an important and basic stage in data processing, in practice, its execution may be rather complex, specially when performed manually and the data lie in a high dimensional manifold.
In this context tools or systems that allow for an automatic visualization of signals play an important role. First, they significantly reduce the overall processing time of data, and  secondly they may reveal information present in the signal which might be barely perceived in a manual analysis.
This section describes a system that can be used for the extraction and visualization of motor unit action potentials from electromyographic signals. It is formed by several basic units which are illustrated in Figure 10. Its input is a raw electromyographic signal, and its output is the visualization of MUAPs or any other information (e.g. noise) present in the signal grouped into logical units (clusters). This visualization allows the researcher to identify outliers, noise and groups of MUAPs. Each of the stages of this system is described below.
Once the EMG signal is detected it is amplified and digitized. An EMG amplifier 2 and data acquisition board 3 are important elements that contribute to the correct collection of data. 2 Typical requirements of the signal conditioner are CMRR > 80 dB, input impedance > 10 15 Ω, noise level < 1. Other relevant units include force sensors for monitoring the level of muscle contraction and biofeedback systems that provide a stimulus and information for subjects about the required task to be executed during the experiment.

Signal detection and acquisition
The stage of signal detection will convert current generated by movement of ions (a physiological process) into current generated by movement of electrons (an electrical process). The choice of the type of electrode to be employed is dependent on the aim of the investigation. For instance, if the interest is to extract and visualize activity of motor unit action potentials from EMG signals, then needle electrodes are traditionally employed, mainly because of the high spatial resolution offered by these sensors. In fact, when studying the state of deep muscles this is the only possible choice. However, many studies [5,11,17,43,49] have shown that surface EMG electrodes with a sufficiently small contact surface may be successfully applied to the detection of MUAP activity in superficial muscles.
The system depicted in Figure 10 has been tested with two types of electrodes for signal detection. They are shown in Figure 11. One is a concentric needle electrode and the other is a customized array of pellet(surface) electrodes whose dimensions follow specifications provided in [43]. Both electrodes are passive and leads connecting them to the amplifier are shielded in order to attenuate the presence of spurious noise activity.
The use of pellet electrodes was a very cheap and simple solution for high spatial resolution signal detection, which yielded outstanding results and consistent detection of EMG signals.

Signal filtering
The main aim of this stage is to reduce the background activity present in electromyographic signals. This is relevant because high levels of background activity (noise) may affect the  detection of useful information (e.g. MUAPs). This stage may be performed by the filtering procedure based on the Empirical Mode Decomposition described above.

Detection of regions of activity
Following the stage of signal filtering the EMG signal is segmented into small windows, so-called regions of activity (RA), that may contain the activity of single MUAPs, MUAP overlaps or noise. A detector of RA was devised for extraction of RA from EMG signals. The block diagram of this is shown in Figure 12 and an example of its application is provided in Figure 13. Its input may be either a raw or filtered EMG signal and the output is a set of regions of activity. Basically, this system will estimate the signal envelope and it will select from it reference points which define the beginnings and ends of RAs. Intuitively, the envelope is the overall shape of the amplitude of the signal and its use instead of the raw signal may simplify, in terms of practical implementation, the solution for this problem, mainly because the signal envelope is positive. For instance, a unique positive threshold may be employed for selection of RAs which are above it.
The envelope of time domain signals can be computed by low pass filtering and curve fitting techniques [21], or by using the Hilbert transform to compute the analytic signal [13]. The signal envelope generated from the first two methods is dependent on the choice of parameters like the filter order and the method for data fitting, whereas in the latter strategy no a priori parameters should be set before calculation of the signal envelope. For this reason this is the method implemented in this system, that is, the absolute value of the complex analytic signal (also known as the complex envelope [13]) is used to estimate the signal envelope.
From the signal envelope, local minima and local maxima points are estimated. This will reduce the number of candidate points to be searched by the peak detector. As a consequence, the processing time will also be reduced, which may be relevant for the analysis of either long or oversampled time-series. Another benefit of this stage is that noise activity may be eliminated.
The role of the peak detector is to search for extrema points which are above or below a threshold. Once a maximum (w o ) is found, the mean of the small window, h o = w o + u, is estimated and only if this mean is above the pre-defined threshold, w o will be selected as the beginning of an RA. In this system u is set to 1 ms. Note that the analysis of the mean of h o may avert the selection of spurious activity (e.g. noise). The end of an RA is detected when a minimum (w f ) is found. In this case, w f is considered to be valid only if the mean of the window h o = w o − u is below the threshold.
The threshold, th, is estimated as th = k × std(W noise ), where W noise is a window of noise directly selected from the analysed signal, std(W noise ) is the standard deviation of W noise , and k is a user-defined constant which controls the threshold level. A typical value for k is 5. It is supported by the Chebyshev theorem [26]. This theorem states that, for any distribution, if k > 1 then the fraction of observations that fall within a range of ±k × std(W noise ) around the mean is at least F = 1 − (1/k 2 ). For instance, F = 0.96 for k = 5, which means in practice that any sample within the interval [mean(W noise ) − k × std(W noise ), mean(W noise ) + k × std(W noise )] has a 96% chance of really being noise.

Feature selection
In the feature selection stage, features will be selected for use in data clustering from a window within the region of activity. Figure 14 shows reference points in time (w o , p o , p, p f , w f ) for definition of this window.
w o and w f are respectively the time when the region of activity starts and ends. These points are estimated by the detector of regions of activity. p is the time when the highest peak in the envelope of the RA occurs. This point is also the point where the variation in amplitude of the signal is maximal within the RA. This is illustrated in Figure 13.
p o and p f are respectively the beginning and end of the window to be selected for analysis. These points are defined as follows: with t o set to 2 ms. Note that the width of the window defined by p o and p f may vary for different RAs. Therefore, interpolation (splines) was employed for selection of 40 samples from each window defined in the interval [p o , p f ]. This means that after feature selection each pattern is represented by a 40D vector with samples obtained via an interpolation procedure.

Data visualization and clustering
In order to ease the application of the sequence of steps detailed in Figure 10 a graphical user interface (GUI) was devised. The main GUI is shown in Figure 15. The system is capable of importing EMG data organized in columns of a text file and storing them in user-defined variables which are available in a list box. The main interface is organized into four logical sections (tabs) that should be accessed sequentially. Figure 16 shows the module which allows the user to filter the EMG signal. The filtering procedure based on the Empirical Mode Decomposition is available here. The result of the automatic detection of regions activity is given in the interface shown in Figure 17.
The results of the data clustering and visualization step are presented in the GUI shown in Figure 18. At this stage patterns are clustered by means of Generative Topographic Mapping (GTM) and data visualization is performed with the GTM grid [3,9]. For generation of the GTM grid, a GTM model with 25 Gaussian functions and 16 basis functions with a width of 1 is fitted to the data. The data can also be projected onto the two-dimensional space so that the user can visualize the distances of distinct groups of MUAPs (see Figure 19).

The application of EMG decomposition in the treatment of neuromuscular disorders triggered by stroke
Stroke, or cerebrovascular accident (CVA), affects a great number of individuals, and is the leading cause of disability among adults [10,24]. After the event, most individuals must deal with severe reduction of motor functionality [31].    Hemiplegia and spasticity are among the most common post-stroke motor deficits. In general, hemiplegia is characterized by an initial flaccid stage, with motor and sensory losses, in which the patient finds himself unable to sustain or move the affected limb. In many cases, the motor sequel evolves into spasticity, a stage characterized by muscle hypertonia. Hemiplegia and spasticity are closely related to the disuse of the affected limb and to secondary changes in muscles, such as: selective atrophy of fast twitch type II muscle fibers, abnormal recruitment of motor units, muscle contractures, and decreased cortical representation, due to disuse of the affected limb [28].
Several rehabilitation therapies have been used in an attempt to recover motor functionality, especially for upper limbs. The majority is based on the premise that the neural system can be retrained [20,34]. The ability of the neural system to adapt to a new structural and functional condition, as well as its response to a traumatic destructive injury or to subtle changes resulting from the processes of learning and memory, is called 'neuroplasticity" [1,12].
Recent studies have shown that behavioral experiments, important tools used in rehabilitation strategies based biofeedback techniques, have a strong impact on the motor cortical representation post-stroke. In this sense, we can also infer that it is possible to use biofeedback strategies for the modulation of neural plasticity, seeking the recovery of motor skills in rehabilitation protocols. The question is: how can we generate the appropriate information for feedback in a situation where the encephalic damage is manifested by an uncontrolled recruitment, or the lack of recruitment, of muscle fibers -leading to involuntary muscle hypotonia or hypertonia? Looking for a solution to this problem, a novel strategy, based in the assessment of the recruitment rate of motor units, is under test (see [39]), to generate control information for a multimodal biofeedback system. In this approach, instead of simply evaluating the process of muscle contraction, the researches decided to focus on the neural control over the muscles. However, the information regarding the recruitment of motor units are not readily available through standard surface EMG. Hence, the biofeedback relies heavily on a robust EMG decomposition system.
The strategy addressed by [39] is based on the discrimination and feature extraction of EMG signals (Motor Unit Action Potential) in order to control a virtual device. The biofeedback protocol immerses the patient in a virtual reality environment in which a representation of the affected limb will be presented. This virtual member is controlled according to the pattern of motor unit recruitment. Thus, although the spastic or hemiplegic patient lacks of proper voluntary control, it is expected that the system will be able to capture small voluntary changes in motor recruitment, as a result of his desire to do so. These changes can then be used to guide the movements of the virtual member. In so doing, the virtual feedback operates as a guide, indicating that the current mental strategy (neuromotor control) is correct and should be encouraged, reinforcing the process of neural reorganization (neuroplasticity).

Future directions
The widespread use of automatic tools for decomposing EMG signals is dependent on how easy is to detect MUAPs from the surface of the skin, on the identification of new applications, either online or offline, which might employ the results of the EMG decomposition, and also on the sharing and availability of developed tools to clinicians and researchers.
Ideally the sensors used for detecting MUAPs from the surface of the skin should be easy of applying and using. As a number of applications require the use of multiple sensors organized on an array the solution to this issue becomes more complex. Therefore, further research on this area, with the aim of developing sensors that allow for the reproduction of experiments is required.
Most applications found in the area of EMG decomposition are solely focused on the development of the automatic tool for decomposing the EMG signal, or in the classification and discrimination of MUAPs. Therefore, the identification of new useful applications are required in order to disseminate the relevance of the technique to other areas. For instance, the use of motor unit information could be employed in robotics, myofeedback, and human-machine interface development.
A major limitation of the results in the area of EMG decomposition is that the developed tools are not shared by researchers. The sharing of these tools together with EMG signal databases could be beneficial to the widespread use of the tools. Furthermore, this would allow researchers to objectively compare distinct solutions.
Considering the application of artefact removal from biomedical signals, the development of algorithms for reducing the influence of noise over signals in real time, without a priori knowledge about the origins and characteristics of the noise, are required. In this way the application of adaptive techniques such as Empirical Mode Decomposition should be further exploited.
It is also expected that the use of these filters can be part of the solution of the problem of EMG decomposition by directly decomposing the raw EMG signal into its motor unit action potential trains. If such filters are developed, EMG decomposition tools could be embedded in hardware speeding up the results from multichannel data.