Broadband noise suppression and feature identification of ECG waveforms using mathematical morphology and embedding theorem

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 1 1 2 ( 2 0 1 3 ) 466–480

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

Broadband noise suppression and feature identification of ECG waveforms using mathematical morphology and embedding theorem T.Y. Ji a , Q.H. Wu a,b,∗ a b

School of Electric Power Engineering, South China University of Technology, Guangzhou 510641, China Department of Electrical Engineering and Electronics, The University of Liverpool, Liverpool L69 3GJ, UK

a r t i c l e

i n f o

a b s t r a c t

Article history:

The paper presents an adaptive morphological filter developed using multiscale mathe-

Received 3 March 2013

matical morphology (MM) to reject broadband noise from ECG signals without affecting the

Received in revised form

feature waveforms. As a pre-processing procedure, the adaptive morphological filter cleans

8 August 2013

an ECG signal to prepare it for further analysis. The noiseless ECG signal is embedded within

Accepted 14 August 2013

a two-dimensional phase space to form a binary image and the identification of the feature

Keywords:

of the feature waveforms is implemented by an adaptive clustering technique according to

waveforms is carried out based on the information presented by the image. The classification ECG signal processing

the geometric information represented by the image in the phase space. Simulation studies

Feature identification

on ECG records from the MIT-BIH and BIDMC databases have demonstrated the effectiveness

Multiscale mathematical

and accuracy of the proposed methods.

morphology

© 2013 Elsevier Ireland Ltd. All rights reserved.

Adaptive filter Embedding theorem

1.

Introduction

The ECG is a signal acquired from the body surface which characterizes the electrical activity of the human heart. A normal heart is reflected on the ECG by three distinctive feature waveforms: the contraction of the atria, which results in the so-called P wave, the QRS complex, which is produced by the contraction of the ventricles, and a T wave, which results from the subsequent return of the ventricular mass to a rest state [1].

The ECG is a quasi-periodic signal where each elementary beat is repeated over time with certain variability on the distance between contiguous beats. Fig. 1 shows a heartbeat and its respective waveform labels. The QRS complexes represent the ventricular activity of the heart; the R waves have the highest amplitude, and an RR interval can be used to define a cycle of the ECG signal. An ECG signal is often disturbed by broadband noise, mainly composed of high-frequency interferences due to electromagnetism and the grounding. Therefore, noise removal is usually the first step of any robust and reliable ECG analysis

∗ Corresponding author at: Department of Electrical Engineering and Electronics, The University of Liverpool, Liverpool L69 3GJ, UK. Tel.: +44 1517944535. E-mail addresses: [email protected], [email protected] (Q.H. Wu). 0169-2607/$ – see front matter © 2013 Elsevier Ireland Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cmpb.2013.08.006

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 1 1 2 ( 2 0 1 3 ) 466–480

Fig. 1 – Feature waveforms of an ECG signal.

[2]. ECG noise is often characterized by the same magnitude and frequency of the feature waveforms; it is thus important to minimize the distortion in the feature waveforms so as to keep those features that would be of most interest in terms of clinical diagnosis, whilst at the same time removing the noise. In our previous work a multi-resolution morphological filter was proposed for use in removing impulsive noise [3]. In Ref. [4] a multistage multiscale mathematical morphology filtering algorithm was used to suppress impulsive noise. A method recently proposed was to preserve QRS complexes whilst removing noise; however, the computation burden was considerably heavier, as it used time-consuming empirical mode decomposition and discrete wavelet transform [5]. Chen et al. proposed a progressive umbra-filling (PUF) procedure to adaptively adjust the filtering scales, which was useful for broadband noise removal [6]. The accurate extraction of feature waveforms is of great concern. Some algorithms detected QRS complexes only, based on their inherent characteristics such as high frequency content and large amplitude [7,8]. Conversely, other algorithms focused to obtain all the characteristic points of the ECG, i.e. the onset, peak and end of the P waves, QRS complexes and T waves [9–11]. Our previous works [3] and [12] proposed morphological approaches for baseline normalization and feature waveform detection. However, they required too much computation time. The detection scheme proposed in this paper involves phase space embedding and MM. It is not the first attempt that has been undertaken to engage the embedding theorem in ECG waveform detection. Paper [13] was a preliminary study that discussed the feasibility of applying phase space embedding to ECG signal processing, and indicated that the P waves, QRS complexes and T waves formed distinctive loops in the phase space. In another work, ECG signal was processed in the phase space for P wave detection [14]. In this study, the ECG signal was embedded in an incremental phase space, and the morphology distance of each ECG point to a reference point was subsequently computed; finally, the point with the maximum distance was detected as the P wave

467

onset. The algorithm was reported to be accurate. However, the reference point needed to be selected manually for each and every heart-beat cycle, and no criteria was provided in respect of how to choose a proper reference point. The classification of the feature waveforms has also been well investigated. A large amount of literature has been devoted to this field. In Ref. [15], several ANN architectures were designed to classify the ECG signals integrating the most common features: arrhythmia, myocardial ischemia, chronic alterations. In Ref. [16], a hidden Markov model-based approach was proposed for automatic beat segmentation and classification. A template cross-correlation technique was applied for ECG waveform identification in [17]. An unsupervised method based on relevance analysis was described to improve ECG heartbeat clustering in [18]. Some literature focused on the extraction of P waves and their morphology classification, and several clinical studies connected certain P-wave properties, such as its width and morphology, with anomalies in the electrical atrial conduction and atrial pathology [9]. Whereas [19] used a cluster analysis method to analyze ECG waveforms to diagnose cardiac arrhythmias, another arrhythmia classification algorithm was based on the density distribution of ECG signals in the phase space [20]. Instead of analyzing the feature waveforms separately, the algorithm took an ECG signal of 5 s duration as a whole, and calculated three density distribution criteria for arrhythmia classification.

2.

Noise removal

2.1.

Drawbacks of conventional filtering techniques

Noise removal is the first step of almost all the signal processing tasks. As the spectrum of pure ECG signals overlaps with that of the noise, it is difficult to separate them in the frequency domain. On the other hand, MM is a technique that processes signals in the time domain. Appendix A.1 introduces the basic concepts of MM. However, it is difficult to select a proper structuring element (SE) of a basic morphological filter, such as opening and closing, to remove noise while preserving signal features. If the SE is too short, the noise cannot be effectively removed; if the SE is too large, the feature waveforms may be filtered as well. The example given in Fig. 2 illustrates the phenomenon: for clarity, fragments of the PR segment and of the QRS complex have been enlarged. The ECG signal is filtered by two openings with a flat SE of size 3 and size 11, respectively. When the SE is short, the resulting signal still vibrates at a high frequency, which is the evidence that there exists a large amount of noise. On the other hand, when the SE is too long, although the UP segment is smoothed, the R wave is impaired at the same time. Therefore, it is necessary to design an adaptive filter that can distinguish feature waves and use SEs of different sizes accordingly.

2.2.

The PUF procedure

The PUF procedure [6] engages the basic morphological filters of opening and closing, and reconstructs a signal using a

468

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 1 1 2 ( 2 0 1 3 ) 466–480

985

1500

1400

980

1300

Magnitude

Magnitude

975

970

1200

1100

965 1000 960

955 130

900

140

150 160 170 (a) Sample

180

800 330

340

350 360 370 (b) Sample

380

Fig. 2 – Filtering results by openings with SEs of different sizes. Dotted line: input signal. Dashed line: the size of SE is 3. Solid line: the size of SE is 11.

series of SEs which are of the same shape but different sizes. The SEs have the form of {rk|r ∈ R, 0 ≤ r ≤ r}, where k is a convex function, r is the size of the SE, and r is the largest size allowed, depending on the application. Applying the SEs in the decreasing order of r, opening and closing result in a coarse-tofine representation of the signal. Based on this property, the PUF procedure makes sure that at each position of a signal, local waveform feature is considered to adaptively select the filtering scale. Besides, the signals obtained from opening are below the original one, while signals resulted from closing are above the original signal. Therefore, the average of opening and closing makes a proper filtering result.

The adaptive morphological filter

1150

In this paper, the adaptive morphological filter is proposed for broadband noise removal, which is a simplification and improvement of the PUF procedure. In order to fit the waveforms of ECG signals, the SE takes a circular function rather than a flat one. The SE is defined as: g(s) =

• Step 1: Denote the input signal by f0 . Define the radius of g0 , r0 , in the first iteration. Set i : =1. • Step 2: Calculate the opening of fi−1 by (16): fi = fi−1 ◦ gi−1 . • Step 3: Calculate the opening residue: Ri = fi−1 − fi . If |Ri (t)| < TR , where T is a pre-set threshold with a small value, the residue of Ri (t) is considered as noise. Otherwise, a decision-making process is performed to determine if it is noise or feature.



r2 − s2 ,

−r ≤ s ≤ r

(1)

where r is the radius of the circular function and the length of g is 2r + 1. The advantage of using a circular SE over a flat one is demonstrated in Fig. 3. Here, the length of the SEs is 21 and the origins of the SE are their geometry centers. As it can be seen, the circular SE generates a curve that fits the waveform of the input signal better than the flat SE does. The adaptive opening performs in a recursive way. In the first iteration, an SE long enough to cover the width of P waves is used; and in each iteration afterwards, the length of the

1100

Magnitude

2.3.

SE decreases at a fixed rate. The adaptive opening works by executing the following steps [21].

1050

1000

950

40

60

80

100

120 140 Sample

160

180

200

220

Fig. 3 – Filtering results by openings with SEs of different shapes. Dotted line: input signal. Dashed line: using a flat SE. Solid line: using a circular SE.

469

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 1 1 2 ( 2 0 1 3 ) 466–480

1600 Magnitude

1400 1200 1000 800

0

500

1000

1500

(a) Sample

Magnitude

990

1160 1140

980

1120 970

1100

960

1080 1060

950 120

140 160 (b) Sample

180

440

450 460 (c) Sample

470

Fig. 4 – (a) Filtering result by the adaptive filter. (b) The enlargement of the first UP segment. (c) The enlargement of the top of the second T wave. Dotted line: input signal. Dashed line: result of the adaptive opening. Dash-dot line: result of the adaptive closing. Solid line: average of adaptive opening and closing. 0 −0.1 lg C(r)

−0.2 −0.3 −0.4 −0.5 −0.6 0.8

1

1.2

1.4 lg r

1.6

1.8

2

2.2

Fig. 5 – lg C(r) against lg r for a de-noised test ECG signal when initial conditions are dE = 2 and  = 1.

- Step 3.1: For each segment of Ri (t) ≥ T, t1 ≤ t ≤ t2 , calculate its width, w = t2 − t1 , and its maximum magnitude, v = max{Ri (t)}. t

- Step 3.3: Otherwise, the segment is considered to be noise. • Step 4: The segments that are considered to be features are added back to fi according to their position:

- Step 3.2: If either of the two indicators exceeds a predefined threshold, that is w > Tw or v > Tv , the segment is considered to be a feature that should not be removed. Denote these Ri (t) by a set R. In order to make the width indicator adaptive, the threshold Tw changes with r.

 f˜i (t) =

fi (t) + Ri (t),

ifRi (t) ∈ R,

fi (t),

otherwise.

(2)

0.65

1

lg CdE(τ)

0.6 0.55 0.5 0.45 0.4 0.35 0

5

10

15 τ

20

25

Fig. 6 – lgC1dE () against  of a de-noised test ECG signal.

30

470

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 1 1 2 ( 2 0 1 3 ) 466–480

• Step 5: Set fi+1 := f˜i . Set ri+1 = ri − rstep and generate gi+1 . Set i : = i + 1. Go to Step 2. • Step 6: The process ends when ri = 1 or when |Ri (t)| < TR for all t. The final output is denoted by fo . The adaptive closing can be performed in the same manner. However, in Step 3, Ri = fi − fi−1 should be used to calculate the closing residue, due to the fact that fi = fi−1 • gi−1 > fi−1 . The final output of adaptive closing is denoted by fc . As the filtering result by opening is always smaller than the input signal and that by closing is always larger, a more practical way is to take their average as the final filtering result. Therefore, the filtering result of the adaptive filter is: ı(f ) =

fo + fc . 2

(3)

An example of the performance of the adaptive filter is given in Fig. 4. The parameters of the adaptive morphological filter are set to be r0 = 25, rstep = 4, TR = 1, Twi = ri /5 and Tv = 20. In order to make it clearer, the figure is enlarged at the first UP segment and the second T wave. As it can be seen, the noise is smoothed out, the feature waveforms are retained, and the detail information, such as the tops of P waves, R waves and T waves and the valleys of Q waves and S waves, are well preserved.

Fig. 7 – Embedding of a de-noised test ECG signal when dE = 2 and  = 15.

3.

Identification of feature waveforms

baseline wander, the magnitude difference between every  samples are calculated:

3.1.

ECG signal embedding

dl (t) = f (t) − f (t + ).

To embed the ECG signals to the phase space, the embedding dimension and the delay constant are determined according to the principles given in Appendix A.2. Here, dE = 2 and  = 1 are selected as the initial condition for the calculation of the real correlation dimension, dc . For a noiseless test ECG signal of length 1000, the relationship between lg r and lg C(r) is plotted in Fig. 5. In this case, dc = 0.4390; hence, dE = 2. For the same test ECG signals of length 3000, the embedding dimension is dE = 3. Thus, the choice of dE depends on the length of the signal being processed. Another factor that influences the embedding result is the delay constant. For the noiseless test ECG signal of length 1000, the relationship between lgCd1E () and  is given in Fig. 6 and the minimum of lgCd1E occurs at  = 15. Fig. 7 shows the embedding of the noiseless test ECG signal when dE = 2 and  = 15. The white pixels in the figure are enlarged to give a clearer view. As the magnitude of the test signal is digitized, the magnitude of the noiseless test signal is also regulated to integers. Thus, the phase space can be presented by an image, in which black pixels form the embedded signal, and white pixels form the background of the image. The dense cluster corresponds to samples belonging to the segments, and the three orbits correspond to the P waves, T waves and QRS complexes, respectively, according to their size.

3.2.

Detection of QRS complexes

The QRS complexes have the distinctive shape of high gradient value, and this characteristic is used to detect the complexes. In order to perform the detection with the tolerance of the

(4)

The difference signal of dl is referred to as the left-hand difference. Local minima of dl with an absolute magnitude value over a preset threshold are recorded as the Q points. The threshold is set at half of the difference between the local maximum and minimum within a window that contains a number of cycles. Similarly, the S points can be located from the right-hand difference of f: dr (t) = f (t) − f (t − ).

(5)

Obtaining a pair of R and S points, the Q point is located between the R and S points, where the maximum gradient value with respect to the R or S point is achieved. Fig. 8 demonstrates the detection of the QRS complexes of a test signal during normal beats. Fig. 9 shows the detection result during a three-beat ventricular tachycardia. The algorithm has been tested on a group of arrhythmia signals and the detection results are listed in Table 1. The ECG signals contain records in the first five minutes and in the whole thirty minutes period respectively. The number of the QRS complexes detected by the algorithm, i.e. the detected beats, is compared with the actual number provided by the database. The errors mainly occur during fusion, which changes the QRS. Take record 208 for example. In the first five minutes, the ECG signal contains 72 fusions, which reduces the accuracy of the detection by 2.7%. On the identification of the R points, an RR interval is considered as a ‘cycle’ of the signal. In order to test its performance on a wider spectrum of ECG leads, the method is also tested on a group of cardiac signals.

471

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 1 1 2 ( 2 0 1 3 ) 466–480

Magnitude

Magnitude

Magnitude

1500

1000 0

500

1000

1500 Sample

2000

2500

3000

0

500

1000

1500

2000

2500

3000

0

500

1000

1500 Sample

2000

2500

3000

500 0 −500

500 0 −500

Magnitude

Fig. 8 – (a) Location of the Q points (in dot), R points (in circle) and S points (in diamond). (b) Left-hand difference. (c) Right-hand difference.

1200 1000 800 6.3

6.35

6.4

6.45

6.5

Sample

4

x 10

Fig. 9 – Location of the Q points (in dot), R points (in circle) and S points (in diamond) during a 3-beat ventricular tachycardia.

Table 1 – Results of QRS detection of arrhythmia records. ECG records

5 min records Total beats

105 109 200 208 212 217 221 228

417 433 433 518 463 363 407 350

Detected beats

30 min records Accuracy rate

Total beats

99.76% 99.77% 97.69% 97.30% 100% 97.52% 99.51% 97.14%

2572 2532 2601 2955 2748 2208 2427 2053

418 432 423 504 463 354 409 340

Detected beats 2530 2516 2536 2908 2740 2144 2461 1995

Accuracy rate 98.37% 99.37% 97.50% 98.41% 99.71% 97.10% 98.60% 97.17%

Table 2 – Results of QRS detection of congestive heart failure records. ECG records

Total beats

Records from lead 1 Detected beats

chf01 chf03 chf05 chf07 chf09 chf11 chf13

331 326 541 418 505 577 491

328 321 529 416 502 563 487

Accuracy rate 99.09% 98.47% 97.78% 99.52% 99.41% 97.57% 99.19%

Records from lead 2 Detected beats 329 327 537 417 497 575 490

Accuracy rate 98.79% 99.69% 99.26% 99.76% 98.42% 99.65% 99.80%

472

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 1 1 2 ( 2 0 1 3 ) 466–480

The group contains 7 ECG signals, and each signal contains records in the first five minutes from two leads, respectively. The detection results are given in Table 2, which shows that the method achieves an average detection accuracy of 98.72% and 99.34% on records from the two leads, respectively. Therefore, it is safe to conclude that the QRS complex detection method performs equally well no matter which lead does the record come from.

3.3.

Detection of P waves and T waves

The detection of P waves and T waves is based on the embedded signal in the phase space. As shown in Fig. 7, samples belonging to the segments form a dense cluster, which can be extracted from the image using morphological reconstruction by dilation. Morphological reconstruction requires two images: a mask image and a marker image. The algorithm of morphological reconstruction by dilation of a mask image f1 from a marker image f2 (f2 ≤ f1 ) is described as follows [22].

respectively. The P and T waves can be classified using the geometric information in the phase space, such as centroid, barycenter and area. To verify its ability of extracting feature waveforms under abnormal morphology conditions, the proposed method is tested on congestive heart failure signals. The statistical results are given in Tables 3 and 4, respectively. Figs. 12 and 13 give two examples of the P and T wave detection results on ECG waveforms with abnormal morphology. As shown in the figures, under moderate conditions the proposed method is able to detect the P and T waves accurately. However, in the case of severe abnormal morphology, as shown in Fig. 12(c) and (d), the detection result is less satisfactory. This is mainly because (i) the location of Q and S points is less accurate, as the gradient value at these points is not the local maximum; (ii) the waveform does not contain clear segment morphology. Thus, it is difficult to identify the dense cluster in the phase space which is supposed to represent the segment. As a result, the objective pixels which are detected belonging to P or T waves may differ from the actual ones.

(0)

• Step 1: Set k : =1. Set f2 = f2 . (k)

(k−1)

⊕ g, where g is the predefined • Step 2: Dilate f2 : f2 = f2 SE. (k) • Step 3: Calculate the point-by-point minimum of f2 and f1 : (k)

ı(k) (t) = min{f2 (t), f1 (t)}. • Step 4: Set k : = k + 1. Go to Step 2. If ı(k) = ı(k+1) , terminate. (k)

For a binary image, the intersection of f2 and f1 is used in step 3. In this paper, the mask image is the embedded signal, X, and the marker image is the erosion of X by an SE of

⎡

⎤

1 1 g = ⎣1 1 0 1

0 1⎦ 1

as shown in Fig. 10(a). The SE can also be set as [0, 1, 0 ; 1, 1, 1 ; 0, 1, 0], [0, 1 ; 1, 1], [1, 1 ; 1, 0] or [1, 0 ; 0, 1] to match the shape of the object. The result of morphological reconstruction by dilation is given in Fig. 10(b). The dense cluster is denoted by Xb . Apparently, Xb ⊆ X. To separate a P or T wave from the segment, the distance from each objective pixel to the set of Xb is calculated. The distance from a pixel x to a set A is defined as the shortest distance from x to a pixel in A: d(x, A) = min{ x − ap }, p

ap ∈ A.

(6)

From a detected Q point leftwards, a string of pixels whose distance to Xb is larger than zero are recorded as the embedding of the P wave, as shown in Fig. 10(c). From the detected S point rightwards, a string of pixels whose distance to Xb is larger than zero are recorded as the embedding of the T wave, as shown in Fig. 10(d). To view it more clearly, the pixels in Fig. 10(c) and (d) are enlarged. Therefore, samples in the time domain that comprise the embedded P wave and T wave are detected to form the feature waveforms of P wave and T wave. The result is given in Fig. 11. In clinic, doctors are interested in the width and peak value of feature waveforms. The measurement of the width and peak value of P waves and T waves of eight arrhythmia signals are given in Tables 3 and 4,

Classification of feature waveforms

4.

In this section, the geometric information of the feature waveforms represented by the embedding in the phase space is used as classification criterion. Consider for example the classification of the P waves: the geometric information includes the length of the P wave in the time domain and the perimeter and area of the embedded P wave in the phase space. Hence, each P wave is coded as a 3-element vector vi ∈ V with V the set of the vectors and i the index of the P waves. The elements of the vector are normalized to [0, 1] for meaningful clustering. The clustering method is based on the distance between vi and a cluster center, denoted by cj ∈ C with C the set of the cluster centers, 1 ≤ j ≤ J the index of the cluster center, and J the number of clusters. For an arbitrary vi , if vi − ck = min{ vi − cj }, 1≤j≤J

1≤k≤J

(7)

where · denotes the Euclidean norm, vi belongs to cluster ck . The cluster centers are adaptively generated and their number is not constrained to a preset number, as the case in [9] where the P waves are always classified into two clusters. The classification procedure is described as follows. • Step 1: Set J : =1. Initialize the first cluster center c1 := v1 . Set i : =1. • Step 2: Calculate the distance dj = vi − cj for all 1 ≤ j ≤ J. If min{dj } exceeds a predefined threshold, , a new cluster j

center is assigned at cJ+1 = vi . Set J : = J + 1. Otherwise, vi is classified to cluster ck if dk = min{dj }. j

• Step 3: Update the cluster centers by averaging the vectors belonging to each cj . • Step 4: Set i : = i + 1. Go to Step 2. Terminate when all the vectors are classified. Fig. 14 shows the classification result of the P waves of the first five minutes of record 106. To make it more

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 1 1 2 ( 2 0 1 3 ) 466–480

473

(a)

(b)

Magnitude

1400 1200 1000 800

800

900 1000 Sample

1100

800

900 1000 Sample

1100

800

900 1000 Sample

1100

(c)

Magnitude

1400 1200 1000 800

(d)

Magnitude

1400 1200 1000 800

Fig. 10 – (a) The marker image. (b) Morphological reconstruction by dilation and the extracted baseline in the time domain. (c) Extracted P wave in the phase domain and time domain, respectively. (d) Extracted T wave in the phase domain and time domain, respectively.

comprehensible, the P waves are normalized taking on zero value at the onset and end samples, as the strategy used in [9], and having the same length. The threshold is set at  = 1 and under this condition, the P waves are classified into two

clusters. As a comparison, when  = 0.8, the P waves are classified into four clusters, as shown in Fig. 15. The same clustering method can also be applied to classify the QRS complexes and T waves.

474

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 1 1 2 ( 2 0 1 3 ) 466–480

1600

Magnitude

1400 1200 1000 800 200

400

600

800 1000 Sample

1200

1400

1600

Fig. 11 – Identification of P waves (in gray) and T waves (in black).

Table 3 – Measurement of the detected P waves. ECG record

Width of P waves (ms) Average

Peak value of P waves (mV)

Standard deviation

Average

Standard deviation

105 109 200 208 212 217 221 228

220.03 330.06 289.72 278.81 224.94 543.02 246.34 218.77

50.86 139.08 124.97 130.54 45.32 168.38 99.21 144.10

0.23 0.23 0.16 0.19 0.28 0.87 0.08 0.20

0.04 0.15 0.16 0.09 0.05 0.47 0.09 0.14

chf01 chf03 chf05 chf07 chf09 chf11 chf13

188.97 220.80 143.68 205.74 147.74 138.19 187.33

72.97 126.97 67.69 92.93 25.00 82.33 106.28

0.23 0.27 0.18 0.25 0.53 0.17 0.23

0.09 0.15 0.08 0.11 0.12 0.10 0.13

Magnitude

500

0 0 −200

−500 1000

Magnitude

200

1100 1200 (a) Sample

1300

−400 1000

500

500

0

0

−500 1900

1950 2000 2050 (c) Sample

2100

−500 1900

1100 1200 (b) Sample

1300

1950 2000 2050 (d) Sample

2100

Fig. 12 – Identification of P waves with abnormal morphology. (a) and (c) Signal segments of chf09 from lead 1. (b) and (d) Signal segments of chf09 from lead 2.

475

Magnitude

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 1 1 2 ( 2 0 1 3 ) 466–480

300

200

200

0

100

−200

0

−400

−100 5300

5400 5500 (a) Sample

5600

Magnitude

100

−600 5300

5400 5500 (b) Sample

5600

2500 2600 (d) Sample

2700

400

0

200

−100 0

−200 −300 2400

2500 2600 (c) Sample

2700

−200 2400

Fig. 13 – Identification of T waves with abnormal morphology. (a) Signal segment of chf01 from lead 1. (b) Signal segment of chf01 from lead 2. (c) Signal segment of chf03 from lead 1. (d) Signal segment of chf03 from lead 2.

5.

Comparison studies

5.1.

Noise removal

It is hard work to present an indicator to evaluate how much the noise is suppressed and how well the feature waveforms are preserved, as there is no such thing as a standard noisefree ECG signal. A comparison study is given in Fig. 16, where the proposed adaptive morphological filter and the MM-based filter proposed in [23] are applied to remove an artificial noise of 12 dB added to ECG record 118. As can be seen from the figure, the adaptive morphological filter is more able to reject the added noise and reveal the original signal. Yet we also

need to address that the computation time of the adaptive morphological filter is 0.26 s while that of the MM-based filter proposed in [23] is only 0.044 s. Both computation time is the average of 10 test runs on a Dell PC with 3.1 GHz CPU and 4.0 GB RAM using MATLAB. For the purpose of quantitative analysis, the following experiment is designed: add artificial noise to original ECG signals, remove the noise, and compare the de-noised signal with the original one using the indicator of spectrum deviation (X, Y): 1 (X, Y) = S S

i=1



xi − (xi + yi )/2 (xi + yi /2)

2

+

yi − (xi + yi /2) (xi + yi /2)

2 (8)

Table 4 – Measurement of the detected T waves. ECG record

Width of T waves (ms) Average

Standard deviation

Peak value of T waves (mV) Average

Standard deviation

105 109 200 208 212 217 221 228

310.89 335.52 398.18 318.53 318.82 586.73 306.30 447.29

64.47 137.90 116.80 106.00 45.52 107.73 81.63 139.63

0.19 0.23 0.34 0.27 0.46 1.03 0.16 0.34

0.06 0.14 0.27 0.18 0.14 0.32 0.07 0.24

chf01 chf03 chf05 chf07 chf09 chf11 chf13

449.09 341.40 260.65 384.00 215.40 236.15 251.58

82.92 166.87 66.88 107.04 89.94 69.34 77.94

0.85 0.15 0.25 0.28 0.26 0.98 0.11

0.26 0.11 0.07 0.09 0.17 0.34 0.06

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 1 1 2 ( 2 0 1 3 ) 466–480

250

250

200

200 Normalised magnitude

Normalised magnitude

476

150 100 50 0

100 50 0 −50

−50 −100 0

150

50

100 150 200 Normalised sample

250

−100 0

50

100 150 200 Normalised sample

250

Fig. 14 – Classification of P waves of an ECG signal ( = 1). Normalized P waves in black are classified to cluster 1 and those in red are classified to cluster 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 15 – Classification of P waves of an ECG signal ( = 0.8). Normalized P waves are classified into four clusters, plotted in black, red, green and blue, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

where X = {x1 , x2 , . . ., xS } and Y = {y1 , y2 , . . ., yS } denote two data vectors and S is their length. If the shapes of X and Y are similar to each other, (X, Y) approaches 0. The algorithm of adaptive morphological filter proposed in this paper is compared with the 3M filter proposed in [4], and simulation studies are carried out on ECG records as mentioned in Tables 1 and 3. A white Gaussian noise with a signal-to-noise ratio (SNR) of

30 dB is added to the original ECG signals. The average spectrum deviation is 0.1488 for the adaptive morphological filter and 0.1640 for the 3M filter. It shows that the filter proposed in this paper is able to remove the noise meanwhile keep the shape of the feature waveforms. The fact that waveform detection is accurately performed also verifies the filter’s ability of shape keeping.

Magnitude

1200 1000 800 600

original ECG de−noised ECG

400 0

500

1000 (a) Sample

1500

2000

Magnitude

1200 1000 800 600

original ECG de−noised ECG

400 0

500

1000 (b) Sample

1500

2000

Fig. 16 – A comparison study between (a) the proposed adaptive morphological filter and (b) the MM-based filter proposed in [23].

477

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 1 1 2 ( 2 0 1 3 ) 466–480

Table 5 – Comparison studies with existing algorithms. Algorithm

FP

FN

Se (%)

+P%

DER

Topological mapping Wavelet transform ANN adaptive filter 3M method Proposed scheme

41 15 10 19 12

4 12 4 7 7

– – – 99.26 99.53

– – – 99.73 99.73

1.75 1.09 0.50 1.01 0.74

5.2.

Waveform detection

A comparison study is carried out to evaluate the performance of the waveform detection schemes proposed in this paper and in other literatures. Performance indexes are derived from [4], which include false negative (FN): failing to detect a true beat, false positive (FP): a false beat detection, sensitivity (Se): Se(%) =

TP TP + FN

(9)

positive prediction (+P): +P(%) =

TP TP + FP

(10)

and detection error (DER): DER(%) =

FT + FN Total number of QRS complexes

(11)

where TP denotes true positive – the total number of the QRS complexes correctly detected by an algorithm. Comparison results are listed in Table 5, where the results of existing algorithms are obtained from [4]. The algorithms are carried out on ECG record 105 and it contains all samples from the thirtyminute recording. As can be seen from the table, the scheme proposed in this paper exceeds most existing waveform detection schemes reported in recent literatures.

6.

Discussion

In the simulation studies, test ECG signals come from MITBIH Arrhythmia Database and BIDMC Congestive Heart Failure Database, respectively. ECG signals from MIT-BIH Arrhythmia Database were digitized at 360 samples per second with 11-bit resolution over ±5 mV range. Sample values thus range from 0 to 2047 inclusive, with a value of 1024 corresponding to 0 mV. ECG recordings from BIDMC Congestive Heart Failure Database were sampled at 250 samples per second with 12-bit resolution over ±10 mV range. Each recording contains two signals obtained from two leads, respectively. The work presented in this paper is based on our previous works [3,24], where impulse noise and baseline drift are removed beforehand. To remove the broadband noise, this paper modifies the procedure of PUF so that the adaptive filter are designed particularly for tackling ECG signals and take less computation time than their original version reported in [6]. The adaptive filter ensures that the feature waveforms are not distorted while the noise is being removed. The parameters of the adaptive morphological filter are set according to the following principles. r0 should be long enough so that g0 can cover the P waves. rstep should be neither too long to miss

any useful information, nor too short to involve redundant computation. TR is set to classify small disturbance as noise. Twi and Tv are used to retain the components that may carry features. Based on these principles and comprehensive trials, the parameters are set to be r0 = 25, rstep = 4, TR = 1, Twi = ri /5 and Tv = 20 in all the simulation studies. Apparently, the setting should vary according to the sampling frequency. The parameters need to be increased/decreased proportionally as the sampling frequency increases/decreases. In order to evaluate how much noise is suppressed and how well the feature waveforms are preserved, the SNR is typically considered. However, in this research, as the original ECG signals contain broadband noise and no standard noisefree ECG signals are available, it is difficult to calculate the SNR. Therefore, a substitute study is performed. The adaptive morphological filter is tested on ECG records from the MITBIH Noise Stress Test Database, which provides records with additive noise of 24, 18, 12 and 6 dB. The SNR of the filtered ECG signals is calculated against the original ECG signals, and the result, which is the average of 30 test runs, is listed in Table 6. As can be seen from the table, if there is quite an amount of additive noise, the adaptive morphological filter is able to reject the noise and enhance the SNR by roughly 60%; whereas if only a small amount of noise is added, the improvement of SNR is not so significant, which is mainly because the broadband noise contained in the original signal is also removed. For the detection of feature waveforms, the paper makes use of the inherent characteristics of the QRS complex to locate the Q, R and S points. Afterwards, the signals are embedded to the phase space according to the embedding theorem and are transformed into binary images. For the original test signal, the embedding parameters are dE = 2 and  = 8. The smaller value of  shows that the signal contains a greater amount of redundant information, which is introduced by the noise, and thus more unnecessary calculation is included. Thanks to the noise removal procedure, the redundant calculation is avoided, and the embedding parameters are dE = 2 and  = 15. Such a setting works for ECG signals of length 500–2000. The windowed data can be read into a buffer and the proposed scheme can be implemented to process the data on line. However, when the length of the ECG signal changes, the values of dE and  should be reset in the same way as described in Appendix A.2.

Table 6 – SNR of the filtered ECG signals. Noisy signals (dB) Filtered signals (dB) Improvement (%)

6 9.93 65.50

12 19.27 60.58

18 25.16 39.78

24 28.57 19.04

478

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 1 1 2 ( 2 0 1 3 ) 466–480

When dE = 2, the embedding signal forms a binary image. This is due to the fact that the ECG signal has integer values. As shown in (20), pixels x1 = (x1 , x1+ ), x2 = (x2 , x2+ ), . . .,xM = (xM , xM+ ) form the object, while the rest pixels form the background. The extraction of the baseline, P waves and T waves is performed based on the embedding in the phase space and can be adjusted to involve one cycle only. Hence, it is robust to the slow changes in the magnitude of the baseline, which is caused by the baseline wander. Consequently, clinically important properties such as the onset, peak and end of feature waveforms can be obtained. Moreover, the geometric information represented by the embedded signal in the phase space is extracted and used for classification. A distance-based clustering technique is adopted and the feature waveforms are classified into a number of clusters according to their morphology. The emergence of a new cluster may correspond to a new clinical phenomenon or even a defect. The medical interpretation of waveform morphology is to be further investigated in future research.

7.

Conclusion

This paper introduces a pre-processing scheme for ECG signals to be pre-treated for medical diagnosis. An adaptive morphological filter has been proposed to remove broadband noise and smooth out the original signals without distorting the feature waveforms or eliminating any useful information. The noise removal procedure also helps impair redundancy for the embedding procedure, which is the second step of the pre-processing scheme. Embedded to the two-dimensional space, the de-noised ECG signal becomes a binary image and the feature waveforms form distinctive shapes. Based on such images, the identification of the feature waveforms is performed automatically, this avoiding the tedious process of manual extraction. The classification of the identified feature waveforms is based on the geometric information represented by the embedded signal in the phase space, and it is carried out with an adaptive approach. By assigning a predefined threshold that indicates data distance, feature waveforms are accordingly classified into a proper number of clusters. Simulation studies have been carried out on ECG records from open databases. The results have shown that the adaptive morphological filter is able to remove broadband noise while keeping the feature waveforms undistorted. Compared to other filters proposed in literature, the adaptive morphological filter has superior performance in both aspects. The accuracy of feature waveform detection has also been demonstrated, and comparison studies have shown that the proposed scheme outperforms a number of methods in literature. For feature waveform classification, the paper presents a classification strategy based on geometric information. As shown in the simulation studies, the P waves can be classified into different numbers of clusters by setting the threshold to different values, which may help study the variability of feature waveforms.

Conflict of interest No conflict of interest.

Acknowledgements This project was supported by National Natural Science Foundation of China (no. 51207058) and Guangdong Innovative Research Team Program (no. 201001N0104744201).

Appendix A. A.1. Mathematical morphology Mathematical morphology (MM) is a relatively new and powerful tool for signal as well as image processing. It involves the interaction between the original signal or image and its corresponding structuring element (SE). The SE corresponds to another signal or image, which is much smaller in size than that of the original, which is being processed. Based on the concept of set theory, MM regards both the signal/image and the SE as sets of samples/pixels. Thus, the interaction comes from set theory, such as translation, reflection, union and intersection. Through the application of appropriate SEs, object features can be extracted via the implementation of MM. Dilation and erosion are the two basic morphological operators, which are denoted by ⊕ and respectively. The dilation and erosion of a signal f(t) by an SE g(s) are defined as follows, respectively: (f ⊕ g)(t) = max{f (t + s) + g(s)}

(12)

(f g)(t) = min{f (t + s) − g(s)}.

(13)

s

s

For image processing, 2-dimensional dilation and erosion can be expressed by: (f ⊕ g)(x) = maxf (x + b)

(14)

(f g)(x) = minf (x + b)

(15)

b∈g

b∈g

where f represents the image and g corresponds to the SE. Such an operation can be interpreted as placing the center of the SE over each sample of the signal, dilation returns the maximum value within the window defined by the SE, while erosion returns the minimum. Two other popular morphological operators are opening and closing, defined respectively as: f ◦ g = (f g) ⊕ g

(16)

f • g = (f ⊕ g) g.

(17)

An important property of these operators is that they can remove sharp peaks and valleys respectively. Opening and closing are morphological filters because they are both

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 1 1 2 ( 2 0 1 3 ) 466–480

increasing and idempotent. The property of increasing can be explained by the following:

 f1 ≤ f2 ⇒

f1 ◦ g ≤ f2 ◦ g f1 • g ≤ f2 • g

f •g•g = f •g

1  (r − |xi − xj |) (i = / j) N→∞ N2



dc = lim

X

⎢ x2 ⎥ ⎢ ⎥ ⎥ ⎣ .. ⎦

=⎢

.

xM

x1+

···

x1+(dE −1)

x2+

···

x2+(dE −1)

.

. . .

..

. . .

xM

xM+

···

⎢ x2 ⎢ ⎣ ..

=⎢

.

⎤ ⎥ ⎥=[x ⎥ 1 ⎦

x2

···

xdE ]

xM+(dE −1) (20)

where dE is the dimension of the phase space,  is the time delay, the column vectors xj (j = 1, · · · , dE ) form the coordinate of each dimension, the row vectors xi (i = 1, · · · , M) represent individual points in the phase space, and the matrix X denotes the embedded signal. Two parameters are required by the delay coordinate embedding: the embedding dimension and the time delay. These parameters should be properly chosen so that the feature hidden within the time series can be presented in the phase space. According to [27], a suggested embedding dimension is dE =2 · boxdim(A) + 1 , where A is the attractor of the dynamical system, boxdim(A) is the system dimension, and x denotes the minimum integer with a value larger than or equal to x. It should be noted that boxdim(A) may be fractional and dE must be integral. In this paper the correlation dimension [26,28] is used as boxdim(A) to calculate dE .

(22)

.

Eq. (21) calculates the number of pairs (xi , xj ) satisfying |xi − xj | < r. To calculate the real correlation dimension, dc , initial values of dE and  need to be selected. Assuming that for small r, C(r) behaves as C(r) ∝ rdc , dc can then be estimated by:

A.2. The embedding theorem A multi-dimensional dynamic system can be expressed by a set of first-order differential equations, for which the number of differential equations corresponds to the number of dimensions of the system. The solutions to these equations are called the states of the system. In order to obtain the characteristics of such a dynamic system, one needs to measure all the possible states. However, this is a challenging task which is often impossible to achieve. Fortunately, it has been proved that the dynamic system can be reconstructed from finite observations of the system [25,26]. Suppose time series x = {x1 , · · ·, xi , · · ·, xN } is observed from a dynamic system. The system can be reconstructed by embedding x to a dE -dimensional phase space. The reconstruction method proposed in [25,26] is the delay coordinate embedding, which is given by:

ifx ≤ 0

0,

1, ifx > 0

(19)

.

⎡ x1

(21)

where xi and xj are two arbitrary points, and (x) is the Heaviside step function: (x) =

Moreover, the opening operator is an anti-extensive filter: f ◦ g < f, while the closing operator is an extensive filter: f • g > f.

⎡ x1 ⎤

N

C(r) = lim

i=1 j=1

The property of idempotent implies that successive applications of openings or closings do not cause any further modifications to the signal: f ◦g◦g = f ◦g

The correlation dimension is determined from the correlation integral, which is defined as: N

(18)

.

479

r→0

lg[C(r)] . lg[r]

(23)

In other words, dc is the slope of the curve of lg [C(r)] as a function of lg [r]. The correlation dimension depends on the value of r. When r is small, the behavior of the correlation dimension is dominated by the characteristics of noise, which has infinite dimensions [29]. In practice, a string of dc are calculated with various values of r. For a range of rL ≤ r ≤ rU , dc (r) is a constant within some tolerance, the correlation dimension is chosen as the average of dc (r) over [rL , rU ] [26]. Finally, the embedding dimension is dE =2dc + 1 . The delay constant should neither be too short to include unnecessary calculation, nor too long so as to miss any useful information. If the value of  is too small, each coordinate is found to be almost the same and the trajectories of the phase space are squeezed along the identity line. On the other hand, if the value of  is too large, in the presence of chaos and noise, the dynamics at a given instant in time become effectively and causally disconnected from the dynamics at a later time instant. Hence, even simple geometric objects look extremely complicated. A natural choice of  is the first minimum of the autocorrelation function such that each coordinate is linearly independent. According to [30],  is equal to the first minimum of lgCd1E , where Cd1E = limCqdE

(24)

q→1

⎛ CdqE = ⎝

ˆl 1

ˆl

⎞(1/q−1) (Pr (xi ()))

q−1

⎠

(25)

i=1

1 (r − |xi () − xj ()|) ˆl ˆl

Pr (xi ()) =

(26)

j=1

ˆl = l − (dE − 1)

(27)

xi () = {xi , xi+ , · · ·, xi+(de −1) }.

(28)

480

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 1 1 2 ( 2 0 1 3 ) 466–480

references

[1] P.S. Addison, Wavelet transforms and the ECG: a review, Physiological Measurement 26 (2005) 155–199. [2] H.H. de Carvalho, R.L. Moreno, T.C. Pimenta, P.C. Crepaldi, E. Cintra, A heart disease recognition embedded system with fuzzy cluster algorithm, Computer Methods and Programs in Biomedicine 110 (3) (2013) 447–454. [3] T.Y. Ji, Z. Lu, Q.H. Wu, Multi-resolution morphological operators for electrocardiogram signal analysis, in: The 26th Chinese Control Conference, vol. 3, Zhangjiajie, China, 26–31 July, 2007, pp. 425–429. [4] F. Zhang, Y. Lian, QRS detection based on multiscale mathematical morphology for wearable ECG devices in body area networks, IEEE Transactions on Biomedical Circuits and Systems 3 (4) (2009) 220–228. [5] A. Kabir, C. Shahnaz, Denoising of ECG signals based on noise reduction algorithms in emd and wavelet domains, Biomedical Signal Processing and Control 7 (September (5)) (2012) 481–489. [6] C.-S. Chen, J.-L. Wu, Y.-P. Hung, Theoretical aspects of vertically invariant gray-level morphological operators and their application on adaptive signal and image filtering, IEEE Transactions on Signal Processing 47 (4) (1999) 1049–1060. [7] S. Suppappola, Y. Sun, Nonlinear transforms of ECG signals for digital QRS detection: a quantitative analysis, IEEE Transaction on Biomedical Engineering 41 (4) (1994) 397–400. [8] Z. Zidelmal, A. Amirou, M. Adnane, A. Belouchrani, QRS detection based on wavelet coefficients, Computer Methods and Programs in Biomedicine 107 (September (3)) (2012) 490–496. [9] A. Herreros, E. Baeyens, R. Johansson, J. Carlson, J.R. Perán, B. Olsson, Analysis of changes in the beat-to-beat P-wave morphology using clustering techniques, Biomedical Signal Processing and Control 4 (October (4)) (2009) 309–316. [10] A. Diery, D. Rowlands, T.R.H. Cutmore, D. James, Automated ECG diagnostic P-wave analysis using wavelets, Computer Methods and Programs in Biomedicine 101 (January (1)) (2011) 33–43. [11] Y. Kutlu, D. Kuntalp, Feature extraction for ECG heartbeats using higher order statistics of WPD coefficients, Computer Methods and Programs in Biomedicine 105 (3) (2012) 257–267. [12] P. Sun, Q.H. Wu, A.M. Weindling, A. Finkelstein, K. Ibrahim, An improved morphological approach to background normalization of ECG signals, IEEE Transactions on Biomedical Engineering 50 (1) (2003) 117–121. [13] M. Richter, T. Schreiber, Phase space embedding of electrocardiograms, Physical Review E 58 (5) (November 1998) 6392–6398. [14] A. Herreros, E. Baeyens, J.R. Perán, R. Johansson, J. Carlson, B. Olsson, An algorithm for phase-space detection of the P characteristic points, in: Proceedings of the 29th Annual

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22] [23]

[24]

[25] [26]

[27] [28]

[29] [30]

International Conference of the IEEE Engineering in Medicine and Biology Society, Lyon, France, 23–26 August, 2007, pp. 2004–2007. R. Silipo, C. marchesi, Artificial neural networks for automatic ECG analysis, IEEE Transactions on Signal Processing 46 (5) (1998) 1417–1425. R.V. Andreao, B. Dorizzi, J. Boudy, ECG signal analysis through hidden Markov models, IEEE Transactions on Biomedical Engineering 53 (8) (2006) 1541–1549. F.J. Chin, Q. Fang, T. Zhang, I. Cosic, A fast critical arrhythmic ECG waveform identification method using cross-correlation and multiple template matching, in: 32nd Annual International Conference of the IEEE EMBS, Buenos Aires, Argentina, 31 August–4 September, 2010, pp. 1922–1925. ˜ D. Cuesta-Frau, J.L. Rodríguez-Sotelo, D. Peluffo-Ordonz, Unsupervised feature relevance analysis applied to improve ECG heartbeat clustering, Computer Methods and Programs in Biomedicine 108 (October (1)) (2012) 250–261. Y.-C. Yeh, C.W. Chiou, H.-J. Lin, Analyzing ECG for cardiac arrhythmia using cluster analysis, Expert Systems with Applications 39 (January (1)) (2012) 1000–1010. A.S. Al-Fahoum, A.M. Qasaimeh, ECG arrhythmia classification using simple reconstructed phase space approach, Computers in Cardiology 33 (2006) 757–760. T.Y. Ji, W.H. Tang, Q.H. Wu, Partial discharge location using a hybrid transformer winding model with morphology-based noise removal, Electric Power Systems Research 101 (August) (2013) 9–16. P. Soille, Morphological Image Analysis, Springer, Berlin, 2003. Y. Sun, K.L. Chan, S.M. Krishnan, ECG signal conditioning by morphological filtering, Computers in Biology and Medicine 32 (November (60)) (2002) 465–479. T.Y. Ji, Z. Lu, Q.H. Wu, Z. Ji, Baseline normalisation of ECG signals using empirical mode decomposition and mathematical morphology, Electronics Letters 44 (January (2)) (2008) 82–83. F. Takens, Detecting strange attractors in turbulence, Lecture Notes in Mathematics 898 (1981) 366–381. A.M. Albano, J. Muench, C. Schwartz, A.I. Mees, P.E. Rapp, Singular-value decomposition and the Grassberger–Procaccia algorithm, Physical Review A 38 (September (6)) (1988) 3017–3026. T. Sauer, J.A. Yorke, M. Casdagli, Embedology, Journal of Statistical Physics 65 (November (3–4)) (1991) 579–616. P. Grassberger, I. Procaccia, Estimation of the Kolmogorov entropy from a chaotic signal, Physical Review A 28 (October (4)) (1983) 2591–2593. J.-P. Eckmann, D. Ruelle, Ergodic theory of chaos and strange attractors, Reviews of Modern Physics 57 (3) (1985) 617–656. W. Liebert, H.G. Schuster, Proper choice of the time delay for the analysis of chaotic time series, Physics Letter A 142 (2–3) (1989) 107–111.