Robust projective filtering of time-warped ECG beats

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Robust projective filtering of time-warped ECG beats M. Kotas ∗ Institute of Electronics, Silesian University of Technology, Akademicka Street 16, 44-100 Gliwice, Poland

a r t i c l e

i n f o

a b s t r a c t

Article history:

In this paper a new version, more immune to noise, of the nonlinear projective filtering

Received 12 January 2008

is presented. The method employs an algorithm of robust principal component analysis

Received in revised form

to signal subspaces construction and as a result it achieves high performance in real elec-

16 May 2008

tromyographic noise environment. Two aspects of the method’s action are investigated: its

Accepted 5 June 2008

ability to suppress noise and its influence on the precision of the QT interval measurement. Then the method influence on evaluation of the beat-to-beat variability of the repolariza-

Keywords:

tion duration is presented. A comparison to the previous versions of the nonlinear projective

Robust principal component

filtering and to the classical linear one is carried out. © 2008 Elsevier Ireland Ltd. All rights reserved.

analysis Projective filtering Dynamic time warping ECG processing

1.

Introduction

The electrocardiographic (ECG) signal, representing the electrical activity of the heart, is often contaminated by noise. Therefore, since the beginnings of the noninvasive cardiology, development of the effective methods of the electrocardiographic noise suppression have stimulated the progress in the field. Great possibilities emerged when the digital filters with linear phase response, allowing for suppression of ECG noise with limited distortions of the desired component [1], were introduced. This type of filters is particularly effective when applied to baseline wander and powerline interference suppression [2,3]. To enable the analysis of the signals contaminated by the wide-band electromyographic (EMG) noise, the synchronized averaging was introduced. Assuming the constant morphology of the ECG beats, the method performs linear time-alignment [4] of the respective cycles and construction of an average one. As a result, the noise is suppressed and the repetitive desired component is enhanced. The assumption of constant morphology of the averaged sig-

∗

nal was relaxed in [5], where the technique of dynamic time warping was employed to the nonlinear alignment of the echocardiographic left ventricular volume cycles. Averaging the aligned cycles allowed to reduce the wide-band noise and the periodical artefact related to heart movements synchronous with respiration. As a result, the method caused reduction of inter-beat variability of the processed signal and improved its quantification [5]. In many applications it is advantageous to suppress noise without the suppression of the variability of the desired signal morphology. In such cases it is necessary to construct a space of possible shapes instead of constructing an average template. This task can be accomplished by application of the principal component analysis [6,7]. However, the intrinsically linear principal component analysis [8] models complicated, possibly nonlinear relationships among different parts of an ECG beat with a limited success [7]. A great progress was achieved when the method of the nonlinear state-space projections (NSSP) was applied to ECG processing [9]. Nonlinear state-space projections allow to suppress noise whose power

Tel.: +48 32 2372019. E-mail address: [email protected]. 0169-2607/$ – see front matter © 2008 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2008.06.007

162

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

spectrum overlaps the spectrum of the desired component [10]. In [11] NSSP was compared to several sophisticated methods of ECG noise suppression and achieved the highest increase of the signal-to-noise ratio in EMG noise environment. The limited number of the method’s applications results from its extremely high computational costs. The modified version of this method, called as projective filtering of timealigned beats, allowed to overcome this problem to some extent [12]. The method was successfully applied to ECG enhancement prior to the measurements of the repolarization duration [13]. The next version, called as projective filtering of time-warped ECG beats (PFTWEB), achieved significantly higher performance [14]. Still, its applicability is limited to the cases of a not particularly high level of noise. The goal of this study is to combine the dynamic time warping with the robust principal component analysis to develop a new version (more immune to noise) of the projective filtering. The filter usefulness will be illustrated with its application to ECG enhancement prior to the analysis of the ventricular repolarization duration (the QT or RT interval [15–17]). The rest of this paper is organized as follows. Section 2 describes the PFTWEB method. In Section 3 the new version of projective filtering is developed, and in Section 4 the numerical experiments are performed. Finally, conclusions are drawn in Section 5.

In [14] it was shown that imposing a few reasonable constraints on the contents of the neighborhoods justifies employing the method of dynamic time warping [18,19] to neighborhoods determination.

2. Projective filtering of time-warped ECG beats

k = {ak + l|l = 0, 1, . . . , Nk },

The method was developed [14] by introduction of a few modifications that raised the performance of the nonlinear state-space projections (NSSP) [9] when applied to ECG signals processing. The fundamental operation is the reconstruction of the state-space representation of the observed noisy signal by application of the Takens embedding operation. A point in the reconstructed space is a vector:

 

x(n) = x n−

m 2

 



−1  , x n−

m 2

 



−2  , . . . , x n+

m  2

 (1)

where x(n) is the processed signal,  the time lag ( = 1 [14]), m is the embedding dimension (here m is even and the m/2 th coordinate of the embedding vector x(n) is the nth sample x(n); such an almost symmetric form of the embedding was applied to simplify the description of the method). In the original NSSP method, for each point x(n) a small neighborhood  (n) is constructed, composed of the points which are close to x(n)  (n) = {k|x(k) − x(n)  < ε},

2.1. Dynamic time warping for neighborhood determination A preprocessing step embraces baseline wander and powerline interference suppression [2,3], QRS complex detection [20], and synchronization of the detected complexes. The last operation produces a set of fiducial marks {rk |k = 1, 2, . . . , K}, corresponding to the same position within the respective detected QRS complexes. To simplify further considerations we assume that the successive complexes belong to the same dominant class. Lets denote: ak = rk − b,

Nk = rk+1 − rk + 2b

(3)

where b ∈ N is a small number introduced to immunize the method against small errors of the fiducial marks determination. The embedding space points are divided into the successive, slightly overlapping sequences of vectors whose time indices belong to the following sets k = 1, 2, . . . K − 1.

(4)

The general concept is to perform the nonlinear alignment of the successive sequences of vectors (k , k = 2, 3, . . . , K − 1) with respect to the first sequence (Fig. 1) and this way to determine the neighborhoods corresponding to the respective positions within 1 . For small b the side points x(a1 ) and x(a1 +N1 ) overlap QRS complexes (see Fig. 1), and thus their neighborhoods can be determined on the basis of the preprocessing step results. The neighborhood of the left side point is a set containing the points that occupy the same position with respect to the successive K − 1 QRS complexes  (a1 ) = {ak |k = 1, 2, . . . , K − 1} and similarly the neighborhood of the right side point is  (a1 +N1 ) = {ak + Nk |k = 1, 2, . . . , K − 1}. Determination of the neighborhoods corresponding to the respective positions within the sequence of vectors from the first set (1 ) is performed in the following way. (1) In the first step each neighborhood of a point from 1 is filled with this point only, and so the time index of this point is included in the following set  (a1 +l) = {a1 + l},

l = 0, 1, . . . , N1 .

(5)

(2)

where || · || denotes the Euclidean distance and ε is the assumed maximal distance between the points. Assuming that in the embedding space locally within a given neighborhood the uncontaminated points are located on, or very close to, a low dimensional hyperplane, NSSP performs projection of x(n) on the determined hyperplane (called as the local signal subspace).

These points become the first mass centers of the respective neighborhoods. Then we set k = 2 and go to the next step. (2) The sequences of points that belong to k are time-warped with respect to the sequence of the mass centers. For this purpose we calculate the so called warping paths wk,l , l = 0, 1, . . . , N1 containing the indices of the vectors x(ak +wk,l ) that are aligned with the successive mass centers x¯ (a1 +l) .

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

The warping path is searched for by minimizing

Qk =

N1 

x(ak +wk,l ) − x¯ (a1 +l) ,

(6)

l=0

while preserving the border conditions, the monotonicity conditions and the condition restricting the number of successive points from k that may be omitted in the warping path [14]. The minimization is performed with the use of the dynamic programming [19]. (3) After the kth sequence k has been time warped with respect to the sequence of the mass centers, the neighborhoods are supplemented with the points whose indices are stored in the warping path  (a1 +l) =  (a1 +l) ∪ {ak + wk,l }

(7)

(4) The mass centers of the respective neighborhoods are updated and k is incremented. If k < K − 1 we go to step 2, otherwise we end the algorithm. After determining the respective neighborhoods we construct the corresponding signal subspaces and then we perform projective filtering.

2.2.

Construction of signal subspaces

Within the respective neighborhoods the local signal subspaces are constructed by application of principal component analysis. Each subspace is spanned by a set of q orthonormal principal directions (axes) ei (i = 1, 2, . . . , q) which have the following properties: the first principal axis e1 assures the maximal variance of the neighborhood points projected on it, the second principal axis e2 assures the greatest variance of the points projected on it in the subspace perpendicular to e1 , etc. If the noise free points of a neighborhood lie in a low dimensional linear manifold or if they are close to it, and if the noise deviations are of moderate energy then we can expect a few principal directions to compose a signal subspace. If, however, the level of noise is too high or some abrupt high energy noise appears then the principal eigenvectors can be stretched in the

163

noise direction. It results from the lack of robustness against outliers of the classical PCA. To immunize PFTWEB against such phenomena, the assumed percentage of the most distant points of the constructed neighborhoods are rejected before application of PCA to signal subspaces construction [14].

2.3.

Projective filtering

For each point under projection (x(n) ) the corresponding sig(o)

(o)

(o)

nal subspace (Eq = [e1 , . . . , eq ]) is searched for [14](q is the assumed dimension of the subspace, o is its number—the position within k ). The projection operation is defined as x

(n)

(o) (o)

= x¯ (o) + Eq Eq

(x(n) − x¯ (o) )

(8)

where x¯ (o) is the neighborhood mass center. Since each point is composed of m successive samples of the processed one-dimensional signal (1), the nth sample of the signal occurs in m trajectory points. Projecting these points into the corresponding signal subspaces results in m possibly different corrections of this sample. The resultant correction can be calculated as the average of those m suggested ones. Detailed description of the operation can be found in [14].

3. Robust projective filtering of time-warped ECG beats (RPFTWEB) Three modifications were introduced to raise the immunity of PFTWEB to high level EMG noise. • Modification 1 (M1). The lack of robustness against outliers of the classical PCA results from the fact that the method is based on the nonrobust statistical estimators: the sample mean and the sample covariance matrix, which are sensitive to outliers themselves. In this study the algorithm of robust PCA, presented in [21], was applied. In this algorithm, the sample mean is replaced ˆ n ) and the by the robust spatial L1-median [22](denoted as m covariance is replaced by the Qn estimator [23] as a measure of dispersion.

Fig. 1 – A section of ECG signal containing two successive QRS complexes. Two delay windows are presented, containing one-dimensional representations of the points x(a1 ) and x(a1 +N1 ) . They are respectively the first and the last point with time indices that belong to the set 1 = {a1 + l|l = 0, 1, . . . , N1 }; see (3) for the definitions of ak and Nk .

164

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

Lets consider application of this algorithm to the analysis of a sample of points whose time indices are contained in a set  : {x(j) |j ∈  }. The Qn estimator is used to find the robust principal directions that maximize the dispersion of the values z(j) = a x(j) , j ∈  , obtained by projecting the points x(j) on these directions. The search for a direction that maximizes an estimator of dispersion is a time consuming problem. A simplification of this problem, proposed in [21], is based on restricting the space of possible solutions to the set  of ˆ n (the spatial L1-median of the directions passing through m the considered sample of points) and through the respective points ˆ n |j ∈  }  = {x(j) − m

(9)

Thus the first “robust eigenvector” based on the Qn estimator is defined as eQn ,1 = arg maxQn ({a x(j) |j ∈  })

(10)

a ∈ An,1 ()

where An,1 () = {(d(j) /d(j) )|d(j) ∈ }. The successive eigenvectors are searched for within the sets of the directions that are projected into the orthogonal compliments of the subspaces spanned by the previous eigenvectors [21]. • Modification 2 (M2). The applied version of robust principal component analysis has the following drawback. When the signals are contaminated by a stationary white noise, the neighborhood points are also contaminated. Thus the set  (9) contains noisy candidates for principal directions. To overcome this problem, before calculating the principal directions the set  undergoes the following operation. For each d ∈  we search for the nearest neighbors a, b ∈  (a is the nearest one). Then we set c = 1 and we check if d − a < ca − b. If the condition is satisfied, we calculate d = ((a + b + c)/3) and we replace d. If it is not satisfied, we exclude d from . If too many vectors are rejected from , we set c = 2c and we repeat the operation on the initial set  (as defined in (9)) until enough vectors remain. This modification allows to reject the vectors which are far from their nearest neighbors (usually as a result of being contaminated) and preserves the less noisy ones with additional smoothing as a result of calculating the mean. • Modification 3 (M3). The third modification is a result of the following reasoning. The original NSSP method needs a high embedding dimension to prevent large errors of neighborhoods determination (e.g. including to a neighborhood the points that occupy completely different positions within the respective ECG beats). Applying dynamic time warping to this operation we introduce the constraints which prevent large errors. Therefore it seems possible to obtain good results of neighborhoods determination with lower embedding dimensions. To check this possibility, the method was modified so that it

allows to apply different dimensions during determination of neighborhoods, and during construction of signal subspaces and ultimate projective filtering. Thus for this method two parameters should be established instead of the dimension m: mnd —the dimension applied for neighborhoods determination, mpf —the dimension applied for the ultimate projective filtering. Owing to the symmetric definition of the embedding operation (1) the time indices of the points contained in the established neighborhoods can without any changes be employed during construction of the signal subspaces. Since the parameters: m, mnd and mpf depend on the sampling frequency, it is advantageous to establish the corresponding time intervals: m , nd and pf (which can be called as embedding or delay windows) and to calculate the corresponding dimensions for a particular sampling frequency applied.

4.

Numerical experiments and discussion

In this section two aspects of ECG signal processing are investigated: the projective filters ability to suppress noise and their influence on the dynamical measurements of the QT interval. Finally, their influence on evaluation of the beat-to-beat variability of the repolarization duration is presented.

4.1.

The projective filters performance

This issue had previously been tested in [7,12,14], and a similar experiment was performed in this study to compare the described approaches to projective filtering of ECG signals. The tests were based on the signals from the MIT-BIH database, provided by the Massachusetts Institute of Technology and Boston’s Beth Israel Hospital. The database consists of several smaller databases for various test purposes. Among others are the ECG Arrhythmia Database (MITDB) and the Noise Stress Test Database (NSTDB) utilized in the experiments that will be presented. The five signals of high quality that had previously been employed in [7,12,14](from the Arrhythmia Database) were used to simulate the desired ECG component. The signals were contaminated artificially with either white Gaussian or real electromyographic noise (from NSTDB) of a few different levels. The first one allows to estimate the methods’ performance in stationary noise of the assumed level [4,6,7,9,12,14,15]. In real ECG signals the noise is usually of different properties. The EMG record from the NSTDB database is an extreme example of it. It is highly nonstationary and highly correlated. Both the types of noise span a kind of range of the cases that can be encountered in real signals. Thus the methods that achieve high performance in such different conditions may be expected to be effective in real life applications. Before constructing the artificially contaminated signals, the desired ECGs and the noise records had been preprocessed separately (like as in [7,12,14]) to suppress the baseline

165

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

Fig. 2 – The average noise reduction factor obtained for SNR = 10 dB in either EMG or white noise environment. For RPFTWEB the horizontal axis corresponds to  pf ;  nd =  pf − 100 ms.

wander and the powerline interference. The contaminated signals underwent processing with the use of the investigated versions of the projective filters. The results of signal enhancement were estimated by application of the following noise reduction factor:

  2  (x(n) − s(n))   n , NRF =   2 (x (n) − s(n))

(11)

n

where s(n) is the preprocessed ECG used to simulate the desired component, x(n) − s(n) = w(n) is the preprocessed noise (added to contaminate the desired ECG), x (n) is the result of the projective filtering, and x (n) − s(n) is the residual noise. Different parameters of the projective filters were applied: a few embedding dimensions in combination with a few different dimensions of signal subspaces. For each set of parameters in each conditions tested, an average noise reduction ratio was calculated. The filters performance for SNR = 10 dB in both types of noise is presented in Fig. 2. We can notice that the original method of projective filtering (NSSP) performed best for q = 1 and application of higher dimensions of local signal subspaces caused deterioration of the results. PFTAB performed best for q = 2 and m = 200 ms. This method can achieve higher NRF when applied with variable dimensions of signal subspaces [12] but this will not be investigated in

this study. Application of dynamic time warping to neighborhood determination (PFTWEB) allowed a distinct increase of the projective filtering performance for q = 2 and q = 3. With such higher dimensions of signal subspaces, the method is more flexible and can preserve more details of the processed signals morphology. The further increase of NRF was achieved by RPFTWEB. For this filter different combinations of nd and pf were tested. The presented case, when nd = pf − 100 ms, appeared particularly advantageous. The above remarks are confirmed by the values gathered in Table 1. The upper six rows contain the results achieved by the compared methods when applied with the parameters assuring their highest performance for SNR = 10 dB (in white Gaussian or EMG noise environment like as in the previous experiment). For reference the best results achieved by principal component analysis of the whole synchronized ECG beats [7] were presented. All the developed versions of projective filtering, including the original method of nonlinear state-space projections, achieved significantly higher performance than this classical approach. However, it was most often RPFTWEB that performed best in the respective conditions tested. The method was particularly effective in EMG environment. Analyzing the results achieved by two versions of robust projective filtering: RPFTWEB*— developed by introduction of modification M1 only, and RPFTWEB— including all modifications, we can compare the significance of those modifications. M1 caused a distinct increase of projective filtering performance in EMG noise environment, but a slight decrease in stationary

Table 1 – The average NRF obtained for all desired ECGs contaminated with either white (WN) or real EMG noise Method

WN SNR = 5 dB

PCA, q = 3 NSSP, m = 500 ms, q = 1 PFTAB, m = 200 ms, q = 2 PFTWEB, m = 400 ms, q = 2 RPFTWEB*, pf = 500, ns = 500 ms, q = 3 RPFTWEB, pf = 500, ns = 400 ms, q = 3 PFTAB, m = 200 ms, q = 3 RPFTWEB, pf = 500, ns = 400 ms, q = 4

3.10 4.22 3.56 4.53 4.32 4.60 3.25 4.53

SNR = 10 dB 2.56 3.28 3.01 3.47 3.10 3.41 2.94 3.51

EMG SNR = 20 dB 1.18 1.33 1.38 1.44 1.39 1.52 1.68 1.70

SNR = 5 dB 1.88 2.34 2.18 2.54 3.22 3.42 1.97 3.24

RPFTWEB* denotes the method developed by introduction of modification M1 (without M2 and M3).

SNR = 10 dB

SNR = 20 dB

1.77 2.28 2.05 2.44 2.89 3.03 1.90 2.91

1.10 1.28 1.26 1.34 1.48 1.53 1.40 1.66

166

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

Fig. 3 – Results of processing the signals of high SNR resampled with the frequency of 120 Hz. For reference the results of the classical low-pass filtering are presented.

white Gaussian noise (if compared to PFTWEB). The next two modifications caused the further increase of NRF in all conditions tested. The lowest two rows of Table 1 contain the results achieved by PFTAB and RPFTWEB when applied with higher dimensions of signal subspaces. The applied dimensions allow to reconstruct the desired ECG with greater precision (which is signalled by the higher NRF values for SNR = 20 dB). Although for these parameters the filters performance in higher level of noise (SNR = 5 dB) is lower, they are applied in the experiments presented in the next subsections. This choice results from the necessity to reconstruct the desired ECG with very high precision to allow accomplishment of a very demanding task of evaluating the beat-to-beat variability of the repolarization duration.

4.1.1.

Projective filtering in low sampling rate

The NSSP method was reported [9] to perform well the signals stored with a low sampling rate. This property was retained by RPFTWEB. It is illustrated in Fig. 3 which shows results of processing the signals of high SNR (chosen among the “desired” ECGs that were employed in the previous experiments). While the left signal contains high and steep QRS complexes, in the right one the complexes amplitude is highly modulated by respiration and cardiac axis movements. RPFTWEB copes very well with high slopes of the QRS complexes which is evident if we compare the difference between the original and the filtered signal after RPFTWEB with the corresponding differ-

ence after the low-pass Hanning filter (H(z) = (z + 2 + z−1 )/4). For the Hanning filter, the difference near QRS complexes is approximately ten times higher. For the signal on the right side of the figure, the method is less effective. Although it copes well with the respiratory modulation of the signal amplitude when the cardiac axis rotation caused a distinct change of the QRS complex shape, the difference s(n) − s (n) increased (compare the third complex before and after RPFTWEB). The QRS complex is rather wide in the discussed case, and the difference s(n) − s (n) caused by the Hanning filter near the respective complexes is of lower amplitude than in the previous case. Nevertheless, except for the third QRS complex, it is comparable or even higher than the difference caused by RPFTWEB. Fig. 4 illustrates RPFTWEB performance after the discussed signals were contaminated with EMG noise. We can notice that the level of the added EMG noise was considerably reduced, and only locally, in the discussed case of the QRS with changed morphology, the introduced signal deformation is of higher energy than the suppressed noise. Before applying the projective filtering we have to suppress the baseline wander and the powerline interference. The significance of signal prefiltering for baseline wander suppression is illustrated in Fig. 5. When the low frequency noise of a high level was simulated and RPFTWEB was not preceded by any filtering operation, the method’s performance was considerably disturbed (the upper case: subplots A–C). Not only the low frequency noise was not suppressed, but reconstruction

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

167

Fig. 4 – Results of processing the signals from Fig. 4 contaminated additively with EMG noise.

of the desired ECG was very poor (see the residual noise in subplot C). When the level of the noise was lower, RPFTWEB managed to reconstruct the desired ECG with higher precision, but again it did not suppress the baseline wander. This confirms the necessity to perform suppression of narrow-band noise before application of projective filtering.

4.2. The influence of projective filtering on the precision of the QT limits determination In the next experiment, as in [13], six signals of high quality from the QT database [24] were exploited to simulate the desired ECG. Each “desired” ECG was contaminated with six different nonoverlapping segments of the electromyographic noise from the MIT-BIH database (before this operation the EMG signal, originally sampled with the frequency of 360 Hz, had been resampled with the frequency applied in the QT database—250 Hz). As in the previous experiment, before adding the desired component and the noise they were prefiltered to suppress the baseline wander and the powerline interference. The simulated noisy signals were used to test the precision of the QRS onset (QRSo ), the T wave peak (Tp ) and the T wave end (Te ) determination. The classical method developed by Laguna et al. [25], slightly modified in [13], was applied to accomplish the measurements. In [13] it was shown that enhancing the ECG signals by PFTAB leads to more accurate measurements. In this study we compare the new method developed (RPFTWEB, ns = 400 ms, pf = 500 ms, q = 4) with the method tested in [13](PFTAB, m = 200 ms, q = 3). For reference we also show the results achieved by application of a linear filter. The Butterworth low-pass filter

of the fourth order and the cutoff frequency of 15 Hz was chosen for this purpose (to avoid phase distortions double filtering was applied—in forward and reverse direction). The measurements were performed on the basis of the desired ECG signals, the signals contaminated with EMG noise, and the contaminated signals enhanced by the tested filters. The root mean-squared error (RMSE) of the measurements as a function of the signal-to-noise ratio is presented in Fig. 6. Each RMSE value in Fig. 6 is an average of 36 values obtained for the respective test signals (the six desired ECGs contaminated by the six nonoverlapping intervals of the EMG noise). The upper part of the figure presents the curves obtained when all errors were taken into account during calculations of the RMSE values. When large errors (greater than 10 ms for QRSo and then 30 ms for Tp and Te ) were rejected before the calculations, the curves presented in the lower subplots were obtained. The percentage of the rejected large errors is presented in Table 2. Analyzing the values presented in Table 2, in the columns corresponding to Tp and Te localization, we can make the following remarks.

• While determining Tp and Te by application of the original procedure [25,13], large errors occurred for SNR = 20 dB, and for lower ratio their number highly raised. • When linear filtering preceded the measurements, the first large errors occurred for 20 dB as well, but their number was much lower. • PFTAB prevented the measuring procedure from large errors for SNR of 20 dB, and RPFTWEB managed to do this even for 15 dB.

168

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

only the low frequency components of the ECG signals but of the higher frequency parts, such as the Q wave, as well.

4.3. Projective filtering for evaluation of repolarization duration variability Precise determination of the Q wave onset and the T wave end in the EMG noise is relatively difficult therefore new variables were defined to quantify the repolarization duration in noisy signals [15–17]. It was shown [15,17] that the most precise measures in such noisy environment can be obtained if the left limit of the repolarization interval is defined as the position of the R wave peak and the right limit as the position of the T wave peak. The RTp interval is highly correlated with the actual QT interval [16]. Therefore this interval was used to compare the investigated methods influence on the analysis of the repolarization duration variability. The second reason to replace the QRS onset by the R wave peak as the left limit was the necessity to allow investigation of the linear filtering (in Fig. 6 we can notice that the applied low-pass filtering makes practically impossible the QRS onset determination). The RTp series obtained after measurements performed on the “clean” and the contaminated ECG signals are presented in Fig. 7. The plots are supplemented by the indices of variability: standard deviation (SD) and root mean square of the successive differences (RMSSD)

   N−1 2  (x(n + 1) − x(n)) RMSSDx =

Fig. 5 – Results of RPFTWEB applied without prefiltering for baseline wander suppression in two cases of different noise level (A–C and D–F, respectively).

• For RPFTWEB the large errors occurred for 10 and for 5 dB but their number was rather low. The RMSE(SNR) curves presented in the upper part of Fig. 6 confirm the above statements. Application of linear filtering resulted in an approximate shift of the curves obtained for Tp and Te by about 5 dB (to the left with respect to the curves obtained when the measurements were not preceded by any type of prefiltering). PFTAB resulted in the shift of about 5–10 dB, and RPFTWEB of about 15 dB. In the lower part of Fig. 6 the shifts are not so great, but comparing the respective plots we should take into account that only a very low number of the errors were rejected for RPFTWEB, slightly more for PFTAB, significantly more for the linear filter, and still much more for the original measuring procedure not preceded by any of the considered filters. We may, thus, conclude that the projective filtering, particularly the new version of the method (RPFTWEB), allows not only to prevent the measuring procedure from large errors but significantly raises the measurements precision as well. Moreover, contrary to the applied low-pass filter, projective filtering improves the analysis of not

n=1

N−1

.

(12)

where x(n), n = 1, 2, . . . , N is the analyzed time series. The series obtained on the basis of the “clean” signal are very similar irrespective of the applied filtering. Moreover, the indices of their variability also remained almost the same. To achieve it, for both versions of the projective filtering, the relatively high dimensions of signal subspaces were applied (those presented in the lower part of Table 1). When the desired signal was contaminated with EMG noise, the original measuring procedure gave poor results. Even for a relatively low level of noise (SNR = 20 dB) the series is rather jagged which is signalled by the higher value of RMSSD. With the growth of the noise level, the series is more and more jagged. Linear filtering improved the measurements a little. Application of PFTAB allowed to obtain similar RTp series down to about 15 dB. Application of RPFTWEB was the most effective. For SNR = 10 dB the series was very similar to that obtained for the clean signal and even for 5 dB there were no large measurement errors. The small arrows over the lowest subplots in Fig. 7 indicate one of the RT intervals measured after the respective types of preprocessing. The corresponding analyzed signal segments are presented in Fig. 8. We can notice that the large error that occurred when no filtering preceded the measurements was caused by high energy noise far apart from the T wave. Linear low-pass filtering prevented the measuring procedure from this error; however, the precision of the measurement was very low. PFTAB decreased the error but it was RPFTWEB that allowed the highest precision of the T wave location. Compar-

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

169

Fig. 6 – The root mean-squared error (RMSE) of the QRSo , Tp and Te determination as a function of the SNR. The respective curves show the results obtained when the measurements were preceded with no filtering (NF), linear filtering (LF) and projective filtering (PFTAB or RPFTWEB). The upper part of the figure presents the curves obtained when all errors were taken into account during calculations of the RMSE values. When large errors were rejected before the calculations, the curves presented in the lower subplots were obtained.

Fig. 7 – The RTp series obtained on the basis of the “desired” ECG signal (the uppermost subplots) and on the basis of the signals contaminated with EMG noise when different filters preceded the measurements. Each series is characterized by two indices of variability (RMSSD and SD). The cases indicated by arrows over the lowest subplots are illustrated in Fig. 8.

170

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

Fig. 8 – The signal segments corresponding to the beat number indicated by the arrows in Fig. 7. The long vertical line shows the Tp position obtained for s(n), the shorter ones the positions obtained for x(n), or x (n) after the respective types of preprocessing. The errors are marked with arrows.

ing the ECG cycle obtained after linear filtering with the ones after projective filtering (either PFTAB or RPFTWEB) we can see that the linear operation caused very significant distortions of the desired ECG morphology. None of the applied versions of the projective filtering caused such disadvantageous effects. The methods preserved not only low but also high frequency components of the processed signal.

4.3.1.

Spectral analysis of RT intervals variability

The influence of RPFTWEB on spectral analysis of the RT intervals variability was illustrated in Fig. 9. The time series presented on the left side of the figure were obtained on the basis of a signal of relatively high SNR. To estimate the power spectral density, each series was resampled with the frequency of 1 Hz (by application of first order interpolation). Then they were filtered to remove the low frequency trend

(the Butterworth high-pass filter of the fourth order and a cutoff frequency of 0.02 Hz was employed—the filtering operation was performed in forward and reverse direction). After such preprocessing the power spectral density estimation was performed by application of the autoregressive modelling [27](employing the Levinson–Durbin algorithm). The same model order equal to 10 was chosen to estimate both PSDRR and PSDRT (although according to the Akaike criterion [28] a lower order was sufficient for PSDRR ). The RR and RT intervals series often contain [26] the very low frequency trend (below 0.03 Hz), and two oscillatory components: of a low frequency—close to 0.1 Hz (LF) and a high frequency—close to 0.25 Hz (HF). The oscillations are related [26] with the slow waves of the arterial pressure (LF) and with respiration (HF). After suppression of low frequency trend the two distinct peaks, corresponding to the LF and HF compo-

Table 2 – The percentage of large measurement errors (greater than 10 ms for QRSo and than 30 ms for Tp and Te ) SNR

No filtering QRSo

25 20 15 10 5

0.1 1.1 4.7 11.8 19.7

LF

PFTAB

Tp

Te

Tp

Te

0 1.0 3.8 9.1 18.7

0 1.0 3.8 9.4 19.9

0 0.1 1.3 4.9 11.9

0 0.1 1.3 5.1 12.4

RPFTWEB

QRSo

Tp

Te

QRSo

Tp

Te

0 0 0.2 1.1 5.8

0 0 0.3 1.4 5.2

0 0 0.3 1.4 5.2

0 0 0 0.1 1.1

0 0 0 0.2 0.9

0 0 0 0.2 1.1

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

171

Fig. 9 – Spectral analysis of RR and RT intervals variability. On the left (from the top): the RR series; the RT series obtained without projective filtering of the analyzed signal; the RT series obtained after RPFTWEB; the RT series obtained after RPFTWEB when the analyzed signal was contaminated with EMG noise. On the right—the corresponding power spectral densities.

nents, were obtained in PSDRR (the uppermost on the right side of the figure). All PSDRT densities presented below contain a distinct LF peak and also a discernible HF component. Whereas the RT series in subplot B was obtained without projective filtering of the analyzed ECG signal, determination of the series in subplot C was preceded by RPFTWEB. We can notice that projective filtering did not change discernibly the distribution of the obtained spectral components. However, the method slightly decreased their amplitude (compare subplots F and G). This can be explained mostly by two factors: a decrease of errors of T maximum determination, and possibly a slight decrease of the true variability of the analyzed signal morphology. However, owing to RPFTWEB, adding the EMG noise of a generally moderate, but highly variable level (the average SNR = 10 dB) did not spoil the RT series determination and the final results of the spectral estimation (compare subplots F–H). We can conclude that the robust projective filtering enables the spectral analysis of the RT series even in cases of rather noisy ECG signals.

5.

Conclusions

The nonlinear projective filtering is a method that allows suppression of noise whose power spectrum overlaps the spectrum of the desired component. Therefore the method can effectively be employed to the enhancement of the ECG signal contaminated by the wide-band electromyographic

noise. However, this type of noise is highly nonstationary. Very often after signal segments of relatively high quality, the noise level grows significantly. Such noisy segments of the processed signal can spoil construction of the signal subspaces by application of the classical principal component analysis. Application of the modified version of the projection pursuit based robust PCA allowed to overcome this problem to a certain extent. As a result, the projective filtering performance in real EMG environment was raised. The developed version of the method was successfully applied to ECG enhancement prior to the analysis of the ventricular repolarization. Compared to the classical linear filtering, the method was much more effective. It was also significantly better than PFTAB, the previous version of the projective filtering. RPFTWEB helped to obtain the repeatable measurements of the repolarization duration in a wide range of the signal-to-noise ratio. Even for rather noisy ECG signals, the achieved precision of the measurements allowed to accomplish successfully one of the most demanding problems of evaluating the beat-to-beat variability of the RT series. The high accuracy of the QT interval measurement proves the projective filtering ability to reconstruct precisely the desired ECG signal embedded in noise. Thus the method can also be helpful when the dynamical measurements of other indices quantifying the ventricular repolarization are considered. This, however, requires the further study.

172

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 2 ( 2 0 0 8 ) 161–172

Acknowledgment The author would like to thank the anonymous reviewers for their valuable comments and suggestions.

references

[1] O. Pahlm, L. Sörnmo, Data processing of exercise ECG’s, IEEE Trans. Biomed. Eng. 34 (1987) 158–165. [2] J.A. van Alsté, W. van Eck, O.E. Herrmann, ECG baseline wander reduction using linear phase filters, Comp. & Biomed. Res. 19 (1986) 417–427. [3] R. Wariar, C. Eswaran, Integer coefficient bandpass filter for the simultaneous removal of baseline wander, 50 and 100 Hz interference from the ECG, Med. Biol. Eng. Comput. 29 (1991) 333–336. [4] R. Jane, H. Rix, P. Caminal, P. Laguna, Alignment methods for averaging of high resolution cardiac signals: a comparative study of performance, IEEE Trans. Biomed. Eng. 38 (1991) 571–579. [5] E.G. Caiani, A. Porta, G. Baselli, M. Turiel, S. Muzzupappa, M. Pagiani, A. Malliani, S. Cerrutti, Analysis of cardiac left-ventricular volume based on time warping averaging, Med. Biol. Eng. Comput. 40 (2002) 225–233. [6] S. Olmos, J. Garcia, R. Jane, P. Laguna, ECG signal compression plus noise filtering with truncated orthogonal expansions, Signal Process. 79 (1999) 97–115. [7] M. Kotas, Application of robust projection pursuit based principal component analysis to ECG enhancement, Biomed. Signal Process. Control 1 (2006) 289–298. [8] I.T. Jollife, Principal Component Analysis, Springer Verlag, New York, 2002. [9] T. Schreiber, D. Kaplan, Nonlinear noise reduction for electrocardiograms, Chaos 6 (1996) 87–92. [10] H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, University Press, Cambridge, 2004. [11] X. Hu, V. Nenov, A single-lead ECG enhancement algorithm using a regularized data-driven filter, IEEE Trans. Biomed. Eng. 53 (2006) 347–351. [12] M. Kotas, Projective filtering of the time-aligned ECG beats, IEEE Trans. Biomed. Eng. 51 (2004) 1129–1139. [13] M. Kotas, Projective filtering of time-aligned ECG beats for repolarization duration measurement, Comp. Meth. Prog. Biomed. 85 (February 2007) 115–123.

[14] M. Kotas, Projective filtering of the time warped ECG beats, Comp. Biol. Med. 38 (January 2008) 127–137. [15] P.E. Tikkanen, L.C. Sellin, H.O. Kinnunen, H.V. Huikuri, Using simulated noise to define optimal QT intervals for computer analysis of ambulatory ECG, Med. Eng. Physics 21 (1999) 15–25. [16] M. Merri, M. Alberti, A.J. Moss, Dynamic analysis of ventricular repolarization duration from 24-hour Holter recordings, IEEE Trans. Biomed. Eng. 40 (1993) 1219–1225. [17] A. Porta, G. Baselli, F. Lombardi, S. Cerutti, R. Antolini, M. Del Greco, F. Ravelli, G. Nollo, Performance assessment of standard algorithms for dynamic R-T interval measurement: comparison between R-Tapex and R-Tend approach, Med. Biol. Eng. Comput. 36 (1998) 35–42. [18] L. Gupta, D.L. Molfese, R. Tammana, P.G. Simos, Nonlinear Alignment and Averaging for Estimating the Evoked Potential, IEEE Trans. Biomed. Eng. 43 (1996) 348–355. [19] H. Sakoe, S. Chiba, Dynamic programming optimization for spoken word recognition, IEEE Trans ASSP ASSP-26 (1978) 43–49. [20] P.S. Hamilton, W.J. Tompkins, Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database, IEEE Trans. Biomed. Eng. 33 (1986) 1157–1165. [21] C. Croux, A. Ruiz-Gazen, A fast algorithm for robust principal components based on projection pursuit, in: Compstat: Proceedings in Computational Statistics, 1996, pp. 211–217. [22] O. Hossjer, C. Croux, Generalizing univariate signed rank statistics for testing and estimating a multivariate location parameter, Nonparamet. Stat. 4 (1995) 293–308. [23] P.J. Rousseeuw, C. Croux, Alternatives to the median absolute deviation, J. Am. Stat. Assess. 88 (1993) 1273–1283. [24] P. Laguna, R.G. Mark, A. Goldberg, G.B. Moody, A database for evaluation of algorithms for measurement of QT and other waveform intervals in the ECG, Comp. Cardiol. 24 (1997) 673–676. [25] P. Laguna, N.V. Thakor, P. Caminal, R. Jane, H. Yoon, New algorithm for QT interval analysis in 24 h Holter ECG: performance and applications, Med. Biol. Eng. Comput. 28 (1990) 67–73. [26] F. Lombardi, G. Sandrone, A. Porta, D. Torzillo, G. Terranova, G. Baselli, S. Cerutti, A. Malliani, Spectral analysis of short term R-Tapex interval variability during sinus rhythm and fixed atrial rate, Eur. Heart J. 17 (1996) 769–778. [27] J.S. Lim, A.V. Oppenheim, Advanced Topics in Signal Processing, Prentice Hall, 1988. [28] H. Akaike, A new look at the statistical model identification, IEEE Trans. Autom. Contr. 19 (1974) 716–723.