Electrocardiogram Fiducial Points Detection and Estimation Methodology for Automatic Diagnose

2.1.1. Choice of Curvature Filter Orders

In order to obtain an accurate result in using curvature filters for wave onset and offset estimation, an appropriate order n must be chosen on each case. In this section we argue about how to choose this optimal order. In order to do that, we test the curvature filters on a test example designed for offset estimation. The picture for onsets is completely similar because of symmetry.

We fix a time wave length w_t, in milliseconds, and a corresponding wave length w_s, in number of samples. The estimated sampling frequency is then obtained by the formula F_s ~ 1000w_s/w_t, in Hertzs, where this quantity is rounded to the closest integer. The base signal is composed by a semi ellipse of length w_s samples and height h (in millivolts), followed by w_s zeroes. Consequently, the signal length is 2w_s. The base signal is perturbed with random noise of a given intensity h_n (in millivolts), giving a set of N_t perturbed signals of length 2w_s.

In the next test example we use the following values: w_t=110 ms, h=0.15 mV, h_n=0.01 mV, and N_t=100 signals (according to the P wave in the Table 1). The wave lengths used to test the curvature filters in this example go from w_s=3 samples up to w_s=100 samples, by increments of 1 sample.

Next, for each noisy signal s_j (with j=1,...,N_t), and for each order n (with n=3,...,2w_s-1), we compute the corresponding normalized curvature coefficients. Note that we do that only for sample points such that the window is completely included into the signal range, as we are dealing with a finite signal; that depends on the order n of the curvature filter used. Then, we estimate the offset by selecting the greatest normalized curvature coefficient; thus, its offset is defined as the corresponding sample point k=o_j,n with greatest normalized curvature coefficient. Note that, as we are looking for offsets for upward waves, we might select points with maximum convexity; i.e., greatest positive curvature.

In Figs. (2 and 3) we show the resulting estimated offsets o_j,n for j=1,...,N_t, and n=3,...,2w_s-1, for the wave lengths w_s=15 samples (corresponding to F_s=136 Hz) and w_s=55 samples (corresponding to F_s=500 Hz). We consider as a reference correct offset the sample point k=w_s.

In view of Figs. (2 and 3) and the rest of simulations, several comments arise about the curvature filters performance:

• If the curvature filter order is too small, then the noise has a profound effect in the curvature coefficients. That is, low order curvature filters do not distinguish between signal twists and noise.
• From certain large enough order, curvature filters are very proficient on the task of detecting the offset.
• The accuracy of the measurement slightly deteriorates as the order increases. On average, a trend to shift outwards the wave is detected.
• The latter convergence for very high orders is spurious since is due to the finite character of the signal.

These remarks suggest that curvature filters could be an useful tool on detecting fiducial points in ECG signals, provided that a suitable order is chosen for the specific task. We think this is valid also for features extraction in signal processing, in general.

In order to measure quantitatively the curvature filters performance in the test example, we define the error as = (w_s being fixed), for any j=1,...,N_t, and n=3,...,2w_s-1. Thus, for each order n we compute the mean error

(15)

Then, the mean optimal order is defined as the (smallest) order n=n_m for which is minimum . In other words, the mean optimal order is the window size for which the greatest-curvature-coefficient method (i.e., the previous one) for detecting the offset is more accurate, on average. Note we avoid too large orders in this criteria.

Another important features which are also highly desirable, apart from accuracy, are fiability and robustness under noise presence. In order to take into account these aspects, we also analyse the measurements error dispersion provided by the greatest-curvature-coefficient method. Hence, we consider the error standard deviation

(16)

In Figs. (4 and 5) we see that the magnitude of the errors dramatically falls off after the initial noise-sensitive period, as expected from the above remarks based on Figs. (2 and 3). Coherently, the mean error and the error standard deviation σ_n also have the same behavior. We define the standard deviation optimal order as the first order n=n_std such that σ_n ≤ σ_n+1 and σ_n is lower than a reference value empirically set to 2. That is, we look for a small enough local minimum of the standard deviation.

As it has been remarked before, immediately after the initial noise-sensitive period the greatest-curvature-coefficient method performs very well, providing a quite low mean error almost stagnant in a long period on n. The same could be said for the standard deviation of the error.

As it is illustrated in Fig. (6), it often occurs that the mean optimal order unnecessarily increases to larger values. Due to computational reasons, it is preferable to use low order filters as long as possible. Another factor to bear in mind is that, generally, waves do not occur isolated in an ECG, but surrounded by other ones. Hence, too high orders are undesirable since the presence of neighboring waves could cause interference.

In order to avoid that, we recommend as a preferential choice the standard deviation optimal order n_std to be used in the greatest-curvature-coefficient method.

In Fig. (7) we can see the optimal orders (mean optimal order n_m and standard deviation optimal order n_std) for w_s in the proposed range (that is, wave lengths from w_s=3 samples up to w_s=100 samples, by increments of 1 sample).

In summary, in view of Fig. (8), we recommend to use a curvature filter order between the 15% and the 30% of the wave width w_s. These values should be increased for very low sampling frequencies, as the resulting from Holter monitors. This quantity could be estimated a priori from the estimated wave length w_t (in milliseconds, which could be taken from references as in Table 1), and the signal sampling frequency F_s (in Hertz), which is assumed to be known, by the formula.

(17)

2.2.1. Preprocessing

Generally, the ECG signal is assumed to be contaminated with several noise sources. This noise could be reduced by using some filters which cut off different frequencies, but definitely cannot be completely avoided. Anyway, for this reason and previous to the feature extraction process, a preprocessing is necessary to be applied to the signal. Consequently, the original signal suffers changes which add some uncertainty to the final results. This phenomena is specially remarkable in the case of P and T waves onset and offset location. That is why we speak about feature “estimation”, and not about “measurement”.

First of all, the amplitude of the original signal is normalized to 1:

(18)

The resulting signal has amplitude 1. Next, the direct level is supressed, obtaining a signal given by

(19)

Obviously, one has mean (ecg2)=0.

After that, a high pass filter of type Butterworth of order 2 is applied to filter ecg2, with a cut frequency of 0.5 Hz, to obtain a preprocessed signal . This step is always present in ECG professional applications, in order to remove very low frequencies.

When applying this preprocessing to signals from the MIT database [73], which proceed from Holters, this last step is very important and helps to improve the accuracy in subsequent estimations. On the other hand, when dealing with signals from pattern generating devices, it is usually required an additional filtering of smoothing or averaging type. In general, the required preprocessing highly depends on the measurement system features used to generate the ECG signal, and could include the use of other filters, as interference supressing filters out of 60 Hz. For the signals we have used in this work, this is not the case.

2.2.2. R_m Peaks Location and Heart Rate Analysis

QRS complex detection is considered a fundamental step in the analysis of ECG signals, because of measurement and characterization of different associated parameters rely on the accuracy with which it is performed [2, 19]. On the one hand, from the QRS complex cardiac frequency can be obtained, and subsequently heart rate is also obtained. On the other hand, R_m peaks almost delimitate each heart beat (more clearly than other waves) in such a way that estimation of the other fiducial points is more easily carried out from them.

In the recent decades many QRS complex detection methods have been reported in the literature [1, 2, 4, 17-19, 24], and references therein), including derivatives approximation, digital filters, wavelets, adaptative filters, neural networks, hidden Markov models or mathematical morphology, among other techniques.

In this work, we make use of the techniques developed in the literature [24]. Concretely, we select the pre-processed signal energy maximum, choosing a threshold which is empirically set. An appropriate choice of these maximum allows to obtain the heart rate (in beats per minute). Moreover, we describe two different methods to locate the R_m peaks in the original signal (in number of samples): one can simply extract the maximum peaks, or one can perform a crossed correlation between the filtered signal and the original signal (normalized and prefiltered).

We start with the preprocessed signal ecg3. Previous to energy processes, we filter it by using a band pass filter, in the band 10-25 Hz or 15-20 Hz, obtaining a filtered signal . This step is mandatory in order to remove artifacts, interferences, and the influence of P and T waves. Note that in this moment we are only interested in detecting R_m peaks (Fig. 9).

Next, we consider the energy signal E₄ given by E₄ = . The energy threshold E_th is empirically set to the value

(20)

Then, the energy maximum are extracted with the constraint E₄ > E_th. In the MATLAB code, the energy peak extraction is performed by using the implemented built-in function findpeaks [74]. The resulting energy peaks have to be sifted out, since they usually appear pairs of “peaks” which are unnaturally quite close (presumably corresponding to the same R wave, due to the noise and filtering processes). These related peaks have to be removed, since they produce non physiological values for the cardiac frequency.

Each energy peak corresponds to an R_m peak, but the last ones need to be located, since the band pass filtering provokes a constant gap between the two signals, which shifts the peaks in time.

What we do is to compute the cross-correlation and finding the correlation coefficient between the energy E₄ of the filtered signal and the energy E₃ (given by ) of the original preprocessed signal. Then, the shift is the difference between the signal length L and the lag corresponding to the correlation coefficient K_cr, the maximum of the cross correlation. In that way, we are able to locate the R_m peaks in the original preprocessed signal ecg3. In the code, the cross-correlation analysis is performed by using the MATLAB built-in function xcorr [74].

Once the R_m peaks of ecg3 are located, the RR intervals are known. The widths of RR intervals, or distance between consecutive R_m peaks, allow to compute the approximated cardiac frequency: If L_RR is the width of a given RR interval (in number of samples), and HR is the heart rate (measured in beats per minute bpm), then

(21)

2.2.3. QRS Complex Location

In order to describe the QRS complex, we have to locate the Q_i, Q_m, S_m and J points of each cycle. Here, for exposition purposes, we consider the R_m point as the maximum of the R wave, and the Q_m and S_m points are understood as the minimum of the Q and S waves, respectively. In addition, Q_i is the starting point of the Q wave and the end of the PQ segment, and J is the ending point of the S wave and the start of the ST segment.

For this task, we proceed by using a modified slope analysis method with restrictions, starting from a peak. This algorithm is inspired by Chapter 9 of [2]. The basic idea is to descend/ascend from a maximum/minimum, respectively, until a change in the sign of the slope is found. On a given sample point k (with 1 ≤ k < L), the slope of the signal ecg3 is defined by

(22)

This method is justified because of in the pronounced ascents and falls of the signal, the influence of the noise is quite lower than in the isolectric segments.

In this analysis we remove the first and final cycles, since they could be incomplete or distorted because of the original signal capture process. Then, consider the original preprocessed signal ecg3, and fix some R_m peak (not the first, not the last one).

On the one hand, we move backward from the R_m peak until we find the first local minimum; that is, until we find a point k whose left neighbor has nonpositive slope s_k-1 ≤ 0. This point is defined as the Q_m point of this cycle. From this point, we continue backward looking for the first local maximum; that is, a point k such that its left neighbor has nonnegative slope This point is defined as the s_k-1 ≥ 0. This point is defined as the Q_i point of this cycle.

On the other hand, we proceed analogously to the right. From the R_m peak, we move forward until we find the first local minimum; that is, until we find a point k with nonnegative slope s_k ≥ 0. This point is defined as the s_m point of this cycle. From this point, we continue forward looking for the first local maximum; that is, a point k with nonpositive slope s_k ≤ 0. This point is defined as the J point of this cycle.

In practice, this method needs to be complemented with different criteria about how far the search should go, related to the estimated length of the Q and S waves. This is because of, depending on the presence or absence of Q and S waves in the ECG signal, the pure slope analysis algorithm could produce fiducial points far away from the R_m peak, which would be unnatural. On the other hand, the greatest/smallest-curvature-coefficient method could assist the slope analysis method, as in the following subsection is illustrated, specially in the case of the Q_i and J points. However, the short width of these waves forces to use a quite small order, which could be imprecise, specially for signals with very low sampling frequency.

2.2.4. Estimation of Remaining Fiducial Points

At this point, what remains is to estimate the fiducial points in P and T waves. We proceed by using a modified slope analysis method [2], as before, combined with the use of curvature filters.

First of all, we find the P_m and T_m peaks. Without loss of generality, we assume we deal with upward P and T waves. The case of downward waves is completely analogous, and an algorithm to distinguish between upward and downward cases is easy to implement by computing maximum and minimum and comparing them (Q_i and J points could be taken as isoelectric reference level for the comparison, respectively).

Hence, from the R_m peak, we look for a global maximum in a preestablished backward/forward range of samples of the signal, respectively. This range could be estimated from the reference values for the PQ and QT intervals (Table 1), the sampling frequency and the previously estimated fiducial points. Thus, the left maximum corresponds to P_m, and the right one corresponds to T_m.

From each one of these points, we descend in order to find the corresponding onset and offset points. The idea is to scan the surrounding signal looking for suitable candidates by using the slope analysis algorithm, and then selecting the best choice with the aid of the curvature filters. The process is as follows.

We start from the P_m point (the process for the T wave is completely analogous). On the one hand, we move backward inside a preestablished range looking for local minimum: Points k whose left neighbors have nonpositive slope s_k-1 ≤ 0. At each step, we store them as candidates. When the Range is exhausted, we select the winner by the greatest-curvature-coefficient method: We compute the curvature coefficients of the candidates, and the sample with greatest one among them is designated as the P_i point of the cycle. On the other hand, from the P_m peak we move forward inside a preestablished range (never farther than the corresponding Q_i point) looking for local minimum: points k with nonnegative slope s_k ≥ 0. They are the candidates, and the P_f is selected among them by the greatest-curvature-coefficient method.

Note that a stopping criteria is necessary in this algorithm. As before, the ranges could be estimated from the reference values for the PQ and QT intervals (Table 1), the sampling frequency and the previously estimated fiducial points. Also, it could happen that candidates which are quite close to the P_m/T_m peaks, respectively, might have to be removed, as at the top of the P and T waves they could exist points with high curvature coefficients, specially in cases of nonstandard morphologies. Consequently, a starting criteria could be also necessary. On the other hand, suitable curvature filters orders for the P and T waves need to be chosen for a well performance of the algorithm (subsection 2.1.1.).