Advanced Methods and Tools for ECG Data Analysis - Part 5

5.3 Wavelet Filtering 145 Figure 5.5 The effect of a selection of different wavelets for filtering a section of ECG (using the first approximation only) contaminated by Gaussian pink noise (SNR = 20 dB). From top to bottom; original (clean) ECG, noisy ECG, biorthogonal (8,4) filtered, discrete Meyer filtered, Coiflet filtered, symlet (6,6) filtered, symlet filtered (4,4), Daubechies (4,4) filtered, reverse biorthogonal (3,5), re-verse biorthogonal (4,8), Haar filtered, and biorthogonal (6,2) filtered. The zero-noise clean ECG is created by averaging 1,228 R-peak aligned, 1-second-long segments of the author’sECG. RMS error performance of each filter is listed in Table 5.1. to the length of the highpass filter. Therefore Matlab’s bior4.4 has four vanishing moments3 with 9 LP and 7 HP coefficients (or taps) in each of the filters. Figure 5.5 illustrates the effect of using different mother wavelets to filter a section of clean (zero-noise) ECG, using only the first approximation of each wavelet decomposition. The clean (upper) ECG is created by averaging 1,228 R-peak aligned, 1-second-long segments of the author’s ECG. Gaussian pink noise is then added with a signal-to-noise ratio (SNR) of 20 dB. The root mean square (RMS) error between the filtered waveform and the original clean ECG for each wavelet is given in Table 5.1. Note that the biorthogonal wavelets with J,K ≥ 8, 4, 3. If the Fourier transform of the wavelet is J continuously differentiable, then the wavelet has J vanishing moments. Type waveinfo(bior) at the Matlab prompt for more information. Viewing the filters using [lpdecon, hpdecon,lp econ, hp econ] = wfilters(bior4.4) in Matlab reveals one zero coefficient in each of the LP decomposition and HP reconstruction filters, and three zeros in the LP reconstruction and HP decomposition filters. Note that these zeros are simply padded and do not count when calculating the filter size. 146 Linear Filtering Methods Table 5.1 Signals Displayed in Figure 5.5 (from Top to Bottom) with RMS Error Between Clean and Wavelet Filtered ECG with 20-dB Additive Gaussian Pink Noise Wavelet Family Original ECG ECG with pink noise Biorthogonal ‘bior’ Discrete Meyer ‘dmey’ Coiflets ‘coif’ Symlets ‘sym’ Symlets ‘sym’ Daubechies ‘db’ Reverse biorthogonal ‘rbio’ Reverse biorthogonal ‘rbio’ Haar ‘haar’ Biorthogonal ‘bior’ N/A indicates not applicable. Family Member N/A N/A bior3.3 dmey coif2 sym3 sym2 db2 rbio3.3 rbio2.2 haar bior1.3 RMS Error 0 0.3190 0.0296 0.0296 0.0297 0.0312 0.0312 0.0312 0.0322 0.0356 0.0462 0.0472 the discrete Meyer wavelet and the Coiflets appear to produce the best filtering performance in this circumstance. The RMS results agree with visual inspection, where significant morphological distortions can be seen for the other filtered sig-nals. In general, increasing the number of taps in the filter produces a lower error filter. The wavelet transform can be considered either as a spectral filter applied over many time scales, or viewed as a linear time filter [(t −τ)/a] centered at a time τ with scale a that is convolved with the time series x(t). Therefore, convolving the filters with a shape more commensurate with that of the ECG produces a better filter. Figure 5.4 illustrates this point. Note that as we increase the number of taps in the filter, the mother wavelet begins to resemble the ECG’s P-QRS-T morphol-ogy more closely. The biorthogonal wavelet family members are FIR filters and, therefore, possess a linear phase response, which is an important characteristic for signal and image reconstruction. In general, biorthogonal spline wavelets allow ex-act reconstruction of the decomposed signal. This is not possible using orthogonal wavelets (except for the Haar wavelet). Therefore, bior 3.3 is a good choice for a general ECG filter. It should be noted that the filtering performance of each wavelet will be different for different types of noise, and an adaptive wavelet-switching pro-cedure may be appropriate. As with all filters, the wavelet performance may also be application-specific, and a sensitivity analysis on the ECG feature of interest is appropriate (e.g., QT interval or ST level) before selecting a particular wavelet. As a practical example of comparing different common filtering types to the ECG, observe Figure 5.6. The upper trace illustrates an unfiltered recording of a V5 ECG lead from a 30-year-old healthy adult male undergoing an exercise test. Note the presence of high amplitude 50-Hz (mains) noise. The second subplot illustrates the action of applying a 3-tap IIR notch-filter centered on 50 Hz, to reveal the underlying ECG. Note the presence of baseline wander disturbance from electrodemotionaroundt = 467seconds,andthedifficultyindiscerningthePwave (indicated by a large arrow at the far left). The third trace is a band-pass (0.1 to 45 Hz) FIR filtered version of the upper trace. Note the baseline wander is reduced 5.3 Wavelet Filtering 147 Figure 5.6 Raw ECG with 50 Hz mains noise, IIR 50-Hz notch filtered ECG, 0.1- to 45-Hz band-pass filtered ECG and bior3.3 wavelet filtered ECG. The left-most arrow indicates the low amplitude P wave. Central arrows indicate Gibbs oscillations in the FIR filter causing a distortion larger than the P wave. significantly, but a Gibbs4 ringing phenomena is introduced into the Q and S waves (illustrated by the small arrows), which manifests as distortions with an amplitude larger than the P wave itself. A good demonstration of the Gibbs phenomenon can be found in [9, 10]. This ringing can lead to significant problems for a QRS detector (looking for Q wave onset) or any technique for analyzing at QT intervals or ST changes. The lower trace is the first approximation of a biorthogonal wavelet decomposition (bior3.3) of the notch-filtered ECG. Note that the P wave is now discernible from the background noise and the Gibbs oscillations are not present. As mentioned at the start of this section, the number of articles on ECG analysis that employ wavelets is enormous and an excellent overview of many of the key publications in this arena can be found in Addison [5]. Wavelet filtering is a lossless supervised filtering method where the basis functions are chosen a priori, much like the case of a Fourier-based filter (although some of the wavelets do not have orthogonal basis functions). Unfortunately, it is difficult to remove in-band noise because the CWT and DWT are signal separation methods that effectively occur in 4. The existence of the ripples with amplitudes independent of the filter length. Increasing the filter length narrowsthetransitionwidthbutdoesnotaffecttheripple.Onetechniquetoreducetheripplesistomultiply the impulse response of an ideal filter by a tapered window. 148 Linear Filtering Methods thefrequencydomain5 (ECGsignalandnoisesoftenhaveasignificantoverlapinthe frequency domain). In the next section we will look at techniques that discover the basis functions within data, based either on the statistics of the signal’s distributions or with reference to a known signal model. The basis functions may overlap in the frequency domain, and therefore, we may separate out in-band noise. As a postscript to this section, it should be noted that there has been much discussion of the use of wavelets in HRV analysis (see Chapter 3) since long-range beat-to-beat fluctuations are obviously nonstationary. Unfortunately, very little at-tention has been paid to the unevenly sampled nature of the RR interval time series and this can lead to serious errors (see Chapter 3). Techniques for wavelet analy-sis of unevenly sampled data do exist [11, 12], but it is not clear how a discrete filter bank formulation with up-down sampling could avoid the inherent problems of resampling an unevenly sampled signal. A recently proposed alternative JTFA technique known as the Hilbert-Huang transform (HHT) [13, 14], which is based upon empirical mode decomposition (EMD), has shown promise in the area of non-stationary and nonlinear JFTA (since both the amplitude and frequency terms are a function of time6). Furthermore, there is striking similarity between EMD and the least-squares estimation technique used in calculating the Lomb-Scargle Peri-odogram (LSP) for power spectral density estimation of unevenly sampled signals (see Chapter 3). EMD attempts to find basis functions (such as the sines and cosines in the LSP) by fitting them to the signal and then subtracting them, in much the same manner as in the calculation of the LSP (with the difference being that EMD analyzes the envelope of the signal and does not restrict the basis functions to be-ing sinusoidal). It is therefore logical to extend the HHT technique to fit empirical modes to an unevenly sampled times series such as the RR tachogram. If the fit is optimal in a least-squares sense, then the basis functions will remain orthogonal (as we shall discover in the next section). Of course, the basis functions may not be orthogonal, and other measures for optimal fits may be employed. This concept is explored further in Section 5.4.3.2. 5.4 Data-Determined Basis Functions Sections 5.4.1 to 5.4.3 present a set of transformation techniques for filtering or separating signals without using any prior knowledge of the spectral components of the signals and are based upon a statistical analysis to discover the underlying basis functions of a set of signals. These transformation techniques are principal component analysis7 (PCA), artificial neural networks (ANNs), and independent component analysis (ICA). 5. The wavelet is convolved with the signal. 6. Interestingly, the empirical modes of the HHT are also determined by the data and are therefore a special case where a JTFA technique (the Hilbert transform) is combined with a data-determined empirical mode decompositiontoderiveorthogonalbasisfunctionsthatmayoverlapinthefrequencydomaininanonlinear manner. 7. This is also known as singular value decomposition (SVD), the Hotelling transform or the Karhunen-Loeve transform (KLT). 5.4 Data-Determined Basis Functions 149 Both PCA and ICA attempt to find an independent set of vectors onto which we can transform data. Those data that are projected (or mapped) onto each vector are the independent sources. The basic goal in PCA is to decorrelate the signal by projecting data onto orthogonal axes. However, ICA results in a transformation of data onto a set of axes which are not necessarily orthogonal. Both PCA and ICA can be used to perform lossy or lossless transformations by multiplying the recorded (observation) data by a separation or demixing matrix. Lossless PCA and ICA both involve projecting data onto a set of axes which are determined by the nature of those data, and are therefore methods of blind source separation (BSS). (Blind because the axes of projection and therefore the sources are determined through the application of an internal measure and without the use of any prior knowledge of a signal’s structure.) Once we have discovered the axes of the independent components in a data set and have separated them out by projecting the data set onto these axes, we can then use these techniques to filter the data set. 5.4.1 Principal Component Analysis To determine the principal components (PCs) of a multidimensional signal, we can use the method of singular value decomposition. Consider a real N× M matrix X of observations which may be decomposed as follows: X = USVT (5.8) where S is an N× M nonsquare matrix with zero entries everywhere, except on the leading diagonal with elements si(= Snm, n = m) arranged in descending order of magnitude. Each si is equal to λi, the square root of the eigenvalues of C = XTX. A stem-plot of these values against their index i is known as the singular spectrum. The smaller the eigenvalues are, the less energy along the corresponding eigenvector there is. Therefore, the smallest eigenvalues are often considered to be associated with the noise in the signal. V is an M× M matrix of column vectors which are the eigenvectors of C. U is an N×Nmatrix of projections of X onto the eigenvectors of C [15]. If a truncated SVD of X is performed (i.e. we just retain the most significant p eigenvectors),8 then the truncated SVD is given by Y = USpVT, and the columns of the N× M matrix Y are the noise-reduced signal (see Figure 5.7). SVD is a commonly employed technique to compress and/or filter the ECG. In particular, if we align M heartbeats, each N samples long, in a matrix (of size N × M), we can compress it down (into an N × p) matrix, using only the first p << M PCs. If we then reconstruct the set of heartbeats by inverting the reduced rank matrix, we effectively filter the original ECG. Figure 5.7(a) illustrates a set of 20 heartbeat waveforms which have been cut into 1-second segments (with a sampling frequency Fs = 256 Hz), aligned by their R peaks and placed side by side to form a 256 × 20 matrix. Therefore, the data set is 20-dimensional, and an SVD will lead to 20 eigenvectors. Figure 5.7(b) is 8. In practice choosing the value of p depends on the nature of the data set, but is often taken to be the knee in the eigenspectrum or as the value where i=1 si > α i=1 si and α is some fraction ≈ 0.95. ... - tailieumienphi.vn