Electromyography Parameter Variations with Electrocardiography Noise

Electromyograms (EMG signals) may be contaminated by electrocardiographic (ECG) signals that cannot be easily separated with traditional filters, because both signals have some overlapping spectral components. Therefore, the first challenge encountered in signal processing is to extract the ECG noise from the EMG signal. In this study, the EMG, mixed with different degrees of noise (ECG), is simulated to investigate the variations of the EMG features. Simulated data were derived from the MIT-BIH Noise Stress Test (NSTD) Database. Two EMG and four ECG data were composed with four EMG/ECG SNR to 32 simulated signals. Following Pan-Tompkins R-peak detection, four ECG removal methods were used to remove ECG with different compensation algorithms to obtain the denoised EMG signal. A total of 13 time-domain and four frequency-domain EMG features were calculated from the denoised EMG. In addition, the similarity of denoised EMG features compared to clean EMG was also evaluated. Our results showed that with the ratio EMG/ECG SNR = 10 and 20, the ECG can be almost ignored, and the similarity of EMG features is close to 1. When EMG/ECG SNR = 1 and 2, there is a large variation of EMG features. The results of our simulation study would be beneficial for understanding the variations of EMG features upon the different EMG/ECG SNR.


Introduction
The electromyography (EMG) signal represents the muscle performance of central and peripheral motor control, especially EMG collected from back muscles, which can be used to evaluate pain around the shoulder, neck, and back. EMG can also be used as an indicator to evaluate ergonomic performance and trunk performance. Surface EMG (sEMG) is also used for body motion detection, such as hand motion detection [1], hand gesture recognition [2], and shoulder motion intention detection, with the aid of a pattern recognition process [3]. However, the collection of EMG signals is easily disturbed by many noises, which affect the interpretation of EMG results and the accuracy of diagnosis. The noise sources of EMG include motion noise, heartbeat, and respiration. The distance between the recording electrode of EMG and the heart may affect the strength of heartbeat noise on EMG, which increases the complexity of EMG signal processing. Therefore, noise processing techniques, such as high-pass filters, are required for further EMG signal processing. The EMG signal is a time series with constantly changing frequency and amplitude with no specific waveform, therefore many time domains, frequency domains, and nonlinear parameters were developed to analyze the EMG signal characteristics [4][5][6][7]. Different EMG parameters have their prominent effects under different EMG measurement environments. Time-domain characteristic parameters are usually mainly the average value, maximum value, summation, or variation of amplitude. This time domain is inclusive of the average amplitude absolute value (AAV), mean absolute value (MAV1 type-1), average amplitude change (AAC), the maximum value of the amplitude contains the amplitude of the first burst (AFB), the sum of the amplitudes contains integrated EMG (IEMG), simple square integral (SSI), and waveform length (WL). The amplitude variation parameters include standard deviation (STD), root mean square (RMS), log detector (LOG), median differential value (MDV), and difference absolute standard deviation value (DASDV). Zero crossing (ZC) computes the signal period in the time domain. Frequency-domain parameters are usually converted from EMG signal to frequency domain through Fourier Transform, and then its total frequency power (TTP), median frequency (MDF), max peak frequency (PKF), the amplitude of peak frequency (PKF.amp), and so on, are calculated. In addition to Fourier transform, there are other transformation functions, such as Wavelet transform [8], Hilbert-Huang Transform, etc. [9]. In addition to time-domain and frequency-domain parameters, nonlinear parameters also play an important role in EMG's characteristics [10].
The most common EMG noise is electrocardiography (ECG). There have been many studies trying to separate the ECG signal from the EMG signal, such as the ICA (independent component analysis) separation technique [11]. Abbaspour et al. use the filter combination technique to detect QRS and construct an ECG signal template; drift noise is removed by a low-pass filter, and denoised EMG is obtained by the subtraction of the ECG template. Their method has achieved impressive results [12]. Christov et al. also used a similar filter combination concept to separate ECG and EMG by dynamic filtration [13,14]. The concept of template subtraction is also used in many studies [15,16]. Combining template subtraction and empirical mode decomposition (EMD) is a recent and very innovative approach [17]. Another approach to separate ECG from EMG is to treat ECG as signal and EMG as noise and remove ECG with R-peak detection. Many interesting algorithms have been proposed and the R-peak detection performance is usually examined with public databases, such as MIT-BIH arrhythmia and Fantasia databases. The moving average (MA) filtering method proposed by Modak et al. has achieved 99.82% sensitivity, 99.88% positive predictivity, and 0.31% detection error rates on R-peak detection performance with the MIT-BIH arrhythmia database [18]. The Zhang et al. algorithm achieved a sensitivity of 99.70% and positive predictivity of 99.93% on R-peak detection performance with R-peak threshold [19]. For noise signal processing, MIT-BIH provided a set of noise-stress test databases in 2000 (noise-stress test database, NSTD), with various SNR ECGs with EMG noise, as well as clean EMG and other noise signals [20]. Other methods, such as event-related moving averages [21], and Ensemble Empirical Mode Decomposition achieved great R-peak detection performance [22]. Table 1 organized various ECG-detection algorithms proposed in recent years, applied under NSTD, and their corresponding ECG R-peak detection performance by sensitivity (Se%) and positive predictivity (+P%). From the 1985 Pan-Tompkins algorithm to Rahul et al. in 2021, sensitivity and positive predictivity increased from (74.46%, 93.67%) to (97.58%, 96.04%); The better the R-peak detection performance, the more complex the algorithm is.
Since EMG is an evaluation tool commonly used in many clinical examinations or sports assessment, such as upper back pain and posture adjustment, the EMG electrodes are placed very close to the upper back, for example, EMG is recorded on specific muscles such as the Trapezius muscle. The artifact from the heartbeats is inevitable during activity. It is necessary to extract cardiac artifacts from it and obtain a clean EMG signal. Therefore, the aim of this study was to explore the variations of common time-domain and frequencydomain EMG parameters under different EMG/ECG SNR with conventional Tompkins R-peak detective and removing algorithm.

Simulated EMG and ECG Dataset
Original data are derived from the MIT-BIH Noise Stress Test (NSTD) Database [30]. The NSTD database consists of 30 min ECG, EMG, and electrode motion artefact recording, respectively. In this study, the first minute of two ECG records (118e00 and 119e00) are used. A0 is derived from 118e00, and B0 is derived from 119e00. Clean ECG is derived from two filters. The first step is to remove the DC component. The second step is to pass through a FIR low pass filter (the order = 36, normalized cutoff frequency = 0.1, hamming window). The filter coefficients are implemented by fir1 function in a signal package with R languages. The two-channel ECG signals of A0 form A02C and A03C; the same process of B02C and B03C is derived from B0. EMG is taken from the ma signals in the NSTD database, which are EMG1, and EMG2; the mixed signals are combined according to the four EMG/ECG intensity ratios: 1, 2/1, 10/1 and 20/1. The abbreviations are Dm.1, Dm.2, Dm.10 and Dm.20, with corresponding SNR (signal-to-noise ratio) are 0 dB (1:1), 6 dB (2:1), 20 dB (10:1) and 26 dB (20:1). SNR is calculated using Equation (1). Table 2 is the list of the mixed signals. Figure 1 is an illustration of some mixed EMG signals. The original signal file and subsequent EMG feature estimation code are provided in Appendix A. SNR = 20 * log 10 EMG ECG (1)

ECG R-Peak Detection
ECG R-peak detection is based on the Pan-Tompkins algorithm [31], and implemented by the R rsleep package, with detect rpeaks code to estimate R-peak location [32], and then calculates the median and standard derivation of R-peak time series. The Pan-Tompkins method consists of a band-pass filter, following with derivatives, squaring, and moving window integration. The R-peak position is determined with a threshold. For the clean ECG dataset, A02C, A03C, B02C, and B03C, the numbers of detected R peaks are used as the reference to compare the calculated R peaks of the corresponding mixed EMG signals. Due to the different SNR, some R peaks will be missing, or some EMG bursts will be misjudged as R peaks, so the following indices are used to evaluate the R-peak detection ability.
True positive (TP): the number of QRS correctly identified. ( False-positive (FP): peaks that do not correspond to the QRS.
False-negative (FN): the number of QRS incorrectly identified.

ECG Deletion
After detecting the R-peak locations in the EMG/ECG simulated signal, the following four methods are used to remove the components of the ECG signal to obtain a denoised EMG signal. Figure 2 depicts the four ECG deletion methods implemented in this study.
F0 Method: No signal processing at all. F1 method: The detected R peak is centered on the window w2, which is covered by a 0.2 sec window width. The signal within window 2 is set as 0.
F2 method: The detected R peak is centered on the window w2, then the signal within the w2 is replaced as the average of w1 and w3. The duration of w1, w2, and w3 are 0.2 s.
F3 method: After EMD decomposes the mixed EMG signals, it takes the first two layers of IMF (intrinsic mode function), IMF1, and IMF2 signals, and recombines them. It is like a high-pass filtering approach to remove ECG contents.

EMG Parameters
There are thirteen time-domain and four frequency-domain EMG features investigated in the study. The list of EMG features and the corresponding formula are as shown in Table 3. In order to understand the degree of EMG parameter variation, the similarity index (SI) is defined as below: SI = EMG feature derived from simulated mixed EMG/EMG feature derived from raw clean EMG.

Statistics
Mean and standard deviation of the following time series were calculated: the numbers of R peaks were detected in the simulated signal, time-domain and frequency-domain of EMG features, and the ratio of its value corresponding to the clean EMG. The paired t-test was used to examine the difference of the EMG features among the four ECG deletion processing methods (F0-F1, F0-F2, F0-F3, F1-F2, F1-F3, and F2-F3). The p-value was set to 0.05. The experiment flowchart of this study is shown in Figure 3.

R-Peaks Detection Performance
The R-peaks detection performance of four simulated signal is listed in Table 4. For ECG A02C, the true R-peak number is 73. The numbers of R peaks detected from the simulated EMG: D1. 1

Similarity for All Features in Different EMG/ECG SNR
Regardless of the ECG deletion methods, the time-domain and frequency-domain EMG features on 16 simulated EMG signals are averaged. The selected features are listed in Table 5. Table 5 shows that, compared with the true EMG (EMG1), the order of the time-domain features are Dm.1 > Dm.2 > Dm.10 > Dm.20 ≥ true EMG. To compare the similarity of EMG features between simulated EMG and true EMG, the ratio of EMG feature to the true EMG is defined as a similarity index (SI). If this ratio is close to 1, it means the EMG feature of the simulated mixed signal is exactly the same as that of the raw true EMG. The mean and standard deviation of the ratios of all time-domain parameters are calculated and listed in Table 6. With a single-sample t-test to test whether there is a difference with ratio = 1, there is no statistical difference with ratio = 1 for F2 at Dm.10, and the mean and standard deviation of the ratio are 0.997 and 0.035, respectively. This means that by using F2 ECG deletion method under SNR = 20dB (Dm.10) condition, the EMG features have the highest similarity with features calculated from true EMG. To understand the influence of the four ECG deletion methods on the EMG features, the paired t-test was used to test the SI difference among F0, F1, F2, and F3 on Dm.1 to Dm.20. The results are shown in Table 7. The findings reveal that F0 has significant differences from F1, F2, and F3 in Dm.1, Dm.2, Dm.10, and Dm.20, and the SI of F0 is also higher than that of F1, F2 and F3, which means that there is the highest derivation on F0 than the other three methods, since there is no removal of ECG components. The F1 method is to delete the signal value within the window length centered at the R peak. F2 is to replace the signal within the window with average EMG around the window. There are significant differences among F0, F1, and F2. F3 is a high-pass filter, and the result of F3 is close to F1. There is no statistical difference between F1 and F3 in Dm.1, Dm.10, and Dm.20. From Table 6, the SI from Dm.1 to Dm.2 are 1.986 (4.927) to 1.299 (1.164) for F1, but it is 1.287 (0.652) to 0.669 (0.373) for F3. The SI for F3 is almost always <1, after Dm.1. Therefore, the disadvantage of the F3 ECG deletion method is that it will greatly suppress the original EMG signal parameters' performance. In Table 6, the first row of ALL signal is used to demonstrate the SI from Dm.1 to Dm.20. Paired t-tests among Dm.1 to Dm.20 of ALL signals are all statistical differences. In other words, different EMG/ECG ratios do affect the degree of variation in EMG features.  Table 8 contains the details of the SI distribution on each time-domain feature for F0, F1, F2, and F3. Combined with Table 6, the variations of SI on each EMG feature among F0 to F3 with different SNR were clearly shown. As shown in Table 8A In other words, these two time-domain EMG features are particularly affected by ECG noise interference. In addition, the F3 method will filter out some EMG components, the SI are <1, as shown in Table 8D.   Table 9 lists the average PKF value among four SNR conditions. PKF is interfered with by ECG noise, and the PKF ratio even exceeds 100. Therefore, PKF is a feature with great variation under noise. Figure 4 directly lists ratios for three frequency-domain features. Unlike the results of time-domain features, there is no consistent rule for frequency-domain features on the SI which vary greatly on each frequency-domain feature. The performance of TTP is similar to traditional time-domain parameters. For MDF, the SI is close to 0.87, and is similar on Dm.2, Dm.10, and Dm.20, regardless of the four ECG removing methods. The SI is almost = 1 for Dm.1, which is the simulated signal with the highest ECG noise. This shows that MDF is very robust for ECG noise.

Discussion
This study discusses the change in EMG parameters in the presence of mixed heartbeat noise. There were several studies on heartbeat signal removal from EMG, such as the subtraction method [33], adaptive filters [34], independent component analysis (ICA), singular spectrum analysis and segmented-beat modulation method [35]. For example, Barrios used the singular spectrum analysis (SSA) to remove ECG without using the ECG reference signal and it outperformed the other methods [36]. Thandar used wavelet transform to compare 9 SNR of EMG/ECG (−20 to 20 dB, in increments of 5 dB) [37], and removed ECG noise by discrete wavelet transform. Thandar's research is similar to this study. Their results found that discrete wavelet transform is better than high-pass filtering in terms of relative error, correlation coefficient, and absolute mean error. Different thresholds are required to remove ECG noise with different SNRs. It is not practical to apply to real-time systems since it takes too much time to find out optimal thresholds. Therefore, they suggest that future research can develop a method for estimating SNR to speed up the noise-removing signal processing procedure. Abbaspour et al. combined wavelet and artificial neural networks to remove ECG from EMG signals [38]. According to the above studies, if the SNR of EMG is known in advance, it can enhance the estimation algorithm performance. Although this study only compares four common ECG removal methods, with the traditional Pan-Tompkins R-peak detection algorithm. The major contribution is to understand the derivation degrees of the extracted EMG parameters. Subsequent research can develop a system for predicting ECG noise, followed by an optimal ECG filtering algorithm, and it is believed that better EMG estimation results can be obtained. In this study, Pan-Tompkins algorithm was used to detect the R peak. According to the study [39], the R-peak performance (Se%, +P%) with the Pan-Tompkins algorithm on the NSTD Database, with the noisiest two sets of signals of 118e_6 and 119e_6, were (71.47%, 60.50%) and (66.48%, 54.21%), respectively. It is very similar to the Dm.10 results demonstrated in Table 4. This means that the R-peak detection result in this study is close to [39] and is acceptable. The main contribution of this paper is to discuss the EMG parameters' changes under EMG/ECG SNR. For different clinical EMG parameters, the results of this paper can be used to estimate whether clinical EMG features are derived from real EMG or from ECG noise. In addition, although it can be used to estimate EMG/ECG SNR, it is not the focus of this article, but this article gives a good direction for a follow-up study to estimate the EMG/ECG SNR in time, with a further rigorous examination.

Conclusions
This study uses simulated EMG and ECG to demonstrate the changes in the timedomain or frequency-domain EMG features under the different EMG/ECG SNR. When the ECG signal strength is relatively lower than EMG, it is only necessary to (1) detect and delete the ECG, and then (2) compensate with the average value of the left and right windows of EMGs, and then (3) the EMG features close to the original EMG can be obtained; different EMG features have different degrees of change after noise interference. Since different EMG features are suitable for different EMG measurement environments, this study can provide a reference for various signal processing of EMG users.

Conflicts of Interest:
The authors declare no conflict of interest.