Contents lists available atScienceDirect
Medical Engineering & Physics
j o u r n a l h o m e p a g e :w w w . e l s e v i e r . c o m / l o c a t e / m e d e n g p h yThe differential method of phase space matrix for
AF/VF discrimination application
Chien-Sheng Liu
a, Wei-Kung Tseng
b, Jen-Kuang Lee
c, Tze-Chien Hsiao
d, Chii-Wann Lin
e,a,f,∗ aDepartment of Electrical Engineering, National Taiwan University, Taipei, TaiwanbDivision of Cardiology, Department of Internal Medicine, E-DA Hospital, I-Shou University, Kaohsiung, Taiwan cDepartment of Laboratory Medicine, National Taiwan University Hospital, Taipei, Taiwan
dDepartment of Computer Science, National Chiao-Tung University, Hsinchu, Taiwan eInstitute of Biomedical Engineering, National Taiwan University, Taipei, Taiwan fInstitute of Applied Mechanics, National Taiwan University, Taipei, Taiwan
a r t i c l e i n f o
Article history:
Received 18 September 2009
Received in revised form 10 February 2010 Accepted 2 April 2010
Keywords: Electrocardiography Phase space matrix Atrial fibrillation Ventricular fibrillation Real-time analysis
a b s t r a c t
The advances in electrocardiographic (ECG) technology have facilitated the development of numerous successful clinical applications and commercial monitoring products for diagnosing disease and moni-toring health. All of these demand the development of smart algorithms and computational resources for the real-time, early indication of critical cardiac conditions. This study presents the development of a Complex Phase Space Difference (CPSD) algorithm with differential method to analyze spatial and temporal changes in reconstructed phase space matrix, and derives an index for real-time monitoring. We used total of 5306 data segments from MIT-BIH, CU, and SCDH databases and clinical trial data to determine the optimal working parameters and verified the classification capability by using a quan-titative index of this algorithm. With threshold values set to 2.0 and 6.0, this method can successfully differentiate normal sinus rhythm (NSR) signals (1.48± 0.21), low risk of atrial fibrillation (AF) signals (3.71± 0.99) and high risk of ventricular fibrillation (VF) signals (9.38 ± 2.22). It is the first real-time algo-rithm that reports the best performance to distinguish AF and VF with sensitivity of 97.9% and specificity of 98.4%. With self-normalization, the algorithm is not subjected to the inter-variability or sampling size effects. Its computational scheme only requires matrices addition and subtraction, and thus significantly reduces the complexity for real-time implementation. It will be able to adopt in different scenarios of tele-healthcare and implantable applications.
© 2010 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction
Healthcare services globally are rapidly evolving from provider-centric to people-provider-centric, with novel technologies to fulfill the demands of personalized, prevented and home-care services[1]. The goal is to develop better support mechanisms that account for individual differences and satisfy personal needs. The combination of computation with physiological signals can yield an easy-to-use real-time index of health status as well as necessary information to support decision-making in healthcare. Bio-signal processing has been an active field of research in biomedical engineering for many decades. It provides a functional assessment of a target physiological system by evaluating acquired biological signals in either the time or the frequency domain for diagnosis or treatment.
∗ Corresponding author at: Institute of Biomedical Engineering, National Taiwan University, Taipei, Taiwan. Tel.: +886 2 33665272; fax: +886 2 23620586.
E-mail address:cwlinx@ntu.edu.tw(C.-W. Lin).
Commonly, less accurate indices are adopted in monitoring appli-cations, such as four of the primary vital signs – temperature, blood pressure, pulse rate, and respiratory rate. Two of these indices are related to the functions of the heart, which can also provide an electrophysiological signal that can be monitored noninvasively using an electrocardiograph (ECG) and electrodes on the surface of the body for several clinical applications[2,3]. However, medical professionals must typically to interpret the recorded signal traces to diagnose disease and so the application of this approach in the monitoring of health poses significant challenges.
Many studies have addressed the problem from a perspective of time–frequency duality[5,6]. The principle of time–frequency duality is that acquired signals can be viewed in either the time domain or the frequency domain to identify embedded informa-tion using proper linear or nonlinear computainforma-tional algorithms. For example, Heart-Rate Variability (HRV) analysis can be per-formed in either the time domain, using the mean and standard deviation in R-R interval variations or the frequency domain, with reference to the power spectral density ratio of the low-frequency 1350-4533/$ – see front matter © 2010 IPEM. Published by Elsevier Ltd. All rights reserved.
and high-frequency bands[7]. This scheme can be used to analyze sympathetic and parasympathetic neuronal amplitude modula-tion. VF detection is performed using time-domain methods to extract morphological features, such as Threshold Crossing Inter-vals[8], Complexity Measure[9], Phase Space Reconstruction[10] and VF Filter[11]algorithms, and frequency-domain algorithms, such as Spectral Analysis[12], Generalized Discriminant Analysis and Support Vector Machine[13]and High-Order Spectral Tech-nique algorithms[14], which are applied to the power spectrum density. Gang et al. undertook the VT and VF by multifractal anal-ysis [15]. Lu et al. used R-R intervals with probability density function for detecting atrial fibrillation[16]. Recently, many nonlin-ear methods, based on the Multilayer Perpectron Neural Network, Adaptive Neuro-Fuzzy Interference System and the Mixture of Expert algorithms, have exploited Artificial Neural Network (ANN) architecture to distinguish among various cardiac diseases[17]. For example, Kannathal et al. used adaptive neuro-fuzzy technique for cardiac state diagnosis[18]. Sarkar et al. employed the Lorenz plot ofRR intervals to classify AT/AF signals[19]. Lemay et al. adopted phase-rectified signal averaging (PRSA) to estimate the dominant frequency during AF[20]. Climent et al. utilized the Poincare Surface Profile (PSP) to detect automatically and quantify preferential clus-ters of RR intervals[21]. Amann et al. summarized the performance of 10 VF algorithms for AED applications[22]. Übeyli used the multilayer perceptron neural network architectures[23], adaptive neuro-fuzzy inference system model[24]and mixture of experts network structure[25]for classification of ECG beats.Table 1 sum-marizes the performance of above-mentioned algorithms for the quantitative detection of cardiac arrhythmia.
From this table, most of the frequency-domain algorithms require complex Fourier transform (FT) computational processes, and the ANN methods need off-line training and extensive mem-ory resources for feature extraction and prediction, potentially limiting their ranges of real-time applications. However, either sen-sitivity or specificity of time-domain methods is lower than the
other methods. Among the time-domain methods, the phase space reconstruction (PSR) method reports best real-time VF/normal dis-crimination. It is thus rather interesting to further explore its feasibility in the AF/VF classification. All the reported phase space analysis algorithms follow Takens’ theorem[26]to reconstruct a multi-dimensional phase space matrix (PSM) from a digitized ECG signal[10,27–30]. The coordination of the projected points on the PSM is determined with proper time delay constant and varies with amplitude, time and periodicity of the acquired signals. Amann et al. calculated the filled area to classify the VF and non-VF signals with improved performance[10]. Richter et al. discussed the pos-sibility of using Takens’ method for ECG fluctuation[27]. Cvikl et al. used plane geometry equation for a planar non-self-intersecting polygon area calculation for QRS complex detection[28]. Ravelli et al.[29]and Fojt et al.[30]utilized the correlation dimension to express complexity of reconstructed PSM under different cardiac rhythms. However, they did not address how to differentiate AF and VF in their article.
Atrial fibrillation (AF) and ventricular fibrillation (VF) are two major cardiac arrhythmias that may occur during the events in which the rhythm of the heart is disordered, and they have a serious impact on public health. The prevalence of atrial fibrilla-tion increases as age increase, patients with AF are at higher risk for future systemic embolic events, including ischemic stroke and ischemic bowel disease. Ventricular fibrillation is a life-threatening disease that may cause sudden cardiac death and emergent elec-trical cardioversion and resuscitation are mandatory. Therefore, to carefully identify and discriminate the differences of atrial fib-rillation and ventricular fibfib-rillation from normal sinus rhythm are important. Cardiologists can easily identify the differences between AF and VF using a standardized recording file from an ECG or a Holter machine. As technology has improved, more tele-machines or implantable devices have become available for detecting many physiological signals. Current algorithms for heart-rate analyses using ECG signals may make errors in discriminating
Table 1
Cardiac arrhythmia detection algorithms for classification.
Algorithm Description Performance (%)
Threshold Crossing Intervals (TCI) The TCI algorithm sets the threshold value (20% of the maximum value with each 1-s segment) to the digitized ECG signal, and analyzes the pulse count and the average time-interval for signal crossing[8].
SEV= 75.1
SPV= 84.4
Complexity Measure (CPLX) The CPLX algorithm selects an appropriate threshold value to translate the ECG signal to a binary string, and obtains its asymptotic behavior[9].
SEV= 59.2
SPV= 92.0
Phase Space Reconstruction (PSR) The PSR algorithm translates the ECG signal to the multi-dimensional phase space matrix, and analyzes the features of geometrical orbits[10].
SEV= 79.0
SPV= 97.8
VF Filter (VF) The VF algorithm applies a narrow band-stop filter with the mean frequency of the relevant ECG signal, and computes the “VF filter leakage” to identify the VF segment[11].
SEV= 18.8
SPV= 100
Spectral (SPEC) The SPEC algorithm applies fast Fourier transform (FFT) to translate the ECG signal to the frequency domain, and analyzes the energy content of various frequency bands[12].
SEV= 29.1
SPV= 99.9
Generalized Discriminant Analysis and Support Vector Machine
(GDA + SVM)
The GDA algorithm maps ECG signals into feature space and performs Linear Discriminant Analysis (LDA) algorithm to replace original space to select discriminating features. The SVM provides effective classification by machine-learning technique[13].
SEA= 94.6 SPA= 99.7 SEV= 100 SPV= 100 SEVA= 100 SPVA= 100
High-Order Spectral Technique (HOST) The HOST algorithm utilizes the ability of high-order statistic techniques (bispectrum) to analyze the nonlinear activities (phase coupling) of ECG rhythm[14].
SEA= 83.3
SPA= 100
SEV= 91.7
SPV= 83.3
Multilayer Perpectron Neural Network (MLPNN)
The MLPNN uses the Lyapunov exponents, wavelet coefficients and power spectral density as inputs and trained with Levenberg-Marquardt algorithm to detect variability[23].
SEA= 95.6
SPA= 98.8
Adaptive Neuro-Fuzzy Interference System (ANFIS)
The ANFIS trains classifier with back-propagation gradient descent method in modeling nonlinear function by ECG signals and adjusted by the error criterion[24].
SEA= 96.7
SPA= 98.8
Mixture of Expert (ME) The ME network structure uses statistical features that decomposed by discrete wavelet transform from ECG signals as input and implemented for ECG beats classification[25].
SEA= 95.6
SPA= 96.7
SEAand SPArepresent sensitivity and specificity in AF/normal discrimination, while SEVand SPVdenote sensitivity and specificity in VF/normal discrimination. SEVAand SPVA
represent sensitivity and specificity in VF/AF discrimination. The sensitivity was defined with the number of true positives divided by the sum of the number of true positive and false negative. The specificity was defined with the number of true negative divided by the sum of the number of true negative and false positive.
between AF and VF when the patient is in extremely tachycar-dia [4]. If the interpretation is not accurate, the outcome may be disastrous, as inappropriate therapies may be delivered by Implantable Cardioverter Defibrillators (ICDs). Irregular periodic-ity, with the different morphologies associated with AF and VF, can be observed in ECG waveforms. However, improving discrim-ination by processing such signals remains a challenging task for tele-healthcare.
By carefully looking into the different PSM patterns, there are minor differences between cases of interest, e.g. AF and VF, espe-cially in the low visited pixel area. In order to take into account such a minor difference, we propose to accumulate the number of visits of each array element in the PSM and then calculate the difference between the subsequent PSM and the reference template. We also apply self-normalization scheme to minimize the inter-variability among different subjects. The normalized values are then tested for classification hypothesis by selecting proper threshold value for its optimal performance with ROC analysis.
2. Materials and methods
This study presents a time-domain analytical algorithm, called “Complex Phase Space Difference (CPSD)”, which correlates the spatial changes with the temporal changes in the phase space matrix patterns reconstructed with various annotated ECG signals from publicly available databases and clinically measured data, to index ECG signals, especially AF/VF, in real time.
2.1. Preprocessing and grouping of ECG signals
Since the CPSD algorithm focuses on the spatial changes between the projection of PSM from the time-domain waveforms, the algorithm must allow few artifacts (such as the baseline drifting, intra-subject variability and AC power-line noise, among others), potentially influencing the accuracy of the classification. In most cases, digital filters can be utilized to reduce or even eliminate the artifacts – especially baseline drifting and power-line noise. The intra-subject variability is minimized by normalization and a statistical approach. All of the ECG signals were filtered using a digital band-pass and notch filters in LabVIEW 8.6 (National Instrument Inc., USA) with a fourth-order Butterworth digital fil-ter function with designated cut-off frequencies of 1–100 and 60 Hz respectively. The high-pass filter (HPF) was employed to eliminate baseline drift, because of the contact between skin and electrode.
In this study, the ECG arrhythmia databases of MIT-BIH[31,32], CU[33,34]and Sudden Cardiac Death Holter (SCDH)[32,35]were used to evaluate the performance of the proposed CPSD algo-rithm. We also used the ECG records of outpatient-visit of internal medicine department at the E-DA hospital (IRB number EMRP-098-005, I-Shou University, Kaohsiung, Taiwan) for verification. The total number of subjects was 45, with demography data of gender (male: 28; female: 17) and age (57.9± 15.5 years). All the recruited subjects were clinically diagnosed and confirmed by two cardi-ologists to cover the range from NSR (n = 36) to low risks (n = 9) of AF, premature ventricular contraction (PVC), atrial premature contraction (APC), right bundle branch block (RBBB), paroxysmal supraventricular tachycardia (PSVT), respectively.Table 2 summa-rizes the annotated ECG signals that were applied in this work. We used 1170 randomly selected ECG segments from MIT-BIH database to determine the optimal working parameters and then evaluated the CPSD performance with the rest of the 4236 segments listed in Table 2.
Table 2
Statistical summary of ECG segments used in evaluation of CPSD algorithm.
ECG signal source Segments Total data length (min) MIT-BIH Database Normal 1170
a 19:30
AF 1186 19:46
CU Database VF 219 3:39
SCDH Database VF 986 16:26
Clinical Measured NormalAF 1664181 27:443:01
aThese 1170 segments were randomly selected from the MIT-BIH database and
used to determin the optimal parameters in the phase space reconstruction process.
2.2. Definition of working parameters
The performance of proposed CPSD depends on the optimal selection of working parameters with well-defined objective func-tions and theoretical assumpfunc-tions. The proposed CPSD algorithm includes four important parameters, i.e. the appropriate data length (L = N× t, where N is the total number of sampled signal and t is the sampling interval), time delay constant (T, where T = k× t, k is the index of T), the matrix dimension (d) and the matrix size (MS) of the PSM. CPSD transforms the one-dimensional time series signal into a d-dimension matrix to improve the classification. A com-prehensive discussion of these theorems and selection methods can be found in the references[36–39]. Briefly, this CPSD method begins with a segment of a time series of digitalized ECG signal, s(n), with serial index n. A multi-dimensional matrix, S(n) = (s(n), s(n + k), s(n + 2k), . . ., s(n + dk)), up to d dimensions can then be con-structed from s(n) with a suitable delay index, k. The optimal values of d and k are then determined by minimizing the Euclidian dis-tance between nearest neighbors in the False Nearest Neighbor (FNN) method[37,38]and the joint probability density function in the Average Mutual Information (AMI) method [39], respec-tively. More details of the derivations are presented in the following sections.
2.2.1. Time delay constant T or delay index k
The phase space matrix, S(n), can have its maximal classification with d-dimensional bases of orthogonal or independent time serial signals, s(n), s(n + k), s(n + 2k), . . ., s(n + dk). The degree of mutual dependence, Mutual Information (MI), of any two sequences can be calculated by their probability density functions (pdf). For the sim-plest case of a two-dimensional S(n), its MI value can be given by (1), where P(s(n), s(n + k)) denotes the joint pdf of s(n) and s(n + k), and P(s(n)) and P(s(n + k)) represent the pdfs of s(n) and s(n + k), respectively. If s(n) and s(n + k) are completely independent of each other, then P(s(n), s(n + k)) = P(s(n))× P(s(n + k)), such that the MI value is zero. The AMI value, I(k), is calculated by(2)and gives the mean amount of dependencies of s(n) and s(n + k) for various delay indexes, k. The minimal value of I(k) is corresponding to the highest independency of two bases and is used to determine the optimal delay index. MutualInformation = log2 P(s(n), s(n + k)) P(s(n)) × P(s(n + k)) (1) I(k) =
s(n),s(n+k) P(s(n), s(n + k)) × MutualInformation (2)2.2.2. Dimension d and data length L
The number of dimensions can often be increased, such that d→ d + 1, to distinguish between overlapping functions in lower dimensions, d. However, this process depends on the optimization using quantitative measures. The Euclidian distance, R, between
any two points in a space can be used as an index to evaluate the separation or classification capability under different dimensions of d and d + 1 as given by(3) and (4) [11]. The FNN value is cal-culated by(5)and compared with a threshold, ε, to evaluate the degree of unfolding that gives the optimal d value(6). Generally, the threshold was advised to be 10 in the previous related stud-ies[18,37]. When the computed Euclidian distance is less than the threshold value ε, the state is set to “False” and the FNN is set to one (the initial FNN value is set to zero). Finally, accumulated FNN is divided by the total number (N) of acquired signal length (L) to yield the FNN(%) by(7). R2d(n, r) = d−1
m=0 [s(n + mk) − s(r)(n + mk)]2 (3) R2 d+1(n, r) = R2d(n, r) + [s(n + dk) − s(r)(n + dk)] 2 (4) R2 d+1(n, r) − R 2 d(n, r) R2 d(n, r) 1/2 = |s(n + dk) − s(r)(n + dk)| Rd(n, r) (5) FNN =⎧
⎪
⎨
⎪
⎩
s(n + dk) − s(r)(n + dk)
Rd(n, r) > ε, FNN = FNN
s(n + dk) − s(r)(n + dk)
Rd(n, r) ≤ ε, FNN = FNN + 1
⎫
⎪
⎬
⎪
⎭
, for n = 1 to N (6) FNN(%) =FNN N × 100% (7)Generally, lower FNN(%) means larger unfolding effect. Thus, one can use the minimal FNN(%) value to optimally determine the dimensions of PSM and the data length of each ECG segment for the CPSD.
2.2.3. Matrix size, MS
Given the optimal matrix dimension determined by above-mentioned FNN method, the size of the matrix, MS, is another important parameter which will affect the performance of CPSD algorithm and subject to further elucidation. It will not only deter-mine the required physical memory space and computational time to calculate the final CPSD value but also affect the mapping res-olution between the analyzed ECG segment and its corresponding location on the PSM. The larger MS value indicated better resolu-tion of the ECG amplitude and yet required larger computaresolu-tional resources. Therefore, it is a trade-off to determine the optimal MS value for practical implementation. The mapping process starts with a maximal normalization of the absolute value of the analyzed signal segment,|s(n)|, and then shifts to the center of the PSM by multiplying a scaling factor of one half of MS to account for the bipolar output below baseline, which is given in(8).
sN(n) =12
s(n) max{|s(n)|}+ 1 × MS (8)
where sN(n) is the normalized s(n).
2.3. Reconstruction of phase space matrix
Based on the above-mentioned method, a two-dimensional PSM, SM, can be constructed from the normalized ECG signal
seg-ment, sN(n), with delay index, k, according to (9). The number
of visited-times on the SMwith coordinate of (sN(n), sN(n + k)) is
counted by(10)from 0 to N− k.
SM= (sN(n), sN(n + k)) (9)
SM[sN(nx), sN(nx+ k)] = SM[sN(nx), sN(nx+ k)] + 1,
for nx= 0 to N − k (10)
where SM[ ] represents the intensity matrix of accumulated
visited-number on the SM.
Following the presentation of the basic method for construct-ing PSM, the constructed reference template, Reference Phase Space Matrix (RPSM), the experimental segments, and Experimen-tal Phase Space Matrix (EPSM) are now used to calculate the CPSD value in the following section.
2.4. Compute CPSD value
Because the proposed CPSD algorithm focuses on the spatial-temporal changes of PSM, we use the absolute values of the differential matrix, DPSM, between EPSM and RPSM to calculate the complexity value (CV) according to(11) and (12).
SDPSM[x, y] = |SEPSM[x, y] − SRPSM[x, y]|, for x, y = 0 to MS
(11)
The complexity value (CV) of DPSM is then defined as the num-ber of all array elements that exceed zero.
CV =
CV + 1, SDPSM[x, y] /= 0 CV, otherwise , for x, y = 0 to MS (12) The first CV of each ECG segment is adopted as a normaliza-tion factor in determining the subsequent CV to obtain the final CPSD values(13)that minimize the possible intra- and inter-patient variations. Self-normalization of each individual’s normal ECG pat-terns enables the establishment of a normal range of baseline CPSD values. We used 1170 randomly selected NSR segments from the MIT-BIH database to determine the effects of inter-subject vari-ability on ECG signals. This hypothesis is tested by curve fitting Gaussian normal distribution function to CPSD values.CPSDn= CVn
CV1,
where n ≥ 1 (13)
2.5. Classification flowchart
The important concepts and the computational algorithms described in the previous sections can now be combined to evaluate the remaining 4236 ECG signals.Fig. 1presents the comprehen-sive flowchart for the classification of ECG by the CPSD algorithm, which comprises two important phases. The first phase is to choose the RPSM and EPSM of NSR to obtain the first CV value as a self-normalization factor. The second phase is the classification of AF and VF using validated thresholds. The moving window scheme is applied to extract an ECG segment with a fixed data length (L). In each moving step (1 s), the CPSD algorithm constructs a new EPSM and calculates its corresponding CPSD value according to(13). The first L + 1 s of all recorded ECG signals is employed to construct one RPSM (0− L s) and one EPSM (1 − (L + 1) s). Other than the CV value, we also used the PSR method[10]to calculate a reference value, DA, which is given by the total number of visited pixels divided
by the total number of pixels of PSM, at the same time to ensure the classification can start with normal RPSM and EPSM. The cor-responding RPSM/EPSM is determined to be normal if the DAfalls
within the published normal range. When these two values meet required conditions, the CPSD algorithm moves to second phase to analyze the subsequent ECG signals. Otherwise, it will return to the first step and load-in another segment for continuous evaluation. The threshold values (TH1and TH2) are utilized to classify AF and
Fig. 1. Flowchart of CPSD algorithm for classification of normal, AF, and VF ECG signals under properly selected conditions on TH1and TH2for AF and VF.
results from the remaining 4236 annotated databases and our IRB data.
3. Results and discussion
Before showing the analysis results by the above-mentioned method and flowchart, we will try to illustrate the example PSM distributions and their discrimination by traditional method and proposed CPSD algorithm with three typical ECG segments of NSR, AF, and VF as shown in theFig. 2. As shown in theFig. 2(a)–(c), three distinct patterns result from the ECG waveforms of increasing level of risks on reconstructed phase space. According to PSR approach [10], the numbers of visited pixels in these three cases are 141, 212, and 328, respectively. Notably, many visits are made to the same pixel because of the data segment exhibits multiple beats with peri-odicity, especially in normal ECG signals. To account for this fact, the accumulated numbers of visited pixels are plotted as an inten-sity map on the PSM and calculated the difference between two PSMs. This approach further amplifies the differences between the
highly ordered distributed pattern of a normal ECG segment and an irregular one.Fig. 2(d) shows the results of a DPSM under nor-mal condition and its reduced total number of visited pixels now to 76. However, in AF and VF cases, the increases in irregularity widen the spreading of the visited patterns as shown inFig. 2(e) and (f). The total numbers of visited pixels for AF and VF are 212 and 348, respectively. This example demonstrates that CPSD can improve the discriminatory performance by simple subtraction of EPSM from RPSM.Table 3shows the information of selected ECG segments used inFig. 2.
In this work, we used 1170 data segments from the MIT-BIH arrhythmia database to determine the optimal working parameters and used the remaining 4236 ECG segments to classify the AF and VF with proper thresholds of CPSD values.
3.1. Optimal working parameters
For the determination of time delay constant in the interval of 0–0.9 s, the minimum value of I(k) is found in between 0.1 and 0.4 s. Since the time delay constant should be shorter than the
Fig. 2. Reconstructed phase space matrices of (a) normal, (b) AF, and (c) VF cases from database based on traditional Takens’ theory, which indicates that the irregularity of
ECG signals increases complexity of the distribution pattern. The absolute difference of EPSM and RPSM in the CPSD algorithm yield the resultant differential patterns for (d) normal, (e) AF, and (f) VF cases.
Table 3
The information of ECG segments used inFig. 2.
ECG signals Record Time portion (min) Number of beats Normal MIT-BIH database
(116.dat) A: 00:01 to 00:08 9 R:00:00 to 00:07 9 AF MIT-BIH database (201.dat) A: 08:05 to 08:12 12 R: 06:52 to 06:59 5 VF CU database (cu17.dat) A: 06:23 to 06:30 28 R: 06:12 to 06:19 10 The A and R represents the analyzed and reference ECG segment used inFig. 2, respectively.
interval between heart beats, the time delay constant T is set to 0.2 s, which corresponds to k = 72 at a sampling rate of 360 Hz. If T exceeds 0.5 s, the variation in I(k) will be larger than otherwise because of overlapping with the subsequent ECG cycle. Since a nor-mal heart rhythm can sometimes reach as high as 180 bpm, one can reasonably select T = 0.2, which is well below the average R-R inter-vals of normal ECG segments, without worry about the problem of overlap.
From the results of FNN(%) calculation, it shows that increasing the dimension, especially from one to two, substantially reducing FNN(%). Although d = 3 yields a lower FNN(%) than d = 2 and sup-posedly has better unfolding capability, it is under the price of exponential increases in the matrix size and computational loads. Therefore, the CPSD algorithm was implemented with d = 2 for the reconstructed PSM in the subsequent calculations. This method is also applied to determine the optimal data length for s(n) with d = 2. Generally, for practical applications of Takens’ method, the ECG seg-ment with a selected data length must include at least one heart beat. According to our experimental results, the data length used in the CPSD algorithm is selected based on a local minimum of FNN(%) at 7 s.
Finally, the optimal matrix size is set to be 40× 40, which gives enough amplitude resolution with reasonable computational resources for real-time requirement comparing to size of 10, 20, and 80.
3.2. Evaluation of CPSD performance
Given the procedures for executing the proposed CPSD algo-rithm and optimally selecting the working parameters, a good reference template must be found to begin the computations. Even
Fig. 3. Histogram of CV values calculated from NSR segments. It is best fitted using
a Gaussian curve with its center at 46.2 and an STD of 8.3 (R2= 0.91).
though it is a kind of “template matching” method, the CPSD algo-rithm uses a 7 s segment of an ECG signal, which typically has 5–11 QRS complexes under normal conditions. The computation of the degree of complexity does not depend on an exact match between each individual template to classify the signals. The degree of complexity within 7 s provides sufficient statistical variances to avoid the “template matching” dilemma. However, the criteria for selecting the reference segment must be justified statistically. Fig. 3displays a histogram of CV values calculated from 1170 ran-domly selected normal ECG segments. It can be best fitted by using a Gaussian distribution curve with its center (mean value) at 46.2 and a standard deviation (STD) of 8.3 (R2= 0.91). We also
calcu-lated the DAvalues by using PSR method and found the values was
0.15± 0.043 (R2= 0.93), which was highly agreed with published
data.
We then used the rest of 4236 the annotated segments inTable 2 to calculate the CPSD values and to validate the proper threshold values in NSR, AF and VF cases. Various irregular ECG waveforms represent higher risks and thus correspond to higher CPSD values. Fig. 4shows the histogram of the resultant CPSD values for NSR, AF and VF. By curve fitting the Gaussian distribution function, we could identify the CPSD values as 1.48± 0.21 for NSR, 3.71 ± 0.99 for AF, and 9.38± 2.22 for VF, respectively. The ECG morphology and temporal variations strongly affected the complexity of the
Fig. 4. Histogram of normal, AF, and VF CPSD values, obtained from 4236 collected ECG segments, and separated by setting thresholds of 2.0 and 6.0 for AF/normal and VF/AF
Fig. 5. The ROC curves of AF/normal, VF/normal and VF/AF, respectively. The inset
was the enlarge area of dash circle.
reconstructed phase space. The improved algorithm can magnify the differences among these signals by taking the absolute value of DPSM. For normal cases, the rhythmic signals yield a “concentrated” pattern, a “narrow” distribution and low CPSD values. In the event of a disordered heart rhythm, as in AF, the waveforms exhibit some features of regular morphology (P wave blocked) but with various periodicities. This increase in complexity increases the CPSD values and broadens the distribution range from 2.0 to 6.0. Completely unpredictable VF signals are associated with irregular morphology and periodicity, which contribute to the highest CPSD values and the widest range.
The optimal threshold values for the classification of these three levels of risk were selected with ROC method.Fig. 5showed the receiver operating characteristic (ROC) curves calculated by changing the threshold values for AF/normal, VF/normal and VF/AF discrimination, respectively. According to this result, the optimal threshold for AF/normal and VF/normal was set to be TH1= 2.0 and
TH2= 6.0 for VF/AF.
To evaluate the performance of the proposed algorithm, the true positive (TP), true negative (TN), false positive (FP) and false negative (FN) were applied to determine the sensitivity (SE) and specificity (SP), respectively. All the classifications made by the CPSD algorithm are then compared to the annotated ECG databases and IRB data for the calculation of SE and SP.
Table 2lists the numbers of normal, AF and VF segments, which are 1664, 1367 and 1210, respectively.Table 4presented the clas-sification results for normal, AF and VF cases with a threshold of two to distinguish abnormal ECG signals (AF and VF) from normal signals and a threshold of 6.0 to distinguish VF signals from AF sig-nals. In both AF and VF cases, all of the signals were accurately identified as true positive, yielding 100% sensitivity. However, 20 normal cases were classified as false positives, resulting in slightly lower specificity of 98.8%. The CPSD algorithm herein was adopted to classify these two abnormal ECG signals using a threshold (TH2)
of 6.0. As presented in the third row ofTable 4, the resultant sen-sitivity was 97.9% and the specificity was 98.4%. A comparison
with the algorithms listed inTable 1reveals that this VF classifi-cation performance is only slightly worse than that of GDA + SVM, but the proposed method has benefit of easily implementation for real-time application. This algorithm also is the best for all cited algorithms in AF analysis.
Despite the successful demonstration of CPSD performance with above-mentioned results, it is also critical to consider several common problems of intra-subject variability for practical clini-cal applications. Among these issues, baseline wandering, rhythmic and QRS amplitude changes are often seen in normal ECG signals. These may cause by the movement artifacts, changes in electrode impedance and movement of the cardiac axis caused by respira-tion. Some of these phenomena have been investigated to ensure the applicability of the CPSD algorithm under these conditions. Fig. 6(a) shows a trace of ECG waveforms with accelerating heart rate from 80 bpm to almost 180 bpm. Its resultant CPSD values (in Fig. 6(d) with open circles) increase from 1.1 to 1.7, which is still well below 2.0 and remains in the NSR group. It is rather interest-ing to further explore the applicability of this algorithm for excise ECG in the future.Fig. 6(b) shows a stable baseline with chang-ing amplitude of QRS complex and its CPSD values in theFig. 6 (d) (filled triangle) remains in the range from 1.1 to 1.5 with ref-erence to a normal range of 1.1–1.4 (filled circle).Fig. 6(c) shows a noise trace of ECG waveforms with QRS amplitude changes and baseline drifting. It is rather interesting to note that the resultant CPSD values remain in the normal range, too. These three com-mon ECG variations have minor effects on the CPSD distribution toward higher values as shown in Fig. 6(d). However, all these changes are well below the threshold value of 2.0 and so remain “normal” with acceptable effects. Again, many of these interfer-ences can be further reduced by using advanced noise reduction methods, including time-varying digital filter[40], adaptive noise reduction [41], empirical mode decomposition [42] or wavelet transform-related algorithms [43]for better input data to begin with.
Even though, our CPSD method can be used in the above-mentioned situations, we did come across one particular interesting case under the extreme condition of VF that might have different interpretation results. As revealed by the previous VF example, the progress of effortless contraction degraded the cardiac muscles and was therefore accompanied by a reduction in ampli-tude. This would in turn reduce the corresponding CPSD values, resulting in so-called “VF attenuation”. Once heart entered the VF stage, it was unlikely to return to the normal (or AF) state without aid such as drug delivery or electrical shock. These reduced CPSD values were identified as VF cases and should be ignored during statistical analysis.
A continuous ECG signal was evaluated by the CPSD algorithm using an extracted VF segment from the record of CU database (cu22.dat).Fig. 7presents the development of continuous ECG sig-nals into VF. The first 8 s of ECGs were initially used to calculate RPSM and EPSM. The first calculated CV value, CV1= 53, satisfied the
specified criteria for the first step in the CPSD algorithm, enabling the CPSD values to be calculated in the next stage. With a step size of 1 s and a moving window with a length of 7 s, the corresponding complexity and CPSD values of the normal segments before 13 s are 1.05, 1.21, 1.21 and 1.35 (Fig. 7(b)), respectively. The CPSD values begin to increase at 12 s because of the presence of VF.
Table 4
Results of analysis of ECG signals (unit: segments).
TP FN TN FP SE (%) SP (%) TH
AF/Normal 1367 0 1644 20 100 98.8 2.0
VF/Normal 1210 0 1644 20 100 98.8 2.0
Fig. 6. Effects of intra-variability on CPSD results for three common normal ECG signals. (a) Accelerating heart-rate signal, (b) change in QRS possibly caused by respiration, (c) baseline drifting with change in QRS for normal ECG signals, and (d) corresponding CPSD histograms in these three cases, with comparison to normal ECG signals.
Fig. 7. Continuous evaluation of VF signals by CPSD algorithm yields updated every second. (a) Normalized ECG signals became VF with envelope modulation and stimulation
artifact. (b) CPSD values that successfully identified VF in 2 s and remained above the VF threshold with minor effects of both artifacts.
4. Conclusions
The proposed CPSD algorithm uses a two-dimensional recon-structed phase space matrix (PSM)[26–30]from each time segment of a digitized ECG signal to compute changes in spatial distribution. Like the published method by Amann et al.[10], two-dimensional delayed-time coordinate is used to construct the PSM and to trace out the projection of time sequence from ECG signals. The dis-tributed pattern of the reconstructed phase space varies with the amplitude, time and periodicity of the acquired signals. Instead of using the trail of projected orbit[28] and concentrated posi-tion [29,30] to distinguish ECG signals, the PSR algorithm[10] calculates the filled area to classify the VF and non-VF signals with improved performance. However, all these methods did not address how to differentiate between AF and VF in this article. By carefully looking into the different PSM patterns, there are minor differences between cases of interest, e.g. AF and VF, espe-cially in the low visited pixel area. In order to take into account such a minor difference, we propose to accumulate the number of visits to each array element in the PSM and then calculate the difference between the subsequent PSM and the reference tem-plate. We also apply self-normalization scheme to minimize the inter-variability among different subjects[44]. The obtained nor-malized CPSD values are then tested for classification hypothesis by selecting proper threshold value for its optimal performance with ROC analysis. In particular, it is the first real-time algo-rithm that reports the best performance to distinguish AF and VF with sensitivity of 97.9% and specificity of 98.4%. According to the experiences, the intra-variability will increase “broader” dis-tribution but still remain in normal range, while the calculated CPSD value of “VF attenuation” phenomena should be ignored. In this proposed algorithm, the computational scheme only requires matrices addition and subtraction, and thus significantly reduces the complexity for real-time implementation and is able to adopt in different scenarios of tele-healthcare and implantable applica-tions.
Acknowledgments
The authors would like to thank the Ministry of Economic Affairs of the Republic of China, Taiwan, for financially supporting this research under contract no. MOEA 96-EC-17-A-05-S1-017. Ted Knoy is appreciated for his editorial assistance.
Appendix A. Supplementary data
Supplementary data associated with this article can be found, in the online version, atdoi:10.1016/j.medengphy.2010.04.001.
Conflict of interest
The authors declare that they do not have any conflict of interest.
References
[1] Abidi SSR. An intelligent tele-healthcare environment offering person-centric and wellness-maintenance services. J Med Syst 2001;25(3):147–65. [2] Callaway CW, Sherman LD, Mosesso VN, Dietrich JTJ, Holt E, Clarkson MC.
Scaling exponent predicts defibrillation success for out-of-hospital ventricular fibrillation cardiac arrest. Circulation 2001;103:1656–61.
[3] Gillberg J. Detection of cardiac tachyarrhythmias in implantable devices. J Elec-trocardiol 2007;40:S123–8.
[4] Israel CW. How to avoid inappropriate therapy. Curr Opin Cardiol 2008;23:65–71.
[5] Husser D, Stridh M, Sörnmo L, Toepffer I, Klein HU, Olsson SB, et al. Elec-troatriography – time–frequency analysis of atrial fibrillation from modified 12-lead ECG configurations for improved diagnosis and therapy. Med Hypothe-ses 2007;68:568–73.
[6] Karrakchou M, Vibe-Rheymer K, Vesin J-M, Pruvot E, Kunt M. Improving car-diovascular monitoring through modern techniques. IEEE Eng Med Biol Mag 1996;15(5):68–78.
[7] Malik M, Bigger JT, Camm AJ, Kleiger RE, Malliani A, Moss AJ, et al. Heart rate variability: standards of measurement, physiological interpretation, and clini-cal use. Eur Heart J 1996;17(3):354–81.
[8] Thakor NV, Zhu Y-S, Pan K-Y. Ventricular tachycardia and fibrillation detec-tion by a sequential hypothesis testing algorithm. IEEE Trans Biomed Eng 1990;37(9):837–43.
[9] Zhang X-S, Zhu Y-S, Thakor NV, Wang Z-Z. Detecting ventricular tachy-cardia and fibrillation by complexity measure. IEEE Trans Biomed Eng 1999;46(5):548–55.
[10] Amann A, Tratnig R, Unterkofler K. Detecting ventricular fibrillation by time-delay methods. IEEE Trans Biomed Eng 2007;54(1):174–7.
[11] Abarbanel HDI, Kennel MB. Local false nearest neighbors and dynamical dimen-sions from observed chaotic data. Phys Rev E 1993;47(5):3057–68.
[12] S B, R R, D C, J M. Algorithm sequential decision-making in the frequency domain for life threatening ventricular arrhythmias and imitative artifacts: a diagnostic system. J Biomed Eng 1989;11(4):320–8.
[13] Asl BM, Setarehdan SK, Mohebbi M. Support vector machine-based arrhythmia classification using reduced features of heart rate variability signal. Artif Intell Med 2008;44:51–64.
[14] Khadra L, Al-Fahoum AS, Binajjaj S. A quantitative analysis approach for cardiac arrhythmia classification using higher order spectral techniques. IEEE Trans Biomed Eng 2005;52(11):1840–5.
[15] Wang G, Huang H, Xie H, Wang Z, Hu X. Multifractal analysis of ventricular fibrillation and ventricular tachycardia. Med Eng Phys 2007;29:375–9. [16] Hong-Wei L, Ying S, Min L, Pi-Ding L, Zheng Z. A probability density
func-tion method for detecting atrial fibrillafunc-tion using R-R intervals. Med Eng Phys 2009;31:116–23.
[17] Acharya UR, Sankaranarayanan M, Nayak J, Xiang C, Tamura T. Automatic iden-tification of cardiac health using modeling techniques: a comparative study. Inf Sci 2008;178:4571–82.
[18] Kannathal N, Lim CM, Acharya UR, Sadasivan PK. Cardiac state diagnosis using adaptive neuro-fuzzy technique. Med Eng Phys 2006;28:809–15.
[19] Sarkar S, Ritscher D, Mehra R. A detector for a chronic implantable atrial tach-yarrhythmia monitor. IEEE Trans Biomed Eng 2008;55(3):1219–24. [20] Lemay M, Prudat Y, Jacquemet V, Vesin J-M. Phase-rectified signal
averag-ing used to estimate the dominant frequencies in ECG signals duraverag-ing atrial fibrillation. IEEE Trans Biomed Eng 2008;55(11):2538–47.
[21] Climent AM, Guillem MıdlS, Husser D, Castells F, Millet J, Bollmann A. Poinca´re surface profiles of RR intervals: a novel noninvasive method for the evaluation of preferential AV nodal conduction during atrial fibrillation. IEEE Trans Biomed Eng 2009;56(2):433–42.
[22] Amann A, Tratnig R, Unterkofler K. Reliability of old and new ventricular fibril-lation detection algorithms for automated external defibrillators. Biomed Eng OnLine 2005;4(60):15.
[23] Übeyli ED. Statistics over features of ECG signals. Expert Syst Appl 2009;36:8758–67.
[24] Übeyli ED. Adaptive neuro-fuzzy inference system for classification of ECG signals using Lyapunov exponents. Comput Meth Programs Biomed 2009;93:313–21.
[25] Übeyli ED. Implementing wavelet transform/mixture of experts network for analysis of electrocardiogram beats. Expert Syst 2008;25(2):150–62.
[26] Takens F. Detecting strange attractors in turbulence. In: Lecture notes in math-ematics. Springer Berlin/Heidelberg, 1981, 10.1007/BFb0091903.
[27] Richter M, Schreiber T. Phase space embedding of electrocardiograms. Phys Rev E 1998;58(5):6392–9.
[28] Cvikl M, Jager F, Zemva A. Hardware implementation of a modified delay-coordinate mapping-based QRS complex detection algorithm. EURASIP J Adv Signal Process 2007;2007(1):104.
[29] Ravelli F, Antolini R. Complex dynamics underlying the human electrocardio-gram. Biol Cybern 1992;67:57–65.
[30] Fojt O, Holcik J. Applying nonlinear dynamics to ECG signal processing. IEEE Eng Med Biol 1998;17:96–101.
[31] Massachusetts Institute of Technology. MIT-BIH arrhythmia database [Online]. Available:http://www.physionet.org/physiobank/database/mitdb.
[32] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, et al. PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation 2000;101(23):e215–20. [33] Massachusetts Institute of Technology. CU database [Online]. Available:
http://www.physionet.org/physiobank/database/cudb.
[34] Nolle FM, Badura FK, Catlett JM, Bowser RW, Sketch MH. CREI-GARD, a new concept in computerized arrhythmia monitoring systems. Comput Cardiol 1986;13:515–8.
[35] Massachusetts Institute of Technology. SCDH database [Online].
http://www.physionet.org/pn3/sddb.
[36] Porta A, Rienzo MD, Wessel N, Kurths J. Addressing the complexity of cardio-vascular regulation. Philos Trans R Soc A: Math Phys Eng Sci 2009;367:1215–8. [37] Hegger R, Kantz H. Improved false nearest neighbor method to detect
deter-minism in time series data. Phys Rev E 1999;60(4):4970–3.
[38] Rhodes C, Morari M. False-nearest-neighbors algorithm and noise-corrupted time series. Phys Rev E 1997;55(5):6162–70.
[39] Fraser AM, Swinney HL. Independent coordinates for strange attractors from mutual information. Phys Rev A 1986;33:1134–40.
[40] Sörnmo L. Time-varying digital filtering of ECG baseline wander. Med Biol Eng Comput 1993;51(5):503–8.
[41] Wu Y, Rangayyan RM, Zhou Y, Ng S-C. Filtering electrocardiographic signals using an unbiased and normalized adaptive noise reduction system. Med Eng Phys 2009;31:17–26.
[42] Blanco-Velascoa M, Wengb B, Barner KE. ECG signal denoising and baseline wander correction based on the empirical mode decomposition. Comput Biol Med 2008;38:1–13.
[43] Sayadi O, Shamsollahi MB. Model-based fiducial points extraction for baseline wandered electrocardiograms. IEEE Trans Biomed Eng 2008;55(1):347–51. [44] Liu CS, Lin YC, Chuang YH, Hsiao TC, Lin CW. Chaotic Phase Space Difference
(CPSD) algorithm for real-time detection of VF, VT and PVC ECG signals. In: 4th European conference of the international federation for medical and biological engineering. Springer Berlin/Heidelberg; 2009, 10.1007/978-3-540r-r89208-3.
