High Sensitivity and Noise Immune Method to Detect Impedance Cardiography Characteristic Points Using Wavelet Transform

INTRODUCTION

Patients with chronic heart failures have periodic episodes of clinical decompensation that do not impair the quality of life only but also consume substantial health care resources such as hospitalization costs. Heart failure management programs have been developed to reduce frequency and severity of these clinical events, but their effectiveness may be limited by physician difficulty in identifying patients at imminent risk. Reliable prediction of these events may avoid the clinicians to intervene aggressively and may minimize potentially the need of hospitalization. It has been hypothesized that the check of hemodynamic variable values might identify critical patients (Packer et al., 2006).

Frequent evaluation and use of these parameters is limited because of the invasiveness of the most accurate investigation techniques. Today, the introduction of impedance cardiography (ICG) has overcome this limitation. It has been widely adopted because it is noninvasive, easy to use and suitable for long-term and continuous monitoring of hemodynamic function (Jensen et al., 1995).

ICG principle is based on measures of the thorax resistance variations caused both by air flow through lungs and blood flow from left ventricle to aorta. The value of these variations makes the recreation of cardiac cycle process possible. Moreover, unexpected increasing or decreasing of impedance can be related to pathological causes such as heart decompensation, heart and valve failures, hypertension and other chronic diseases (Treister et al., 2005).

Wang et al. (1995) have been developed both to study physiological mechanisms for understanding origin and meaning of ICG signals and to improve effectiveness and applicability of ICG diagnostic test adopting advanced signal processing techniques.

In particular, many efforts have been done to implement automatic detection of reference points in biological signals. However, existing peak detection algorithms are difficult to automate for generic use because either they rely on a number of parameters that need to be customized for a particular application or they use reference informations that are highly specialized for a particular application.

Most of the proposed methods make use of filtering technique (band pass filtering and temporal filtering) (Leski and Tkacz, 1992; Pan and Tompkins, 1984), or adaptive thresholding technique (Sun et al., 2005; Suppappola and Sun, 1994). All the previous techniques exhibit limitations when real signal are adopted (Sun et al., 2005). In fact, the first drawback of filtering-based approach is that frequency variations in the signal under test (due to different causes such as, for instance, cardiac frequency changes) may adversely affect the method performance. For instance, the frequency band of some biological signal, such as ECG, differs for different subject and can change for the same subject due to particular events. The second problem in the filter based algorithms is the frequency band overlapping of noise and some biological signals. Therefore, the choice of a suitable bandwidth is a trade off between noise reduction and loss of high frequency details while the duration of the sliding window is a trade off between false and missed detections.

Whereas, the main problem of thresholding techniques are their sensitivity to baseline variations and signal intensity. This high noise sensitivity can be a problem for some types of signals having low signal to noise (S/N) ratio.

An extensive overview of the various algorithm for signal peak detention can be found in (Kohler et al., 2002) which includes approaches based on neural networks, adaptative filters, Hidden Markov models, Hilbert transform, too.

In this study, an improved signal processing technique suited for the implementation in design tools, is presented. The new technique implements a fast parallelized algorithm based on singularity detector adopting a multiscale wavelet transform. The obtained simulation results show the method independence from noise and artefacts.

ICG TECHNIQUE

Impedance Cardiography is the study of the cardiac function by means of thorax electrical impedance measurements. High frequency (20-100 kHz), low intensity current (1-5 mA rms) is injected through the thorax by some electrodes and the impedance change is sensed by measuring a voltage across other electrodes. No risk of physiological effects is present because various tissues of human body are not excitable at this frequency and at this low current level (Patterson, 1989). The impedance variation can be used for diagnostic informations and for Stroke Volume (SV) estimation by using blood flow appropriate model. The term SV indicates the amount of blood pumped by the heart left ventricle in one contraction.

Figure 1 shows a typical impedance waveform obtained from electrodes in which the characteristic points are indicated.

Steady-state basal impedance comprises the major component of total thoracic impedance and two oscillating components: the respiratory and the cardiac component. Since air is a poor electrical conductor, the impedance (Z) increases during inspiration and decreases during expiration, resulting in an oscillation of Z equal to the respiratory rate.

Pulsating blood flow through the thoracic aorta causes shifts in the thoracic impedance as a function of changes in blood volume. In fact, since blood is the most conductive tissue between the electrodes, the rise of blood volume causes the reduction of Z value. Moreover, as blood flow increases, Z decreases since the blood speed causes alignment of erythrocytes to produce a more conductive path.

Fig. 1:	Typical impedance waveforms from the thorax of a human subject, where, A: Atrial activity, B: Synchronous wave with first heart sound, C: Largest decrease in impedance during systole, X: Aortic valve closure, Y: Pulmunaary valve closure and O: Largest decrease in impedance during diastolw close in time to mitral valve opening

This oscillating component of total thoracic impedance is displayed as pulsating variation in impedance that can be expressed adopting the derivative (dZ/dt). Measurements of thoracic impedance changes (dZ/dt waveform) during the cardiac cycle are used to calculate SV. There are several ways to calculate SV (Bernstein, 1986; Kubicek et al., 1974; Sramek, 1982) generally all the equations take into account position and value of the C-point related to the B-point and the X-point.

In this study, a fast parallelized real time method for the exact localization of C points is presented.

THE WAVELET TRANSFORM

Biomedical signals are generally non stationary with significant events characterized by both time location and frequency content. Therefore, the frequency analysis of such signals through Fourier technique is unsatisfactory, because it is efficient in providing frequency discrimination but poor in time localization. Vice versa the wavelet transform provides temporal and spectral information simultaneously, so it is suitable for determining characteristic points of non stationary signals, such as the ICG signal. It decomposes a time variant signal into several components having various scales or resolutions. A suitable time and frequency limited wavelet is chosen as the mother. Scaling and shifting the mother wavelet, a family of functions called daughter wavelets is generated. For small values of the scale factor, wavelet is constructed in the time domain and gives informations about signal fine details. Therefore, a global view of signal is obtained by large value of scale factor (Haddad, 2003). The wavelet transform of a time signal at any scale is the convolution of the signal and a time-scaled daughter wavelet.

There are essentially two types of wavelet decompositions: the redundant ones (continuous wavelet transform) and the nonredundant ones (orthogonal, semi-orthogonal, or biorthogonal wavelet bases) (Unser and Aldroubi, 1996). The first type is preferable for feature extraction because it provides a description that is truly shift-invariant. The second type is preferable for data reduction or when the orthogonality of the representation is an important factor. However, the choice between these types of decompositions has to take into account computational considerations, too. A decomposition in terms of wavelet bases using Mallat fast algorithm is typically orders of magnitude faster than a redundant analysis, even if the fastest available algorithms are used (Rioul and Duhamel, 1992; Unser, 1994).

As the aim of this paper is the implementation of a fast parallelized algorithm, a nonredundant wavelet decomposition has been chosen. To determine the best wavelet function to be used, the ICG signal properties have been studied, such as the shape and the time localization of characteristic events. Since the temporal signal shape is an important parameter, orthogonal wavelets are unsuitable to be used. In fact they are unable to provide symmetry in the time domain as a non-linear phase shift is introduced during the analysis. Signal shape is maintained if phase shift is linear. Thus, the wavelet function to be adopted should be a symmetrical function (Dinh, 2001). The Spline wavelets have properties satisfying the previous requirements. The higher order of the Spline wavelet results in the sharper frequency response of the equivalent FIR filter, that is always desirable. But the FIR equivalent filter of the higher order Spline wavelet has longer coefficient series leading to larger computational time consumption. Therefore, the Cubic Spline wavelet is assumed to have an order high enough for this application.

The traditional wavelet theory considers the decomposition algorithm with an iterative structure (in particular an asymmetrical tree structure) that does not efficiently merge with the novel computational techniques, such as parallel processing, concurrent programming and implementation in design tools (Cohen and Kovacevic, 1996). In this study, the a’ trous and the Mallat algorithms for parallelized filter bank design are used (Yang et al., 2005). The algorithm generates a set of parallelized perfect-reconstruction filter banks for an arbitrary number of end-nodes of a traditional tree structure (Koh and Rodriguez-Marek, 2003).

Method description: The method presented in this section processes the first derivative of the impedance cardiography signal and allows to determine the C peak positions in time domain (Fig. 1).

The tested ICG signal has been sampled at 250 Hz (Fig. 2).

For the method implementation, named C_point_DETECTOR, no external trigger source is necessary; therefore the ICG signal is the only input.

Figure 3 shows the adopted algorithm model in which the obtained results are indicated as:

The method is carried out in two steps (Fig. 4): Wavelet decomposition and C_Point_Detection.

Fig. 2:	ICG signal

Fig. 3:	Algorithm model realized with the software tool MATLAB Simulink®

Fig. 4:	‘C point_DETECTOR’ System model

Fig. 5:	Wavelet ‘bior3.3’. (a) Decomposition scaling function phi, (b) Decomposition wavelet function psi, (c) Decomposition low-pass filter, (d) Decomposition high-pass filter, (e) Reconstruction scaling function phi, (f) Reconstruction wavelet function psi, (g) Reconstruction low-pass filter and (h) Reconstruction high-pass filter

Wavelet decomposition: A fit mother wavelet is fundamental for high performance of the proposed ICG detection technique. A wide variety of functions can be chosen as mother wavelet but wavelet bior 3.3 has been adopted because it allows the perfect signal reconstruction keeping phase shift linear. Wavelet bior 3.3 is a cubic spline which is characterized by high frequency splitting and low band coherence (Fig. 5a-h).

To locate C points, the method decomposes the ICG signal into six dyadic scales so to reduce noise sensitivity, significantly (Fig. 6). It is evident that various frequency components of input signal are presented in the different scales. Lower frequency components are in high level scales while higher frequency components are in low level scales.

According to the power spectra of ICG signal, noise and artefacts, it is evident that the larger contribute of the true signal is located in scales 4, 5 and 6, while the scales 1, 2 and 3 are mostly affected by noise (Xu and Liu, 2005).

Adopting a soft thresholding technique applied to levels 1, 2, 3, noise is reduced and then signal is reconstructed in time domain for further processing.

ICG signal representation in wavelet domain is particularly indicated for noisy signals. In fact for searching characteristic points, processing of noise free scales only, such as scales 4-5-6, is sufficient. The method greater advantage is evident in comparison with frequency analysis which needs informations of the overall bandwidths for the C point localizations.

The method uses an evolution of the classical Mallat decomposition, called a’ trous algorithm. The a’ trous algorithm for non-orthogonal wavelets uses the same filter bank structure as the Mallat algorithm (Mallat, 1989), but differs for high pass and low pass FIR filters. It has been demonstrated that after the application of wavelet filters for j-times, the precision of a’ trous algorithm is 2j time higher then Mallat algorithm (Table 1) (Wang et al., 1995; Shensa, 1992).

Fig. 6:	Decomposition of ICG signal over six diadic scales using a’ trous algorithm

Fig. 7:	Method procedure for too closed peaks management

Table 1:	Precision of the Mallat and the a’ trous algorithms varying decomposition levels

Usually, wavelet decomposition algorithms make use of filters in a tree structure. This is unsuitable for both the parallel computing and the implementation in design tools. To overcome these limits, equivalent parallel filter banks are used. However, as it is well known, the output signal realignment is necessary to equalize the delay introduced by each filter. This structure makes the algorithm attractive for a hardware implementation (Koh and Rodriguez-Marek, 2003).

C_point_Detection: At the second step the method first spots the positions of all ICG signal peaks, then the C-point locations.

A point of maximum value is present in component signals (scales 4, 5, 6) in the same location of each singularity in ICG signal. Adopting a parallel algorithm, the proposed method looks for local maximum points inside scale 4, scale 5 and scale 6 using a thresholding technique: only maximum points higher than the adopted threshold value are taken into account. Occurrence of too close peaks is solved replacing them with one new virtual peak located in a middle position between them. Figure 7 shows the spatial representation of local maximum points after the replacement of too close maximums. For instance, peaks P and Q of scale 4 selected by the method as too close points, are replaced by the point R in a new spatial representation named scale 4a.

Each processed scale gives information about the same signal, so points located inside each scale are related to each others. For this reason, the method executes a tracing across scales to locate overlapping singular points. However, as the perfect coincidence of two or more points is almost impossible, a confidence interval of occurrence is associated to each point.

For software implementation, a window having a suitable size is constructed around each singular point considering one scale at a time.

The generic window in a generic scale is analyzed to determine if peaks of other scales are located inside it. The research is iterated for both all the windows of a scale and for all the scales. The method defines as valid peaks, those that are present in two out of three scales at least. The best estimate of the valid peak position (the true C-point position) is obtained averaging peak positions inside the different scales.

Figure 8a shows the characterization of a new spatial representation named scale 4b for peaks of scale 4a consequently to the definition of a confidence interval for each point. For instance, the method defines a window for peak S and verifies if peaks of other scales are located inside this window. In this example, point T of scale 6a is inside the window. Therefore point U is defined as valid value in the new representation named scale 4b. U point location is in the middle between S and T positions. Figure 8b and c show the same analysis for points belonging to scales 5a and 6a, respectively.

Points of scales 4b, 5b and 6b are grouped in one spatial representation indicated in Fig. 9 with TEMP. Points that are too closed in TEMP are replaced with one peak having the middle location as shown in representation named OUTPUT.

Fig. 8:	(a-c) Parallel procedure for peak validation

Fig. 9:	Procedure last step for C point location

In conclusion, the following steps characterize the new method: