Electrocardiogram (ECG) is the P-QRS-T wave indicating the electrical activity of the heart. It originates at the sinoatrial (SA) node and the ECG signals are extracted by placing sensors at the limb extremities of the subject. The shape of the ECG signal indicates the cardiac health of the subject [1]. It is very difficult to decipher the subtle changes in these signals by our eye alone. These bio-signals are non-stationary and non-linear in nature.

Heart rate variability reflects the autonomic control of the cardiovascular system [2-3]. It is a simple, noninvasive technique which provides an indicator of the dynamic interaction and balance between the sympathetic nervous system and the parasympathetic nervous system. These signals are not linear in nature and hence, analysis using nonlinear methods can unveil the hidden information in the signal.

The HRV can be analyzed using (i) time domain analysis, (ii) frequency domain and (iii) non-linear methods. Time domain methods of HRV are simplest to use as they are based on common statistical measures. Time domain methods are the statistical methods and based on mean and standard deviation of the RR intervals [4].

A typical power spectrum of the heart rate signal has three main frequency regions [2]. The high frequency power (HF: 0.15-0.5 Hz) band reflects respiratory sinus arrhythmia (RSA) and cardiac vagal activity. Low frequency (LF: 0.04-0.15 Hz) power is related to baroreceptor control and is mediated by vagal and sympathetic systems. Very-low-frequency (VLF: 0.0033-0.04 Hz) power is related to the thermoregulatory, vascular mechanisms, and renin-angio tension systems.

Cardiac arrhythmia detection algorithms have been proposed using non-linear techniques [5-8]. Dingfie et al. have classified cardiac arrhythmia into six classes using autoregressive modeling [9].

The power spectrum provides the signal’s power (energy per unit time) falling within its frequency components [10-12]. But, it does not give any information about the phase relations between the frequency components. HOS can be used to analyze the nonlinear signals (VHF amplitudes), which may involve cross phase relations [13].

The Bispectrum is the Fourier Transform of the third order correlation of the signal which indicates cross correlation between frequency components in a two-dimensional (2-D) frequency plot [14-15]. It gives information about the phase coupling between the frequency components at f₁, f₂ , and f₁ + f₂ . The bicoherence is the normalized bispectrum.

Higher order spectra (HOS) are spectral representations of moments and cumulants and can be defined for deterministic signals and random processes. They have been used to detect deviations from Gaussianity and identify non-linear systems [14]. HOS based features can be formulated to be rotation, translation and scaling invariant when applied to one and two dimensional pattern recognition [16]. Research in the application of HOS includes detection of photomontage [17], mine detection [18], study of machine faults [19], recognition of viruses from electron microscopic images [20], analysis of biosignals like ECG [21] and EEG [22] etc.

The bicoherence approach was used to analyze the unexpected sudden very high frequency (VHF) peaks in the heart transplant subjects in [13]. They have indicated two types of peaks in this type of patients. The first type of VHF peaks reflects the cardiac-respiratory coupling, imposed by non-sinusoidal breathing. The second type was not induced by the respiratory system.

HOS have been widely used in the analysis of the EEG signals and recently, it has been applied to the HRV signal. The coupling between cardiac and respiratory activity was also studied using the bivariate data and cross-bispectrum computation between the ECG and respiratory signals [23]. They have absorbed the nonlinear coupling during both spontaneous and paced respiration. A technique for time-variant analysis of quadratic phase coupling (QPC) of the heart rate data was studied [24]. They have observed a phase co-ordination between the 10-second-rhythm and respiratory sinus arrhythmia (RSA) as well as a non-linear coupling (amplitude modulation) between these HRV components in neonatal heart rate signals.

Recently Chua et al., have analyzed the cardiac arrhythmia using heart rate signals using bispectrum invariant features and phase entropies [25]. They have proposed different bispectrum and bicoherence plots. In this paper, we derived features from the HOS and used them for classification. We were able to classify normal HRV signals and four other classes of arrhythmia with an average accuracy of above 85%.

The layout of the paper is as follows: Section (2) presents the data acquisition process and preprocessing of the raw cardiac signals. Section (3) deals with the HOS analysis namely the bispectrum, and the features of HOS, statistical analysis and SVM classifier used, Section (4) presents the results of the study. The discussion on our data analysis is presented in section (5). Finally the paper concludes in section (6).

The data used for analysis in this work was obtained by recording ECG signals (of 10 to 15 minutes duration) of over 330 patients in reclining condition. The permission from the clinicians and staff of the Kasturba Medical College Hospital, Manipal, India has been obtained to collect the data for this study. The original analog signal recorded in magnetic tapes (Holter recording) was subsequently converted to digital equivalent by discrete time sampling and quantization. A sampling rate of 320 samples/sec was chosen due to the fact that almost all the useful frequency components in the ECG signal were below 100 Hz. The numbers of dataset in each class is shown below in Table 1.

For the purpose of this study, the cardiac data that include both the normal and other cardiac abnormalities are classified into five categories or classes, namely,

1. Normal sinus Rhythm (NSR)
2. Pre-ventricular Contraction (PVC)
3. Complete Heart Block (CHB)
4. Sick Sinus Syndrome (SSS)
5. Congestive Heart Failure (CHF)

A brief description of the different cardiac classes is given below.

2.1. Preprocessing

The preprocessing of the ECG signals runs through the following steps.

• The data is passed through a low pass filter with a cut-off frequency of 35Hz to remove unwanted high frequencies present in the ECG signal.
• A high pass filter with cut-off frequency 0.3Hz is applied to suppress the baseline wander that is present in the signal.
• A band-stop filter of cut-off frequencies 50 or 60Hz is used to suppress the power-line interference noise.
• Finally a median filter is used to extract baseline wander of the processed ECG signal, and then it is subtracted from the processed ECG signal to effectively remove all the baseline wander.

In order to detect the R peaks of the ECG, Tompkins algorithm [26-27] is employed on the ECG data.

The interval between two successive QRS complexes is defined as the RR interval (t_r-r seconds) and the heart rate (beats per minute) is given as:

(1)

The HRV signal is analyzed using different higher order spectra (also known as polyspectra) that are spectral representations of higher order moments or cumulants of a signal. In particular, this paper studies features related to the third order statistics of the signal, namely the bispectrum. The Bispectrum is the Fourier Transform of the third order correlation of the signal and is given by

(2)

where X(f) is the Fourier transform of the signal x(nT) and E[.] stands for the expectation operation. In practice, the expectation operation is replaced by an estimate that is an average over an ensemble of realizations of a random signal. For deterministic signals, the relationship holds without an expectation operation with the third order correlation being a time-average. For deterministic sampled signals, X(f) is the discrete-time Fourier transform and in practice is computed as the discrete Fourier transform (DFT) at frequency samples using the FFT algorithm. The frequency f may be normalized by the Nyquist frequency to be between 0 and 1.

3.1. Higher Order Spectral Features

In our earlier study, we proposed the general patterns for different classes of arrhythmia. Bispectral entropies [25] were derived to characterize the regularity or irregularity of the HRV from bispectrum plots. The formulae for these bispectral entropies are as follows:

Normalized Bispectral Entropy (BE 1):

(3)

where Ω = the region as in Fig. (1).

Fig. (1). Non-redundant region of computation of the bispectrum for real signals. Features are calculated from this region. Frequencies are shown normalized by the Nyquist frequency.

(4)

where Ω = the region as in Fig. (1).

The normalization in the equations above ensures that entropy is calculated for a parameter that lies between 0 and 1 (as required of a probability) and hence the entropies (P₁ and P₂) computed are also between 0 and 1.

In this study we also make use of features related to moments [28] and the weighted centre of bispectrum (WCOB) [29] to characterize these plots. The features related the moments of the plot are:

The sum of logarithmic amplitudes of the bispectrum:

The sum of logarithmic amplitudes of diagonal elements in the bispectrum:

The first-order spectral moment of amplitudes of diagonal elements in the bispectrum:

These features (H₁-H₃) were used by Zhou et al. [28] to classify mental tasks from EEG signals.

The definition of WCOB [29] is given by:

where i, j are the frequency bin index in the non-redundant region.

Blocks of 1024 samples, corresponding to 256 seconds at the re-sampled rate of 4 samples/sec were used for computing the bispectrum. These blocks were taken from each HRV data record with an overlap of 512 point (i.e. 50%).

3.3. Support Vector Machine (SVM) Classifier

In this study a kernel-based classifier is adopted for classification of the cardiac abnormalities. Herein, the attribute vector is mapped to some new space. Despite the fact that classification is accomplished in a higher dimension space, any dot product between vectors involved in the optimization process can be implicitly computed in the low dimensional space [30].

For a training set of instance-label pairs (x_i , y_i ), i = 0,…, l where x_i ∈ Rⁿ and y_i = {-1, 1} and if ø(.) be a non-linear operator mapping the attribute vector x to a higher dimensional space, the optimization problem for the new points ø(x) becomes

subject to the constraints

(6)

where C>0 is the penalty parameter for the error term and ξ_i are a set of slack variables that are introduced when the training data is not completely separable by a hyper plane. The SVM finds a linear separating hyper plane with the maximal margin in this higher dimensional space. As in the linear case, the mapping appears in terms of the kernel function K(x_i , x_j )= ø(x_i)^Tø(x_j). Despite the fact there are several kernels, typical choice for kernels are radial basis functions. The RBF kernel non-linearly maps samples into a higher dimensional space. There are several methods that can be used to extend a binary class SVM to multi-class SVM. In this work, we used the One against all method SVM to classify the five classes of HRV data [31]. We perform an initial grid search to obtain a suitable regularization constant and it was set to 190 for our experiment [32].

Table 2 shows the range of values for all the seven features for the five classes. The distribution of various features for the fives classes of arrhythmia are shown in Fig. (2a-g). The result of ANOVA with features obtained from HOS for various kinds of cardiac diseases is listed in Table 2.

Fig. (2a-d). Distributions of various features extracted from the bispectrum (a) bispectrum entropy 1 (P1), (b) bispectrum entropy 2 (P2), (c) sum of the logarithmic amplitudes of the bispectrum (H1), (d) sum of logarithmic amplitudes of diagonal elements (H2).

Fig. (2e-g). Distributions of various features extracted from the bispectrum (e) first-order spectral moment of amplitudes of diagonal elements in the bispectrum (H3), (f) the weighted centre of bispectrum frequency index 1 (f1m), (g) the weighted centre of bispectrum frequency index 2 (f2m).

For Normal cases, heart rate varied continuously between 60bpm-80 bpm. The bispectrum entropies (P₁ and P₂) appear to be high due to higher variation in the heart rate. The mean value of P₁ is 0.719 while that of P₂ is 0.43. The mean values of moments H_1, H_2, H₃ are 2.81e5, 1.29e3 and 1.42e5 respectively. The WCOB mean values for f_1m, f_2m are 60 and 22.32 respectively. It may be that these values are related to the rate of breathing and its harmonics. And there may be a modulating effect on the heart rate variability due to the breathing pattern.

In the Pre-Ventricular Contraction, there will be ectopic beats in the normal ECG signals. The mean entropies (P₁ and P₂) indicated in the Table 2 are correspondingly higher than the normal case due to higher variation. The mean values of moments H_1, H_2, H₃ are 2.81e5, 1.29e3 and 1.42e5 respectively. The WCOB mean values for f_1m, f_2m are 126.7 and 56.35 respectively.

In CompleteHeartBlock (CHB) the Atrio-ventricular node is unable to send electrical signals rhythmically to the ventricles and as a result the heart rate remains low. The bispectrum entropies (P₁ and P₂) indicated in the Table 2 are lower as compared to the normal subject due to the reduced beat to beat variation. The mean values of moments H_1, H_2, H₃ are 1.79e5, 8.94e2 and 8.94e4 respectively. The WCOB mean values for f_1m, f_2m are 41.95 and 12.91 respectively.

In SSS -III there is a continuous variation of heart rate between Bradycardia and Tachycardia. The bispectrum entropies (P₁ and P₂) indicated in the Table 2 are comparable with the normal case due to the higher variation in the beat to beat. The mean values of moments H_1, H_2, H₃ are 4.64e5, 1.93e3, 2.38e5 and 8.94e4 respectively. The WCOB mean values for f_1m, f_2m are 62.85 and 31.05 respectively.

In the case of Congestive Heart Failure (CHF) the heart is unable to pump the blood (supply enough oxygen) to the different parts of the body. The mean bispectrum entropies (P₁ and P₂) are lower than the normal case due to the reduced beat to beat variation. The mean values of moments H_1, H_2, H₃ are 2.02e5, 9.74e2, 1.02e5 and 8.94e4 respectively. The WCOB mean values for f_1m, f_2m are 33.71 and 10.50 respectively.

Fig. (2a-d) indicates the distribution of P₁, P₂, H₁ and H₂ features. And the distribution of H₃, f_1m and f_2m are shown in Fig. (2e-g). Table 3 shows the confusion matrix obtained by the SVM classifier. The result of classification efficiency is shown in Table 4. Our results show that, our proposed method can classify the unknown cardiac class with an efficiency of about 85%, sensitivity and specificity of 90% and 87.93% respectively (Table 5).

Different non-linear methods have been used to classify the cardiac classes using the heart rate signals [6, 33-34]. In all these studies, different non-linear parameters namely correlation dimension, Lyapunov exponent, approximate entropy, fractal dimension, Hurst exponent and detrended fluctuation analysis have been used to identify the unknown class of the disease.

In this work, we have applied higher order spectra as a non-linear tool to analyze cardiac signals. Table 5 shows the results of the application of HOS for cardiac signals. We have used the SVM and bispectral features to diagnose the different cardiac arrhythmia.

Bispectrum and bicoherence helps in the reduction of Gaussian noise reduction and detection of nonlinearities (phase coupling) in the signals [14, 35]. The real and spurious VHF (very high frequency) peaks were detected using bicoherence technique in about 50% of heart transplant(HT) subjects using the heart rate signals [13]. They have observed two types of peaks in HT subjects. First peak reflects the nonlinear response to the cardiovascular system to the respiration. Second peak indicates the nonlinear response to the nonrespiratory spectral component of unknown origin.

Different ranges of bispectrum entropies and bispectrum invariants have been proposed for different cardiac arrhythmia [25]. They have also proposed unique bispectrum and bicoherence plots for different cardiac diseases.

One of the major challenges in non-linear biosignal processing is the presence of intra-class variation. Another challenge is that there are overlaps among the derived features for various arrhythmias. Hence in our present work, we have used two bispectrum entropies and three features related to the moments and two weighted centre of bispectrum as descriptors to differentiate different arrhythmia. These features were then fed to the SVM classifier for automated classification. We achieve about 85% of classification accuracy with the current set of features. The accuracy may be further increased by extracting better features and taking more diverse training data.

A comparison with other related work of arrhythmia classification based of heart rate signals and features based on non-linear methods is shown in Table 6.

Acharya et al. pioneered the work on HRV classification using non-linear parameters such as spectral entropy, Poincare plot geometry and Largest Lyapunov exponent. In this initial attempt, they were able to achieve an average accuracy of 95% for four classes [36] and 85.36% for 8 classes [37]. Recently, Kanalthal et al. have [33] improved this work using adaptive neuro-fuzzy fuzzy inference system (ANFIS) classifier and they were able to achieve the efficacy to an average value of 94.09%. Recently, Acharya et al., have classified the cardiac states using using fast Fourier transform (FFT), auto regressive (AR), auto regressive moving average (ARMA) and moving average (MA) modeling techniques into 9 classes [38]. They have used first three peak amplitudes and corresponding peak frequencies as the features into the neural network classifier. And they were able to classify correctly up to 83.83% using ARMA modeling technique.

This is an initial work with HOS features and with further enhancement and more training data, it is possible that the efficacy could be further improved.

The Heart rate variability signal can be used as reliable indicator of some cardiac diseases. In this work, we have extracted different HOS features from heart rate signals for automated classification. We have evaluated the effectiveness of different bispectrum entropies, moments and weighted centre of bispectrum as features for the classification of various cardiac abnormalities. Our proposed system that utilizes a combination of the different features with a Support Vector Machine classifier is able to identify the unknown cardiac class with a sensitivity and specificity of 90% and 87.9% respectively. The accuracy of our proposed system depends on several factors, such as the size and quality of the training set, the rigor of the training imparted, and also the extracted HOS parameters.