Averaging of nonlinearly aligned signal cycles for noise suppression

Biomedical Signal Processing and Control 21 (2015) 157–168

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Averaging of nonlinearly aligned signal cycles for noise suppression Marian Kotas ∗ , Tomasz Pander, Jacek M. Leski Institute of Electronics, Silesian University of Technology, Akademicka 16, 44-100 Gliwice, Poland

a r t i c l e

i n f o

Article history: Received 28 March 2015 Received in revised form 19 May 2015 Accepted 8 June 2015 Available online 27 June 2015 Keywords: Nonlinear alignment Dynamic time warping Event-related potentials Evoked potentials ECG signal QT interval Noise suppression

a b s t r a c t Averaging of nonlinearly aligned (time-warped) signal cycles is an important method for suppressing noise of quasi-periodical or event related signals. However, in the paper we show that the operation of time warping introduces unfavorable violation of the requirements that should be satisfied for effective averaging and, as a result, it causes poor suppression of noise. To limit these effects, we redefine the matrix of the alignment costs. To improve results of averaging in cases of variable energy noise, we apply weighting of the summed signal samples. The derived formula gives smaller weights for more noisy signal cycles and this way limits their influence on the constructed template. The proposed modifications caused a significant increase of the Noise Reduction Factor (NRF) in the experiments on the simulated evoked potentials. Whereas the greatest NRF obtained by the reference methods in nonstationary white noise environment was equal to 1.55, for the new method proposed we achieved a value of 4.44. For non-stationary colored noise, the corresponding values were 1.44 and 2.99. Moreover, application of the developed method to ECG signal processing, prior to the measurements of the QT interval, significantly improved the measurements immunity to noise. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Dynamic time warping (DTW) is a technique that uses dynamic programming [3] to determine the best alignment of two time series (or sequences of vectors) and to provide a measure of their morphological similarity. This similarity measure was developed for spoken words recognition [29,22] and later applied in a variety of pattern recognition problems, for instance to ECG beats clustering [27], to knee joint sound classification [11], to sleep stage classification [18] and in [33] more generally to different biomedical time series clustering. In this paper, we focus on application of the DTW technique to noise suppression. Analysis of biomedical signals is often possible only after effective suppression of noise they are embedded in. Therefore development of methods dealing with this problem has for many decades been an objective of unceasing efforts. When the noise component frequency-band overlaps that of the desired component, suppression of noise is particularly difficult; however, when the processed signals are of repetitive shape, a relatively simple method of coherent averaging appears advantageous [28]. Since the early years of computerized electrocardiography, this method has successfully been applied to ECG signals enhancement [26] and even to fetal ECG extraction from the maternal

∗ Corresponding author. Tel.: +48 32 2372019; fax: +48 32 2372225. E-mail address: [email protected] (M. Kotas). http://dx.doi.org/10.1016/j.bspc.2015.06.003 1746-8094/© 2015 Elsevier Ltd. All rights reserved.

abdominal bio-electrical signals [10]. Particularly important is its application in exercise electrocardiography where measurements of the signal features (amplitudes and slopes of the respective ECG waves, as well as their time lengths) require prior suppression of high level muscular noise. To accomplish this requirement effectively, the method has for many years been modified and improved [2,17,21]. However, with the progress in the analysis and interpretation of electrocardiogram, the method intrinsic constraint of forming an average template of an ECG cycle and loosing the information on its variability appeared unacceptable. It stimulated the efforts to develop new methods of noise suppression, preserving the shapes of the individual cycles. Their reconstruction on the basis of orthogonal functions derived with the use of principal component analysis appeared rather advantageous [12] as well as application of a slightly similar method of projective filtering of time-aligned ECG beats [13]. However, these methods are less effective in cases of high QT interval variability, associated with a changing position within an ECG cycle of a relatively low T wave. Reconstruction of such signals with the use of projective filtering appeared much more precise and immune to noise when a kind of nonlinear alignment (time warping) of ECG cycles preceded the projections [14]. A novel approach, based on two-dimensional warping of the ECG signals, both along the time and the amplitude axes, was proposed in [30]. A comparative study of the approaches allowing for averaging of signals with different types of shape variability was performed in [4].

158

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

approach, the warping path is subject to the following constraints [29]: • Boundary conditions: i1 = j1 = 1, iK = Ny , jK = Nx , which force the warping path to start with the pair x(1), y(1) and to end with x(Nx ), y(Ny ) . • Continuity conditions: ik+1 − ik ≤ 1, jk+1 − jk ≤ 1, which prevent the elements of both signals from being omitted in the warping path. • Monotonicity conditions: ik+1 ≥ ik , jk+1 ≥ jk , which force the elements of both signals to occur in the warping path in a non-decreasing order.

Fig. 1. Graphical representation (A) of the allowed step directions corresponding to (3) and an example of a warping path {(i, j)} = {(1, 1), (2, 2), (2, 3), (3, 4), (3, 5), (4, 6), (5, 6), (6, 7)} (D) relating two different signals: y(n) (B) and x(n) (C).

A similar progress can be observed in the development of methods for processing of so-called evoked potentials. Such signals respond to a stimulus with a series of positive and negative deflections from the baseline. The time between a stimulus and a particular deflection is called as latency. The latencies convey important information on physiological mechanisms evolving in the brain. Unfortunately the evoked potentials are mixed with spontaneous EEG signal and can be of low signal-to-noise ratio. Different types of filtering techniques can be applied to raise this ratio [5,1,25,36]; in cases of a very high level of noise, however, they can rarely produce signals of sufficient quality to detect the important deflections and measure their amplitudes and latencies. Thus of high importance are the methods of averaging [6,31] and their modifications [37,20,7]. Application of nonlinear alignment of the evoked potentials, prior to their averaging, performed with the use of the technique of dynamic time warping, was proposed in [24,9]. However, the procedure of dynamic time warping does not distinguish the signal and the noise components of the aligned signals and therefore noise can cause their erroneous alignment [24]. To prevent this effect, trilinear modeling of evoked potentials was applied in [34] to filter single responses before aligning and averaging them. In this study, we present more details showing how unfavorable can be the influence of noise on the results of time warping and averaging. Our aim is to propose an approach limiting these undesired effects. The rest of the paper is organized as follows. Section 2 presents the classical approach to dynamic time warping, a few methods of nonlinearly aligned cycles averaging known from the literature, an approach limiting the unfavorable influence of DTW on averaging, and a new method of weighted averaging of nonlinearly aligned cycles. Section 3 shows the experiments on the simulated evoked potentials and on real ECG signals. Finally, conclusions are drawn in Section 4.

With the above constraints, the length of the warping path satisfies the conditions max {Ny , Nx } ≤ K < Ny + Nx . We search for the path that minimizes the total cost of the alignment Q =

K 

dik ,jk .

(2)

k=1

It can be found with the use of dynamic programming [3]. y ,j=Nx is conFirst a matrix of cumulative costs G = [gi,j ]i=N i,j=1 structed, according to the recursive definition [29] gi,j = di,j + min{gi−1,j−1 , gi−1,j , gi,j−1 }

(3)

where g1,1 = d1,1 , and g0,j =∞, gi,0 =∞ for all i, j (the pictorial representation of this definition, which shows the allowed stepdirections of the warping path, is presented in Fig. 1(A)). Then, starting from the upper right corner of G (i = Ny , j = Nx ), the warping path is searched for by backtracking along the allowed step-directions through the minima of the cumulative costs in G, until i = 1 and j = 1. While backtracking, the allowed step-directions are opposite to those presented in Fig. 1(A). The back tracking can easily be performed without any conditional operations if, during determination of the cumulative costs matrix, we store the directions taken to reach each (i, j) point of this matrix in an auxiliary matrix. Then starting from the upper right corner of this auxiliary matrix, we simply move back along the stored directions. Applied to time warping of the signals x(n) and y(n), this classical approach produced the warping path presented in Fig. 1(D). We can notice that to achieve nonlinear alignment of both signals, two elements of y(n) and one element of x(n) had to occur two times in the warping path. 2.2. Existing methods of averaging of nonlinearly aligned signal cycles Let us denote the signal cycles to be averaged as xl (n),

n = 1, 2, . . ., Nl , l = 1, 2, . . ., L,

(4)

2. Methods 2.1. Dynamic time warping Let us consider two time signals of possibly different length: x(n), n = 1, 2, . . ., Nx and y(n), n = 1, 2, . . ., Ny . To perform their noni=N ,j=Nx linear alignment, we calculate the matrix of costs D= [di,j ]i,j=1y whose classical definition is as follows: 2

di,j = (y(i) − x(j)) .

(1)

Each di,j corresponds to the alignment of y(i) and x(j). The warping path P of the signals y(n) and x(n) consists of the ordered pairs of time indices: P = {(ik , jk )|k = 1, 2, . . ., K}, which indicate the elements of the successive aligned pairs: y(ik ), x(jk ). In the classical

where L denotes the number of cycles; the lth cycle is of length Nl . A few methods have already been proposed for constructing a mean cycle (also called as a template) by averaging of the respective time-warped signal cycles [9,23]. 2.2.1. Averaging of successive nonlinearly aligned cycles In this method, proposed by Gupta et al. in [9], we take the first cycle x1 (n) as the first template t1 (n) and we align the template nonlinearly with the subsequent cycles; each time after the alignment, the template is updated. By warping the current (l − 1)th template with the next lth cycle, we obtain the lth warping path P(l) , consisting of the succes(l) (l) sive pairs: P(l) = {(ik , jk )| k = 1, 2, . . ., Kl }, whose number Kl (the

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

length of the warping path) becomes the length of the new updated template (l)

tl (k) =

(l)

(l − 1)tl−1 (ik ) + xl (jk ) l

,

k = 1, 2, . . ., Kl .

(5)

The method was abbreviated by the authors as NLAAF2, which stands for Nonlinear Alignment Averaging Filter 2 (as it was the second of the two NLAAF versions developed in [9]). In the course of adding the successive warped cycles to the template, it is gradually growing in length. Also the warping paths relating the respective cycles with the template change after each time-alignment. In our work, we not only create a template but primarily we use it to reconstruct the individual signal cycles. However, in order to be able to reconstruct them on the basis of the created template, we have to remember which positions within a template are associated with the individual samples of the respective cycles. To this end we have to gradually update this (z,l) information after each addition of a cycle. Let us denote as pk the time index of a sample from the zth cycle that is associated with the kth position within the lth template. To easily interpret this variable, we should take into account that the kth element of the new updated (lth) template is a mean of different samples from the (z,l) averaged l cycles (z = 1, . . ., l); pk is the number (index) of a sample from the zth cycle. In the first template t1 (k) = x1 (k), k = 1, 2, . . ., (1,1) K1 (K1 = N1 ), thus simply pk = k, k = 1, 2, . . ., K1 . After alignment of the first template with the second cycle, we obtain the warping (2) (2) path P(2) = {(ik , jk )|k = 1, 2, . . ., K2 } and according to (5) we form the second template t2 (k). Time indices of the samples from the 2nd cycle that are associated with the successive positions within (2,2) (2) this template are given by pk = jk , k = 1, 2, . . ., K2 . Since, however, the template length could have grown, the relations between the time indices of the samples from the first cycle (z = 1) and the positions within the new template (l = 2) must be updated in the (1,2) (1,1) following way: pk = p (2) , k = 1, 2, . . ., K2 . i

k

Generally, after warping the (l − 1)th template with the lth cycle we obtain the lth warping path P(l) and we form the lth template tl (k). Then the indices of the samples from the lth cycle are associated with the positions within this template in the following way: (l,l)

pk

(l)

= jk ,

k = 1, 2, . . ., Kl

(6)

and the time indices of the samples from the previous cycles are updated as follows: (z,l)

pk

(z,l−1)

= p (l) i

, k = 1, 2, . . ., Kl ;

1 ≤ z < l.

(7)

k

After adding the last Lth cycle and using (6) and (7), we obtain the (l,L) final information pk , k = 1, 2, . . ., KL , l = 1, 2, . . ., L, on the association between the samples from the respective cycles and the positions within the final (Lth) constructed template. To reconstruct the nth sample of the lth cycle, we find the positions within the final Lth template, it is associated with: l (n) = (l,L) {k|pk = n} and we calculate the mean of the corresponding values xl (n) =

1 |l (n)|



tL (k),

(8)

k∈l (n)

where | · | denotes cardinality of a set. 2.2.2. Pairwise addition of nonlinearly aligned cycles This method was also proposed by Gupta et al. in [9]. It is based on pairwise averaging of the time-warped signal cycles. In this approach it is assumed that the number L of the cycles is a power of 2. In the first averaging stage, the cycles are grouped into L/2 of the successive pairs {x1 (n), x2 (n)}, {x3 (n), x4 (n)}, . . ., {xL−1 (n), xL (n)}.

159

The paired cycles undergo time warping and averaging, and L/2 mean cycles are obtained. Then the obtained mean cycles are again grouped into L/4 of the successive pairs and undergo the same operations (time warping and averaging). This operation is continued until one mean cycle has been constructed. Each time the warping path has been determined, the information on the positions of the respective samples of the warped cycles is updated and, finally, after one mean cycle has been formed, we reconstruct the individual signal cycles, similarly as in the previous method. The method was abbreviated by the authors as NLAAF1, as it was the first version of Nonlinear Alignment Averaging Filter proposed in [9]). 2.2.3. Barrycenter averaging of nonlinearly aligned cycles In this method, proposed in [23], we randomly select one of the L cycles as the initial form of the constructed template t(n) and then we align nonlinearly all cycles with this template, obtaining L warping paths. Then the operation called by the authors as the DTW Barycenter Averaging (DBA) is performed: each sample of t(n) is updated as the average of those samples of the respective L cycles that were aligned with the previous version of t(n). The operation is repeated the assumed number of times and then the individual cycles are reconstructed on the basis of t(n) and the final warping paths. 2.3. Weighted averaging of nonlinearly aligned cycles (WANAC) In order to improve suppression of nonstationary noise of variable energy, we propose to apply weighted averaging of nonlinearly aligned signal cycles. The method is slightly similar to the method of barrycenter averaging. First we create an initial form of a template t(n). Then we align it nonlinearly with all cycles to be averaged, we calculate the weights that depend on the costs of the respective cycles alignment with the template, and we perform weighted averaging to update the template. These operations are performed repeatedly until the specified criterion has been satisfied. Applying a vector notation, let w = [w1 , w2 , . . ., wL ] represent a vector of real nonnegative weights associated with the respective L signal cycles and let t = [t(1), t(2), . . ., t(Nt )] denote the constructed weighted average of the nonlinearly aligned signal cycles. Since the signal cycles can be of different length, for them we preserve the previous notation: xl (n), n = 1, 2, . . ., Nl , l = 1, 2, . . ., L. The set of warping paths relating t with the respective cycles is denoted as P = {P(l) | l = 1, 2, . . ., L}. We define a criterion function that is based on a weighted sum of the cumulative alignment costs J(w, t, P) =

L 

(wl )m G(t, xl ),

(9)

l=1

where m is a weighting exponent m ∈ (1, ∞), G(· , ·) is a cumulative cost of nonlinear alignment

G(t, xl ) =

Kl  

(l)

(l)

2

t(ik ) − xl (jk )

k=1

=

Kl 

dik ,jk .

(10)

k=1

Criterion (9) can be interpreted as a measure of the total dissimilarity between t and the set of cycles: xl , l = 1, 2, . . ., L. Determination of an optimal weighting vector w* , an optimal “warped mean” t* and an optimal set of warping paths P∗ can be achieved by minimization of the defined criterion function J(w∗ , t∗ , P∗ ) = min J(w, t, P). w,t,P

(11)

160

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168 (l)

(l)

Since this minimization is very complex, we propose to solve iteratively a few simpler problems, performing a kind of alternated optimization.

of k for which ik = n in sets l (n) = {k|ik = n}, we can rewrite (20) as

• For a fixed template t and any fixed vector of weights w, criterion function (9) reaches its minimum when for each cycle the cumulative cost of its alignment with t is minimal. These minimal cumulative costs can be calculated separately for all cycles with the use of dynamic time warping. This way we obtain L warping paths, i.e. set P and L cumulative costs of alignment

J(w, t, P) =

gl = G(t, xl ),

l = 1, 2, . . ., L.

(12)

• Having the cumulative costs (12) computed and assuming that they are fixed, we can proceed with the determination of the optimal weights w* . However, we can see that criterion function (9) descends to zero when all weights are equal to zero. This is a trivial solution. To exclude it, we should minimize (9) preserving the linear constraint L 

wl = c,

c > 0.

(13)

l=1

To derive the appropriate formula, we define the Lagrangian of (9) using (12) and constraint (13) LJ (w, ) =

L 

(wl )m gl − 

 L 

l=1

 wl − c

.

(14)

l=1

Setting the Lagrangian gradient to zero, we obtain

∂LJ (w, )  wl − c = 0 = ∂ L

(15)

l=1

and

∂LJ (w, ) = m(ws )m−1 gs −  = 0, ∂ws

s = 1, 2, . . ., L

(16)

From (16) we get ws =

  1/(m−1) m

m

=

(17)

L

c

1/(1−m) g l=1 l

.

(18)

Combining (17) and (18) we get ws =

c

L  gl 1/(1−m) . l=1

(19)

gs

• Having the weights computed according to (19) and assuming that they are fixed, we can derive formula for an updated mean vector t by differentiating (9) with respect to the elements of t. Let us rewrite the definition of the criterion function J(w, t, P) =

L  l=1

(wl )m

Nt  

(l)

2

t(n) − xl (jk )

.

(21)

n=1 k∈l (n)

l=1

This way we have obtained the criterion which is an explicit function of t(n). We can now set the gradient of (21) with respect to t(n) to zero and we get L    ∂J(w, t, P)  (l) (wl )m = 2 t(n) − xl (jk ) = 0. ∂t(n) l=1

(22)

k∈l (n)

From (22) we get the formula

L t(n) =

l=1

(wl )m

L

l=1



(l) x (j ) k∈l (n) l k

(wl )m |l (n)|

,

(23)

needed to update the elements of the template (as we can see, because of the denominator that acts as a normalizing factor, constant c from formula (19) does not influence (23) – it is reduced and therefore its value is immaterial; during computations we can set c = 1). The initial form of a template t(0) is created in the following way. Like in NLAAF1, we group the successive cycles in pairs (if the number of cycles is odd, we neglect the last one). Then we perform nonlinear alignment and averaging within each pair and we choose the pair for which the total cost of the alignment was lowest. The “warped mean’ of the chosen pair of cycles is used as the initial template. On the basis of the above considerations, we obtain the following iterative algorithm: 1. Initialize t(0) . Fix m > 1, c > 0,  > 0, rmax »1. Set the iteration index r = 1. 2. Align nonlinearly t(r−1) with respect to the individual signal cycles: determine L warping paths minimizing (10) and L cumu(r) lative costs of alignment gl , l = 1, 2, . . ., L, defined by (12). (r)

1/(1−m) gs .

From (17) and (15) we get

  1/(m−1)

L 

(wl )m

Kl  

(l)

(l)

2

t(ik ) − xl (jk )

.

(20)

k=1

We know that each sample t(n) appears in this criterion at least L times (once for each cycle, the template is aligned with); some samples appear more times. It depends on the results of (l) time-warping: t(n) occurs for all ik = n. If we gather the values

3. Calculate elements ws of the weights vector w(r) using (19) and (r) the calculated alignment costs gl .

4. Update elements of template t(r) using (23) and weights vector w(r)

.



5. If t(r) − t(r−1) / t(r−1) >  and r < rmax then r ← r + 1 and go to (2) else stop.

After the final form of the template has been obtained, we can reconstruct the individual cycles using the final warping paths relating the template with the respective cycles. To this end for each sample of each cycle we find the samples (of the template) it was associated with, and we calculate their mean. Remarks According to the described algorithm, the constructed template is nonlinearly aligned with each individual signal cycle at most rmax times. It means that the proposed method computational costs are increased approximately rmax times (at most) with respect to the NLAAF2 method, which requires only a single alignment of the template with each signal cycle. In the experimental section, we use rmax = 20. With c = 1 Eq. (19) can be rewritten as

 ws =

L    gs 1/(m−1) l=1

gl

(−1) .

(24)

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

We can notice that

lim ws =

m→1+

⎧ 1, gs ≤ gl ⎪ ⎪ ⎨

for all l;

(25)

⎪ ⎪ ⎩ 0, otherwise

what means that a trivial solution is obtained with only one nonzero weight, corresponding to the cycle for which gs is the smallest. For m tending to infinity all weights tend to the same value: ws = 1/L. In such a case the WANAC method would operate similarly to the method of barycenter averaging. Thus, with the increase of m the weighted averaging is less and less selective. In our experiments, we use m = 2. 2.4. Redefinition of a cost matrix

2.4.1. Conditions for effective averaging of linearly aligned signal cycles Various conditions should be satisfied to assure effective suppression of noise and enhancement of the desired component as a result of averaging. Among the most important are the following two. Firstly, all the cycles that undergo averaging must contain the same deterministic desired component (xl (n) = s(n) + l (n), l = 1, 2, . . ., L). Secondly, the noise l (n) must be of zero mean (its expected value equal to zero: E(l (n)) = 0, where E( · ) denotes the expectation operator). Fulfillment of these requirements assures that by averaging the aligned cycles, we obtain an unbiased estimate of the desired component:



E



1 xl (n) L L

l=1



=E





1 (s(n) + l (n)) L L

l=1



1 = E s(n) + l (n) L L

the fact that these samples are matched by minimization of the squared differences between them). By averaging of the matched pairs of samples, the noise components of the same polarity will not be suppressed much. Thus the template created by averaging the first two nonlinearly aligned cycles will be corrupted at the positions of the considered corrupted pairs of samples. Then alignment with the third and the subsequent cycles is likely to give similar effects. As a result, we can expect that after the nonlinear alignment of the cycles under averaging, the assumption concerning the zero mean of noise of the samples to be averaged will not be preserved. Let us now focus on the desired components of the samples to be averaged. We assume that the desired components of the processed signals were formed by different deformation of the time axis of the same deterministic signal s(): xl () = s(tl ()) + l ()

Averaging of time aligned signal cycles has long been applied to increase SNR of quasi-periodical signals [26,28]. In this technique, we first detect so-called fiducial points within all cycles to be averaged, then we perform their time-alignment with respect to these points, and finally their averaging. Contrary to the operation proposed in this paper, linear time-alignment is performed: the time axis of each signal cycle is preserved without any deformation, the cycles are only shifted in time.

= s(n).

(26)

l=1

Detectors of fiducial points are developed so as to match to the shape of the desired component and are usually immune to noise. Therefore location of the detected points is scarcely influenced by noise, and the operation of time-alignment does not violate much the specified conditions (apart from causing the well known filtering effects [28] in cases of synchronization errors). More severe can be the effects introduced by the operation of nonlinear alignment, described in Section 2. The ability to repeat any sample of the time-warped signals the requisite number of times, to obtain good matching between two nonlinearly aligned signal cycles, unavoidably leads to violation of the two crucial conditions. 2.4.2. The influence of nonlinear alignment on the conditions for effective averaging Let us first consider if after the nonlinear alignment of the cycles we can still assume the zero mean of the noise components of the samples to be averaged. For instance in the NLAAF2 method: if a sample of the first cycle is heavily corrupted with noise, it is likely to be aligned with one or more samples of the second cycle that are also corrupted with noise of the same polarity (which results from

161

(27)

where  denotes a continuous time variable, tl () is a function deforming time axis of the lth signal (cycle). Unfortunately, we do not know those deforming tl (·) functions and we only try to compensate their effects with the use of dynamic time warping. However, by minimizing the squared differences between the corrupted signal samples, we do not assure that their desired components are the same. Thus for noisy signals even this important condition is not satisfied. Let us analyze the effects of time warping of two noisy cycles. We select pairs (ik , jk ) minimizing the sum of the squared differ2 ences (y(ik ) − x(jk )) for k = 1, 2, . . ., K, preserving all the conditions specified in Section 2. If for some k we obtain y(ik ) = x(jk ), we cannot infer that the desired components sy (ik ) and sx (jk ) are equal. We only know that sy (ik ) + y (ik ) = sx (jk ) + x (jk ) and thus sy (ik ) − sx (jk ) = x (jk ) − y (ik ). For noise of a considerable level, the sum on the right is very rarely small; we can thus conclude that the desired components: sy (ik ) and sx (jk ) are rarely even approximately equal. What worse, averaging these equal samples, we do not change their values and consequently we leave the noise components unchanged. This shows that the operation of time warping, aimed to minimize (2) with the costs defined by (1) and the conditions specified in Section 2 can lead to a poor suppression of noise. 2.4.3. The modified costs of alignment To prevent the described unfavorable effects of time warping, we propose to apply the following definition of cost matrix elements:





 di,j = x(i) − y(j)

(28)

where · denotes the Euclidean norm, and vectors x(n) and y(n) are defined as follows: x(n) = [x(n − v), x(n − v + 1), . . ., x(n), . . ., x(n + v)]T

(29)

which means that they are composed of 2v + 1 successive signal samples (with the nth sample being the central one). By such a redefinition of a cost matrix, we assure that longer signal intervals are matched to each other instead of single samples of warped signals. With this definition, however, we have to overcome some difficulties incurred. As specified in Section 2, the operation of dynamic time warping is intended to align two time signals of possibly different length: x(n), n = 1, 2, . . ., Nx and y(n), n = 1, 2, . . ., Ny . Applying (28) to calculate d1,1 , we have to use x(n) and y(n), −v < n ≤ v + 1, and to calculate dNy ,Nx we need to know x(n) for Nx − v ≤ n ≤ Nx + v and y(n) for Ny − v ≤ n ≤ Ny + v, thus we have to know v samples preceding and v samples succeeding the time-warped signal cycles. Fortunately, many practical applications of high significance are not excluded by this requirement. For

162

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

instance such signals as evoked potentials or electrocardiograms are recorded continuously and the signal cycles are selected from these longer records on the basis of some additional information, like the moments when the potentials were evoked. The second difficulty is that the created template should also be longer (it should also contain v samples preceding and v succeeding its main part). To construct such longer templates, we assume that for each cycle under averaging these beginning v and ending v samples are preserved without any time warping. This way the beginning part of the constructed template is simply the mean (or the weighted average in the case of the new method proposed) of these beginning parts of the respective cycles, and similarly is it with the ending parts. When the WANAC method is applied with the modified definition of a cost matrix, this definition should also replace the classical one in (10). 2.4.4. The alternative definition of the modified alignment costs Costs (28) can be obtained from costs (1) with the use of a kind of filtering

  z=v 



2  (i) (j)

di,j = y − x = (y(i + z) − x(j + z))

In the first stage of ECG signal processing, we perform QRS complexes detection. This way we obtain so-called fiducial points rl , l = 1, 2, . . ., L, corresponding to the central position within the respective detected L complexes. Then we cut the signal into the time segments beginning before the detected complexes and ending before the successive ones: rl −  ≤ n < rl+1 − , l = 1, 2, . . ., L − 1, where  denotes the assumed shift of the segments onsets with respect to the determined fiducial points (for the sampling frequency of 1000 Hz, we used  = 25). These signal segments undergo time warping, averaging and reconstruction with the use of the described methods. 3. Numerical experiments and discussion In this section, we will investigate the presented methods performance when applied to the simulated evoked potentials and to real ECG signals. 3.1. Evaluation of the results Applying the methods to process a set of simulated evoked potentials, we can asses the overall results with the use of the general noise reduction factor

z=−v

  z=v  = di+z,j+z ,

2.5. Application to ECG processing

(30)

z=−v

i.e. by adding 2v + 1 elements of a cost matrix based on (1) and finally taking the square root of the sum. This filtering is performed along the direction parallel to the main diagonal of this matrix (i = j). In fact, if we calculated cost matrix based on definition (1) for extended signal cycles y(n), −v < n ≤ Ny + v and x(n), −v < n ≤ Nx + v, we could obtain the cost matrix defined by (28) with the use of such filtering. Remarks It should be emphasized that the described filtering of the cost matrix is not equivalent to the low-pass filtering of the signals before the alignment costs determination. The latter approach could only be applied in cases of smooth signals. Application of definition (28) is not limited to such cases only. What can affect its accuracy of the warping path determination is the irregularity of the true path to be estimated. With the use of (28) such path cannot be recovered precisely (and, actually, using (1) we can achieve this goal in noiseless cases only). Since, however, there are not any physiological basis of such irregularity for most biomedical event related signals, the requirement of the warping path smoothness does not seem to limit much the area of practical applications of the alignment costs proposed. To determine each entry of the cost matrix using the norm function (28), we have to calculate 2v + 1 differences, 2v + 1 operations to square these differences, 2v sums and 1 square root. The number of these operations is approximately 2v times larger than the number required if the classical definition of alignment costs is applied. However, determination of the cumulative costs matrix and the operation of backtracking are unmodified. Since these operation are much more time consuming, particularly in the MATLAB environment, the overall growth of the computational complexity of the operation of time warping is less severe. In our experiments, the execution times of the unoptimized MATLAB function performing modified costs based DTW increased up to about 4 times (for 2v = 200, Nx = Ny = 400) with respect to the classical approach.

gnrf

  l=L n=Nl 2  (xl (n) − sl (n)) =   l=1 n=1  2 l=L n=Nl  n=1

l=1

(31)

xl (n) − sl (n)

where xl (n) denotes the lth signal cycle (of lenght Nl ), sl (n) is its desired component and xl (n), the reconstructed cycle. When more detailed analysis of the results is required, an individual noise reduction factor can be calculated for each cycle

(l) nrf

  n=Nl 2  (xl (n) − sl (n)) =   n=1  . 2 n=Nl 

(32)

xl (n) − sl (n)

n=1

Applying the methods to ECG signals, we will investigate their influence on the precision of the QT interval measurement. Since the measurement errors increase the QT series variability, which is of diagnostic significance [32,19], we will calculate two basic indices of this variability: the standard deviation of the QT time series, and the root mean square of the successive differences (RMSSD [32]), defined as follows:

 RMSSDx =

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

N−1

.

(33)

where {x(n) | n = 1, 2, . . ., N} is the analyzed time series. 3.2. Averaging of nonlinearly aligned evoked potentials Evoked potentials can be regarded as a kind of signal cycles. In this section, we make an experiment on the data simulating such potentials. 3.2.1. Test signals To test the methods described, we have created a kind of a small database. It consists of 5 sets of simulated evoked potentials that can be used to find the proper values of the methods parameters and of 20 sets of signals that can be exploited in the final tests.

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

163

Fig. 2. Visualization of the way the test signals are constructed: (A) the shape function defined by (34); (B, C) examples of the functions (35) used to deform the time axis; (D,E) the discrete time signals (a samples corresponding to time variable 0 < t ≤ 1) constructed according to (36).

To simulate the desired components of the constructed signals, we used the shape function proposed in [35]: 2

2

2

f (z) = 15e−20(z−0.7) − 0.5e−50(z−0.45) + 0.6e−100(z−0.3) 2

2

− 0.6e−150(z−0.2) + 0.5e−200(z−0.15) ,

(34)

and to deform the time axis, we applied the second order polynomial function also proposed in [35]: 2

tb (z) = z + bz − bz ,

(35)

which satisfies the conditions t(0) = 0 and t(1) = 1; thus, by using different values of parameter b, we can obtain more or less severe deformation of a time axis in the interval [0, 1], preserving its border values. The shape function is presented in Fig. 2(A). For independent variable t in the range [0, 1] it resembles visual evoked potentials following the target stimuli [25], with the largest peak similar to the so-called P300 wave. Beyond interval [0, 1] the function quickly descends to zero. The desired components of the signals were generated according to the formula s(n) = f

  n  tb

a

,

n = −99, −98, . . ., Ns + 100,

(36)

where parameter a specifies the number of time samples within the interval (0, 1] and Ns > a is the number of samples of the main part of the signal (which undergoes time warping). With the assumed sampling rate of 500 Hz and Ns = 400, we obtained signals of 0.8 s length (plus the preceding and succeeding parts of 0.2 s). The way the signals were constructed is illustrated in Fig. 2. For different values of parameters a and b different deformations of the time axis were caused and different final forms of these signals were obtained. To show precisely relation between shape function (34) and the constructed signals, only their parts corresponding to t ∈ (0, 1] are presented, but obviously, according to formula (36), the generated signals extend beyond this time interval. Part A of the created database consists of 5 sets of 16 simulated evoked potentials and corresponding noise components (white or colored autoregressive, simulating the spontaneous EEG signal). To create the desired signals, random values of parameters a ∈ N and b ∈ R were generated, with 275 ≤ a ≤ 325 and −0.2 ≤ b ≤ 0.2, both with uniform distribution (and a rounded afterwards). With such values of these parameters, the latency of the simulated P300 wave varied from about 400 to about 500 ms, which is typical for visual evoked potentials triggered by target stimuli [25]. This part of the database was used to select the proper value of parameter v,

deciding on the length (2v + 1) of the vectors used for determination of the modified costs of alignment (28). Part B of the database contains 20 sets of simulated evoked potentials, generated in a slightly different way. First, for each set of signals additional variables ra and rb , deciding on the range of a and b, were randomly established, with 0 ≤ ra ≤ 100 and 0 ≤ b ≤ 0.5. Then for each signal of the created set, the values of a and b were generated, with 250 ≤ a ≤ 250 + ra and −rb ≤ b ≤ rb . Finally the signals were created according to formula (36). By this approach we obtained sets of different properties, depending on the values of ra and rb . For large values of these parameters, we simulated evoked potentials with great differences of their time evolution within a set; for their small values the signals were more similar. 3.2.2. Selection of the signal vectors length For this purpose, we used part A of the database. The noisy test signals (evoked potentials) were formed in the following way: xl (n) = sl (n) + cl cel (n)

(37)

where sl (n) denotes the desired component, generated according to (36); el (n) is the noise component (either white gaussian or colored autoregressive). Noise components el (n), l = 1, 2, . . ., L had the same RMS value and thus the coefficients c and cl decided on the obtained signal to noise ratio. We created the signals with 4 types of noise: stationary white gaussian (SWN), nonstationary white gaussian (NWN), stationary colored (SCN) and nonstationary colored (NCN). In stationary types of noise cl = 1 and the value of c assures the general SNR in the created set of evoked potentials. To obtain nonstationary types of noise (with various RMS in the respective signal cycles), we generated random values of coefficients cl , with uniform distribution in the range 0.1 ≤ cl ≤ 1. We applied the NLAAF2m method (subscript m tells that the method was based on the modified alignment costs) to suppress noise in all sets of signals and for each set the general noise reduction factor ( gnrf ) was computed. For each type of noise, we computed the mean of the  gnrf values obtained for the respective five test sets. These values are presented in Fig. 3 as functions of v. Comparing the results obtained for different types of noise, we can notice that colored noise was more difficult to be suppressed than the white one. It is caused by more disadvantageous influence of this type of noise on the determination of warping paths. For SNR = 5 dB  gnrf grows with the growth of v, thus we can notice that application of longer signal vectors makes the method

164

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

Fig. 3. Dependance of the general noise reduction factor on the applied length (2v + 1) of the signal vectors used by NLAAF2m for cost matrix determination: results obtained for different types and levels of noise and the mean of all 12 cases are presented.

more immune to high level of noise. On the other hand, for lower level of noise (SNR = 20 dB) accurate reconstruction of the signal cycles is more crucial and smaller values of v are more favorable. The mean  gnrf reaches the highest value for v = 55 and it seems to be a good compromise between these contradictory requirements (although in different noise environment different values of this parameter can be more favorable). Therefore, we apply this value in the final tests, for all methods based on the modified version of a cost matrix. 3.2.3. Comparison of the methods This experiment was based on part B of the created database. The test signals were formed as in the previous subsection. For different methods and different noise conditions a mean of  gnrf values obtained in 20 test sets was calculated (see Table 1). The most striking difference is between the left three and the right three columns. It shows high significance of the introduced modification of the cost matrix definition. Only after this modification, nonlinear alignment did not cause such disastrous violation of the conditions for effective averaging, and the mean  gnrf was significantly increased. Comparing the results obtained for different types of noise, we can notice that generally colored noise is more troublesome to deal with. Higher correlation between the noise components of the signals to be aligned results in a greater impact of these components on the results of alignment. The important difference between the methods developed by Gupta et al. [9] (NLAAF1 and NLAAF2) and the other two is that these latter perform nonlinear alignment and construction of the template repeatedly. After creating the first form of this template, the operation is repeated many times and the template is refined. When the classical definition of a cost matrix is applied, this causes unexpected effects. Instead of achieving better suppression of noise, the procedure is less effective. To explain this effect, we must again refer to the unfavorable results of nonlinear alignment; repeating this operation, we introduce more and more precise matching of the template to the individual noisy signal cycles and more severe violation of the conditions for effective averaging. As a result, suppression of noise is less effective (observe the results obtained for SWN or SCN when classical time warping was applied). Only for nonstationary types of noise, weighted averaging slightly compensated these unfavorable effects, and for NWN the WANACm method achieved the best results. When the modified costs of alignment were applied, these unfavorable effects of repeated averaging were limited and did

not cause such deterioration of the DBAm and WANACm methods performance. Moreover, when nonstationary types of noise were processed, the repeated weighted averaging (WANACm ) caused significant increase of the general noise reduction factor  gnrf . In Fig. 4, we have visually presented exemplary results of averaging. To this end, we chose the signals with nonstationary colored noise of high level (NCN, SNR = 5 dB). According to Table 1, among the three reference methods based on the traditional costs of alignment (NLAAF1c , NLAAF2c and DBAc ), the NLAAF2c method appeared most effective for this type and level of noise. For this method, we selected the set of signals in which it achieved the highest value of  gnrf . Then we computed the individual noise reduction (l)

factor (nrf ) for each evoked potential and on this basis, we selected the EP with the best results of signal enhancement. Similarly we selected the EP with the best results obtained with NLAAF2m , and the corresponding EP for the WANACm method. We can notice that all three methods achieved the best signal (10) enhancement for the same EP (x10 (n)). The obtained values of nrf emphasize the advantageous results of the introduced modifications of the averaging procedure. Modified costs of alignment raised (10) nrf achieved by NLAAF2 from 2.37 (NLAAF2c ) to 6.10 (NLAAF2m ), and with repeated weighted averaging (WANACm ) the value 9.89 was obtained (and a very effective reconstruction of the desired component of the evoked potential, indeed).

3.3. Application to ECG processing for more reliable assessment of QT interval variability The QT interval reflects the time between depolarization and repolarization of the heart ventricles. It is an important electrocardiographic parameter, often used to quantify the duration of ventricular repolarization [19]. Assessment of the repolarization duration variability requires very precise measurements of the onset and the offset of this interval in a number of the subsequent ECG beats. In this subsection, we want to show how the investigated methods can help to raise the accuracy of these measurements. The experiment was based on ECG signals selected in [13] from the QT database [15] of the Physionet resources [8]. Since the signals were aimed to simulate the desired ECG, they had to be of very high SNR. To this end, we had carefully performed [13] visual inspection of the beginning 80 s segments of the first channel of many records. This way we selected a few signals (of the assumed length of 80 s, sufficient to implement the investigated methods effectively) of very high quality with a negligible level of noise.

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

165

Table 1 Results of signals enhancement (mean  gnrf ) in part B of the database. NLAAF stands for Nonlinear Alignment Averaging Filter, DBA for DTW Barrycenter Averaging and WANAC for Weighted Averaging of Nonlinearly Aligned Cycles. Method

Time warping Classical cost

Modified cost

SNR [dB] 5

10

20

5

10

20

NLAAF2 NLAAF1 DBA WANAC

Stationary white noise (SWN) 1.31 1.32 1.34 1.34 1.24 1.24 1.12 1.19

1.32 1.33 1.25 1.24

2.64 2.70 2.67 2.65

2.61 2.65 2.61 2.61

2.02 2.10 2.13 2.09

NLAAF2 NLAAF1 DBA WANAC

Non-stationary white noise (NWN) 1.55 1.55 1.55 1.54 1.37 1.37 1.86 1.83

1.51 1.50 1.37 1.84

2.78 2.81 2.79 4.44

2.68 2.75 2.71 3.99

2.07 2.14 2.13 2.54

NLAAF2 NLAAF1 DBA WANAC

Stationary colored noise (SCN) 1.32 1.26 1.31 1.25 1.10 1.09 1.20 1.13

1.16 1.16 1.06 1.08

2.18 2.28 2.30 2.16

2.14 2.20 2.24 2.11

1.70 1.79 1.81 1.74

NLAAF2 NLAAF1 DBA WANAC

Non-stationary colored noise (NCN) 1.44 1.34 1.42 1.33 1.13 1.12 1.43 1.35

1.21 1.20 1.08 1.18

2.25 2.36 2.36 2.99

2.15 2.24 2.28 2.66

1.70 1.78 1.83 1.93

They included traces of healthy: sel16483, sel16786, sel16795 as well as of diseased subjects: sel808 (supraventricular arrhythmia), sele0122 (coronary artery disease) and sele0170 (coronary artery disease). These ‘desired’ signals were additively contaminated with the EMG noise from the MIT-BIH database. Since signals from the QT database are sampled with the frequency of 250 Hz and the EMG noise record with the frequency of 360 Hz, the latter was resampled with the frequency of 250 Hz. The desired signals and the noise component were separately filtered to suppress the baseline wander and the powerline interference. Then, before addition to the desired components, the noise was multiplied by a coefficient to assure the assumed average SNR of the created noisy signals. It should, however, be emphasized that, since the level of the EMG noise is highly variable, the SNR of many individual ECG beats differed much from this assumed, average one. Later these signals were resampled with the frequency of 1000 Hz and then processed

with the use of the tested averaging methods (each signal individually, according to the description provided in Section 2.5). The enhanced signals were analyzed for QRS onset (QRSo ), T wave peak (Tp ) and T wave end (Te ) determination. The classical method developed by Laguna et al. [16], slightly modified in [13], was applied to accomplish the measurements. For reference the measurements were also performed on the ‘desired’ ECG components. To assess the influence of the proposed modifications on the averaging performance, we chose three methods to be compared: NLAAF2c , NLAAF2m and WANACm . WANACm is the method whose performance is to be assessed (weighted averaging of the nonlinearly aligned cycles, based on the modified costs of alignment). NLAAF2c was the most effective among the reference methods based on the classical costs of alignment in the experiment on the simulated evoked potentials (similar was the performance of NLAAF1c , but application of this method is limited to the cases

(l)

Fig. 4. The best results: nrf , of individual evoked potentials enhancement with the use of NLAAF2c , NLAAF2m and WANACm , achieved by the respective methods for the same 10th evoked potential (l = 10) of the processed set of potentials. In each subplot there are, from the top: the noisy EP, the reconstructed EP and its desired component; for better presentation the signals are vertically shifted; a.u. stands for arbitrary units.

166

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168 10

[mV]

6

2 1 0 −1

[mV]

8

2 1 0 −1

2

desired signal

200

E

F

400

600

time [ms]

800

0.05 1 time [s]

2

0 0

Fig. 5. A noisy ECG signal of SNR = 5 dB (A), the results of its enhancement with the use of: NLAAF2c (B), NLAAF2m (C) and WANACm (D), and the desired component of the processed signal (E). Double-sided arrow F shows an extent of a single segment that undergoes time warping and averaging.

QRS

10

20

30

60

70

80 90 segment number

T

p

e

25

25

20

20

No ANAC NLAAF2

15

15

NLAAF2

10

10

5

5

5

50

These visual results are confirmed in Fig. 7 where we have presented the root mean squared errors (RMSE) of the respective characteristic points measurements in the noisy and the enhanced signals. We can see that very high signal to noise ratio is required to determine the respective positions precisely. For Tp and Te the measurements not preceded by any of the investigated methods were accurate only for SNR greater than 20 dB. For 20 dB the measurements contained large errors, which significantly increased the RMSEs. Applying the NLAAF2c method, based on the classical costs of alignment, we slightly improved the results and the curve RMSE(SNR) was approximately shifted to the left by about 5 dB. Applying the modified costs (NLAAF2m ), we obtained even more significant improvement: the curve was shifted by additional 10dB. With the use of weighted averaging (WANACm ), we obtained relatively accurate measurements even for SNR=5dB. When QRS onset determination is considered, the NLAAF2c method appears almost useless. Only when the modified costs of alignment were applied (NLAAF2m , WANACm ) the precision of measurements was significantly raised. The QT series determined without signal enhancement, or with application of either NLAAF2c or WANACm are presented in Fig. 8. While the former approach produced visually acceptable measurements, i.e. the series without large errors only for SNR greater or equal to 20 dB, and application of NLAAF2c helped only a little, the WANACm method helped to obtain acceptable results even for SNR = 5 dB. T

o

6

40

Fig. 6. (a) Results of ECG signal enhancement with the use of the WANACm method, and the QRSo , Tp and Te positions determined on the basis of the desired ECG (the arrows directed downwards) and on the basis of the noisy and the enhanced signal (the arrows directed upwards); (b) bar plot of the weights calculated by WANACm , corresponding to the respective nonlinearly aligned signal segments.

when the cycles number is an integer power of 2). Comparison of NLAAF2m with NLAAF2c will additionally reveal the significance of the alignment costs modification. Visual results of ECG enhancement with the use of these methods are presented in Fig. 5. As we can see, for the first, very noisy QRST segment, application of the classical costs of alignment (NLAAF2c ) produced rather poor results. With the modified costs, the NLAAF2 method appeared much more effective (plot C). It was able to reconstruct the desired ECG with much lower level of the EMG noise. However, when weighted averaging was applied (WANACm ), the noise vanished almost completely (plot D). For the second, less contaminated QRST segment, all the methods performed much better, but again modification of the alignment costs and introduction of weighting discernibly improved the results of ECG enhancement. In Fig. 6 we have illustrated the proposed method (WANACm ) influence on the precision of the QT interval measurement. Uppermost is the desired ECG signal, below is the noisy one, and below, the signal obtained with the use of the WANACm method. We can observe that, like in Fig. 5, also in this case the method caused significant improvement of the ECG signal quality. This favorable operation of the method resulted among others from the proper estimation of the weights corresponding to the respective averaged signal segments. The weights, presented in part b of the figure, limited the influence of the noisy signal segments and thus improved the results of averaging. As a consequence, large errors of Tp and Te determination, which occurred when the measurements were performed on the noisy signal, were precluded after the signal enhancement.

RMSE [ms]

e

weights

b) D

T

enhanced signal (WANACm)

2 1 0 −1 0

0 −2 0

Tp

noisy signal

4

[mV]

[mV]

QRSo

a) A B C

4 3

c

m

2

WANAC

m

1 0 0

10 20 SNR [dB]

30

0

10

20 SNR [dB]

30

0

10

20 SNR [dB]

30

Fig. 7. Illustration of the investigated methods influence on the precision of QRSo , Tp and Te determination in the noisy ECG signals. The RMS errors are presented as functions of SNR. No ANAC denotes the measurements not preceded by averaging of nonlinearly aligned cycles.

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

QT [ms] QT [ms]

340 320 300

QT [ms]

340 320 300

QT [ms]

340 320 300

QT [ms]

SNR=30dB: RMSSD=3.35, SD=3.80 340 320 300

340 320 300 0

SNR=30dB: RMSSD=3.32, SD=3.95 340 320 300

20dB: 8.56, 6.71

20dB: 5.00, 4.58

20dB: 2.59, 3.50 340 320 300

15dB: 8.84, 6.78

15dB: 2.84, 3.61

340 320 300 10dB: 48.01, 38.23

340 320 300 10dB: 26.25, 20.45

10dB: 3.99, 3.94

340 320 300 5dB: 90.12, 66.38

50 beat number

SNR=30dB: RMSSD=2.48, SD=3.35 340 320 300

340 320 300 15dB: 32.92, 23.05

340 320 300 5dB: 63.51, 45.94

100

340 320 300 0

167

50 beat number

5dB: 6.39, 5.17

100

340 320 300 0

50 beat number

100

Fig. 8. Results of QT series determination (with QRSo and Te as the interval onset and offset, respectively) in the noisy signals (on the left side) and in the signals enhanced with NLAAF2c (in the middle) or WANACm (on the right side). The SNR of the analyzed signal and the variability indices of the obtained QT series are given over each subplot. For the series obtained on the basis of the ‘desired’ ECG: RMSSD = 2.77 and SD = 3.61.

Comparing the indices of QT variability that were obtained for the ‘desired’ ECG (see the caption in Fig. 8) with those presented over the respective subplots, we can notice that for SNR of 30 dB the WANACm method slightly decreased this variability. However, these are the errors of the QT interval limits determination that can falsely increase the series variability, and we can observe that the measurements on the original, not enhanced signals can cause such effects even for a very low level of noise (in the leftmost column, for SNR = 30dB, RMSSD = 3.35 instead of 2.77). Such an increase of the series variability for such a low average level of noise shows how difficult it is to achieve reliable assessment of the QT interval variability in real signals where unexpected, momentary bursts of noise can happen. As it is visible in Fig. 5, the EMG noise record from the MIT-BIH database contains such phenomena. In spite of that, the WANACm method helped to obtain the similar values of RMSSD for SNR in the range of 15–30 dB and of SD even down to 10 dB. These results show high advantageous can be the application of the new method proposed to ECG processing. 4. Conclusions We have presented three methods of nonlinearly aligned cycles averaging, known from the literature, and we have described the reasons of their ineffective processing of noisy signals. The fundamental reason is that nonlinear alignment of signals under averaging violates the basic requirements that should be fulfilled to make this averaging effective. We have proposed a modification of a cost matrix that can diminish these unfavorable results of the procedure of dynamic time warping. Application of this modification has significantly improved the results of evoked potentials enhancement with the use of all the three reference methods. To exploit the fact that the noise components of the processed signals are often nonstationary, we have developed a method of weighted averaging of nonlinearly aligned signal cycles. With this method, applied with the modified alignment costs, we have achieved significant improvement of evoked potentials enhancement in nonstationary noise environment. When the most effective among the reference methods and the newly developed one were applied to ECG signal enhancement prior to the measurements of the QT interval, these findings were confirmed. Whereas the reference method based on the classical alignment costs improved the measurements only insignificantly, its application with a modified costs resulted in much better signal enhancement and subsequently in much more precise measurements of the repolarization duration. However, it was the weighted

averaging of nonlinearly aligned cycles, based on the modified alignment costs, that allowed for most noise immune determination of the QT time series and for most reliable assessment of its variability. Acknowledgement This research was supported by statutory funds (BK-2015) of the Institute of Electronics, Silesian University of Technology. The work was performed using the infrastructure supported by POIG.02.03.01-24-099/13 grant: GeCONiI–Upper Silesian Center for Computational Science and Engineering. References [1] E.A. Bartnik, K.J. Blinowska, P.J. Durka, Single evoked potential reconstruction by means of wavelet transform, Biol. Cybern. 67 (1992) 175–181. [2] E. Bataillou, E. Thierry, H. Rix, O. Meste, Weighted averaging using adaptive estimation of the weights, Signal Process. 44 (1) (1995) 51–66. [3] R. Bellman, S. Dreyfus, Applied Dynamic Programming, Princeton Univ. Press, New Jersey, 1962. [4] S. Boudaoud, H. Rix, O. Meste, Integral shape averaging and structural average estimation: a comparative study, IEEE Trans. Signal Proc. 53 (10) (2005) 3644–3650. [5] S. Cerutti, G. Baselli, D. Liberati, G. Avesi, Single sweep analysis of visual evoked potentials through a model of parametric identification, Biol. Cybern. 56 (1987) 111–120. [6] G.D. Dawson, A summation technique for the detection of small evoked potentials, Electroenceph. Clin. Neurophysiol. 6 (1954) 65–84. [7] H. Gibbons, J. Stahl, Response-time corrected averaging of event-related potentials, Clin. Neurophysiol. 118 (2007) 197–208. [8] A.L. Goldberger, L.A. Amaral, L. Glass, J.M. Hausdorff, P.C. Ivanov, R.G. Mark, J.E. Mietus, G.B. Moody, C.K. Peng, H.E. Stanley, PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals, Circulation 101 (23) (2000) E215–E220. [9] L. Gupta, D.L. Molfese, R. Tammana, P.G. Simos, Nonlinear alignment and averaging for estimating the evoked potential, IEEE Trans. Biomed. Eng. 43 (4) (1996) 348–356. [10] E.H. Hon, S.T. Lee, Averaging techniques in fetal electrocardiography, Med. Electron. Biol. Eng. 2 (1) (1964) 71–76. [11] K.S. Kim, J.H. Seo, J.U. Kang, C.G. Song, An enhanced algorithm for knee joint sound classification using feature extraction based on time-frequency analysis, Comp. Methods Programs Biomed. 94 (2009) 198–206. [12] M. Kotas, Application of projection pursuit based robust principal component analysis to ECG enhancement, Biomed. Signal Process. Control 1 (4) (2006) 289–298. [13] M. Kotas, Projective filtering of time-aligned ECG beats for repolarization duration measurement, Comp. Methods Programs Biomed. 85 (2) (2007) 115–123. [14] M. Kotas, Robust projective filtering of time-warped ECG beats, Comp. Methods Programs Biomed. 92 (2) (2008) 161–172. [15] 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, Comput. Cardiol. 24 (1997) 673–676.

168

M. Kotas et al. / Biomedical Signal Processing and Control 21 (2015) 157–168

[16] P. Laguna, N.V. Thakor, P. Caminal, R. Jane, H. Yoon, New algorithm for QT interval analysis in 24h holter ECG: performance and applications, Med. Biol. Eng. Comput. 28 (1990) 67–73. [17] J.M. Leski, Robust weighted averaging, IEEE Trans. Biomed. Eng. 49 (8) (2002) 796–804. [18] X. Long, J. Foussier, P. Fonseca, R. Haakma, R.M. Aarts, Analyzing respiratory effort amplitude for automated sleep stage classification, Biomed. Signal Process. Control 14 (2014) 197–205. [19] M. Malik, V. Batchvarov, Measurement, interpretation and clinical potential of QT dispersion, J. Am. Coll. Cardiol. 36 (2000) 1749–1766. [20] C.D. McGillem, J.I. Aunon, Measurements of signal components in single visually evoked brain response, IEEE Trans. Biomed. Eng. 24 (1977) 232–241. [21] A. Momot, Methods of weighted averaging of ECG signals using Bayesian inference and criterion function minimization, Biomed. Signal Process. Control 4 (2) (2009) 162–169. [22] C. Myers, L.R. Rabiner, A.E. Rosenberg, Performance tradeoffs in dynamic time warping algorithms for isolated word recognition, IEEE Trans. Acoust. Speech Signal Process. 28 (6) (1980) 623–635. [23] F. Petitjean, A. Ketterlin, P. Gancarski, A global averaging method for Dynamic Time Warping, with applications to clustering, Patt. Recog. 44 (2011) 678–693. [24] T.W. Picton, O.G. Lins, M. Scherg, The recording and analysis of event-related potentials, in: F. Boller, J. Grafman (Eds.), in: Handbook of Neuropsychology, Vol.10, Elsevier, New York, 1995, pp. 3–73. [25] R. Quian Quiroga, Obtaining single stimulus evoked potentials with wavelet denoising, Physica D 145 (2000) 278–292. [26] P. Rautaharju, H. Blackburn, The exercise electrocardiogram: experience in analysis of “noisy” cardiograms with a small computer, Am. Heart J. 69 (4) (1965) 515–520.

[27] J.L. Rodriguez-Sotelo, D. Peluffo-Ordon, D. Cuesta-Frau, G. CastellanosDominguez, Unsupervised feature relevance analysis applied to improve ECG heartbeat clustering, Comp. Methods Programs Biomed. 108 (2012) 250–261. [28] O. Rompelman, H.H. Ros, Coherent averaging technique: a tutorial review, J. Biomed. Eng. 8 (1986) 24–35. [29] H. Sakoe, S. Chiba, Dynamic programming algorithm optimization for spoken word recognition, IEEE Trans. Acoust. Speech Signal Process. 26 (1) (1978) 43–49. [30] M. Schmidt, M. Baumert, A. Porta, H. Malberg, S. Zaunseder, Two-dimensional warping for one-dimensional signals-conceptual framework and application to ECG processing, IEEE Trans. Signal Proc. 62 (21) (2014) 5577–5588. [31] W.W. Surwillo, Cortical evoked potentials in monozygotic twins and unrelated subjects: comparisons of exogenous and endogenous components, Behav. Genet. 10 (1980) 201–209. [32] 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. Phys. 21 (1999) 15–25. [33] J. Wang, P. Liu, M.F.H. She, S. Nahavandi, A. Kouzani, Bag-of-words representation for biomedical time series classification, Biomed. Signal Process. Control 8 (2013) 634–644. [34] K. Wang, H. Begleiter, B. Porjesz, Warp-averaging event-related potentials, Clinical Neurophysiology 112 (2001) 1917–1924. [35] K. Wang, T. Gasser, Synchronising sample curves nonparametrically, Ann. Stat. 27 (1999) 439–460. [36] T. Wang, L. Lin, A. Zhang, X. Peng, C.A. Zhan, EMD-based EEG signal enhancement for auditory evoked potential recovery under high stimulus-rate paradigm, Biomed. Signal Process. Control 8 (2013) 858–868. [37] C.R. Woody, Characterisation of an adaptive filter for the analysis of variable latency neuroelectric signals, Med. Biol. Eng. 5 (1967) 539–553.