**3. 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.

#### **3.1. Conventional methods for EMG signal noise removal**

#### *3.1.1. Low-pass differential filter*

4 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

part of the data set (normally the first 3 to 5 seconds), and then the non-classified data set is

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

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

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

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

grouped into one of the classes defined in the first step.

Computational Intelligence in Electromyography Analysis –

264

that it makes no assumptions about the data set to be grouped.

data points (i.e that some MUAPs are missing).

measure of the accuracy in real data decomposition.

**3. Artefact removal from EMG signals**

Since its introduction, the low-pass differential filter (LPD) [44] has been widely employed in EMG signal processing [15, 19, 32, 40–42]. This filter is implemented in the time-domain as:

$$y\_k = \sum\_{n=1}^{N} (\mathbf{x}\_{k+n} - \mathbf{x}\_{k-n})\_\prime \tag{1}$$

where *xk* is the discrete input time-series and *yk* 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:

$$y[k] = \sum\_{n=1}^{N} x[k+n] - \sum\_{n=1}^{N} x[k-n]. \tag{2}$$

Its representation in the frequency domain may be obtained via its Z transform as follows,

$$\begin{aligned} Z(y[k]) &= Z\left(\sum\_{n=1}^{N} x[k+n]\right) - Z\left(\sum\_{n=1}^{N} x[k-n]\right) \\ Z(y[k]) &= \sum\_{n=1}^{N} Z(x[k+n]) - \sum\_{n=1}^{N} Z(x[k-n]) \\ Z(y[k]) = Y(z) &= \sum\_{n=1}^{N} z^n X(z) - \sum\_{n=1}^{N} z^{-n} X(z) \\ Y(z) &= X(z) \left(\sum\_{n=1}^{N} (z^n - z^{-n})\right) \Big|\_{z = e^{\frac{j2\pi f}{f\sigma}}} \end{aligned}$$

where the ratio *Y*(*z*)/*X*(*z*) is the filter transfer function *H*(*z*), *f* is the frequency in Hz and *fsr* is the sampling frequency in Hz.

Figure 1 presents the results of the estimate of the filter transfer function with different sizes of windows, *N* = 40 and *N* = 20. Note how the cut-off frequency is shifted to a higher frequency when the window size is reduced.

**Figure 2.** Adaptive noise canceller. *x*[*n*] is the noise corrupted signal, *s*[*n*] is the noise free signal, *y*[*n*] the

EMG Decomposition and Artefact Removal 267

noise and *r*[*n*] the reference noise. *y*ˆ[*n*] is the estimated noise, whereas *s*ˆ[*n*] the estimated signal.

**Figure 3.** Example of application of the adaptive noise cancelling on a sample of EMG signal. The corrupted signal (top) is contrasted with the cleaned signal. The power spectra of the signals are shown

be present in the filtered signal and depending on the level of noise false spikes may be

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

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

at the bottom and illustrate the effect of the application of the filter on the 50 Hz component.

cancellation based on the least mean square (LMS) algorithm [2].

generated.

*3.1.3. Adaptive 50/60 Hz noise cancellation*

**Figure 1.** Frequency response of an LPD filter with window width 20 and 40.

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]:


#### *3.1.2. 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:

$$y\_k = \sum\_{n=1}^{N} w(n)(\mathbf{x}\_{k+n} - \mathbf{x}\_{k-n})\_\prime \tag{3}$$

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

**Figure 2.** Adaptive noise canceller. *x*[*n*] is the noise corrupted signal, *s*[*n*] is the noise free signal, *y*[*n*] the noise and *r*[*n*] the reference noise. *y*ˆ[*n*] is the estimated noise, whereas *s*ˆ[*n*] the estimated signal.

**Figure 3.** Example of application of the adaptive noise cancelling on a sample of EMG signal. The corrupted signal (top) is contrasted with the cleaned signal. The power spectra of the signals are shown at the bottom and illustrate the effect of the application of the filter on the 50 Hz component.

be present in the filtered signal and depending on the level of noise false spikes may be generated.

#### *3.1.3. Adaptive 50/60 Hz noise cancellation*

6 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

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

• Under conditions of low signal-to-noise ratio the relatively strong high frequency noise

• 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

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

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

*<sup>w</sup>*(*n*)(*xk*<sup>+</sup>*<sup>n</sup>* − *xk*−*n*), (3)

**Figure 1.** Frequency response of an LPD filter with window width 20 and 40.

*yk* =

*N* ∑ *n*=1

(background activity) may be accentuated;

Computational Intelligence in Electromyography Analysis –

noise components will pass through the filter.

*3.1.2. Weighted low-pass differential filter*

discussed in [46]:

266

time-domain as:

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.

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

EMG Decomposition and Artefact Removal 269

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 Matlab Wavelet Toolbox<sup>1</sup> is used for data processing (see Figure 4). The EMG signal is decomposed into three levels (*N* = 5) and the *db*5 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

The results of the signal decomposition are depicted in Figure 4. A window of the input EMG

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

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

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 *gi* in a time-series is a local maximum if *gi > gi*−<sup>1</sup> and *gi > gi*+1, and a sample *qi* is a local minimum if *qi < qi*−<sup>1</sup> and *qi < qi*+1, where *<sup>i</sup>* is a

of this procedure is illustrated in the following example.

the signal and at the same time to preserve its shape.

signal together with its filtered version are presented in Figure 5.

an alternative to the approach based on wavelets described above.

<sup>1</sup> This toolbox can be invoked by typing *wavemenu* in the Matlab command environment.

*3.2.1. The Empirical Mode Decomposition method*

functions that satisfy two conditions [21]:

discrete time.

**3.2. EMG signal filtering based on Empirical Mode Decomposition**

*3.1.5. Example of application*

MUAPs [18, 37, 48].

```
function [EstimatedNoise,filteredSig50Hz] = NoiseRemoval(x,fs)
%% Function name....: NoiseRemoval
% Description......:
% This function estimates a filtered signal based on a
% 50Hz noise optimal filter (LMS algorithm)
% Parameters.......:
% x ......-> input time-series
% fs......-> sampling frequency (Hz)
% Return...........:
% EstimatedNoise... -> estimated 50 Hz noise
% filteredSig50Hz.. -> filtered signal (50 Hz noise free signal)
mu = 0.001; %constant of convergence of the LMS algorithm
taps = 10; %filter order
t = [0:1:length(x)-1]/fs;
d = sin(2*pi*t*50); % reference noise (component to be removed from the input signal)
ha = adaptfilt.lms(taps,mu);
[y,e] = filter(ha,d,y);
[EstimatedNoise,filteredSig50Hz] = filter(ha,d,e);
```
#### *3.1.4. 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, 33]:


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, *tN*, is selected and soft-thresholding is applied to detail coefficients as shown in Equation 4.

$$Y\_N = \text{sign}(D\_N)(|D\_N| - t\_N)\_+ \tag{4}$$

where *YN* is the de-noised version of the *Nth* detail coefficients and the function (*x*)+ is defined as

$$\mathbf{u}(\mathbf{x})\_{+} = \begin{cases} \mathbf{0}, & \mathbf{x} < \mathbf{0} \\ \mathbf{x}, & \mathbf{x} > \mathbf{0}. \end{cases} \tag{5}$$

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.

#### *3.1.5. Example of application*

8 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

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.

%% Function name....: NoiseRemoval

% Description......:

268

% Parameters.......:

% Return...........:

taps = 10; %filter order t = [0:1:length(x)-1]/fs;

applied [14, 16, 33]:

as shown in Equation 4.

defined as

ha = adaptfilt.lms(taps,mu); [y,e] = filter(ha,d,y);

*3.1.4. Signal filtering based on wavelets*

1. Signal decomposition;

3. Signal reconstruction.

2. Detail coefficients thresholding;

function [EstimatedNoise,filteredSig50Hz] = NoiseRemoval(x,fs)

% 50Hz noise optimal filter (LMS algorithm)

% EstimatedNoise... -> estimated 50 Hz noise

mu = 0.001; %constant of convergence of the LMS algorithm

% x ......-> input time-series % fs......-> sampling frequency (Hz)

Computational Intelligence in Electromyography Analysis –

[EstimatedNoise,filteredSig50Hz] = filter(ha,d,e);

% This function estimates a filtered signal based on a

% filteredSig50Hz.. -> filtered signal (50 Hz noise free signal)

d = sin(2\*pi\*t\*50); % reference noise (component to be removed from the input signal)

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

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, *tN*, is selected and soft-thresholding is applied to detail coefficients

where *YN* is the de-noised version of the *Nth* detail coefficients and the function (*x*)+ is

0, *x <* 0

(*x*)+ =

*YN* = *sign*(*DN*)(|*DN*| − *tN*)+ (4)

*<sup>x</sup>*, *<sup>x</sup> <sup>&</sup>gt;*<sup>=</sup> 0. (5)

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<sup>1</sup> is used for data processing (see Figure 4). The EMG signal is decomposed into three levels (*N* = 5) and the *db*5 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.

## **3.2. 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.

#### *3.2.1. 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 *gi* in a time-series is a local maximum if *gi > gi*−<sup>1</sup> and *gi > gi*+1, and a sample *qi* is a local minimum if *qi < qi*−<sup>1</sup> and *qi < qi*+1, where *<sup>i</sup>* is a discrete time.

<sup>1</sup> This toolbox can be invoked by typing *wavemenu* in the Matlab command environment.


**Figure 4.** Graphical user interface of the Wavelet Matlab Toolbox used for signal denoising. The EMG signal and its denoised version are shown. The wavelet components are individually denoised by using a thresholding scheme.

**Figure 5.** Segment of the EMG signal shown in Figure 4. The filtered signal and the residue is also

procedure was proposed for EMG signal filtering [3, 7]:

filtered signal is obtained as a linear summation of thresholded IMFs.

1. Decompose the signal into IMFs; 2. Threshold the estimated IMFs;

3. Reconstruct the signal.

applied to individual IMFs.

based on wavelets may also be applied to intrinsic mode functions. Thus the following

EMG Decomposition and Artefact Removal 271

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,..., *IMFN*, where *N* is the number of IMFs. These IMFs are then soft-thresholded, yielding *tIMF*1,..., *tIMFN*, which are thresholded versions of the original components. The

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, *tn*, *n* = 1, . . . , *N*, is selected and soft-thresholding is

The threshold *tn* 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

presented.

2. At any point in the time-series, the mean value of the envelopes, one defined by the local maxima (upper envelope) and the other by the local minima (lower envelope), is zero.

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.

#### *3.2.2. 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

**Figure 5.** Segment of the EMG signal shown in Figure 4. The filtered signal and the residue is also presented.

based on wavelets may also be applied to intrinsic mode functions. Thus the following procedure was proposed for EMG signal filtering [3, 7]:


10 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

Computational Intelligence in Electromyography Analysis –

**Figure 4.** Graphical user interface of the Wavelet Matlab Toolbox used for signal denoising. The EMG signal and its denoised version are shown. The wavelet components are individually denoised by using

2. At any point in the time-series, the mean value of the envelopes, one defined by the local maxima (upper envelope) and the other by the local minima (lower envelope), is zero.

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

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

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

on the available data, whereas the other normally uses a set of pre-fixed filters.

*3.2.2. The algorithm for signal filtering based on the EMD*

a thresholding scheme.

270

to attenuate its presence.

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,..., *IMFN*, where *N* is the number of IMFs. These IMFs are then soft-thresholded, yielding *tIMF*1,..., *tIMFN*, 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, *tn*, *n* = 1, . . . , *N*, is selected and soft-thresholding is applied to individual IMFs.

The threshold *tn* 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 272 Computational Intelligence in Electromyography Analysis – A Perspective on Current Applications and Future Challenges EMG Decomposition and Artefact Removal <sup>13</sup>

12 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

**Figure 6.** The EMD method is employed for decomposing the input signal into IMFs (*IMF*1,..., *IMFN*, where *N* is the number of IMFs). These IMFs are soft-thresholded, yielding *tIMF*1,..., *tIMFN*, which are thresholded versions of the original components. The filtered signal is obtained as a linear summation of thresholded IMFs.

**Figure 8.** EMG signal and its filtered version. The residue, which is the difference between those signals

EMG Decomposition and Artefact Removal 273

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 (*IMF*1, . . . , *IMF*8) are obtained

The subsequent step is to threshold the components *IMF*1, . . . , *IMF*8. 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

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

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

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

of this component and the attenuation of the noise in the EMG signal is apparent.

**4. A system for extraction and visualization of MUAPs**

is also shown.

manifold.

and they are presented in Figure 9(a).

(0.03 s and 0.07 s) indicated in Figure 8 were selected.

**Figure 7.** Illustrative example showing the estimate of thresholds *tn*. First, a window of noise (i.e. samples within the interval [a,b]) is selected from the arbitrary time-series. The interval [a,b] is employed for selection of noise in the IMFs, and then from each IMF, thresholds *t*1,..., *tN* are defined. Note that *std* is the standard deviation.

these are regarded as the required thresholds (*t*1, ..., *tN*). Figure 7 illustrates the procedure for estimating *tn*.

#### *3.2.3. 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.

12 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

Computational Intelligence in Electromyography Analysis –

**Figure 6.** The EMD method is employed for decomposing the input signal into IMFs (*IMF*1,..., *IMFN*, where *N* is the number of IMFs). These IMFs are soft-thresholded, yielding *tIMF*1,..., *tIMFN*, which are thresholded versions of the original components. The filtered signal is obtained as a linear summation of

**Figure 7.** Illustrative example showing the estimate of thresholds *tn*. First, a window of noise (i.e. samples within the interval [a,b]) is selected from the arbitrary time-series. The interval [a,b] is employed for selection of noise in the IMFs, and then from each IMF, thresholds *t*1,..., *tN* are defined.

*3.2.3. An example of application of the filter based on the EMD method*

these are regarded as the required thresholds (*t*1, ..., *tN*). Figure 7 illustrates the procedure for

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

thresholded IMFs.

272

Note that *std* is the standard deviation.

signal shown in Figure 8 is used.

estimating *tn*.

**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 (*IMF*1, . . . , *IMF*8) are obtained and they are presented in Figure 9(a).

The subsequent step is to threshold the components *IMF*1, . . . , *IMF*8. 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.

#### **4. 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 Computational Intelligence in Electromyography Analysis –

14 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges 274 A Perspective on Current Applications and Future Challenges EMG Decomposition and Artefact Removal <sup>15</sup>

**Figure 11.** Types of electrodes for EMG signal detection. (A) A customized array of pellet electrodes. (B)

EMG Decomposition and Artefact Removal 275

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

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

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

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.

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

successfully applied to the detection of MUAP activity in superficial muscles.

shielded in order to attenuate the presence of spurious noise activity.

Example of a single pellet electrode. (C) Concentric needle electrode.

task to be executed during the experiment.

**4.1. System description**

*4.1.2. Signal filtering*

*4.1.1. Signal detection and acquisition*

**Figure 9.** (a) Intrinsic mode functions (*IMF*1, . . . , *IMF*8) obtained from the EMG signal presented in Figure 8. (b) De-noised intrinsic mode functions.

**Figure 10.** Block diagram of the system for extraction and visualization of MUAPs. After detection (1) and acquisition (2) the EMG signal is filtered (3), and signal windows, designated as regions of activity (RA), are extracted from it (4). In (5) features are selected from RA and then employed for their clustering and visualization.

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 amplifier2 and data acquisition board3 are important elements that contribute to the correct collection of data.

<sup>2</sup> Typical requirements of the signal conditioner are CMRR *>* 80 dB, input impedance *>* 10<sup>15</sup> Ω, noise level *<* 1.2 *μV*.

<sup>3</sup> A possible model would be PC-card DAS 16/16 Measurement Computing

#### A Perspective on Current Applications and Future Challenges EMG Decomposition and Artefact Removal <sup>15</sup> EMG Decomposition and Artefact Removal 275

**Figure 11.** Types of electrodes for EMG signal detection. (A) A customized array of pellet electrodes. (B) Example of a single pellet electrode. (C) Concentric needle electrode.

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.

## **4.1. System description**

14 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

**(a) (b)**

**Figure 9.** (a) Intrinsic mode functions (*IMF*1, . . . , *IMF*8) obtained from the EMG signal presented in

**Figure 10.** Block diagram of the system for extraction and visualization of MUAPs. After detection (1) and acquisition (2) the EMG signal is filtered (3), and signal windows, designated as regions of activity (RA), are extracted from it (4). In (5) features are selected from RA and then employed for their

secondly they may reveal information present in the signal which might be barely perceived

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

Once the EMG signal is detected it is amplified and digitized. An EMG amplifier2 and data acquisition board3 are important elements that contribute to the correct collection of data. <sup>2</sup> Typical requirements of the signal conditioner are CMRR *>* 80 dB, input impedance *>* 10<sup>15</sup> Ω, noise level *<* 1.2 *μV*.

and groups of MUAPs. Each of the stages of this system is described below.

<sup>3</sup> A possible model would be PC-card DAS 16/16 Measurement Computing

Figure 8. (b) De-noised intrinsic mode functions.

Computational Intelligence in Electromyography Analysis –

clustering and visualization.

in a manual analysis.

274

#### *4.1.1. 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.

#### *4.1.2. 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 Computational Intelligence in Electromyography Analysis –

16 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges 276 A Perspective on Current Applications and Future Challenges EMG Decomposition and Artefact Removal <sup>17</sup>

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

EMG Decomposition and Artefact Removal 277

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

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

The role of the peak detector is to search for extrema points which are above or below a threshold. Once a maximum (*wo*) is found, the mean of the small window, *ho* = *wo* + *u*, is estimated and only if this mean is above the pre-defined threshold, *wo* 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 *ho* may avert the selection of spurious activity (e.g. noise). The end of an RA is detected when a minimum (*wf* ) is found. In this case, *wf* is considered to be valid only if the mean of the

The threshold, *th*, is estimated as *th* = *k* × *std*(*Wnoise*), where *Wnoise* is a window of noise directly selected from the analysed signal, *std*(*Wnoise*) is the standard deviation of *Wnoise*, 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*(*Wnoise*) around the mean is at least *<sup>F</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> (1/*k*2). For instance, *<sup>F</sup>* <sup>=</sup> 0.96 for *<sup>k</sup>* <sup>=</sup> 5, which means in practice that any sample within the interval [*mean*(*Wnoise*) − *k* × *std*(*Wnoise*), *mean*(*Wnoise*) + *k* × *std*(*Wnoise*)] has

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 (*wo*, *po*, *p*, *pf* , *wf* ) for

*wo* and *wf* 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

of RAs which are above it.

envelope.

eliminated.

window *ho* = *wo* − *u* is below the threshold.

a 96% chance of really being noise.

*4.1.4. Feature selection*

definition of this window.

**Figure 12.** Block diagram of the detector of regions of activity.

**Figure 13.** Example of detection of regions of activity. The experimental surface EMG signal and its complex envelope are shown. *wo* and *wf* are respectively the beginning and end of a region of activity. The peak (*p*) within each RA is marked with ∗ and it is used as a reference point for feature selection. The inset shows a window of noise with the pre-defined threshold level.

detection of useful information (e.g. MUAPs). This stage may be performed by the filtering procedure based on the Empirical Mode Decomposition described above.

#### *4.1.3. 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 (*wo*) is found, the mean of the small window, *ho* = *wo* + *u*, is estimated and only if this mean is above the pre-defined threshold, *wo* 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 *ho* may avert the selection of spurious activity (e.g. noise). The end of an RA is detected when a minimum (*wf* ) is found. In this case, *wf* is considered to be valid only if the mean of the window *ho* = *wo* − *u* is below the threshold.

The threshold, *th*, is estimated as *th* = *k* × *std*(*Wnoise*), where *Wnoise* is a window of noise directly selected from the analysed signal, *std*(*Wnoise*) is the standard deviation of *Wnoise*, 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*(*Wnoise*) around the mean is at least *<sup>F</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> (1/*k*2). For instance, *<sup>F</sup>* <sup>=</sup> 0.96 for *<sup>k</sup>* <sup>=</sup> 5, which means in practice that any sample within the interval [*mean*(*Wnoise*) − *k* × *std*(*Wnoise*), *mean*(*Wnoise*) + *k* × *std*(*Wnoise*)] has a 96% chance of really being noise.

#### *4.1.4. Feature selection*

16 Computational Intelligence in Electromyography Analysis: A Perspective on Current Applications and Future Challenges

**Figure 13.** Example of detection of regions of activity. The experimental surface EMG signal and its complex envelope are shown. *wo* and *wf* are respectively the beginning and end of a region of activity. The peak (*p*) within each RA is marked with ∗ and it is used as a reference point for feature selection.

detection of useful information (e.g. MUAPs). This stage may be performed by the filtering

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

The inset shows a window of noise with the pre-defined threshold level.

*4.1.3. Detection of regions of activity*

regions of activity.

procedure based on the Empirical Mode Decomposition described above.

**Figure 12.** Block diagram of the detector of regions of activity.

Computational Intelligence in Electromyography Analysis –

276

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 (*wo*, *po*, *p*, *pf* , *wf* ) for definition of this window.

*wo* and *wf* 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

**Figure 14.** Identification of reference points in time for feature selection. *wo* and *wf* indicate the beginning and end of the region of activity. *p* is the time when the highest peak in the envelope and within the interval [*wo*, *wf* ] occurs, and *po* and *pf* define the boundaries of a rectangular window for feature selection.

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.

*po* and *pf* are respectively the beginning and end of the window to be selected for analysis. These points are defined as follows: *po* = *p* − *min*{*to*, *p* − *wo*} and *pf* = *p* + *min*{*to*, *wf* − *p*}, with *to* set to 2 ms. Note that the width of the window defined by *po* and *pf* may vary for different RAs. Therefore, interpolation (splines) was employed for selection of 40 samples from each window defined in the interval [*po*, *pf* ]. This means that after feature selection each pattern is represented by a 40D vector with samples obtained via an interpolation procedure.

#### *4.1.5. 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 15.** Main graphical user interface of a software that implements the sequence of steps detailed in

EMG Decomposition and Artefact Removal 279

**Figure 16.** Graphical user interface which allows the user to filter the input signal.

Figure 10.

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).
