Comparison of computation time for estimation of dominant frequency of atrial electrograms : Fast fourier transform , blackman tukey , autoregressive and multiple signal classification

Dominant frequency (DF) of electrophysiological data is an effective approach to estimate the activation rate during Atrial Fibrillation (AF) and it is important to understand the pathophysiology of AF and to help select candidate sites for ablation. Frequency analysis is used to find and track DF. It is important to minimize the catheter insertion time in the atria as it contributes to the risk for the patients during this procedure, so DF estimation needs to be obtained as quickly as possible. A comparison of computation times taken for spectrum estimation analysis is presented in this paper. Fast Fourier Transform (FFT), Blackman-Tukey (BT), Autoregressive (AR) and Multiple Signal Classification (MUSIC) methods are used to obtain the frequency spectrum of the signals. The time to produce DF was measured for each method. The method which takes the shortest time for analysis is selected for real time application purpose.


INTRODUCTION
During AF, atrial electrical activity is described as chaotic and random [1].Frequency domain analysis can help to interpret such activity.DF analysis has been used in several fields such as acoustic emission [2], speech perception [3] and cardiac arrhythmias [4].DF is defined as the frequency of the signal at the point where the frequency spectrum has the maximum value.
The current available treatment of AF consists of medication, electrical cardioversion and ablation.Typically, antiarrhythmic medications such as amiodarone, propafenone and flecainide are used for sinus rhythm control [5] but they are not effective to prevent irregular rhythms.
If medication is not suitable to control the sinus rhythms, electrical cardioversion is considered [6].The procedure is called defibrillation.It can be performed in 2 different ways: at the hospital (electrical cardioversion) or using an implantable cardioverter-defibrillator device (ICD).
When medication and electrical cardioversion control fail, ablation is usually the next step to be considered.Catheter-based radio frequency ablation therapies offer a chance to cure AF.The chances for people survival were 98% and 95% in year 1 and year 2 respectively [5].The ablation performed in early paroxysmal AF has the highest chance for success compared to its use for persistent and permanent AF [7].The success rate of ablation at DF site is shown by significant prolongation of the atrial fibrillation cycle length (AFCL) compared to the site with non-dominant frequency [8,9].This result supports the use of DF mapping to identify suitable ablation targets.
Catheter ablation of AF is a complex interventional electrophysiologic procedure.Its risk is higher than that of ablation for other cardiac arrhythmias [10].The worldwide survey of AF ablation reported that 6.0% complication was recorded (524/8745) [11].Univariate analysis stated that sex, age, insertion difficulty, length of the inserted catheter, type of catheter and, term of insertion were among contributors for catheter related blood stream infection [12].The risks include requiring prolonged hospitalization, long-term disability or death [13].
Since the duration of the ablation and overall proce-dure is one of the contributors to blood stream infection, computation time is critical for any technique to be used during any surgical procedure.Shortening the analysis time can significantly influence the risk of overall surgical procedure.Our aim in this study is to compare the computation time for four different spectrum analysis techniques and compare the different techniques in their ability for producing accurate results.With current technology the data collection capability for non-contact mapping systems consists of catheter mounted multielectrode array allowing the recording up to 3000 points of virtual electrophysiological data at 1200 Hz.With this number of data points, it is important to identify an approach that could be used to produce DF with minimum computation time.
This will in turn help electrophysiologist to identify the dangerous frequencies as target ablation sites and to apply ablation effectively in an expeditious manner.Frequencies between 6-15 Hz are classified as 'dangerous' as the atria fires in an irregular and chaotic fashion.
The main objective of this study is to reduce the risk during catheter procedure by selecting the technique for estimation of DF with the shortest computation.

Fast Fourier Transform (FFT) Technique
The application of the Fast Fourier Transform (FFT) is important in signal processing computations using specially in spectral analysis [13].The discrete Fourier transform of   x nT is given by The FFT algorithm is the main choice for obtaining the Discrete Fourier Transform (DFT) because of its speed.
Welch technique based on FFT algorithm is used to compute the periodogram.This technique uses signal segmentation and averaging to improve statistical properties of the spectral estimates [14].The equation for computing the averaged periodogram is: where is the periodogram of the m th segment of data.

  ˆm k P f
The ends of the data sequence are tapered to zero when data is windowed.DF was obtained from the Welch/ FFT analysis.There are a few considerations to the use of DF analysis especially for the electrograms [4]: a) The estimation of the power spectrum uses the FFT because this is an efficient and widely available method.
b) The maximum frequency is inversely proportional to the sampling interval t  as t  2 1 .
The frequency resolution in the power spectrum is inversely proportional to the length of the signal where T is the length of the signal collected   This has implications for the interpretation of the DF obtained.

Blackman-Tukey (BT) Technique
This method was proposed by Blackman and Tukey (1958) [13].The estimator is based on the Wiener-Khinchin theorem, which states that the power spectrum density is the Fourier Transform of the autocorrelation function of the series [15,16].
If fewer data points are used, the variance of the autocorrelation estimate increases and the estimate become less reliable [13].The averaging associated with windowing a series decreases the resolution of the method, from the frequency intervals of 1/N, to a windowed frequency interval of about 1/M, where M is the length of the window.As a result, wider windows yield higher spectral resolution, and vice versa.The drawbacks of Blackman-Tukey are suppression of weak signal main-lobe responses by strong signal sidelobe and frequency resolution limited by the available data record duration.
This method involves 3 steps [13]: 1) The autocorrelation sequence of the data is estimated using the formula; where: x n is the n-th value of the data series, N is the total number of points in the data series m is the autocorrelation lag ^ signifies estimated value.
2) Windowing the autocorrelation sequence.The windowing of the autocorrelation sequence is normally carried out using a hamming window that has the following characteristics: 3) Finding the Fourier transform to yield the following PSD estimate [15]: where: w (m) is a window function of length 2(M-1)-1, which is zero for |m|

M-1 
The PSD estimate is determined over the frequency where f s is the sampling fre-quency.

Autoregressive (AR) Technique
Autoregressive technique involves selection of model order and the estimation of model parameters from the collected data.This technique is widely used because it has a rational transfer function and the algorithm to estimate the parameters results in a system of linear equations [15].In this model, the data x n are considered to be the output of a system with a random input u n .The relation between input and output are described by the following difference equation which represents the auto-regressive moving average ARMA (pole-zero) model: where: p,q are model orders for the AR and the MA parts, respectively a,b are model parameter (a pk being the k th parameter of the p th -order model) The autoregressive (AR) part of this general equation is given by: Knowing p past values of the series we can estimate the next output as: Equation ( 8) defines the AR with all-pole model.While, the MA model is given by: The prediction error e pn is defined as the difference between the actual sample value x n and its predicted value ˆn x .
The power spectrum of an AR process is: where 2 e  is the variance of e pn .

Multiple Signal Classification (MUSIC)
The multiple signal classification is the improvement of Pisarenko harmonic decomposition [17].This is an eigen-based subspace decomposition method and is best used to estimate the frequencies of complex sinusoids in additive white noise [18].
Assume that y(n) is a random process that consists of p complex exponentials in white noise, The autocorrelation matrix of a noisy signal can be written as: where: xx is the autocorrelation matrix of the signal is the autocorrelation matrix of the noise.
Eigen-analysis is used to partition the eigenvectors and the eigenvalues of the autocorrelation matrix of a noisy signal into 2 subspaces, signal subspace and noise subspace.The eigenvector can be divided by two groups, the p signal eigenvector corresponding to the p largest eigenvalues and M -p noise eigenvector corresponding to the smallest eigenvalues [19].
Ideally, p of these roots will lie in the unit circle on the frequencies of the complex exponentials.The smallest eigenvalue determines the noise variance and the corresponding eigenvector is the prediction polynomial for the clean signal [20].
In MUSIC, the effects of spurious peaks are reduced by averaging.The power spectrum estimate is defined as where N is the dimension of the eigenvectors v k is the k-th eigenvector of the correlation matrix p is the dimension of the signal subspace.Since

 
xx P f has its zeros at the frequencies of the sinusoids, hence the reciprocal of   xx P f has its poles at these frequencies [18].Then, the MUSIC spectrum can be written as

METHOD
The spectrum estimation analysis was tested using sinusoid signals in noise to identify their DF using the 4 different approaches.In FFT analysis, Welch method estimates the frequency spectrum of the signals.The Fourier transform was calculated separately for each of the segments.The analysis was performed by taking the average of the segments, each with 2048 points, using a 4096-point FFT and an overlap of 50% with a hamming window.
In Blackman-Tukey method, the autocorrelation was determined using 4096 points.This autocorrelation function was windowed with the same hamming window used for the FFT approach (4096-point).
The autoregressive spectrum was determined based on the finite linear aggregate of the previous value of the process and the current value of white noise.The model order used was p = 8.As for MUSIC, two complex exponentials with the same window are used to identify the spectrum.

RESULTS
The spectrum estimation was determined using several test frequencies: 8.0 Hz, 7.7 Hz, 7.0 Hz, 6.3 Hz, 5.5 Hz and 5.0 Hz. Figure 1 shows the graphical results of FFT, Blackman-Tukey, Autoregressive and MUSIC in which tested using sinusoidal wave of 8 Hz.DF is being consistently estimated as 8Hz.
Table 1 shows the DF value comparison for all techniques against the base Frequency.
The results show consistent value of DF produced by all techniques.
Table 2 shows that, the percentage error is of the same order of magnitude for these 4 techniques.For calculation purpose we are using 3000 data series to calculate total computation time.As mentioned before 3000 points is the state-of-art for mapping atrial activation.
Table 3 shows a significant difference in processing time between the 4 techniques tested.In summary, FFT takes 22.8 seconds to process 3000 signals compared to 340.9 s for BT, 206.4 s for AR and 10851.8s for MUSIC.As discussed before, all techniques produce the same estimation for DF.Therefore, the time taken to produce the estimation is the main factor for the choice of which technique should be used.

CONCLUSIONS
The study shows significant difference in computation time between the techniques while the estimated value for DF was the same for all four techniques.The choice of a particular technique contributes to the overall surgical procedure time.As discussed before the duration of catheter deployment is a factor that can contribute to blood stream infection.Therefore, the selection of the fastest technique will reduce the risk of the procedure.
As shown by the result of this study, FFT is the best technique as it produces accurate estimates of DF in a speed compatible for quasi real-time analysis.

Figure 1 .
Figure 1.Comparison of Fast Fourier Transform, Blackman-Tukey, Autoregressive and Multiple Signal Classification spectral estimates.

Table 1 .
DF value comparison for all technique against the base Frequency.

Table 2 .
Calculation of error for FFT, Blackman-Tukey, Autoregressive and Multiple Signal Classification spectral estimates.

Table 3 .
Computation time for 3000 signals.