A DENOISING METHOD OF DIAPHRAGM ELECTROMYOGRAM SIGNALS BASED ON DUAL-THRESHOLD FILTER

1. Introduction

In December 2019, COVID-19 had appeared in Wuhan and spread rapidly throughout the country within few months. COVID-19 is highly infectious, whose main routes of transmission are through respiratory droplets and close contact. The virus can lead to acute pneumonia, multiple organ failure, and other diseases, which seriously threaten the life of the patients.^1,2 Previous research has shown that diaphragmatic electromyography (EMGdi) signal is an effective indicator for breath monitoring, which can well reflect the inspiratory and expiratory times of respiratory center and the degree of respiratory effort.³ Moreover, the diaphragm is a vital auxiliary organ for human breathing, and the respiratory signals transmitted during its operation are not affected by lung disease and airway resistance. During treatment, EMGdi signals can be used to monitor the patient’s respiratory system. In addition, the EMGdi signals play a significant role in triggering the ventilator and adjusting it to cooperate with human breathing.^4,5

The frequency range of EMGdi signals is mainly from 0.5Hz to 250Hz, where the dominant power spectrum ranges from 125Hz to 150Hz.⁶ As additive interference, the frequency range of electrocardiogram (ECG) signals is from 0Hz to 75Hz. Therefore, the frequency contents of ECG and EMGdi signals overlap and it is difficult to remove ECG signals from EMGdi signals using the conventional filter banks. As reported in pioneering studies, the 30-Hz high-pass fourth-order Butterworth filter and frequency-domain filter were proposed.⁷ They eliminate information between 10Hz and 30Hz to filter the QRS complex out, but dampen the EMG signals’ amplitude in the process. Considering that the frequency contents of ECG and EMGdi signals overlap, Christov et al.⁸ and Kahl and Hofmann⁹ applied the template subtraction method to reduce the impact of cardiac activity by subtracting a template from the EMG signals. The template is made of the QRS complex by averaging multiple QRS complexes. However, this method will cause a partial loss of EMG signals. In addition, the processing methods in the time domain include mathematical morphology, which applies the idea of image processing to signals processing and removes the waveforms with ECG signal characteristics.^10,11 However, the processing methods in the time domain are not perfect. With the development of mathematical analysis, wavelet transform has become an advanced technique which combines time and frequency aspects to analyze continuous signal. With multiple resolutions, wavelet transform can describe the nonstationary characteristics of signals well, and has the advantages of flexible selection of basic functions.^12,13,14 Abbaspour et al.¹⁵ proposed a method based on the combination of adaptive neuro-fuzzy inference system (ANFIS) and wavelet transform to eliminate the ECG interference in the surface EMG signals. The low-frequency component containing part of the noise can be removed by combining the wavelet transform. Luo and Yang¹⁶ proposed a stationary wavelet transform algorithm to remove the ECG interference. They constructed the threshold by Donoho contraction algorithm to eliminate the wavelet coefficients of ECG. Finally, the denoised signals are reconstructed by the processed wavelet coefficients. Wu et al.¹³ and Taelman et al.¹⁷ combined wavelet transform with an independent component method to remove ECG signals and artifacts in myoelectricity. Besides, wavelet transform is widely used in the processing of nonlinear and nonstationary signals. Chauhan et al.^18,19,20 and van Leuteren et al.²¹ proposed evolutionary algorithm with adaptive wavelet mutation strategy to identify the bearing defect. Although the above-mentioned contributions were effective in normal ECG scenarios, the case of irregular ECG or weak EMGdi signals continues to be a bottleneck for the ECG cancellation from EMGdi signals.

In order to overcome this issue, we propose a denoising method for ECG cancellation from the EMGdi signals. Specifically, a dual-threshold filtering algorithm is applied, which can effectively remove the irregular ECG interference and better perform for the case of weak EMGdi signals. Based on experimental results using clinical data, it is shown that the proposed denoising scheme performs superior to the state-of-the-art methods.

2. Materials and Methods

2.1. Data acquisition

All EMGdi signals data used in this paper are obtained from pioneering clinical studies^22,23 by the Guangzhou Institute of Respiratory Diseases. The investigative protocol was approved by the institutional ethics committee of Chinese State Key Laboratory of Respiratory Disease, and informed consents were obtained from all the subjects who participated. The study was registered in ClinicalTrials.gov (NCT01782768) on April 29, 2013. The signals acquisition system mainly consists of five-lead esophageal electrodes and the PowerLab medical physiological signals acquisition system of ADInstruments, Australia, with a sampling frequency of 2000Hz. The signals have been preprocessed in the acquisition process as follows: 1000 times of magnification, 1-kHz low-pass filter, 15-Hz high-pass filter, and 50-Hz stopband. This preprocessing step filters out most of the high- and low-frequency noise interferences and 50-Hz power frequency interference.^24,25 We included the data of 10 patients for analysis, with an average age of ($67.5\pm 10.1$) years. Figure 1 shows the samples taken from three representative clinical EMGdi signals, which are interfered by different types of ECG signals. In Fig. 1(a), the EMGdi signals interfered by ECG can be obviously seen. We defined such a type of signals as $A$-type signals. Clinically, however, the EMGdi signals are weak in comparison with ECG, and patients could be suffering from cardiac disease, leading to the abnormal ECG interference. We plot two typical examples as shown in Figs. 1(b) and 1(c). As a comparison with that of $A$-type signals, we define such clinical signals as $B$-type signals.

2.2. Decomposition of ECG and EMGdi signals using wavelet transform

The continuous wavelet transform can be defined as

where the reasonable term ${\varphi }_{a,b}(t)$ is included to make the function analysis controllable in both the time domain and the frequency domain. The wavelet function is given by where $a$ is a scale factor used for restoration later, $a$ is multiplied by 2, and the window width is doubled; $b$ is a translation factor; and $\omega (t)$ is a wavelet basis function, and different basis functions can be selected according to different signals. The wavelet transform is different from the short-time Fourier transform (STFT). The STFT adds windows to the signals and performs the Fourier transform. The wavelet transform uses a finite-length attenuated wavelet basis as a window function for wavelet decomposition.

Figure 2 shows the decomposition of EMGdi signals from ECG interference using a five-scale Haar wavelet transform, where cd1, cd2, cd3, cd4, and cd5 are high-frequency components of different scales, and ca5 is the low-frequency component on the fifth scale. Since the sampling frequency of the experimental data is 2000Hz, the frequency range of ca5 is from 0Hz to 62.5Hz, which is similar to the main frequency band of the ECG signals. According to the wavelet transform, ECG interference can be removed from the EMGdi signals.

2.3. ECG positioning using adaptive wavelet threshold

From Fig. 2, although most of the ECG interference can be decomposed from EMGdi signals using the fifth-scale wavelet coefficients, i.e., ca5, some residual components of ECG signals are still mixed with the EMGdi signals, as shown in cd1. In order to remove the residual interference of ECG, we propose an adaptive wavelet-aided ECG positioning scheme. Specifically, the wavelet coefficients of ca5 are employed to detect ECG signals, and then, the position of detected ECG is used to remove the residual interference from EMGdi signals in cd1. In the QRS complex, the peak of the $R$-wave is the largest,²⁶ so the main task of QRS complex positioning is how to detect the position of the $R$-wave peak. In previous studies, QRS complex detection methods are mainly aimed at regular ECG signals. The period and $R$-wave amplitude of normal ECG signals are stable, and the detection is not difficult. However, abnormal ECG signals have problems such as unstable periods and large changes in $R$-wave peak, which bring great troubles in the removal of ECG interference from EMGdi signals. To solve the above problems, we propose the wavelet adaptive threshold method, which is the first “threshold” of the double-threshold filtering algorithm. The process is as follows:

(1)	The original signals were decomposed into five different dimensions with five-scale Haar wavelet. The five dimensions are represented by $\mathrm{\text{Ca}}(5,k)$, $\mathrm{\text{Cd}}(5,k)$, $\mathrm{\text{Cd}}(4,k)$, $\mathrm{\text{Cd}}(3,k)$, $\mathrm{\text{Cd}}(2,k)$, and $\mathrm{\text{Cd}}(1,k)$. $\mathrm{\text{Ca}}(j,k)$ and $\mathrm{\text{Cd}}(j,k)$ represent the low-frequency coefficients and high-frequency coefficients, respectively, where $j$ represents the scale of the decomposition, and $k$ is the position information of the wavelet coefficients.
(2)	The frequency range of the low-frequency coefficient $\mathrm{\text{Ca}}(5,k)$ is from 0Hz to 62.5Hz, on which the ECG signals are mainly concentrating. So the low-frequency coefficient $\mathrm{\text{Ca}}(5,k)$ is used to detect the QRS complex. Record the maximum value and position information between every two positive-to-negative mutation points on the wavelet coefficient. Since humans will have 1–1.6 heartbeat records per second,26 take the maximum value per second in the first five seconds as th(1), th(2), th(3), th(4), and th(5), and use the mean value of these five numbers as the initial threshold for the first five ECG positions. Starting from the sixth threshold, it is automatically updated by the unequal-weight weighted average method, and the expression can be written as $$\mathrm{\text{th}}(n)=\left\{\begin{array}{cc}\frac{1}{5}\sum _{x=1}^5\mathrm{\text{th}}(x),\hfill & 0<n\le 5,\hfill \\ k_1\mathrm{\text{th}}(n-2)+\frac{1}{5}{k}_2\sum _{x=n-6}^{n-1}\mathrm{\text{th}}(x),\hfill & n>5,\hfill \end{array}\right.$$(3) where $k_1$ and $k_2$ are the threshold coefficients, which meet $k_1+k_2=1$. They represent the proportion of the old threshold and the new threshold. The larger the $k_2$ is, the faster the old thresholds are forgotten. The position corresponding to the wavelet coefficients greater than the current threshold is the position of $R$-wave peak in ECG interference, and the position of $R$-wave peak in the QRS complex is denoted as $P_n$. Every time the algorithm detects the peak position of an $R$-wave, the threshold is updated immediately to detect the next peak position. Therefore, the algorithm can always effectively predict the change of the ECG signals and then adjust the new threshold so that ECG signals with different amplitudes can be detected. The automatic update of the threshold is completed before the current data is processed, so as to improve the running rate of the algorithm.
(3)	Use Peakindex to express the detected peak position information collection. The Peakindex can be expressed as follows : $$\mathrm{\text{Peakindex}}=\{P_1,P_2,P_3,\dots ,P_n\}.$$(4) In order to verify the effectiveness of the QRS detection algorithm, we consider the irregular ECG data from MIT-BIH Arrhythmia database27 (nos. 103 and 221) as examples to test the performance of the so-designed QRS complex detection. Figures 3(a) and 4(a) plot the employed ECG signals with 1-dB additive white noise. By using the so-designed filter and QRS complex detection, it can be seen from Figs. 3(b) and 4(b) that the $R$-wave peak can be effectively determined even in a low signal-to-noise ratio (SNR) condition (i.e., 1dB).

2.4. ECG elimination by average energy threshold

In previous studies, the processing of the interference interval is usually by truncating or shrinking the “larger” wavelet coefficients, and the wavelet coefficients less than or equal to the threshold are not processed. Zhan et al.²⁸ applied the “larger-coefficient-truncated” method to deal with interference, and the shrinkage functions are defined as

Luo and Yang¹⁶ also adopted this method. However, this method will damage part of the EMGdi signals, because these signals and the ECG signals linearly overlap in time domain. In this paper, we propose a method for processing interference intervals with average energy threshold. First, the average energy of interference neighborhood is used as the basis of interference interval threshold. Second, the wavelet energy coefficients greater than neighborhood average energy threshold are replaced by this threshold. Finally, the wavelet energy coefficients less than this threshold will not be processed. Since the neighborhood energy of each QRS complex interference changes, the average energy threshold is also constantly changing. The steps are as follows:

(1)	The wavelet coefficients decomposed by the five-scale Haar wavelet are squared to obtain the wavelet energy coefficients of the signals. The expression can be written as $${\mathrm{\text{PW}}}_x(j,k)=|C_x(j,k){|}^2,$$(6) where PW${}_x(j,k)$ represents the wavelet energy coefficients with the $j$th scale and sampling point $k$.
(2)	To reduce the loss of the EMGdi signals, we set a new threshold for the interference interval to smooth the interference interval. The duration of the QRS complex in healthy adults is between 0.07s and 0.10s.26 The $T$ value is selected as the interference duration in an appropriate range according to different patient conditions. According to the characteristic of the QRS complex, the range of ECG interference is determined in the ratio of 1:1 before and after the $R$-wave peak. The $T$ can be expressed as $$T=[P_n-0.5t\times f,P_n+0.5t\times f],$$(7) where $f$ is the frequency of different layers of the wavelet decomposition. Let $t_1=P_n-0.5t\times f$, $t_2=P_n+0.5t\times f$, where $t_1$ is the starting position of the interference interval, and $t_2$ is the ending position of the interference interval.
(3)	Take 0.06-s wavelet energy coefficients on both sides of the interference interval $T$ to find the average energy of the neighborhood of the interference interval. The average energies of the left and right neighborhoods of the $P_n$th interference interval can be expressed as $$A_r(P_n)=\sum _{{k=P}_n-t_1-0.06\times f}^{P_n-t_1-1}{\mathrm{\text{PW}}}_x(j,k)∕0.06\times f,$$(8) $$A_l(P_n)=\sum _{k=P_n+t_2+1}^{P_n+t_2+0.06\times f}{\mathrm{\text{PW}}}_x(j,k)∕0.06\times f.$$(9) The average energy threshold is obtained by adding the average energies of the left and right neighborhoods. The threshold was calculated using the equation given by $${\mathrm{\text{TH}}}_j(P_n)=A_r(P_n)+A_l(P_n).$$(10) After the threshold value is obtained, the interference interval $T$ is denoised according to the following formula : $${\mathrm{\text{newPW}}}_x(j,k)=\left\{\begin{array}{cc}{\mathrm{\text{PW}}}_x(j,k),\hfill & {\mathrm{\text{PW}}}_x(j,k)\le {\mathrm{\text{TH}}}_j(p_n),\hfill \\ {\mathrm{\text{TH}}}_j(p_n),\hfill & {\mathrm{\text{PW}}}_x(j,k)>{\mathrm{\text{TH}}}_j(p_n).\hfill \end{array}\right.$$(11) According to Eq. (11), the wavelet energy coefficients after denoising are obtained. Since the wavelet energy coefficients are obtained through the square value of the wavelet coefficients, it is necessary to find the square root and give the original positive and negative signs to obtain the denoised wavelet coefficients. After wavelet reconstruction, the denoised EMGdi signals are obtained.

3. Results and Discussion

3.1. Performance verification

We evaluated the use of our algorithm for the removal of ECG contamination from EMGdi recordings in the clinical signals. And its performance (ECG removal as well as EMGdi preservation) was compared to two previously reported algorithms, which are the wavelet-based adaptive filtering algorithm (WAF)²⁸ and stationary wavelet transform algorithm (SWT).¹⁶ The clinical EMGdi signals with ECG interference are shown in Figs. 5(a), 6(a), and 7(a), and the processing results of wavelet-based adaptive filtering algorithm, stationary wavelet transform algorithm, and dual-threshold filtering algorithm are shown in Figs. 5(b), 6(b), 7(b); 5(c), 6(c), 7(c); and 5(d), 6(d) and 7(d), respectively. Figure 5 shows that the three methods have good performance in removing the ECG interference. Remarkably, residuals of ECG peaks often remained and the EMG signals’ amplitude dampened after applying WAF and SWT. However, that is not conducive to the subsequent use of EMGdi signals.

3.2. Results analysis

For a clinical signal, we cannot accurately obtain its real signal, thus, SRN is not useful to evaluate the efficiency of algorithms. Power spectral density (PSD) is a commonly used method to describe EMG signals,⁶ which can compare the relative contributions of ECG and EMGdi. In order to analyze the ECG interference removal performance of our algorithm comprehensively, two commonly used characteristic values of the power spectrum are used in the experiments,^3,29 high frequency-to-low frequency ratio ($H∕L$) and median frequency ($F_m$). Here, $H∕L$ is the ratio of power in the range of 125–150-Hz frequency band to the power in the range of 25–50-Hz frequency band. The $H∕L$ can be expressed by

The dominant power of the ECG signals lies below 50Hz, whereas that of the EMGdi signals lies between 125Hz and 150Hz. Therefore, $H∕L$ can effectively reflect the relative contributions of ECG and EMGdi.

The median frequency is defined as

where $f$ is the frequency, and $P(f)$ is the power spectrum of the corresponding frequency. $F_m$ can reflect the changes of ECG and EMGdi contents of the signals before and after processing. To evaluate the performances of the three algorithms, we present the PSDs before and after processing the $A$-type signals and the $B$-type signals in different methods.

According to the experimental results shown in Table 1, after being processed by the three algorithms, respectively, the $H∕L$ ratios are higher than those of the original signals, indicating that most of the ECG signals in the range of 25–50-Hz frequency band have been removed. Especially, the $H∕L$ ratio of the signals processed by the dual-threshold filtering algorithm is significantly higher than those of the other two algorithms. These results indicate that the dual-threshold filtering algorithm performs better in removing ECG interference. A comparison is also shown in Table 2. It can be seen that, even when the EMGdi signals are weak and the ECG signals are abnormal, the algorithm in this paper still has better performance in removing ECG interference.

From the experimental results shown in Table 3, we can see that the median frequency of the signals increases significantly after processing by the three algorithms. That means the proportion of EMGdi content of the signals has gone up. Remarkably, the experimental results of the dual-threshold filtering algorithm show the most obvious increase in the median frequency of the original signals. We can also learn from Table 4 that the algorithm we proposed still has better performance in retaining the EMGdi signals than the other two algorithms, even when the EMGdi signals are weak and the ECG signals are unsteady and irregular.

4. Conclusion and Further Work

In order to remove ECG interference from EMGdi signals, this paper proposed a dual-threshold filtering algorithm, which performed robustly under the condition of irregular ECG interference or weak EMGdi signals. For validations, two types of clinical irregular ECG signals were considered in experiments. Experimental results showed that the proposed denoising scheme yielded a superior solution compared to the existing methods, especially for the practical scenarios of weak EMGdi signals.

However, this study has some limitations. Specifically, by taking advantage of the adaptive threshold filtering, the detection of QRS complex entails relatively high computations, thus failing to provide a real-time process for conventional central processing unit. In future work, light-weight QRS detection scheme will be studied.

Ethical Compliance

All EMGdi signals data used in this paper are obtained from the study “Comparative Effects of Different Noninvasive Ventilation Mode on Neural Respiratory Drive in Recovering AECOPD Patients”^22,23 at the Guangzhou Institute of Respiratory Diseases. The investigative protocol was approved by the institutional ethics committee of Chinese State Key Laboratory of Respiratory Disease, and informed consents were obtained from all the subjects who participated. The study was registered in ClinicalTrials.gov (NCT01782768) on April 29, 2013.

Conflict of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

This work was supported by Science and Technology Program of Guangzhou (Nos. 202002030353 and 2019050001), Science and Technology Planning Project of Guangdong Province (No. 2019A1515011940), and Industry-Academia-Research Innovation Project of Blue-Fire of Ministry of Education (No. CXZJHZ201803).