Application of Stockwell Transform and Shannon Energy for Pace Pulses Detection in a Single-Lead ECG Corrupted by EMG Artifacts

: Electrocardiogram (ECG) analysis is important for the detection of pace pulse artifacts, since their existence indicates the presence of a pacemaker. ECG gives information on the proper functionality of the device and could help to evaluate the reaction of the heart. Beyond the challenges related to the diversity of ECG arrhythmias and pace pulses, the existence of electromyogram (EMG) noise could cause serious problems for the correct detection of pace pulses. This study reveals the potential of a methodology based on Stockwell transformation (S-transform), subsequent Shannon energy calculation and a threshold-based rule for pace artifact detection in a single-lead ECG corrupted with EMG noise. The design, validation and test are performed on a large, publicly available artiﬁcial database acquired with high amplitude and time resolution. It includes various combinations of ECG arrhythmias and pace pulses with di ﬀ erent amplitudes, rising edges and total pulse durations, as well as timing that corresponds to di ﬀ erent pacemaker modes. The training was done over 312 (ECG + EMG) signals. The method was validated on 390 clean ECGs and independently tested on 312 (ECG + EMG) and 390 clean ECGs. The achieved accuracy over the test dataset was Se = 100%, PPV = 98.0% for ECG corrupted by EMG artifacts and Se = 99.9%, PPV = 98.3% for clean ECG signals. This shows that, despite EMG artifacts, the S-transform could distinctly localize the pace pulse positions and, together with the applied ShE, could provide precise pace pulses detection in the time domain.


Introduction
Pacemakers are devices that deliver electrical impulses to the heart chambers in order to provoke regular cardiac contractions. Electrocardiogram (ECG) analysis is important for the detection of pace pulse artifacts, since their existence indicates the presence of a pacemaker. ECG gives information on the proper functionality of the device and could help to evaluate the reaction of the heart. Statistics on the occurrence of pace pulses can be used to measure the evolution of arrhythmias and can eventually fine-tune the electrical or pharmacological therapy. However, not all pacemakers are equipped with event recorders, and some of them provide only concise statistics about pacing, but not information about the timing of the pace pulses [1].
One more reason for the need to detect pace pulses is the necessity to reject artifacts before automatic ECG analysis. As recently reported, the presence of pacing artifacts may impede the diagnosis of ventricular fibrillation during cardiac arrest [2].
In the best scenario, the ECG of a patient with an implanted pacemaker consists of the following: • natural ECG signal with amplitude up to 2 mV and frequency range between 0 Hz and 150 Hz; • pace pulses with a wide amplitude range from 100 µV to as high as 100 mV, pulse width from 100 µs to 2 ms, measured on the patient's skin, and a fast rising edge that could be as narrow as 10 µs [3].
Depending on the underlying heart disease, the pacemaker could be programmed to pace the upper chambers, the lower chambers or both at a different rate.
According to the number of active leads, the pacemaker types include the following: • single-chamber pacemaker, which uses one electrode placed either in the right atrium (when the sinus node is generating pulses that are too slow or irregular) or in the right ventricle (when the electrical flow is slowed or blocked in the region of the AV node, and the normal impulses from the atria cannot reach the ventricle); • dual-chamber pacemaker, which uses one electrode in the upper chamber and one electrode in the lower chamber of the heart, thus regulating the heart's electrical activity in both chambers; • bi-ventricular pacemaker, which uses three leads, first into the right atrium, second into the right ventricle and third into the left ventricle, applied for cardiac resynchronization therapy.
According to programming, the pacemakers could be as follows: • fixed-rate pacemaker, which produces pace pulses at a steady rate, regardless of the heart's own electrical activity; • pacemaker "on demand", which monitors the heart rhythm and generates electrical pulses only if the heart is beating too slowly or if it misses a beat; • rate-responsive pacemaker, which speeds up or slows down the heart rate depending on how active the person is.
Depending on the pacemaker type, pace artifacts could appear at different places in the PQRST segment, e.g., preceding the P wave or the QRS complex for single-chamber pacemakers, before both P wave and QRS complex in the case of dual-chamber devices, and two close in time artifacts preceding the QRS generated by the bi-ventricular pacemaker. According to the device programming, the pulses could be fixed-rate, on demand or rate-responsive [4].
In real environments, the ECG of a patient usually contains noise components such as: • low-frequency noise (0.05-1) Hz, which is generated by the patient's movements and respiration; • power-line interference (50 or 60 Hz); • electromyographic (EMG) noise caused by muscle electrical activity, which is reported to reach up to 1000 Hz [5].
There are two general approaches for pace pulses detection: hardware and software; these are outlined as follows: • Hardware approach-analyzes the analogue ECG signal in the process of ECG acquisition. This type of pace pulses detection is implemented by Texas Instruments relying on attenuation of the QRS complexes and consequent steep slope detection [6]. • Software approach-analyzes the digital ECG signal, which should be digitized at high sampling rates that preserve the frequency content of the pacing pulses [7][8][9][10][11]. They involve: • band-pass filtering that has been specially designed to distinct the pace pulses from the real ECG components and low-frequency noises [12,13]; • search for steep edges in a single lead [14,15] or multi-lead ECG [3]; • data transformations such as the first derivative of the ECG, power envelopes in the spectral range of (500 ÷ 2200) Hz or moving statistical window followed by peak selection based on the statistical properties of the transformed signal [16]; • application of a vectorcardiogram for precise detection of biventricular pacing [17].
Recently, a hybrid approach for pace pulses detection operating at low sampling frequency has been proposed [18]. The method utilizes simple analog processing to extract the pace pulses energy via filters tuned to capture spectral contents outside ECG frequencies of interest, followed by straightforward digital pace pulse detection at a low sampling rate.
The requirements of different standards, related to the parameters of the pace pulses that should be detected, are not very strong. ANSI/AAMI EC11 [19] have set the following ranges: (0.1 ÷ 2) ms for pulses duration; (2 ÷ 250) mV for amplitude range, up to 100 impulses per minute pacing rate and less than 100 µs rising edge duration. According to the IEC60601-2-27 standard [20], the lower limit of the pulse duration range is 0.5 ms, while the upper limit of the amplitude range is as high as 700 mV. However, modern pacemakers could generate pace pulses falling below these limits [4]. The different pace pulse morphologies combined with the changes in the paced patient population presenting with a variety of arrhythmias require a significant increase in the complexity of pace pulse detection algorithms [21].
The variety of ECG rhythms combined with pace pulses that differ in amplitude, duration and appearance in the ECG is not the only challenge for the pace pulse detection algorithm. Additional difficulties could be posed by the presence of EMG noise, which overlaps with pace pulses and could prevent their correct detection. On the one hand, low-amplitude pace pulses that could hardly be distinguished from environmental noise reduce the sensitivity (Se) and could lead to erroneous alarms for pacemaker damage. On the other hand, false-positive detections that decrease the positive predictive value (PPV) could lead to inaccurate interpretation of pacemaker rhythms by automated algorithms. In our previous study [22] we reported a drop in the positive predictive value between 10% and 22% when the algorithm described in [14] was applied on ECG signals with EMG artifacts.
Recently, a number of methods for EMG suppression in ECG recordings have been proposed, based on locally adaptive, low-pass Myriad filters [23], dynamic procedure for separation of ECG and EMG [24,25], adaptive periodic segment matrix and singular value decomposition [26], electrocardiogram reconstruction by a Segmented-Beat Modulation method [27], etc. Despite that all these methods address the preservation of QRS complexes and all other components of the diagnostic ECG, a technology that does not influence the steep pace pulses outside of the frequency range of the diagnostic ECG has yet to be announced.
In the last decade, Shannon energy and the Stockwell transform (called S-transform, [28]), which analyzes the ECG in the time-frequency domain, were applied separately [29,30] or in combination [31] for heartbeat detection.
The aims of this study are as follows: • To assess the potential of the Stockwell transform to highlight pace pulses on the background of intensive EMG artifacts in the time-frequency domain and to involve it in a brand new software approach for pace pulses detection.

•
To evaluate the ability of the elaborated software methodology to provide high accuracy for clean ECG signals and ECGs corrupted with EMG noise on an independent test dataset that has not been involved in the training. For that purpose, the digital ECG analysis is applied on a high-resolution, publicly available ECG database with preserved frequency content [32] that contains various arrhythmia types combined with EMG noise and pace pulse artifacts.

•
To compare the performance of the presented software approach for pace pulses detection to the performances of a time domain algorithm previously developed by us, as well as to a recently published hybrid approach for pace pulses detection.

ECG Database
The ECG signals used for training and testing are taken from the publicly available 'ArtPacedECGdb' database, which is described in detail in [32]. It includes ECG recordings in lead II with 10 s duration that represent different ECG rhythms generated by a Heidelberger Praxisklinik (HKP) simulator, including: All these rhythms are combined with artificially superimposed pace pulses that cover a wide range of rising edge (<10 µs ÷ 100 µs), total pulse durations (100 µs ÷ 2 ms) and amplitudes (100 µV ÷ 3 mV), and they correspond to various pacemaker modes and programming (i.e., atrial, ventricular pacing or both; fixed rate or on demand; bi-ventricular pacing).
The database involves a total of 1404 recordings with annotated positions of the pace pulses, and it has been previously applied for training and test purposes in [18,22]. The total number of clean ECGs with pace pulses is 780, while the paced ECGs that contain EMG noise are 624. The signals are recorded with 9.81 µV/LSB amplitude resolution at 128 kHz sampling rate, which preserves the steep raising and trailing edges of the pace pulses and the high-frequency waves of the EMG artifacts. The signals are down-sampled to 16 kHz in order to develop and test a method for pace pulses detection that operates with ECG signals acquired at a reasonable sampling rate.
In this study we divided the 'ArtPacedECGdb' database into three non-overlapping subsets: • training subset including 312 ECG signals corrupted by tremor (further referred as ECG + EMG) for development of the method; • validation subset including 390 clean ECG signals intended to verify the correct operation of the method on clean ECGs with pace pulses. This subset has been previously used for the training of the method described in [14], which is used here for comparison; • test subset including 312 (ECG + EMG) and 390 clean ECGs for independent testing of the designed method.

Methods
A block diagram of the proposed method for pace pulses detection is presented in Figure 1. The designed method operates over 10 s ECG buffer, which is converted in the time-frequency domain via application of the S-transform for localization of the pace pulses positions. The task is "returned" in the time domain via subsequent calculation of the Shannon energy (ShE) over the time-frequency representation of the 10 s ECG buffer. The calculated ShE is subjected to a threshold rule for detection of the pace pulse positions. A detailed description of the applied processing procedures and threshold-based rule is presented below.

S-Transform
S-transform, introduced by Stockwell et al. [28], is a time-frequency analysis method, which is deduced from short-time Fourier transform and continuous wavelet transform. S-transform is given by the equation where ECGT(t) is the ECG signal in the time domain, τ denotes the time of spectral localization and f is the frequency in the Fourier transform. The Fourier transform of ECGT(t) is The S-transform could be represented as operations on the Fourier spectrum ECGF(f) of the time signal ECGT(t): where f ≠ 0 and α is a frequency variable.

S-Transform
S-transform, introduced by Stockwell et al. [28], is a time-frequency analysis method, which is deduced from short-time Fourier transform and continuous wavelet transform. S-transform is given by the equation where ECG T (t) is the ECG signal in the time domain, τ denotes the time of spectral localization and f is the frequency in the Fourier transform. The Fourier transform of ECG T (t) is The S-transform could be represented as operations on the Fourier spectrum ECG F (f) of the time signal ECG T (t): where f 0 and α is a frequency variable. The S-transform of a discrete time series, for which f = n/NT, τ = jT and α = m/NT, is given by The output of S-transform is a matrix of complex values with rows indicating the frequency of the analyzed ECG signal and columns indicating the time.

Shannon Energy
In order to detect the pace artifacts highlighted in the time-frequency plane via the S-transform, we calculated the Shannon energy over S in the frequency range of (1000-2000) Hz, where we expect negligible influence of the EMG components. The result is the energy of the local spectrum obtained with S-transform for each sample j of the original ECG. ShE of each row of the extracted S-matrix in the selected frequency band is where n0 and n1 correspond, respectively, to f 0 = 1000 Hz and f 1 = 2000 Hz.

Threshold-Based Rule for Pace Artifact Detection
Considering the negligible influence of the EMG and ECG components in the frequency range for which ShE is calculated (i.e., 1000-2000 Hz) a threshold-based rule for detection of the pace pulses positions in the ShE signal should be applicable on the one hand and easy to adjust on the other hand. In order to prevent the influence of steep artifacts with unknown origin on the detection accuracy, we decided to design a threshold-based rule (presented via Equation (6)) involving the mean value of the rectified ShE signal.
abs(ShE( j)) > K * mean(abs(ShE)), The value of the coefficient K is empirically adjusted as follows: • The pace pulses detection accuracy (defined in the next Section 3.4) is estimated for K values in the range (2, 4, 6, . . . 16, 18) over the training dataset (312 ECG+EMG recordings), in order to highlight K values that provide high accuracy on ECG+EMG signals.

•
For the highlighted K values, the pace pulses detection accuracy is also estimated over the validation dataset (390 clean ECGs), in order to select a single K value that simultaneously provides high accuracy for ECG+EMG and clean ECG signals.

Accuracy Assessment
The potential of the designed threshold-based rule to correctly detect pace artifacts is assessed via sensitivity (Se) and positive predictive value (PPV), calculated as follows: where TP stands for the true positive detections, i.e., the number of pace pulses that are correctly detected; FN stands for the false negatives, i.e., the number of pace pulses that are not detected; FP stands for the false positives, i.e., the number of erroneously detected pace pulses. The influence of the signal-to-noise ratio (SNR) is also estimated considering the pace pulses amplitudes as the signal amplitude and the mean value of the absolute EMG amplitude in zero lines as the noise amplitude.

Results
The potential of the S-transform to highlight pace artifacts in an ECG recording corrupted with EMG is illustrated in Figure 2. Although the amplitude of the pace artifacts is comparable to the EMG amplitude, their positions are distinctly localized in Figure 2b. Figure 2c shows the calculated Shannon energy, which is suitable for subsequent analysis via the designed simple threshold-based rule.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 14 EMG amplitude, their positions are distinctly localized in Figure 2b. Figure 2c shows the calculated Shannon energy, which is suitable for subsequent analysis via the designed simple threshold-based rule. The results achieved during the adjustment of the amplitude threshold K are presented via the Se and PPV curves in Figure 3. Three values of K (8, 10, 12) lead to Se > 99% and PPV > 99% over the training dataset. The analysis over the validation dataset highlighted K = 10 as the value that assures the most balanced accuracy (see Table 1).   The results achieved during the adjustment of the amplitude threshold K are presented via the Se and PPV curves in Figure 3. Three values of K (8, 10, 12) lead to Se > 99% and PPV > 99% over the training dataset. The analysis over the validation dataset highlighted K = 10 as the value that assures the most balanced accuracy (see Table 1).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 14 EMG amplitude, their positions are distinctly localized in Figure 2b. Figure 2c shows the calculated Shannon energy, which is suitable for subsequent analysis via the designed simple threshold-based rule. The results achieved during the adjustment of the amplitude threshold K are presented via the Se and PPV curves in Figure 3. Three values of K (8, 10, 12) lead to Se > 99% and PPV > 99% over the training dataset. The analysis over the validation dataset highlighted K = 10 as the value that assures the most balanced accuracy (see Table 1).    The results achieved via the designed threshold-based rule with K = 10 over the test dataset (for both clean ECG and ECG+EMG) are presented in Table 2, once in total for all analyzed ECG recordings and a second time as mean ± standard deviation, median value and min-max range. Additionally, the Se and PPV values over the entire test dataset are illustrated in Figure 4. The results achieved via the designed threshold-based rule with K = 10 over the test dataset (for both clean ECG and ECG+EMG) are presented in Table 2, once in total for all analyzed ECG recordings and a second time as mean ± standard deviation, median value and min-max range. Additionally, the Se and PPV values over the entire test dataset are illustrated in Figure 4.   The influence of SNR on the accuracy is illustrated in Figure 5 in terms of PPV. As could be seen from Figure 4a, the SNR does not affect the Se, which is 100% (i.e., all pace pulses in the ECG+EMG test dataset are detected, despite the SNR).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 14 The influence of SNR on the accuracy is illustrated in Figure 5 in terms of PPV. As could be seen from Figure 4a, the SNR does not affect the Se, which is 100% (i.e., all pace pulses in the ECG+EMG test dataset are detected, despite the SNR). Examples of ECG+EMG signals with correctly detected and not detected pace artifacts, when K is set to 12, and erroneously detected pace artifacts, when K is set to 10, are presented in Figures 6−8, respectively. Examples of ECG+EMG signals with correctly detected and not detected pace artifacts, when K is set to 12, and erroneously detected pace artifacts, when K is set to 10, are presented in Figures 6-8 The influence of SNR on the accuracy is illustrated in Figure 5 in terms of PPV. As could be seen from Figure 4a, the SNR does not affect the Se, which is 100% (i.e., all pace pulses in the ECG+EMG test dataset are detected, despite the SNR). Examples of ECG+EMG signals with correctly detected and not detected pace artifacts, when K is set to 12, and erroneously detected pace artifacts, when K is set to 10, are presented in Figures 6−8, respectively. . The pace pulses are annotated with red asterisks '*'. The amplitude of ShE is limited due to the low amplitude of the pace pulse combined with the not very steep slope. Seventeen from the total number of 18 pace artifacts in the 10 s ECG segment (1st subplot) are slightly above the threshold and are correctly identified. The pace artifact marked with a red ellipse is slightly below the threshold of K*mean (abs(ShE)), when K = 12. The second subplot represents a detailed view of the low-amplitude pace artifact, which is not detected.  Figure 7. Example of false negative (FN) error just after the 11x10 4 sample. ECG (blue trace), the calculated ShE (black trace), the ShE threshold (red line). The pace pulses are annotated with red asterisks '*'. The amplitude of ShE is limited due to the low amplitude of the pace pulse combined with the not very steep slope. Seventeen from the total number of 18 pace artifacts in the 10 s ECG segment (1st subplot) are slightly above the threshold and are correctly identified. The pace artifact marked with a red ellipse is slightly below the threshold of K*mean (abs(ShE)), when K = 12. The second subplot represents a detailed view of the low-amplitude pace artifact, which is not detected.  . The pace pulses are annotated with red asterisks '*' on the 1st subplot. The false detection is due to the sharp EMG wave, which results in ShE exceeding the threshold of K*mean (abs(ShE)), when K = 10. The first subplot is identical to Figure 3. The second plot represents a detailed view of the EMG wave, which is detected as a pace artifact.

Discussion
Generally, a limited number of studies has addressed the problem for pace pulses detection, which might be due to the fact that there is no publicly available ECG database containing recordings from paced patients. This study reveals the potential of a methodology based on Stockwell transformation and subsequent Shannon energy calculation for pace artifacts detection in a single-lead ECG corrupted with EMG noise. The design, validation and test are performed on a large, publicly available artificial database [32] recorded at a high sampling rate of 128 kHz, thus preserving even the steepest pulses with rising edges about 10 µs. The benefits of using the 'ArtPacedECGdb' database, especially for training and validation, are as follows: • 'ArtPacedECGdb' database contains a great variety of ECG rhythms combined with pacing pulses that cover a wide range of amplitudes, rising edges and total pulse durations. The pace pulse positions are consistent with single-chamber, dual-chamber and bi-ventricular pacemaker types, as well as fixed-rate and 'on-demand' pacing, thus mimicking a variety of pacemaker types and modes. • 'ArtPacedECGdb' database provides paced ECGs corrupted by EMG artifacts, which is the only physiological noise that could compromise the methods for pace pulses detection. • 'ArtPacedECGdb' database is accompanied by reliable annotation of the positions of the pace pulses.
The benefit of using the ArtPacedECGdb database for test purposes is its public accessibility, which gives the opportunity for fair comparison.
The proposed methodology for pace pulses detection is further advanced by ECG analysis in the time-frequency domain, which provides the opportunity to correctly localize even very low amplitude pace pulses (~100 µV) on the background of significant EMG artifacts based on their spectral components. While the time domain does not contain information about the frequency content of the signals, and the Fourier transform cannot depict how the ECG spectra change with time, the S-transform applied in this study provides an accurate presentation of the ECG frequency content changes over time (represented in the second subplot of Figure 2). The calculated Shannon energy over the S-transform of the ECG under analysis provides information about the energy concentration along the frequency axis at a given time instant and successfully localizes the positions of even low-amplitude pace pulses in the time domain.
The proposed algorithm handles equally well (mean Se~100%, mean PPV > 98%) with paced ECGs corrupted or not by EMG noise. This fact also confirms the potential of the S-transform to reliably represent the ECG signal components generated by different sources. The analysis of the results presented in Table 2 and Figure 4 shows that the accuracy of the presented method is mostly influenced by the simulated pacemaker mode, instead of the EMG presence in the tested ECG. The drop downs in Figure 4 appear particularly for ECG signals with simulated bi-ventricular pacing, for which the minimal value of PPV in Table 2 is observed. SNR does not influence the observed high Se = 100%, which means that the EMG artifacts have not impeded the correct detection of the pace pulses. The decreased PPV (mean and minimal values) for SNR in the range 0.2 ÷ 0.75 (see Figure 5) is due to the increased number of false positive detections, generally for ECG recordings with pace pulse amplitude <200µV. This could be explained by the decreased value of mean (abs(ShE)) in Equation (6), which is crossed also by some EMG artifacts.
Comparisons of the performance of the method presented in this study to the accuracies of methods for pace pulses detection reported in the literature are presented in Table 3. It includes an algorithm previously designed by us [22], which was tested over the ECG signals from the test dataset described in section 'ECG database', as well as to a recently reported algorithm [18], which was trained and tested also on 'ArtPacedECGdb' database. Table 3. Accuracy of the designed threshold-based algorithm for pace pulses detection over the test dataset (390 clean ECGs, 312 ECG+EMG) for K = 10. Se and PPV achieved in our previous study over the same datasets [22], as well as Se and PPV reported by Nallathambi et al. [18].

This Study
Previous Study [ In a previous work [22] we have shown that the EMG artifacts significantly reduce the PPV of an algorithm designed to detect pace pulses in the time domain (see Table 3). The Stockwell transformation applied in this study distinctly localizes the pace pulse positions, despite the EMG artefacts, which have a higher amplitude (as evident from the time-frequency contour in the second subplot of Figure 2). The calculation of Shannon energy over S successfully "returns the task" in the time domain, retaining the benefit to highlight even small-amplitude pace pulses.
The results reported by Nallathambi et al. [18] are comparable to those achieved in this study, although the methodology is quite different. In [18] the authors applied the final detection procedure over an ECG sampled at a relatively low sampling frequency (250 Hz). This is done at the cost of specific hardware preprocessing that suppresses everything except the pace pulses in the analogue ECG. In our study we relied on digital processing of high-frequency ECG (16 kHz).
The results presented in this study show that the proposed method is precise for pace pulses detection in ECG corrupted by EMG artifacts (Se = 100%, PPV = 98.0%) as well as for clean ECG signals (Se = 99.9%, PPV = 98.3%). It provides significant PPV increase (about 9% points) for ECG+EMG signals, when compared to an algorithm operating in the time domain [22].

Limitations
Although the presented database contains a variety of ECG arrhythmias combined with different pacing pulses, we consider it a limitation that the presented algorithm for pace pulses detection is not tested on real signals from patients with pacemakers. Unfortunately, such an ECG database is still not publicly available.

Conclusions and Future Work
In this work, we pointed out the advantage of using S-transform to access to the frequency content of the ECG signals and localize the existent pace pulses, as well as the benefit of using Shannon energy to correctly detect even low-amplitude pace pulses. The performance of the designed algorithm is high (Se~100% and PPV~98%), when evaluated over a large database that includes various combinations of ECG arrhythmias and pace pulses with different amplitudes, rising edge and total pulse durations, as well as timing that corresponds to different pacemaker modes.
Future work would be related to further investigations of the frequency range in which the Shannon energy is calculated (selected to be between 1000 and 2000 Hz in this study), and application of the designed algorithm over real ECG signals with pace artifacts.

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