A simple SSA-based de-noising technique to remove ECG interference in EMG signals

Biomedical Signal Processing and Control 30 (2016) 117–126

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

A simple SSA-based de-noising technique to remove ECG interference in EMG signals Jorge Barrios-Muriel a,∗ , Francisco Romero a , Francisco Javier Alonso a , Kostas Gianikellis b a b

Dept. of Mechanical, Energy and Materials Engineering, School of Industrial Engineering, University of Extremadura, Spain Dept. of Didactics of the Musical, Plastic and Corporal Expression, Faculty of Sport Sciences, University of Extremadura, Spain

a r t i c l e

i n f o

Article history: Received 7 August 2015 Received in revised form 31 May 2016 Accepted 3 June 2016 Keywords: Electromyography Biomedical signal processing Singular Spectrum Analysis ECG cancellation

a b s t r a c t Electromyography (EMG) signals provide significant information of muscle activity that may be used, among others, to estimate the activation stages during a certain activity or to predict fatigue. Heart activity or electrocardiogram (ECG) is one of the main contamination sources, especially in trunk muscles. This paper proposes a novel method based on Singular Spectrum Analysis (SSA) and frequency analysis to separate both signals present in the raw data. The performance of the method has been compared in time and frequency domains with traditional high-pass filtering or novel techniques such as Complete Ensemble Empirical Mode Decomposition or Wavelets analysis. The results show that for both time and frequency domains the proposed approach outperforms the other methods. Thus, the proposed SSA approach is a valid method to remove the ECG artifact from the contaminated EMG signals without using an ECG reference signal. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Surface electromyography (EMG) is a technique concerned with the development, acquisition and analysis of myoelectric signals. The main problems in EMG recordings arise from the environmental noise, as well as from other electrical signals within the body, such as the electrocardiogram (ECG) contamination, frequent in EMG recordings from trunk muscles [1,2]. The ECG signal influences both the amplitude and frequency content of the EMG signal [3,4]. These effects are reflected in the frequency domain by altering the lower frequencies of the power spectrum, and also by the distortion of the EMG signal amplitude in the time domain. Thus, the removal of ECG from EMG presents a major challenge in the subsequent extraction of accurate and useful information [5]. The ECG interference is characterized by a large amplitude and the overlapping of the EMG frequency causing a distortion of its frequency content. Conventional ECG removal procedures include high-pass filtering (HPF), usually employing finite impulse response, or Butterworth filters with a cut-off frequency of 30 Hz [6,7]. Nevertheless, these HPF techniques can fail since the ECG has a frequency spectrum that overlaps markedly with the EMG recording.

∗ Corresponding author at: Dept. of Mechanical, Energy and Materials Engineering, School of Industrial Engineering, University of Extremadura, Avda de Elvas S/N, 06006 Badajoz, Spain. E-mail address: [email protected] (J. Barrios-Muriel). http://dx.doi.org/10.1016/j.bspc.2016.06.001 1746-8094/© 2016 Elsevier Ltd. All rights reserved.

Hof [8] suggested a simple method that consists in the simultaneously acquisition of EMG and ECG signals. By combining a reference ECG measurement with a transfer function, the artifact signal can be removed from EMG. Alternatively, subtraction methods [9–12] and adaptive filters [4,5,13,14] have been proposed on the assumption that the ECG is independent of the EMG signal. Unfortunately, the interference between ECG and EMG signals depends on the location of the muscles and its position with respect to the heart [15]. Moreover, ECG removal methods based on the subtraction of a template signal require the identification of the time intervals corresponding to the ECG artifact in the EMG recordings. Thus, an automatic procedure based on this principle is not feasible, and semi-automatic ECG detection algorithms are required. Recently, new approaches based on advanced statistical vectorial decomposition have been presented to remove the ECG from the electromyogram, namely independent component analysis (ICA) [16–18], wavelet transform [19,20] and combination of ICA with the previous one [21]. However, the major drawback of these methods are the selection of the appropriate wavelet mother and the corresponding decision thresholding. At the same time, ICAbased methods require significant time processing, as well as an adequate number of input signals from multiple muscles to successfully isolate the ECG component [22]. A relatively new decomposition technique, the empirical mode decomposition (EMD) has been found to be very advantageous to

118

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

decompose non-stationary signals [26]. Its application to remove artifacts has been applied to de-noise ECG signals [23] and also combined with ICA [24]. The quality of the separation in noisy signals has even been improved with latest versions, first, the ensemble EMD (EEMD) [25] and then, the complete EEMD with adaptive noise (CEEMDAN) [26]. On the other hand, the use of non-parametric techniques such as Singular Spectrum Analysis (SSA) can be an effective and powerful tool to decompose the signal into a set of additive time series from which identify the interest signal from noise [27]. This methodology has been applied successfully to muscle activity onset detection and filtering of EMG signals in the time domain [28,29]. Nevertheless, this technique has not been applied to the frequency domain in ECG artifact removal, and may provide a powerful tool to automatically identify the frequency content of the ECG present in a contaminated EMG recording (EMGc ). In this paper we first describe the different advanced methods available in the literature. Afterwards, we present a novel method for ECG removal based on SSA in the frequency domain. The results, are then compared with the methodologies described. 2. Methods There are different methods in the literature to separate the ECG artifact from the EMG signal. Here, the current methodologies are presented and discussed. These are compared and contrasted against the SSA-based method proposed in this work that allows us to identify and separate the interest signal, i.e., EMG, from noise (ECG). 2.1. High pass filtering Traditional ECG removal procedures are based on high-pass filtering (HPF), as most of the ECG’s spectral power is located below 30 Hz. Thus, the recommendation to remove ECG artifact is to use HPF with a cut-off frequency of 30 Hz. In fact, it is the most common used method compared to all other methods [5]. In this study, we use a 4th order bi-directional high-pass Butterworth filter (HPF) to obtain a de-noised EMG signal (EMGu ).

and n = 1,. . .,N where N is the number of samples. Then, the final oscillatory modes are obtained by averaging: 1 IMFrk (n) R R

IMFk (n) =

This procedure improves the quality of the separation, and the mode mixing alleviation, but EEMD created some new problems, e.g., the decomposition is not complete and different realizations of signal plus noise may produce different number of modes. A variation of the EEMD algorithm has been proposed to correct these problems [32]. This method, namely “Complete Ensemble Empirical Mode Decomposition with Adaptive Noise” (CEEMDAN) adds a particular noise at each stage, and achieves a complete decomposition with no reconstruction error. The CEEMDAN procedure ensures the completeness property of the proposed decomposition and thus, provides an exact reconstruction of the original data, as well as a better spectral separation of the modes with a lower computational cost [32]. The CEEMDAN computes the first decomposition by averaging the R IMFs and a unique residue (res) is retained: res1 (n) = EMGc (n) − IMF1 (n)

(3)



where IMF1 is the same as the first IMF obtained by EMMD. Then EMD is performed over a set of signals obtained by adding different noise realization to res1 (n), so the next IMF (labeled IMF2  (n)) can be found by averaging the corresponding set of IMFs. The corresponding residue is, then, res2 (n) = res1 (n) − IMF2  , and so on, until the stop criterion is achieved. Once the IMFk (n) are obtained, x(n) can be expressed as:

x(n) =

K 

IMFk + R(n)

(4)

k=1

where the final residue R(n) satisfies:

R(n) = x(n) −

2.2. Complete ensemble empirical model decomposition

(2)

r=1

K 

IMFk

(5)

k=1

The EMD is a new signal analysis tool which is able to decompose any complex time series into a set of spectrally independent oscillatory modes [25]. These intrinsic mode functions (IMFs) are a finite and often a small number of amplitude and frequency modulated (AM/FM) zero-mean signals [30]. The advantage of EMD is the decomposition of the signal in a natural way where a priori knowledge of the signal of interest embedded in the data series is not needed. Despite its proven usefulness, one of the major drawbacks of the original EMD algorithm is its high sensitivity to noise and, in some cases, a mode-mixing problem due to the presence of very similar oscillations in different modes. To overcome these problems, a new method has been proposed: the Ensemble Empirical Mode Decomposition (EEMD) [31]. It performs the EMD over an ensemble of the signal plus Gaussian white noise. The EEMD algorithm defines the IMF set for an ensemble of trials (EMGr ), where each one was obtained by applying EMD to the interest signal (EMGc ) with independent, identically distributed and zero-mean, added white noise to the signal: EMGr (n) = EMGc (n) + wr (n), r = 1, ..., R

(1)

where wr (n) are different realizations of white noise with variance ε, producing IMFrk (n), R is the number of signals in the ensemble

with K total number of modes. The CEEMDAN version of the EEMD algorithm is used in this study. In order to select the ε value several simulations have been carried out with different values of ε (2, 0.2, 0.02 and 0.002), the best results were obtained for ε = 0.02, which was used in all simulations. To reconstruct EMGu , we used a number of IMFs between 6 and 8, which only depends on the noise of the signal. The number of simulation for CEEMDAN is executed 20 times per each value of Signal to Noise Ratio (SNR) and per each test, in order to obtain statistically meaningful results.

2.3. Wavelets analysis Wavelet analysis (WA) is a sophisticated technique for signal compression and noise reduction, which is very useful in different fields of biomedical signal processing. WA is based on similar ideas as the EMD algorithm. Both methods decompose the signal into components of disjointed spectra. Discrete wavelet transform (DWT) uses filter banks for construction of the multi-resolution analysis with relatively low computation time [33]. With this approach, the EMGc is passed through a low-pass filter and a high-

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

119

SSA

1. Embedding

2. SVD

3. Grouping Factor vector ( )

yes ,

Eigenvalue ( )

yes

∈

> 30

>

(

no

…

)

no ∈

4. Reconstruction

ECG

EMGu

Fig. 1. SSA algorithm and modifications in the grouping step to separate adequately ECG and EMG signals.

Fig. 2. (a) Raw EMG signal for dynamic test. (b) Contaminated EMG signal. (c) Raw EMG signal for isometric test (d) Contaminated EMG signal. SNR = 1 dB ( = 5).

120

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

pass filter (quadrature mirror filter) to obtain the approximate coefficients (cAi ) and the detail coefficients (cDi ):

⎧  j ⎪ cA(j+1) (n) = h[n − 2k]an , j = 1, 2, .., M − 1 ⎪ ⎨ n  (j+1) (n) = ⎪ g[n − 2k]ain , j = 1, 2, .., M − 1 cD ⎪ ⎩

(6)

n

where hj (n) and gj (n) are the low-pass and high-pass interpolating filters, corresponding to the selected wavelet mother function and j = 1, 2, 3,. . ., M is the decomposition level. This process splits the aj frequency information roughly in half, partitioning it into a set of detailed coefficients cD(j+1) and a set of approximation coefficients cA(j+1) . The procedure continues iteratively until the desired M-level decomposition is obtained. The decomposition level is at most [log2 L] where L is the length of the EMG signal. In this study the ECG noise is separated from the EMG by using Daubechies wavelet decomposition [20], [34–36], specifically by using db4 and db7 to obtain the best decomposition of the EMG signal [37]. To reconstruct the EMGu from a given wavelet coefficients, an universal threshold [34] is selected for each decomposition level and an inverse wavelet transform is performed. 2.4. Singular Spectrum Analysis SSA is a non-parametric technique of time series analysis based on principles of multivariate statistics. By applying this method it is possible to obtain a decomposition of the given signal into an additive set of independent time series. The resulting set of series can be interpreted as a trend representing the signal mean at each instant, a set of periodic series, and an aperiodic noise [27]. The procedure is usually divided into the following steps: embedding, singular value decomposition, grouping and reconstruction. In this work, we introduce several modifications in the grouping and reconstruction steps to separate adequately the ECG signal from EMG. Step 1. Embedding Let F = (f0 , f1 , . . ., fN−1 ) be the N-length time series representing the original signal. Let L be the window length, with 1 < L < N and L an integer. It is possible to construct a Hankel matrix, called the trajectory matrix by means of the sliding window that transverses the series: T

Xj = (fj−1 , fj , . . ., fj+L−2 ) , j = 1, 2, . . ., K

f0

f1

···

···

fN−L

⎞

⎜ f ⎟ f2 . . . . . . fN−L+1 ⎟ ⎜ 1 ⎜ ⎟ ⎜ . ⎟ .. .. .. . ⎟ .. . X=⎜ . . . . ⎜ ⎟ ⎜ ⎟ . . ⎜ ⎟ .. fN−1 ⎠ ⎝ fL−2 fL−3 . . fL−1

fL−2

···

···

ˆ i > rms( ˆ 1,  ˆ 2 , ...,  ˆ L) 

(9)

L

ˆ i = i / where   . This procedure ensures the selection of the i=1 i indices representing the highest level of information within the signal. The non-selected indices are included in the EMG subset, as in this method we interpret this signal as noise. The next stage is to analyze the frequency domain of the factor vector corresponding to the retained indices. To do so, given Vi corresponding to any of these indices, we propose, for each of them:

(7)

where K = N − L + 1 is the number of columns, i.e., the number of different possible positions of the said window. The matrix X = (X1 ,X2 ,. . .,XK ) is a Hankel matrix since all elements on the diagonal i + j = constant are equal, i.e., each ascending skew-diagonal from left to right is constant:

⎛

U1 ,U2 , . . ., Ud are the corresponding eigenvectors, and the vectors

Vi , namely the factor vectors, are obtained from Vi = XT · Ui / i for i = 1, 2, . . ., d. Step 3. Grouping The next step consists of partitioning the set of indices {1, ..., d} into m disjoint subsets I1 , . . ., Im . Let I = i1 , . . ., ip be one of these partitions. Then, the trajectory matrix EI corresponding to the set I is defined as EI = Ei1 + Ei2 + . . . + Eip . Once the matrices have been calculated for the partitions established, I1 , . . ., Im , the original time series trajectory matrix can be expressed as the sum of the trajectory matrices corresponding to each partition: X = EI = EI1 + EI2 + . . . + EIp . When dealing with high SNR signals, this step can be simplified by approximating matrix X by the summation of the first r elementary matrices, i.e., X ≈ E1 + E2 + . . . + Er , however, this simplification requires an adequate selection of the r parameter. In this work the following criteria was stablished to select the partitions (see Fig. 1). The aim of this step is to identify the ECG components in order to set two separate this signal from noise, i.e., the EMGu signal. By doing so, the computational cost is reduced considerably as the reconstruction is performed for only one group instead of for all the components. As explained in Ref. [27], if L  K then the factor vector in the eigentriple obtained in the SVD has a greater similarity with the component reconstructed from this eigentriple than the eigenvector. Therefore, the result of the reconstruction from each eigentriple can be predicted with the help of this factor vector. To select the appropriate components to reconstruct the ECG signal, the singular spectrum is studied. From the information given by the eigenvalues, expressed as percentage of the norm retained of matrix S, we first select those whose values are above the rms, i.e.:

1. Obtain the frequency at which the power spectral density (PSD) is maximum (fm ), using Welch’s overlapped averaged periodogram method by averaging 0.5 s, 50% overlapped Hamming-windowed epochs. 2. Include in the ECG subset the indices satisfying fm < 30 Hz as the frequency spectrum of individual QRS complex is found in the range of 0–30 Hz [38,39]. On the other hand, the indices satisfying fm > 30 Hz are included in the EMG subset. At the end of this step the following description is obtained:

(8)

fN

Step 2. Singular value decomposition (SVD) of the trajectory matrix. The trajectory matrix may be expressed as the summation of d rank-one elementary matrices X = E1 + E2 + . . . + Ed [27], where d is the number of non-zero eigenvalues in decreasing order 1 , 2 , . . ., d of the L × L matrix S = X · XT . The elementary matrices are given by Ei = i · Ui · Vi T i = 1, 2, . . ., d, where

X = EI = EI1 + EI2 = EIECG + EIEMG

(10)

Step 4. Reconstruction (diagonal averaging) At this step, each elementary matrix Ei is transformed into a principal component of length N by applying a linear transformation known as diagonal averaging or dehankelization. The elementary matrices are not themselves Hankel matrices. Thus, to reconstruct each principal component the average along the diagonals i + j = constant is calculated. The diagonal averaging algorithm [27] is as follows: let Y be any of the elementary matrices Ei of dimension L × K, the elements of which are yij , 1 ≤ i ≤ L, 1 ≤ j ≤ K.

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

121

Fig. 3. Information preserved after filtering respect the EMGraw . Top: dynamic tests. Bottom: isometric tests.

The time series G (principal component) corresponding to this elementary matrix is given by:

gk =

⎧ k+1 ⎪ 1  ⎪ ⎪ ym, k−m+2 ⎪ ⎪ k+1 ⎪ ⎪ m=1 ⎪ ⎪ L∗ ⎨ 1 ym, k−m+2

for 0 ≤ k < L∗ − 1 for L∗ − 1 ≤ k < K ∗

L∗ ⎪ ⎪ m=1 ⎪ ⎪ ∗ N−K ⎪ +1 ⎪ ⎪ 1 ⎪ ym, k−m+2 ⎪ ⎩ N−k

(11)

for K ∗ ≤ k < N

m=k−K ∗ +2

where L∗ = min(L, K), K ∗ = max(L, K)and N = L + K − 1. The expression in Eq. (11) corresponds to an averaging of the matrix elements over the diagonals i + j = k + 2. By using the modification proposed in the grouping step, only one reconstruction corresponding to the ECG group is performed. Therefore, the EMGu signal can be computed as: r = EMGu = F − G = F − ECG

(12)

where r represents the residual signal, i.e., the EMGu signal, F is the original signal and G the reconstructed component or ECG signal.

3. Simulations & data analysis To validate the effectiveness of the proposed approach, the results obtained using the SSA-based method are compared and contrasted against the procedures recommended in the literature. Namely, the results of HPF, CEEMDAN and multi-resolution wavelet analysis are compared with the proposed SSA-based algorithm. The EMG signal was recorded by means of standard passive surface electrodes (TrignoTM Wireless EMG System, Delsys Inc.) following the recommended standard given by the SENIAM [40]. A single voluntary subject was considered for this methodological study, since signal characteristic of ECG and EMG are similar between subjects [5]. Before placing the measurement electrode, the placement site was identified [9]. The electrodes sites were initially cleaned with isopropyl alcohol to remove oils and surface residues to avoid impedance mismatch and to improve the SNR. The EMG electrode was placed in the right biceps brachii and ECG electrode along the midsagittal line on bony aspects of the manubrium and upper sternum, with a reference electrode over the bony portion of the ipsilateral clavicle. This position was chosen to minimize the influence of muscle activity on the ECG signal. The recorded EMG signal corresponds to both dynamic and isometric contractions and is not influenced by ECG interference signal

122

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

Fig. 4. Power spectrum density of the denoised signal for SNR = 1 dB. Top: dynamic contractions. Bottom: isometric contractions.

[5]. Participant was informed of the purpose of the study before providing written consent in accordance with procedures approved by the University of Extremadura. The participant performed 20 dynamic and isometric (elbow at 90◦ of flexion) contractions test, respectively, during 10 s at 25% of the MVC and with at least 2 min rest between series. The signals were acquired and exported to PC by means of Delsys Base Station using the software provided by Delsys at a sampling frequency of 2000 Hz. The processing of the EMG signals is performed in MATLAB® , running on an Intel® CoreTM i5 CPU at 3.20 GHz. On the other hand, the different levels of contamination of the ECG signal were generated in MATLAB® [41]. Therefore, the contaminated EMG signal, was defined as:

evaluate the performance of the methods to retain the information present in the signal, we define the following index for the time domain:

EMGc (t) = EMGraw (t) +  · ECG (t)

  N 1 RMSE =  (EMGraw − EMGu )2

(13)

where EMGraw (t) is the acquired surface EMG signal, ECG(t) is the acquired electrocardiogram signal, and  is a scale factor that determines the energy ratio between the EMG and ECG component. For this study,  is selected as one of the values (0.5, 1,1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5) corresponding to the range 1–10 dB of input SNR. An example of a clean and a contaminated surface EMG signal for dynamic and isometric contractions is shown in Fig. 2. For each SNR, the de-noising methods are applied to the corresponding EMGc and an estimate of the EMGu is obtained. To

ˇ=

RMS(EMGraw ) RMS(EMGu )

(14)

The ˇ coefficient represent the amount of information preserved after filtering respect the EMGraw , i.e., it provides the amount of EMG signal that is retained in EMGu signals. The results for dynamic and isometric tests are shown in Fig. 3. In the same way, the root mean squared error (RMSE) can also be used to measure the difference between the raw signal and the filtered one in the time domain. The RMSE is defined as:

N

(15)

j=1

Ideal values for the RMSE, in this approach, are zero. The results are listed in Tables 1 and 2 and discussed in the following section. Regarding the frequency domain, the commonly used quantitative method for describing an EMG signal is the PSD, which is routinely performed by the periodogram method. To assess quan-

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

123

Table 1 Performance of the different methods expressed in terms of mean and standard deviation (SD) of RMSE (mV × 10−5 ) for different SNR in dynamic contractions. Best method for each SNR level is highlighted.

SNR level

Method HPF Wavelet “db4” Wavelet “db7” CEEMDAN SSA L=500 SSA L=250 SSA L=100 SSA L=50

1 4,873 (±0,2840) 1,4048 (±0,1201) 1,2607 (±0,1223) 3,2802 (±0,1801) 1,8559 (±0,0335) 0,4832 (±0,0584) 0,7628 (±0,0582) 1,2581 (±0,0781)

2 4,7889 (±0,2893) 1,3867 (±0,1216) 1,2568 (±0,1225) 3,165 (±0,1743) 1,6712 (±0,0336) 0,4684 (±0,0583) 0,7589 (±0,0588) 1,2551 (±0,0779)

3 4,7124 (±0,2954) 1,3703 (±0,1230) 1,2518 (±0,1227) 3,0516 (±0,1861) 1,4867 (±0,0337) 0,4553 (±0,0580) 0,7563 (±0,0601) 1,2528 (±0,0776)

4 4,6438 (±0,3006) 1,3557 (±0,1243) 1,248 (±0,1228) 2,9154 (±0,1539) 1,3027 (±0,034) 0,4445 (±0,0576) 0,7576 (±0,0635) 1,2500 (±0,0771)

5 4,5836 (±0,3054) 1,3429 (±0,1254) 1,2447 (±0,1229) 2,8061 (±0,1734 ) 1,0812 (±0,1041) 0,4387 (±0,0573) 0,7795 (±0,0865) 1,2467 (±0,0767)

6 4,5321 (±0,3098) 1,3319 (±0,1265) 1,2419 (±0,1225) 2,6799 (±0,1621) 0,7792 (±0,1513) 0,4889 (±0,0984) 0,6996 (±0,0631) 1,2433 (±0,0776)

7 4,4895 (±0,3136) 1,3228 (±0,1274) 1,2396 (±0,1221) 2,5533 (±0,1659) 0,5096 (±0,1911) 0,4902 (±0,0474) 0,5856 (±0,0512) 1,2634 (±0,1135)

8 4,4562 (±0,3167) 1,3157 (±0,1281) 1,2378 (±0,1221) 2,3646 (±0,1698) 0,6733 (±0,1743) 0,4098 (±0,0429 ) 0,5086 (±0,0536) 1,2044 (±0,1348)

9 4,4323 (±0,3192) 1,3105 (±0,1287) 1,2364 (±0,1228) 2,1222 (±0,1588) 0,5629 (±0,1312) 0,4041 (±0,1971) 0,4710 (±0,0712) 0,9824 (±0,2454)

10 4,4180 (±0,3210) 1,3074 (±0,1291) 1,2354 (±0,1201) 1,9055 (±0,1395) 1,0178 (±0,2686) 1,8556 (±0,2893) 1,4784 (±0,058) 1,5877 (±0,0582)

Table 2 Performance of the different methods expressed in terms of mean and standard deviation (SD) of RMSE (mV × 10−5 ) for different SNR in isometric contractions. Best method for each SNR level is highlighted.

SNR level

Method HPF Wavelet “db4” Wavelet “db7” CEEMDAN SSA L=500 SSA L=250 SSA L=100 SSA L=50

1 4,7246 (±2,2157) 9.,5069 (±0,8608) 9,1294 (±0,4161) 10,3546 (±1,0196) 2,8634 (±0,2005) 2,2115 (±0,2527) 3,3952 (±0,9279) 6,2115 (±1,4558)

2 4,722 (±2,2179) 9,503 (±0,8611) 9,1291 (±0,4162) 10,3332 (±1,0530) 2,8829 (±0,1962) 2,2847 (±0,3812) 3,7168 (±1,2208) 5,0997 (±0,2305)

3 4,7196 (±2,2199) 9,4995 (±0,8613) 9,1287 (±0,4163) 10,2868 (±1,0358) 3,0645 (±0,1962) 4,2409 (±1,3998) 3,5771 (±1,0033) 5,0657 (±1,0754)

4 4,7176 (±2,2216) 9,4964 (±0,8616) 9,1284 (±0,4164) 10,2438 (±1,0328) 3,3553 (±0,9204) 4,7454 (±1,4878) 3,9896 (±1,5867) 4,6659 (±0,4270)

5 4,7158 (±2,2231) 9,4937 (±0,8618) 9,1282 (±0,4165) 10,212 (±0,9732) 4,2849 (±1,1787) 4,3918 (±1,2270) 5,8756 (±1,9992) 6,114 (±2,1054)

6 4,7143 (±2,2244) 9,4915 (±0,8620) 9,1279 (±0,4166) 10,1801 (±0,8703) 5,293 (±1,4094) 4,7751 (±0,7554) 7,1918 (±1,2986) 7,315 (±2,0788)

titatively the difference between the EMGraw and the EMGu signals, the magnitude squared coherence is defined as [4]:

f 2 =

f 2 f1

f1

(P (EMGraw ))2

(16)

(P(EMGraw )P(EMGu ))

where P(EMGraw ) and P(EMGu ) are the PSD of the acquired and filtered EMG signal, respectively. Specifically, PSD is calculated using Welch’s overlapped averaged periodogram method by averaging 0.5 s, 50% overlapped Hamming-windowed epochs. The ideal value for , in this approach, is one. The results of the coherence analysis are listed in Tables 3 and 4. 4. Results and discussion The performance of the methods to retain the information present in the acquired EMG signal is shown in Fig. 3. Low values of SNR correspond to signals with a high level of contamination. Therefore, major difficulties are expected for the different methods to suppress the noise, and consequently, a major variability is observed in the results. Results show two tendencies: first, the previously commented higher variability for lower SNR. Second, the beta value decreases as the noise decreases (increase in SNR). This fact confirms that for lower levels of noise, the information in the input signal is quite similar to the output signal, and therefore the general trend correspond to a decreasing in the beta index, close to 1 as the noise level is reduced to the minimum. The results show that the SSA approach, using L = 500 and L = 250 window lengths,

7 4,7131 (±2,2255) 9,4896 (±0,8621) 9,1278 (±0,4166) 10,164 (±0,8839) 6,2842 (±1,3002) 5,4167 (±1,5048) 5,7323 (±1,5955) 5,3004 (±1,0415)

8 4,7121 (±2,2263) 9,4881 (±0,8622) 9,1276 (±0,4166) 10,1251 (±1,0909) 5,1555 (±2,2324) 5,7261 (±1,6603) 6,7301 (±1,9037) 3,9884 (±0,5140)

9 4,7114 (±2,2269) 9,4871 (±0,8623) 9,1275 (±0,4166) 10,1171 (±0,8789) 5,407 (±1,0924) 5,6645 (±1,2282) 6,321 (±0,8806) 5,8414 (±1,4787)

10 4,7111 (±2,2273) 9,5865 (±0,8623) 9,1274 (±0,4166) 10,1061 (±1,2707) 3,8062 (±0,5553) 3,3030 (±0,2048) 3,2745 (±0,0409) 4,1517 (±0,0781)

retains more information than other methods. The WA approach provides similar results that CEEMDAN for dynamic contractions. However, a better performance is shown in dynamics contractions. The performance of HPF method is clearly worse than other methods in dynamic tests in most cases, however, it shows better results, in isometric test, than CEEMDAN and Wavelet db4. The behavior of all methods, as expected, improve as the SNR increases, i.e., when the amplitude of the interference (ECG signal) is lower. These results are consistent with previous studies which compare different methods to remove ECG interference in EMG signals [19,42]. Figs. 4 and 5 show the PSD of the EMGraw and the PSD of the EMGu signals for input SNR = 1 dB and SNR = 10 dB respectively. In both figures, the existence of high differences between the denoised and cleaned data in the low frequency can be observed. The proposed SSA approach adequately fit the original spectrum from 20 Hz, as shown in Figs. 4 and 5. WA and CEEMDAN methods removed low frequency components, leading to a loss of information in the frequency content of the EMG signal. Fig. 5 shows that for high SNR values, the HPF and CEEMDAN method do not adequately remove the ECG component, and low frequency components remain in the PSD, thus this method at high levels of contamination is not valid. Low frequencies may be retained due to a poor decomposition of the EMGc signal in both methods, allowing the preservation of some energy of the EMG signal in the ECG artifact and therefore that remains in the reconstruction. The SSA method for high levels of decomposition (L = 500) may not adequately remove the ECG, as there is still energy in the range 0–30 Hz, but from 20 to 450 Hz, which is the frequency range of the EMG

124

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

Table 3 Magnitude-squared coherence () values expressed in terms of the mean and standard deviation (SD) for dynamics contractions. Best method for each SNR level is highlighted.

Method HPF Wavelet “db4” Wavelet “db7” CEEMDAN SSA L=500 SSA L=250 SSA L=100 SSA L=50

Magnitude-squared coherence (σ) 1 0,9130 (±0,1458) 0,3926 (±0,4371) 0,6625 (±0,4187) 0,3255 (±0,4193) 0,9087 (±0,2359) 0,9210 (±0,2769) 0,8709 (±0,2874) 0,7545 (±0,3343)

2 0,9133 (±0,1454) 0,3944 (±0,4376) 0,6650 (±0,4180) 0,3273 (±0,4202) 0,9281 (±0,2261) 0,9227 (±0,2765) 0,8825 (±0,2919) 0,7655 (±0,3257)

3 0,9136 (±0,1451) 0,3965 (±0,4381) 0,6677 (±0,4171) 0,3313 (±0,4232) 0,9685 (±0,2060) 0,9173 (±0,2776) 0,8950 (±0,2788) 0,7798 (±0,3152)

4 0,9139 (±0,1445) 0,3988 (±0,4391) 0,6706 (±0,4162) 0,3345 (±0,4248) 0,9477 (±0,2269) 0,9169 (±0,2776) 0,8507 (±0,2466) 0,7923 (±0,3034)

5 0,9142 (±0,1439) 0,4015 (±0,4395) 0,6738 (±0,4153) 0,3384 (±0,4267) 0,9121 (±0,2395) 0,9220 (±0,2757) 0,8821 (±0,2223) 0,8583 (±0,2501)

6 0,9145 (±0,1431) 0,4046 (±0,4401) 0,6774 (±0,4142) 0,3432 (±0,4279) 0,9423 (±0,2528) 0,9212 (±0,2765) 0,8734 (±0,2204) 0,7678 (±0,3411)

7 0,9149 (±0,1421) 0,4083 (±0,4404) 0,6814 (±0,4129) 0,3462 (±0,4299) 0,9670 (±0,2623) 0,9608 (±0,2594) 0,8530 (±0,2529) 0,7813 (±0,3932)

8 0,9154 (±0,1408) 0,4130 (±0,4403) 0,6860 0,4114) 0,3503 (±0,4322) 0,9781 (±0,2459) 0,9285 (±0,2376) 0,8803 (±0,2326) 0,8117 (±0,3740)

9 0,9160 (±0,1389) 0,4194 (±0,4407) 0,6914 (±0,4095) 0,3533 (±0,4328) 0,9804 (±0,2115) 0,9593 (±0,2518) 0,9011 (±0,2339) 0,9764 (±0,1139)

10 0,9175 (±0,1349) 0,4289 (±0,1458) 0,6982 (±0,4079) 0,3567 (±0,4187) 0,9507 (±0,2076) 0,9768 (±0,1233) 0,9860 (±0,1682) 0,9861 (±0,1068)

Table 4 Magnitude-squared coherence () values expressed in terms of the mean and standard deviation (SD) for isometric contractions. Best method for each SNR level is highlighted.

Method HPF Wavelet “db4” Wavelet “db7” CEEMDAN SSA L=500 SSA L=250 SSA L=100 SSA L=50

Magnitude-squared coherence (σ) 1

2

3

4

5

6

7

8

9

10

0,9308 (±0,1368) 0,4210 (±0,4022) 0,6670 (±0,4023) 0,3514 (±0,3606) 0,9515 (±0,2313) 0,8457 (±0,1322) 0,7880 (±0,3328) 0,6269 (±0,3227)

0,9336 (±0,1347) 0,4229 (±0,4028) 0,6699 (±0,4016) 0,3552 (±0,3648) 0,9527 (±0,2332) 0,8492 (±0,2324) 0,7909 (±0,3324) 0,6190 (±0,3221)

0,9366 (±0,1327) 0,4251 (±0,4045) 0,6730 (±0,4006) 0,3619 (±0,3658) 0,9531 (±0,2335) 0,8532 (±0,2327) 0,7957 (±0,3310) 0,6255 (±0,3177)

0,9387 (±0,1310) 0,4276 (±0,4043) 0,6766 (±0,3995) 0,3649 (±0,3701) 0,9561 (±0,2337) 0,8579 (±0,2329) 0,7041 (±0,3267) 0,6456 (±0,3105)

0,9408 (±0,1294) 0,4305 (±0,4052) 0,6807 (±0,3980) 0,3740 (±0,3727) 0,9607 (±0,2341) 0,8639 (±0,2330) 0,7120 (±0,3255) 0,6512 (±0,3082)

0,9427 (±0,1279) 0,4340 (±0,4061) 0,6856 (±0,3963) 0,3783 (±0,3767) 0,9611 (±0,2342) 0,8745 (±0,2325) 0,7252 (±0,3196) 0,6410 (±0,3080)

0,9444 (±0,1265) 0,4382 (±0,4073) 0,6914 (±0,3942) 0,3869 (±0,3802) 0,9635 (±0,2285) 0,8768 (±0,2320) 0,7632 (±0,3079) 0,6370 (±0,3116)

0,9458 (±0,1249) 0,4436 (±0,4087) 0,6987 (±0,3917) 0,3941 (±0,3876) 0,9606 (±0,2332) 0,8826 (±0,2315) 0,6677 (±0,3983) 0,6997 (±0,3853)

0,9470 (±0,1232) 0,4512 (±0,4106) 0,7082 (±0,3889) 0,4124 (±0,3905) 0,9719 (±0,2322) 0,8012 (±0,2238) 0,6857 (±0,3958) 0,6776 (±0,3989)

0,9482 (±0,1207) 0,4636 (±0,4136) 0,7212 (±0,3851) 0,4302 (±0,3958) 0,9710 (±0,2289) 0,8980 (±0,2307) 0,6514 (±0,3671) 0,7321 (±0,3761)

signal, the SSA method fits adequately the raw signal and therefore, the information relative to the EMG is maintained. This is also proven with values of Magnitude-squared coherence () close to 1 shown in Tables 3 and 4. These results shown as the proposed method is the one that best fits with the original spectrum for all the SNR values. Thus, this approach can be used to remove the ECG artifact in any context. An advantage of the proposed method is that no reference electrode of ECG is necessary for the removal of the ECG artifact. Moreover, in comparison to other methods, in the proposed scheme only one parameter, namely the window length, must be chosen compared to other methods as the reconstruction step has been automated. Tables 1 and 2 show RMSE values for dynamic and isometric contractions, respectively. These results provide information about the behavior of the signal in the time domain. These score provides a relative performance of each method for each SNR level. Higher values of RMSE indicate less ability to remove ECG artifact from the EMGc . In dynamic contractions (Table 1) and overall, the SSA method provides the best results, regardless of the value of the window length (L). It is worth to mention that for some low SNR levels, if L is significantly low, the decomposition may lose resolution and therefore the ability of the method to retain the interest signal as stated in our previous work [29]. Signal/noise grouping and reconstruction are critical steps in SSA that underlie any application of this method. From the preceding description it is apparent that the steps depend on one basic parameter that must be assigned or chosen by the practitioner, i.e., the window length. Clearly the values chosen for L will affect the performance and it is therefore necessary to ensure that the tech-

niques used for the selection of L will lead to strong separability and minimization (in some sense) of the error in the reconstruction step. For this reason, the use of different values of L is recommended in order to determine the behavior of the de-noising as a function of SNR and then to select the window length that provides the best result. In this study, we used different values of L (50, 100, 250, 500) and the best results were obtained between 250 and 500, where the resolution given by the method was the highest. Therefore it is recommended to use high values of L to obtain best results, i.e., at least a quarter of the sampling frequency [29]. Attending to the other methods, CEEMDAN, on the contrary, shows better performance for higher values of SNR and therefore, its use is recommended for small ECG contamination levels. Values in HPF method remain constant for different SNR values and in both situations (dynamic and isometric contractions) as a consequence of the linearity of the method. Regarding isometric contractions, Wavelet and CEEMDAN give the poorest results as both methods were proposed for stationary signals. SSA method provides similar results for isometric and dynamic contractions as no assumption is made on the input signal. Tables 3 and 4 illustrate the Magnitude-squared coherence () of PSDu respect of PSDraw . These results show the behaviour of the different methods in the frequency domain. The SSA method, especially for L = 500, shows the best results for both situations. The mentioned window length provides the highest resolution for the sampling frequency used. Overall, in dynamic contractions, as SNR rises, results improve for any SSA test. Regarding other methods, HPF provides better results than Wavelet and CEEMDAN in isometric contraction, as these methods are not able to separate

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

125

Fig. 5. Power spectrum density of the denoised signal for SNR = 10 dB. Top: dynamic contractions. Bottom: isometric contractions.

Fig. 6. EMG and ECG reconstructions from a contaminated signal for SNR = 1 dB (left) and SNR = 10 dB (right) using the proposed methodology (L = 500). Top: contaminated signal. Middle: reconstructed EMG signal. Bottom: reconstructed ECG signal.

126

J. Barrios-Muriel et al. / Biomedical Signal Processing and Control 30 (2016) 117–126

adequately EMG and ECG in non-stationary tests. At this point it is worth to mention that the proposed methodology provides a solution not only to remove ECG artifact from EMG signals, but also to separate and reconstruct both signals (Fig. 6). The level of ECG contamination in EMG signals is highly dependent on the placement of the EMG electrodes which is often related to the position of the muscle group with respect to the trunk. Such contamination may become more severe in certain muscles such as pectoralis major, rectus abdominis and left trapezius. In this sense, as the proposed SSA approach is not influenced by the SNR level, the performance of the method should be kept constant when applied to a variety of muscle groups. Moreover, the method does not alter the spectral properties of the EMGu , thus the obtained results can be used too in latter time-frequency analysis to gather useful information. Finally, the execution time on an Intel processor 3.20 GHz has been evaluated. HPF proved to be the fastest method while the CEEMDAN and SSA approaches are the slowest. Further research is needed if de-noising EMG signals in real time is desired.

5. Conclusion In this paper, we have proposed a novel method based on SSA to de-noise the EMG signal contaminated with ECG artifacts. The performance of the method has been presented in the results, and has also been compared in time and frequency domains to other available algorithms used to decompose single-channel signals, HPF, WA and CEEMDAN. In general, the performance of the methods depends on the SNR level. The results show that in both time and frequency domains the proposed approach outperforms the other methods, as it introduces the smallest distortion in the EMG signal and its spectrum. Accordingly, we can conclude that the proposed SSA approach is a valid method not only to remove the ECG artifact from the contaminated EMG signals without using an ECG reference signal, but also to separate both signals. The resulting de-noised EMG signals retain the spectral properties required for a subsequent time-frequency analysis.

References [1] C.J. De Luca, The use of surface electromyography in biomechanics, J. Appl. Biomech. 13 (1997) 1–38. [2] D. Farina, C. Cescon, R. Merletti, Influence of anatomical, physical, and detection-system parameters on surface EMG, Biol. Cybern. 86 (6) (2002) 445–456. [3] Y. Hu, J.N. Mak, K.D. Luk, Effect of electrocardiographic contamination on surface electromyography assessment of back muscles, J. Electromyogr. Kinesiol. 19 (1) (2009) 145–156. [4] G. Lu, J.-S. Brittain, P. Holland, J. Yianni, A.L. Green, J.F. Stein, T.Z. Aziz, S. Wang, Removing ECG noise from surface EMG signals using adaptive filtering, Neurosci. Lett. 462 (2009) 14–19. [5] J. Drake, J. Callaghan, Elimination of electrocardiogram contamination from electromyogram signals: an evaluation of currently used removal techniques, J. Electromyogr. Kinesiol. 16 (2) (2006) 175–187. [6] M.S. Redfern, R.E. Hughes, D.B. Chaffin, High-pass filtering to remove electrocardiographic interference from torso EMG recordings, Clin. Biomech. 8 (1) (1993) 44–48. [7] R.G.T. Mello, L.F. Oliveira, J. Nadal, Digital Butterworth filter for substracting noise from low magnitude surface electromyogram, Comput. Methods Programs Biomed. 87 (1) (2007) 28–35. [8] A.L. Hof, A simple method to remove ECG artifacts from trunk muscle EMG signals, J. Electromyogr. Kinesiol. 19 (6) (2009) e554–e555. [9] P.R.B. Barbosa, J. Barbosa-Filho, C.A.M. De Sá, E.C. Barbosa, J. Nadal, Reduction of electromyographic noise in the signal-averaged electrocardiogram by spectral decomposition, IEEE Trans. Biomed. Eng. 50 (1) (2003) 114–117. [10] M. Zivanovic, M. Gonzlez-Izal, Nonstationary harmonic modeling for ECG removal in surface EMG signals, IEEE Trans. Biomed. Eng. 59 (6) (2012) 1633–1640. [11] Y. Deng, W. Wolf, New aspects to event-synchronous cancellation of ECG interference: an application of the method in diaphragmatic EMG signals, Biomed. Eng. IEEE Trans. 47 (9) (2000) 1177–1184.

[12] A. Bartolo, C. Roberts, Analysis of diaphragm EMG signals: comparison of gating vs. subtraction for removal of ECG contamination, J. Appl. Physiol. 80 (6) (1996) 1898–1902. [13] S. Sanei, T.K.M. Lee, V. Abolghasemi, A new adaptive line enhancer based on singular spectrum analysis, IEEE Trans. Biomed. Eng. 59 (2) (2012) 428–434. [14] I. Christov, I. Daskalov, Filtering of electromyogram artifacts from the electrocardiogram, Med. Eng. Phys. 21 (10) (1999) 731–736. [15] M. Nitzken, N. Bajaj, S. Aslan, Local wavelet-based filtering of electromyographic signals to eliminate the electrocardiographic-induced artifacts in patients with spinal cord injury, J. Biomed. Sci. Eng. 6 (7B) (2013) 1–32. [16] B. Azzerboni, Neural-ICA and wavelet transform for artifacts removal in surface EMG, 2004 IEEE Int. Jt. Conf. on Neural Networks 2004. Proceedings 4 (2004) 3223–3228. [17] J.N.F. Mak, Y. Hu, K.D.K. Luk, An automated ECG-artifact removal method for trunk muscle surface EMG recordings, Med. Eng. Phys. 32 (8) (2010) 840–848. [18] X. Ren, Z. Yan, Z. Wang, X. Hu, Noise reduction based on ICA decomposition and wavelet transform for the extraction of motor unit action potentials, J. Neurosci. Methods 158 (2) (2006) 313–322. [19] H.Y. Lin, S.Y. Liang, Y.L. Ho, Y.H. Lin, H.P. Ma, Discrete-wavelet-transform-based noise removal and feature extraction for ECG signals, IRBM 35 (6) (2014) 351–361. [20] V. von Tscharner, B. Eskofier, P. Federolf, Removal of the electrocardiogram signal from surface EMG recordings using non-linearly scaled wavelets, J. Electromyogr. Kinesiol. 21 (4) (2011) 683–688. [21] J. Taelman, Wavelet-independent component analysis to remove electrocardiography contamination in surface electromyography, Eng. Med. Biol. Soc. 2007. EMBS 2007 29th Annu. Int. Conf. IEEE (2007) 682–685. [22] N. Willigenburg, A. Daffertshofer, Removing ECG contamination from EMG recordings: A comparison of ICA-based and other filtering procedures, J. Eectromyogr. Kinesiol. 22 (3) (2012) 485–493. [23] M. Agarwal, R.C. Jain, Ensemble empirical mode decomposition: an adaptive method for noise reduction, IOSR J. Electron. Commun. Eng. 5 (5) (2016) 2278–8735. [24] P. Torres, M.E. Colominas, M.A. Schlotthauer, G. Flandrin, A complete ensemble empirical mode decomposition with adaptive noise, Acoustics Speech and Signal Processing (ICASSP) (2011). [25] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. A Math. Phys. Eng. Sci. 454 (March (1971)) (1998) 903–995. ´ M. De Vos, I. Gligorijevic, ´ J. Taelman, S. Van Huffel, Source [26] B. Mijovic, separation from single-channel recordings by combining empirical-mode decomposition and independent component analysis, IEEE Trans. Biomed. Eng. 57 (9) (2010) 2188–2196. [27] N. Golyandina, V. Nekrutkin, A.A. Zhigljavsky, Analysis of Time Series Structure: SSA and Related Techniques, vol. 90, Chapman & Hall/CRC, 2001. [28] P. Zhou, X. Zhang, A novel technique for muscle onset detection using surface EMG signals without removal of ECG artifacts, Physiol. Meas. 35 (2014) 45–54. [29] F. Romero, F.J. Alonso, J. Cubero, G. Galán-Marín, An automatic SSA-based de-noising and smoothing technique for surface electromyography signals, Biomed. Signal Process. Control 18 (April) (2015) 317–324. [30] L. Cohen, Time-frequency distributions-a review, Proc. IEEE 77 (7) (1989) 941–981. [31] Z. Wu, N.E. Huang, Ensemble empirical mode decomposition: a noise-assisted data analysis method, Adv. Adapt. Data Anal. 1 (01) (2009) 1–41. [32] M.A. Colominas, G. Schlotthauer, M.E. Torres, P. Flandrin, Noise-assisted emd methods in action, Adv. Adapt. Data Anal. 4 (4) (2012) 1250025. [33] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, 1999. [34] D. Moshou, I. Hostens, Wavelets and self-organising maps in electromyogram (EMG) analysis, Eur. Symp. Intell. Tech. (2000) 186191. [35] A. Phinyomark, C. Limsakul, P. Phukpattaranont, The Usefulness of Wavelet Transform to Reduce Noise in the SEMG Signal, INTECH Open Access Publisher, 2012. [36] U. Sahin, F. Sahin, Pattern recognition with surface EMG signal based wavelet transformation, Syst. Man Cybern. (SMC), 2012 IEEE Int. Conf. (2012) 295–300. [37] F.A. Mahdavi, S.A. Ahmad, M.H. Marhaban, The utility of wavelet transform in surface electromyography feature extraction – a comparative study of different mother wavelets, Proc. World Acad. Sci. Eng. Technol. 7 (2) (2013) 7–12. [38] C.-H. Lin, Frequency-domain features for ECG beat discrimination using grey relational analysis-based classifier, Comput. Math. Appl. 55 (4) (2008) 680–690. [39] J.V. Basmajian, C.J. De Luca, Muscles Alive: Their Functions Revealed by Electromyography, Williams and Wilkins Press, Baltimore, USA, 1985. [40] H.J. Hermens, B. Freriks, C. Disselhorst-Klug, G. Rau, Development of recommendations for SEMG sensors and sensor placement procedures, J. Electromyogr. Kinesiol. 10 (5) (2000) 361–374. [41] R. Karthik, ECG Simulation Using Matlab, College of Engineering, Anna University, Chennai, 2007. [42] J. Taelman, B. Mijovic, S. Van Huffel, S. Devuyst, T. Dutoit, ECG artifact removal from surface EMG signals by combining empirical mode decomposition and independent component analysis, Biosignals (2011) 421–424.