Analysis of QRS-T subtraction in unipolar atrial fibrillation electrograms

This paper presents a QRS-T subtraction approach for atrial fibrillation (AF) intracardiac atrial electrograms (AEG). It also presents a comparison between the proposed method and two alternative ventricular subtraction techniques: average beat subtraction (ABS) using a fixed length window and an approach based on flat interpolation for QRS cancellation. Areas of the atrium close to the mitral valve showed stronger ventricular influence on the AEGs when compared with the remaining atrial regions. Ventricular influence affects the spectral power distribution of the AEG and can also affect the estimation of the dominant frequency unless the whole ventricular activity influence (QRS-T) is removed. The average power after QRS-T subtraction is significantly reduced for frequencies above 10 Hz (mostly associated with QRS complexes), as well as for frequencies between 3 and 5.5 Hz, (mostly related to T waves). The results indicate that the proposed approach removes ventricular influence on the AF AEGs better than the QRS cancellation method. Spectral analysis showed that both the ABS and the proposed method do well and no method should be preferred to the other. In the time domain, the proposed approach is matched to the lengths and timings of onset and offset for individual QRS-T segments while the ABS approach uses an arbitrary length around the QRS for the pattern used for QRS-T removal.


Introduction
Atrial fibrillation (AF) affects about 1 % of the worldwide population [7] and its rate increases with the population age to over 9 % for those above 64 years old [5]. AF classification evolves from paroxysmal to persistent, long persistent and permanent [24]. In AF, the atria do not empty well between contractions and this can lead to the formation of blood clots and thromboembolic events, increasing the risk of stroke fivefold [4]. AF is responsible for frequent medical appointments and hospital admissions, contributing to elevated costs of health care [2]. Despite intensive research, the mechanisms of generation and maintenance of AF are still not well understood.
Cardiac signals from both surface and intracardiac recordings have been widely used in AF studies aiming to understand the atrial electrical behaviour to improve efficacy of interventional ablation treatment, but the presence of ventricular activity (VA) on these signals has hampered the analysis with possibility of distorted results. It is important to study atrial activity (AA) without the influence of the VA to improve understanding of AF genesis and maintenance. As AA and VA overlap in time and frequency domains within the relevant AF frequency range (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) Hz) [16], the use of linear filtering solutions is ineffective to separate them and isolation of the AA is difficult.
Averaging consecutive QRS-T segments to generate a template that is then subtracted from the raw signals has been the method of choice for ventricular cancellation [18]. One classical implementation of this approach is the average beat subtraction (ABS) method, which relies on the fact that AF is uncoupled to VA [18,19]. This method is characterised by a single-lead analysis of the ECG and by subtraction of a mean beat template. The template has the length of the shortest RR interval [19]. The alignment is done either with the QRS fiducial point at the centre of the template or at the sample located after 30 % of the beginning of the template, with the length of the template equal to the shortest RR interval and with the T-wave end defined as the last sample of that segment [18]. A limitation of this approach is that the QRS onset and T-wave end are fixed and the 30-70 % ratio used in their study frequently produces truncation of the T wave, but other choices might give better results.
An alternative technique recently published defines the QRS onset at 41.67 ms prior to the QRS fiducial point and the QRS offset at 50 ms after the QRS fiducial point and avoids the problem of identifying a QRS template by replacing the signal over the QRS period with an interpolation, with the authors favouring a flat line interpolator [1]. The three limitations of this approach are: (1) the length of the QRS to be removed is arbitrary and might not correspond to the actual onset or offset of some QRS complexes, resulting QRS residuals on the 'clean' signal; (2) it does not attempt to preserve the AA during the QRS, as that part of the signal is simply replaced with a flat line, for instance, and (3) it does not remove the ventricular repolarisation component.
This study presents an improved VA cancellation method for unipolar noncontact intracardiac atrial electrograms (AEGs). The method does not suffer from some of the limitations of the classical ABS approach or of the 'QRS cancellation' method, as it is based on a segmentation of each individual QRS complex [13] and T-wave end on the ECG [15]. The QRS-T templates are matched both in shape and in length to the individual signals (both onset and offset of the QRS-T complex). It also preserves the AA of the interval corresponding to the QRS-T interval.
The proposed method was used to remove VA for persistent AF (persAF) patients and the results were analysed in the time and frequency domains before and after VA subtraction. The proposed method was also compared to the classical ABS and the 'QRS cancellation' approaches.

Study population
This study enrolled eight male symptomatic, drug refractory patients (age 47 ± 4.2 years) with longstanding persAF (duration 103 ± 27 months) who underwent an electrophysiological (EP) procedure guided by the EnSite Array balloon with noncontact mapping (EnSite 3000, St. Jude Medical Inc.). The noncontact array balloon with 64 electrodes was placed in the left atrium transeptally, and electrograms were collected using inverse solution mathematics and interpolation in an inside-out fashion allowing the reconstruction of up to 3,600 simultaneous unipolar AEGs projected to the atrium's anatomical 3D representation performed at the beginning of the EP procedure. The array was not moved after geometry creation to avoid distortion of the isopotential maps [3,6,8,11,17,21] and the distance between the centre of the balloon and the patient's endocardium wall did not exceed 4 cm (checked with fluoroscopy).
The EP procedure was followed by ablation of the four pulmonary veins' (PVs) ostia, aimed at the electrical isolation between the potential firing triggers localised in the PVs [10] and the left atrium. Ablation was performed either using point-by-point ablation or with multielectrode PV ablation catheter. The study had informed written consent from all patients undergoing AF ablation prior to the EP study. Collection of patient data was performed in accordance with institutional ethics guidance.

Electrograms
For each patient, the ECG and 2,048 simultaneous left atrial AEGs were recorded over 20 s. The minimum and maximum numbers of QRS-T complexes for the recording period were, respectively, 24 and 37 (population average 29.62 ± 4.40). The AEG signals were sampled at 1,200 Hz and high-pass filtered at 1 Hz. The ECG was band-pass filtered between 0.5 and 50 Hz. Both filters are embedded in the EP mapping system.

Proposed QRS-T cancellation method
The general block diagram for the proposed QRS-T subtraction method is presented in Fig. 1. Firstly, the R peak, QRS onset and T-wave end were detected on the ECG (lead I for 6 patients and V5 for 2 patients) using the methods described by Madeiro et al. [13] (for QRS segmentation) and Qinghua et al. [15] (for detection of T-wave end). After segmentation of all QRS-T complexes in the ECG, the corresponding time locations of the fiducial points were used for the AEG signals (Fig. 2). A pattern was obtained using median operator after their QRS-T fiducial points had been aligned. The instant of time with the maximum cross-correlation between the pattern and the local VA in each intracardiac segment was used for time alignment. Then, subtraction was performed between the local VA and the pattern. The steps are detailed further on.
To detect the QRS onset and T-wave end, first of all a 50 Hz notch filter was applied on the ECG to remove mains hum. Then, the timing position of the R-wave peaks and their respective QRS onset-offset were identified [13]. The R-wave detection is divided in two stages: training and analysis. In the training stage, an initial 10 s long interval of the ECG is used. A detector that includes wavelet transform, first-derivative filter, Hilbert transform and squaring function is applied over this training interval (see details in [13]). Then, we detect the earliest QRS fiducial points on the ECG and compute the average and standard deviation of the intervals between beats. After the training stage, we use an adaptive threshold technique to detect each QRS complex over the remaining of the ECG signal. This threshold is defined as [13]  If eventual occurrences of false-positive or false-negative are identified, based on the stored average and standard deviation of the intervals between beats, then the above referred pre-processing techniques are re-applied over the corresponding signal excerpts to emphasise QRS complexes and attenuate artefacts, noise and P/T waves (see details in [13]).
After detection of QRS fiducial points, we compute the QRS envelope. This envelope is defined as the absolute value of the analytical (complex) signal whose real and imaginary parts are the filtered ECG (real part) and the Hilbert transform of the filtered ECG (imaginary part) [13,14]. Considering the QRS complex as a modulated waveform, the beginning and end of the QRS complex envelope calculated using the Hilbert transform coincide with the QRS onset and offset [13,14,20]. For QRS delineation, we search for onset and end point within a time window defined by [R P -300(ms), R P ? 190(ms)], where R P is the location of the detected R peak and, for each window, the filtering process previously described for QRS detection is used. Calling the QRS envelope signal E[n], an area indicator A[n] is computed for a given QRS envelope [9,13] A½n to the offset of each QRS complex [9]. For QRS onset detection, the surface indicator A[n] is defined as reaching its maximum value at the beginning of each QRS envelope. After the QRS segmentation, a baseline correction was applied by means of a cubic spline interpolation anchored on the QRS onset points [12]. The detection of T-wave end was also based on an indicator of the area of the T wave [15] and no additional pre-processing was necessary. Letting Q f and Q i be the samples (indices) related, respectively, to each QRS offset and the subsequent QRS onset, we define a moving sample k, Q f ?w B k B Q i -p-1, where p is the number of samples corresponding to 16 ms and w is the number of samples corresponding to 128 ms [15]. Then, for each k value, the following metrics were calculated If the inequality is satisfied, then the T-wave end in the current QRS offset-onset interval is located at the sample k 0 , if k 0 [ k 00 , or k 00 , if k 0 [ k 0 (biphasic T waves). Otherwise, the T-wave end is located at the sample with the maximum amplitude of |A(k)| (monophasic T waves) [15].
After the QRS-T segmentation (black thick line in Fig. 3a), a pattern using median values was computed. First, the ventricular fiducial points were aligned and then a median operator on each of the matching points was applied. Finally, each local VA delimited in the AEG signals was aligned with the computed pattern and then subtracted (Fig. 3d).
A comparison between three cancellation methods on an AEG signal is illustrated in Fig. 3. In Fig. 3a, QRS-T segments are highlighted in the original electrogram. Figure 3b presents an AEG signal after QRS cancellation [1]. Figure 3c and d show the AEG signal after QRS-T subtraction, respectively, using the ABS and the proposed approach. The patterns used by both QRS-T methods to cancel the VA are presented in Fig. 3e and f, respectively.
To select the sliding window segment to compute a QRS-T pattern, we initially define the first segment as the interval between 0 and 7 s. Then, for the next segment window, we define its end as the location of the next T-wave end and its beginning as 7 s earlier. We proceed with an analogous approach till the last QRS-T complex.
To calculate a QRS-T pattern in a given segment window, we initially compute the mean QT interval length, called Lpat. The detailed steps are described as follows: 1. We define an array of zeros with the same length of Lpat; 2. Then, we select the AEG segment which begins at the time instant related to the first QRS onset and ends at the time instant related to the first T-wave end point; 3. If that segment is larger than Lpat, then we fill out the array with the segment from the QRS onset until the sample related to the duration given by Lpat; 4. Otherwise, if the segment referred in the step (2) is shorter than Lpat, then we fill out the array referred in the step (1) with the corresponding segment and we leave the last array samples with zeros; 5. For each other window in the AEG segment beginning at the next QRS onsets and finishing at the corresponding T-wave ends, we align it with the window referred in the step (2) using cross-correlation; 6. After aligning each QRS-T window with the window referred in the step (2), each one fills out an array of zeros similar to that one referred in the step (1), observing the same rule explained in steps (3) and (4); Fig. 3 Ventricular subtraction in AEG signals: a AEG raw signal, b after QRS cancellation, c after QRS-T subtraction using ABS, d after QRS-T subtraction using the proposed method, e QRS-T patterns for the proposed approach considering the intervals between 0 and 7 s and 3-10 s, respectively, and f computed patterns calculated by the ABS approach for the same intervals 7. We build a matrix (Ns 9 Lpat) with the window defined in the step (2) and the set of aligned windows, where Ns is the number of QRS-T complexes in the corresponding AEG segment; 8. We obtain the median of each column of the matrix and a resultant array is built with the set of median results, which corresponds to the QRS-T pattern associated to that AEG segment.
The VA subtraction in the AEG signals is performed as follows: 1. For each QRS-T complex in the AEG segment, the delay between the corresponding fiducial points from the AEG raw signal and the QRS-T pattern was obtained using cross-correlation. 2. Then, for each QRS-T complex in the AEG segment, we select a window that begins at the sample related to QRS onset corrected by the delay and ends after Lpat. 3. Finally, we subtract the local VA influence identified in the AEG by the QRS-T pattern.

Time and frequency analysis approaches used to compare the cancellation methods
The segmentation of each QRS complex and T wave in the ECG is the first step for QRS-T subtraction in AEGs.
Examples of QRS and QRS-T segmentation using the three different methods are presented in Fig. 2. The first approach (Fig. 2a) was implemented by Ahmad et al. [1] and the QRS onset (marked with squares) and offset (circles) are marked, respectively, as 41.67 ms (50 samples) prior to the QRS fiducial point, and 50 ms (60 samples) after the QRS fiducial point. Figure 2b presents the classical ABS approach: The QRS-T pattern length has the same length as the shortest RR interval [18]. In addition, the onset (marked with squares) and offset (circles) from the QRS-T complex were defined such that 30 % of the corresponding window precedes the QRS fiducial point and 70 % follows it. The proposed method is presented in Fig. 2c. The delineation of each QRS-T complex is done by individual detections of QRS onset (marked with squares) and T-wave end (circles) [13,15]. The corresponding time location of each fiducial point detected in the ECG is marked in the AEG signal (Fig. 2d). The dotted vertical trace on the left corner in Fig. 2 shows that the QRS fiducial points of the ECG (Fig. 2a-c) are aligned with the VA disturbances on the AEG (Fig. 2d). The rectangular solid trace box between around 6.0 and 7.5 s highlights the ventricular segmentation in a single VA on the ECGs proposed by the three methods ( Fig. 2a-c) and the corresponding AEG (Fig. 2d). The bold black solid traces in the AEG (Fig. 2d) are the projections of the QRS-T segments as identified by the proposed method.
The first two vertical traces localised before the VA fiducial point in the AEG are the QRS onset segmentation, respectively, identified by the ABS approach (dashed trace) and the proposed method (dotted trace). In addition, the T-wave end segmentation identified by the ABS (dashed trace) and the proposed method (dotted trace) are highlighted on the trace in Fig. 2d. Spectrum analysis was used to help compare the techniques in the frequency domain and was performed using the fast Fourier transform (FFT) with a hamming window to reduce spectral leakage. Zero padding was introduced to improve spectral peaks identification (frequency steps of 0.05 Hz).
The frequency spectrum impact of each ventricular cancellation technique in different areas of the left atrium was also investigated (Fig. 4). This allows one to observe the ventricular influence on the AEGs in areas of relevance to the arrhythmia (PVs, roof and septum) and to observe the VA impact in areas even far from the mitral valve (MV). Figure 4 presents the average spectrum of four preclassified areas of the left atrium: close to the PVs, roof, septum and close to the MV. For each area, 12 points were selected following their identification. For the PV region, three points were selected close to each of the four PVs. Spectrum analysis allowed identifying the spectrum for each individual AEG over 20 s. The AF physiological range was between 3 and 12 Hz using a window of 4 s, and an overlap of 50 % (moving 2 s at a time). The spectral resolution was 0.25 Hz and zero padding was applied to produce frequency steps of 0.05 Hz. The average spectrum from the 12 points for each 4 of 20 s was calculated (total of 9 windows). The average of all windows represents the average spectrum of each particular left atrium area, making possible to compare spectra for each ventricular far-field cancellation.
The average spectrum of the entire left atrium was also investigated. In this analysis, we performed the average of the spectra of the 2,048 simultaneous AEGs for the different protocols: raw signals, after QRS cancellation and after each of the QRS-T subtraction methods using the spectral analysis strategy previously described. The whole left atrium spectrum allows one to observe the global effect of the cancellation techniques on the AEGs. Figure 4 also hints that a wrong identification of the dominant frequency (DF) is possible if VA is not properly removed: on the spectra corresponding to the areas close to the MV, if one looks at the raw signal, a DF of about 4.8 Hz would be identified. For 'QRS only' removal, a DF of 2.8 Hz would be identified and for QRST-ABS or our QRST segmentation approaches the correct DF, at 6.5 Hz, is identified.
The ratio between the frequency components' energy in the range (5.5-12 Hz) and (3-5.5 Hz), defined here as P ratio , is shown in Tables 1 and 2 for each subtraction method. The  Table 1 are classified by specific regions of the left atrium. In Table 2, this ratio is calculated for all of the 2,048 AEGs (the whole of the left atrium) for each patient. This metric, P ratio , allows identifying the impact of the T-wave subtraction over the AEG's spectrum, where higher values mean better repolarisation subtraction (small residual levels of the T-wave influence on the AEG signals). The equation for the power ratio index is where P ratio is the metric (similar to a SNR) that we use to measure the quality of VA removal and P a-b is the area under the curve for the power spectrum between the frequencies a and b, in Hz. Figure 2 highlights the main differences between the proposed method and the ABS method for segmentation of the QRS-T: as the ABS method is based on a fixed time length around the QRS fiducial point, as we mentioned before, the 30-70 % ratio used in their study frequently produces truncation of the T wave (as in the case shown in Fig. 2b), but other choices might give better results. One justification for using such approach is that the QRS fiducial point is easier to identify reliably than the end of the T wave. Provided that the delimitation of the T wave is done correctly, we argue that the proposed method should be preferred. Figure 2d highlights this difference: Following the proposed approach the segment highlighted in the AEG (segment between the dotted lines) would be processed, while following the ABS approach the segment between the dashed lines would be processed (clearly not cancelling the whole of the T wave, see Fig. 2c).

Frequency domain analysis
The frequency spectra of the AEGs illustrated in Fig. 3(a-d) are shown in Fig. 5(a-c). The frequency spectrum of the AEG signal after the proposed QRS-T subtraction method (black trace in Fig. 5a-c) is compared with the spectra of its respective raw signal (Fig. 5a, grey trace), the electrogram after QRS cancellation (grey trace in Fig. 5b) and the electrogram after the QRS-T subtraction by the ABS method (Fig. 5c, grey trace).
The results of using P ratio for the comparison of the techniques are summarised in Tables 1 and 2 and show that either the ABS method or the technique proposed here do better than QRS only removal technique and, as expected, the 'raw' signal (with no QRS or QRS-T removal). Based on the spectral analysis on AF electrograms, both the ABS method and the proposed method do well and no method should be preferred to the other.

Discussion
In the ABS-based approach for QRS-T subtraction, several QRS-T segments showed an inaccurate delimitation of the QRS onset and T-wave end (see Fig. 2b). In addition, the ABS-based template is often longer when compared with that of the proposed approach (Fig. 3e-f). We observed that the main reasons are (Fig. 2): (1) for the presented case, the 30 % segment that precedes each QRS fiducial point covers a time longer than the 'true' interval between a QRS onset and its fiducial point; (2) the 70 % segment after QRS fiducial point often falls before the T-wave end. Instead of using 30/70 other choices for the ratio such as 20/80, 15/85 are used, that might give better results, but we have not tested that. The proposed approach aims to match the actual length and timings of onset and offset for all individual QRS-T segments and does adaptively using a 7 s window.
For Ahmad's QRS cancellation approach [1], it can be observed that: (1) QRS onset and offset are based on a fixed window, which certainly does not reflect the dynamic changes in the ventricular depolarisation; (2) there is no cancellation of the ventricular repolarisation component; (3) there is no alignment between the QRS fiducial point in the ECG and the local VA detected in each AEG signal; (4) during the subtraction, an interpolation method is applied (spline, flat or linear) to connect the onset and offset points, so any AA that overlaps with the QRS segment is also removed.
Regarding the frequency domain analysis, the AEG spectrum shown in Fig. 3 illustrates that the proposed method attains a considerable attenuation for frequencies above 7 Hz (see also Fig. 5a). Comparing the results of the QRS cancellation method with the proposed approach (see Fig. 5b), we observed that the frequencies between 3 and 5 Hz were reduced suggesting that this range is mostly associated with the T-wave frequency content. For frequencies above 7 Hz, both methods produced similar attenuation which, we argue, is mostly related to the subtraction of the QRS activity. In addition, for frequencies between 5 and 7 Hz, only a slight reduction in power was observed, basically preserving the signals within that frequency range. The spectra for both QRS-T subtraction approaches showed a similar overall behaviour for this illustrative example.
The average spectral behaviour for the four distinct areas of the atrium is shown in Fig. 4. We observed different degrees of VA contribution over the raw signal across the spectrum. Areas close to the MV showed stronger frequency components (higher amplitudes) related to VA when compared with areas such as the roof and around the PVs. Regarding the proposed approach, it was noticed that components of frequencies above 10 Hz, and also the components of frequency range around 3-6 Hz (3-5.5 Hz PVs, 3-5.2 Hz roof; 3-6 Hz septum; and 3-5.7 Hz MV) were attenuated due to, respectively, ventricular depolarisation and repolarisation subtraction. An agreement on these ranges was also observed when the overall average spectrum of the atrium was calculated. When both QRS-T subtraction methods were compared, similar results were identified in the frequency domain.
As an additional analysis, we also investigated the effectiveness of the ventricular repolarisation subtraction. The ratio between the frequency components related to physiological AF range excluding the VA repolarisation (5.5-12 Hz) and the reported frequency components related to VA repolarisation (3-5.5 Hz), P ratio , was computed (Eq. 7; Tables 1, 2). Previous studies have presented a slightly different approach [22,23], using a ratio between 4 and 9 Hz over the total spectrum power with focus on the QRS template cancellation on the ECG. Comparing P ratio of signals which had 'QRS cancellation' with P ratio of the raw signal, in Table 1, we observed that it has decreased for all atrium areas. This suggests that (1) the approach of  Fig. 3a after QRS-T subtraction using the proposed method (black trace) compared with the spectra (grey trace) of the: a raw signal, b after the QRS cancellation, and c after QRS-T subtraction (ABS) QRS segment cancellation using flat interpolation results in a distorted reduction for all frequency components; (2) in particular the power above 10 Hz (corresponding mostly to the QRS complex) had a great reduction when compared with the raw signal; (3) for the T-wave frequency range, we observed a minor reduction, although reduction of T-wave activity is not a target for this cancellation method.
When the signals after both QRS-T subtraction methods were compared with the previous results, an improvement of P ratio for all atrium areas could be observed. This suggests that: (1) frequency components corresponding to both QRS and T wave were dramatically reduced; and/or (2) the remaining frequency content, which is mostly related to the main AA within the physiological AF range (above 5.5 Hz), clearly stands out from the frequency content related to the range (3-5.5 Hz), as consequence of an effective repolarisation subtraction.
When the ratios between both QRS-T subtraction methods were compared, the ABS approach produced a slightly higher P ratio around the MV and in areas close to the PVs (2.45 and 3.38) when compared with the proposed method (2.41 and 3.26). In contrast, for the roof and septum, the proposed method presented slightly higher ratios (5.33 and 3.76) when compared with the ABS (5.07 and 3.71). We extended the analysis considering the average spectrum for the 2,048 AEG signals of the whole atrium for individual patients ( Table 2). As we can see from both tables, our method produces similar P ratio index results to those of the ABS approach, but our templates are matched to the actual length and timings of onset and offset for individual QRS-T segments rather than using an arbitrary length based on the shortest R-R interval (Fig. 3) and positioning the QRS onset fiducial point at either 50 or 30 % of the total time. Considering the capabilities of the proposed method for adapting itself to the existing interand intra-patient variability in QRS durations and QT intervals, which define the length for each QRS-T complex, we argue that it should be preferred to the ABS approach.
Finally, concerning the limitations, it is well known that T-wave end location is a difficult problem, especially when the signal to noise ratio is low [15]. Also the presence of the f-waves related to AF over the ECG adds an extra challenge, e.g. an extra physiological 'noise', to time location of T-wave end. These limitations may directly affect the QRS-T pattern and the accuracy of the QRS-T subtraction in the AEG signals. However, the identification of the average of all QT intervals as the length of QRS-T pattern does reduces the effect of eventual misdetections.
The main clinical objective of identifying DF in AEGs is the location of key candidates for ablation to minimise the amount of burning and maximise the outcome of ablation procedures in persistent AF. Cancellation of ventricular influence on high density unipolar AF electrograms is a prerequisite to avoid 'contamination' of the AEGs that might alias as possible AF triggers. As well as the QRS complex, VA repolarisation also needs to be subtracted and using standard filters is not feasible as the frequency range of VA overlaps with the physiological AF spectrum. It should be performed using individual 'customized' and adaptive QRS-T segmentation and a coherent subtraction strategy. We have presented a novel approach that performs such cancellation and we have compared its performance to two alternative techniques. Evidence has been presented that shows that the new approach is capable of reducing the influence of the VA on the AEGs and that it does not suffer from some of the disadvantages of the alternative techniques it was compared to.