Nonparametric Spectral Estimation Technique to Estimate Dominant Frequency for Atrial Fibrillation Detection ()
1. Introduction
Atrial fibrillation (Afib) is a supra ventricular tachycardia (SVT) in which the pacemaker activity of the sinoatrial node is overwhelmed by disorganized electrical impulses, causing rapid and irregular quivering of the atria. Atrial cells will fire at rates of about 350 - 600 bpm. Reference [1] was amongst the pioneers to propose 3- fold theory to explain Afib mechanism. Multiple wavelet hypotheses were proposed by [2] , and verified by [3] . Reference [4] attained the breakthrough by identifying Afib initiation areas in pulmonary veins. All these models indicated that although Afib appears to be an irregular arrhythmia, it follows a definite pattern of initiation and maintenance [5] . Therefore, researchers preferred frequency domain analysis over time domain analysis to detect: 1) Afib, 2) areas of rapid activations, 3) differences in pathophysiology, and 4) changes in rate as a result of interventions. For frequency domain analysis, the pre-requisite for Afib detection and its further analysis has been the estimation of reliable dominant frequency (DF) using different parametric and nonparametric spectral estimation techniques [6] .
Parametric spectral estimation techniques try to fit a parametric model to signal by minimizing certain cost function. The model order is not known in practice. An underestimation of the model order results in missing of signal poles thus reducing the resolution and causing some harmonic components to be unidentifiable. An overestimation of model order produces poles linked to the noise as well, which makes noise more prominent in spectrum by creating false peaks.
Nonparametric spectral estimation techniques do not make any assumptions on data-generating process or model. Also, the nonparametric Fast Fourier Transform (FFT) based spectral estimation techniques are more time-efficient than the parametric spectral estimation technique for real-time applications (e.g. to minimize the duration of catheter deployment) [7] . Popular nonparametric FFT-based spectral estimation techniques used for DF estimation have been: 1) periodogram, 2) modified periodogram, 3) Bartlett, and 4) Welch.
Periodogram method (without windowing) has been used by [8] -[12] . Modified periodogram method, using Hanning/Hamming windows, has been applied by [13] -[15] . However, the issue of bias and variance with periodogram and modified periodogram methods affects the quality of these spectral estimators [16] . Reference [17] suggested the Bartlett method with Hanning window to minimize the issue of bias and variance. Welch method has been applied by [18] -[20] .
DF is the frequency containing maximum power [21] . Estimation of the local atrial activation rate in AF is the key advantage of DF analysis [22] . For Afib, the DF lies in the range of 3 - 12 Hz [21] -[24] . For periodic activations, e.g. Normal Sinus Rhythm (NSR), DF gives the robust estimation of activation rate. However, for irregular activation like Afib, the estimation may not be as robust. Therefore, DF analysis for Afib is built upon selection of the lowest-noise signal. Regularity index (RI), defined as “ratio of the power spectral area under the DF and its harmonics divided by the total spectral area” has been used to ensure reliability of DF [6] -[20] . Mostly, DF with RI > 0.2 has been considered reliable for Afib detection and its further analysis [25] [26] .
Different researchers have used different nonparametric FFT-based spectral estimation techniques (periodogram, modified periodogram, Bartlett, and Welch) for DF estimation to carry out Afib detection/analysis [6] - [26] . No study has been carried to find out the most appropriate nonparametric FFT-based spectral estimation technique. Also, no details have ever been mentioned regarding calculation of power spectral areas, which are required for RI computation.
In this paper, the objective is to find the most appropriate nonparametric FFT-based spectral estimation technique to estimate reliable DF for Afib detection. As periodogram and modified periodogram are inconsistent spectral estimators and yield fluctuating estimates from successive realizations [27] , Bartlett method [28] using Hanning window, as suggested by [17] , and Welch spectral estimation method [29] are compared. Also, power spectral areas required for RI have been calculated using numerical analysis.
This paper is organized in the following way. The detailed description of the signal processing steps and algorithms is given in Section 2. Section 3 discusses the results. Limitations and conclusion/future recommendations related to this research work are highlighted in Section 4 and Section 5 respectively.
2. Methodology
2.1. Overview
In this work, ICEGMs for Afib and non-Afib (NSR) cases have been acquired from High Right Atrium (HRA) and Superior Vena Cava (SVC) for frequency domain analysis. ICEGMs are recorded using bipolar catheters with minimum inter electrode spacing (IES), which eliminated the need for cancellation of far field ventricular activity [30] . Pre-processing steps similar to [31] with addition of zero padding for better frequency resolution have been used. DF has been calculated using both Bartlett with Hanning window and Welch methods. RI is calculated using Simpson 3/8 and Trapezoidal rules [32] . Afib detection algorithm has been applied to find out Afib and non-Afib (NSR) cases. An overview of methodology is depicted in Figure 1.
2.2. Data Set and Signal Acquisition
The data set consists of ICEGMs obtained from Armed Forces Institute of Cardiology (AFIC), Rawalpindi, Pakistan, and online MIT Physionet database [33] . The study protocol has been approved by local institutional review board.
AFIC data set consists of ICEGMs for both Afib and non-Afib (NSR) cases from de identified patients. The data has been recorded from 10 patients using quadripolar catheter with 2-5-2 mm IES at sampling rate of 2000 hertz. To avoid far field ventricular activity, ICEGMs from HRA Distal catheter have been used for analysis. HRA Distal catheter configuration is shown in Figure 2.
MIT Physionet data consists of ICEGMs for Afib cases from de identified patients. The data has been recorded from the right atria of 8 patients using decapolar catheter with 2-5-2 mm IES at sampling rate of 1000 hertz. To avoid far field ventricular activity, ICEGMs extracted from catheters CS 12 and CS 34 (close to annulus of SVC) have been used for further analysis in this paper. Catheter configuration for MIT Physionet data is depicted in Figure 3.
Figure 2. Catheter configuration-AFIC data.
Figure 3. Catheter configuration-MIT Physionet data.
2.3. Signal Pre-Processing
Pre-processing of ICEGMs is carried out to obtain the signal corresponding to the local depolarization, and to restrict the frequency spectrum to physiological range. These include: 1) band pass (40 - 250 hertz) filtering using 3rd order Butterworth filter, 2) obtaining monophasic waveform by rectification and 3) low pass filtering (15 hertz) using 3rd order Butterworth filter.
2.4. Spectral Estimation Using Bartlett Method with Hanning Window
The Bartlett Method divides the signal segment of length N into K sub segments, with each sub segment having length L = N/K [23] . As suggested by [17] , the modified periodogram method is then applied to each K sub segment. The average of the resulting estimated power spectral density is taken as the estimated power spectrum. Mathematically, it is given as
(1)
In this case, the signal segment (N) of 10 secs duration is defined as an “event”. It is further divided into non overlapping sub segments (K) of 2 secs each. Modified periodogram with Hanning window w(n) of each sub segment (k) is computed. The average of these five spectral estimates gives Bartlett spectral estimate. Frequency resolution is kept to be 0.1221 hertz. Figure 4 depicts the signal segment and sub segments for Bartlett method with Hanning window.
2.5. Spectral Estimation Using Welch Method
The Welch method eliminates the tradeoff between spectral resolution and variance in the Bartlett method by allowing the sub segments to overlap [27] . In the Welch method, the signal segment of length N is divided into K sub segments, with each sub segment having length L = N/K. The K sub segments of length L are overlapped and the modified periodogram are computed from the overlapped K sub segments. Also, the periodogram are normalized by the factor U to compensate for the loss of signal energy owing to the windowing procedure [22] . Mathematically, it is calculated as
(2)
where D is the offset between two consecutive sub segments and U is given as
(3)
The signal segment (N) of 10 secs duration is defined as an “event”. It is divided into overlapping sub segments (K) of 2 secs each. The overlap percentage is kept to be 50% (D). Modified periodogram of each sub segment with Hanning window w(n) is taken. The average of these nine spectral estimates gives Welch spectral estimate. Frequency resolution is kept to be 0.1221 hertz. Figure 5 depicts the signal segment and sub segments for Welch method.
2.6. Dominant Frequency Estimation
The frequency having maximum peak in Power Spectral Density (PSD) estimate within range of 1 - 12 hertz is declared as DF. The maximum peak and the peaks of up to three harmonics (if exists) are found for both methods.
Figure 4. Signal segment and sub segments for Bartlett method with Hanning window.
Figure 5. Signal segment and sub segments for Welch method.
2.7. Calculation of Regularity Index (RI)
To ensure the reliability of DF, RI is computed by dividing the spectral area under DF and its harmonics within approximately1 hertz base by the spectral area for 3 - 12 hertz. The spectral area under curve is computed using Simpson 3/8 and Trapezoidal rules by using algorithm as depicted in Figure 6.
2.8. Afib Detection Algorithm
Afib detection algorithm is based on DF and RI. Figure 7 depicts the detection algorithm for Afib. After ICEGMs pre-processing, spectral estimate of the signal have been found using both Bartlett (with Hanning window) and Welch methods. Maximum peak in the range of 1 to 12 hertz has been found and declared as DF. For DF less than 3 hertz or greater than 12 hertz, RI is not calculated. For DF between 3 to 12 hertz, peaks of up to three harmonics (if exists) are also found. RI is computed to ensure reliability of DF between 3 - 12 hertz.
Afib is said to be detected if DF is in range of 3 - 12 hertz, and RI is greater than 0.2. The event is declared as non-Afib for following: 1) DF is less than 3 hertz or greater than 12 hertz, and 2) DF is between 3 - 12 hertz, but RI is less than 0.2.
3. Results and Discussions
In this work, 2500 events (2000-Afib and 500-non-Afib) have been tested by Afib detection algorithm using Bartlett (with Hanning window) and Welch methods. All these events have been marked as Afib and non-Afib by the doctors for comparison and validation of the algorithm.
3.1. Afib Event-I-Detection by Both Bartlett (Using Hanning Window) and Welch Method
In this Afib event, Afib has been detected by both Bartlett (using Hanning window) and Welch methods. Figure 8 depicts the pre-processing steps of ICEGMs. Figure 8(a) presents the ICEGMs recorded using catheters. These ICEGMs have been filtered using 3rd order band pass Butterworth filter so as to accentuate signal corresponding to local depolarization (see Figure 8(b)). Figure 8(c) depicts the absolute value of band pass filtered signal. The rectification has been performed to transform biphasic waveform into monophasic waveform. Figure 8(d) shows the low pass rectified band pass ICEGMs. Low pass filtering (LPF) with cut off frequency of 15 hertz has been carried out to limit the frequency spectrum to physiological range.
PSD for the signal depicted in Figure 8(d) has been computed using both Bartlett (with Hanning window) and Welch methods. PSD comparison is shown in Figure 9. For both methods, DF (maximum peak) turned out to be 6.1 hertz and RI was 0.50. The Afib detection algorithm has successfully declared the event as Afib for both spectral estimation techniques.
3.2. Afib Event-II-Detection by Welch Method Only
In this Afib event, Welch method has been successful in detecting Afib, however, Bartlett method (using Hanning window) declared it as non-Afib. Same pre-processing steps as of Figure 8 have been applied to ICEGMs. PSD has been computed using both Bartlett (using Hanning window) and Welch methods. On PSD comparison, as depicted in Figure 10, Welch PSD has DF (maximum peak) at 5 hertz with RI = 0.23; however, Bartlett PSD has DF (maximum peak) at 4.64 hertz with RI = 0.18. Although, DF using Bartlett method is within range of 3 - 12 hertz, however, RI is less than 0.2. This indicates that sufficient power is not contained by the peak frequency so as to qualify it as reliable DF. The Afib Detection Algorithm has successfully declared the event as Afib for
Figure 6. Algorithm to find area under curve using numerical techniques.
Figure 8. (a) ICEGMs recorded using catheters; (b) Band passed ICEGMS; (c) Rectified band passed ICEGMS; (d) Low passed rectified band passed ICEGMS.
Welch method only.
3.3. Non-Afib Event-Detection by Both Bartlett (Using Hanning Window) and Welch Method
In this non-Afib event, NSR ICEGMs have been tested using both spectral estimation techniques. Both methods give same DF values which remained within range of 1 - 2 hertz. The Afib detection algorithm declared all such events as non-Afib for both Bartlett (using Hanning window) and Welch methods with 100% accuracy.
Figure 9. Power spectral density comparison.
Figure 10. Power spectral density comparison.
3.4. Analysis
For NSR, the Afib detection algorithm has been 100% accurate for either of the spectral estimation techniques. The reason for high accuracy and precision is the inherent periodicity of NSR. The frequency of regular pattern of NSR waveforms has been accurately estimated by both spectral estimation techniques.
For Afib events, it is found that Welch method has more accurately detected Afib events with 98% accuracy; however, Bartlett method using hanning window has an accuracy of 90%. Table 1 and Table 2 depict the range of DF and RI for both methods. It is evident that for Bartlett method (using Hanning window) DF falls within
Table 1. DF range for Bartlett and Welch methods.
Table 2. RI range for Bartlett and Welch methods.
range of 3.42 - 11.47 hertz with RI range 0.15 - 0.73. However, for Welch method the DF ranges between 3.42 - 9.52 hertz with RI range of 0.18 - 0.70.
The reasons for high accuracy of Welch method are: 1) variance of Welch method is less than that of Bartlett method, 2) averaging nine overlapped spectral estimates for Welch method in comparison of five nonoverlapped spectral estimates of Bartlett method helped in reduction of estimation error, and 3) reduction is estimation error helped to reproduce reliable DF more closer to the local activation rate.
4. Limitation
In this work, ICEGMs for analysis have been acquired from HRA and SVC only. ECG signals, ICEGMs from left atrium, and ICEGMs from other locations of right atrium have not been studied in this paper; however, this limitation does not affect the accuracy of this research work. To address the limitation, the atrial signals can be extracted from ECG/ICEGMs, and same algorithm can be applied on extracted atrial signal for estimation of DF.
5. Conclusion
Clinical applicability of frequency domain analysis for Afib has been shown by a number of researchers. However, no study has been carried out to standardize the nonparametric FFT-based spectral estimation technique for estimation of reliable DF used for Afib detection. Periodogram and modified periodogram spectral estimators are inconsistent, and carry a greater degree of estimation error. In this work, comparison of Bartlett method (with Hanning window) and Welch method is carried out. Welch method is found to be more accurate because of its better statistical properties and can be standardized not only for Afib detection but also for other frequency domain analysis like detection of areas of rapid activations, differences in pathophysiology, and changes in rate as a result of interventions. The approach presented in this paper can be modified and applied to ECG signals, ICEGMs from left atrium, and ICEGMs from other locations of right atrium to extend its applicability.
Acknowledgements
The authors wish to thank AFIC, Rawalpindi, Pakistan for their support during the course of this research, and ICT R & D Fund, Pakistan for providing the funding for this work.