Coronary Artery Disease (CAD) is a leading type of heart disease in the world caused by the gradual build-up of plaque on the walls of the arteries. Due to CAD’s high incidence rate and mortality, it is very harmful to human health. CAD can develop slowly and silently over years without any symptoms. Early diagnose of CAD is one of the most important medical research areas. Diastolic murmurs that occur as additional components in the heart sound signal provide clinicians with valuable diagnostic and prognostic information about the function of heart valves. When coronary arteries become narrowed or blocked, the turbulence appears which is produced by blood moving across the stenotic arteries. During the relatively quiet diastolic period of the cardiac cycle, the murmurs are likely to be loudest when coronary blood flow is maximal. Initial studies show that diastolic murmurs produced by coronary arterial stenosis contain higher frequency components.
The heart sound signal represents the mechanical activity of the cardiohemic system, which is complicated and non-stationary. It contains physiological and pathological information between the heart and the various parts of the body, so it can be used in diagnosis of heart disease. Heart sound has been widely used in diagnosis of heart disease and many methods have been adopted to aid the diagnosis [1, 2]. The heart sound signal generally can be separated into four parts: the 1st heart sound S1, the systolic period, the 2nd heart sound S2 and the diastolic period, shown in figure 1.
Diastolic murmurs occur between S2 and the next S1 when the heart muscle relaxes between beats. Heart murmurs are usually considered pathological. They can be caused by some kinds of heart attacks, such as coronary artery stenosis, aortic regurgitation, etc. Diastolic murmurs can provide clinicians with valuable diagnostic and prognostic information about the function of heart valves.
Short Time Fourier Transform, Wigner-Ville Distribution and Wavelet Transform, etc., have some inherent limitations [3, 4, 5]. Short Time Fourier Transform involves an intrinsic trade-off between time resolution and frequency resolution. In Wigner-Ville distribution, the inherent cross-term interferences often mask the true time-frequency information associated with the signal of interest. The wavelet transform has received considerable attention in recent years. It provides a multi-resolution representation of signals, however, it is not adaptive in nature; once the wavelet mother function is given, one will have to use it to analyse all the data. In addition, the wavelet transform also underlies an uncertainty principle. In 1998, Dr.Norden Huang proposed a novel signal processing algorithm: the Hilbert Huang Transform (HHT) [6, 7]. It has proved to be a powerful tool to analyse non-stationary and nonlinear signals. The key parts of HHT are the Empirical Mode Decomposition (EMD) and Hilbert transform. EMD can decompose adaptively diastolic murmurs into a finite and usually small number of Intrinsic Mode Functions (IMFs) that admit a well-behaved Hilbert transform. The Hilbert transform of IMFs can yield instantaneous frequency and instantaneous amplitude. The local energy and instantaneous frequency derived from the IMFs give the fine-resolution frequency-time distribution of the energy that is designated as the Hilbert spectrum. The three-dimensional distribution can reflect the inherent essential characteristic of the signal.
The paper is organized as follows: section 2 introduces generalized wavelet shrinkage denoising method. In section 3, the Hilbert spectrum based on EMD and marginal spectrum distributions of diastolic murmurs are studied; a new method to restrict the end effect of EMD is proposed in section 4.In section 5, the algorithm based on the Empirical Mode Decomposition (EMD) and Teager Energy Operation (TEO) is proposed as an effective approach for estimating the instantaneous frequency of diastolic murmurs. Finally, some conclusions are given in section 6.
2. Wavelet shrinkage method
We consider the following model of a discrete noisy signal:
The vector x represents noisy signal and is an unknown original clean signal. z is independent identity distribution Gaussian white noise with mean zero and unit variance. is intensity of noise. For simplicity, we assume intensity of noise is one.
The step of wavelet shrinkage is defined as follows:
Apply discrete wavelet transform to observe noisy signals.
Estimate noise and threshold value, thresholding the wavelet coefficients of the observed signal.
Apply the inverse discrete wavelet transform to reconstruct the signal.
The wavelet shrinkage method relies on the basic idea that the energy of signal will often be concentrated in a few coefficients in the wavelet domain while the energy of noise is spread among all coefficients in the wavelet domain. Therefore, the nonlinear shrinkage function in the wavelet domain will tend to keep a few larger coefficients over threshold value that represent the signal, while noise coefficients’ down threshold value will tend to reduce to zero.
In wavelet shrinkage, how to select the threshold function and how to select the threshold value are most crucial. Donohue introduced two kinds of thresholding functions: ‘hard threshold function’ and ‘soft threshold function’.
The hard threshold function (2) results in larger variance and can be unstable because of the discontinuous function. The soft threshold function (3) results in unnecessary bias due to shrinkage of the large coefficients to zero. We build the generalized threshold function:
When m is an even number:
When m is odd number:
When m=1, it is the soft threshold function; when m= , it is the hard threshold function. When m=2 it is Non-Negative Garrote threshold function. We show slope signal as an example, Figure2 illustrates the generalized threshold functions for different m.
It can clearly be seen that when the coefficient is small, the smaller m is, the closer the generalized function is to the soft threshold function; when the coefficient is big, the bigger m is, the closer the generalized function is to the hard threshold function; when m lies between 1 and , the general threshold function achieves a compromise between hard and soft threshold function. With careful selection of m, we can achieve better denoising performance [8, 9].
We derived the exact formula of mean, bias, variance and risk for the generalized threshold function.
And are density and probability function of standard Gaussian random variable respectively. Then:
, ,, are the mean, bias, variance and risk of generalized threshold function. When m is 1, 2, , they are the mean, bias, variance and risk of the risk of soft, Non-Negative Garrote, hard threshold functions, respectively.
The soft threshold function provides smoother results in comparison with the hard threshold function; however, the hard threshold function provides better edge preservation in comparison with the soft threshold function. The hard threshold function is discontinuous and this leads to the oscillation of denoised signal. The soft threshold function tends to have bigger bias because of shrinkage, whereas the hard threshold function tends to have bigger variance because of discontinuity. The Non-Negative Garrote threshold function is the trade-off between the hard and soft threshold function. Firstly, it is continuous; secondly, the shrinkage amplitude is smaller than the soft threshold function.
Stein Unbiased Risk Estimate (SURE)  is an adaptive threshold selection rule which is data driven. The threshold value minimizes an estimate of the risk.If is weakly differentiable, for single coefficient, is true risk. Then
is the unbiased risk estimate of :
Where and by partial integral
is the unbiased risk estimate of true risk .
For the generalized threshold function (5) and single coefficient, when m is even,
The SURE is
When m is odd,
The SURE is
Suppose wavelet coefficients are , the threshold value λ is set to minimize the estimate of the risk for the given data,
For hard threshold function (m is ∞), because H (x) is discontinuity, the SURE is illogical.
The noisy PCG signal is processed using the method mentioned above. For the generalized threshold functions, parameter m is selected as 2 which is simple and provides a good compromise between the hard and soft threshold function. The data-driven SURE threshold value is used. The filtered PCG signal is illustrated as figure 4(a). The phase space diagram of the filtered PCG signal is shown in figure 4(b). From visual inspection of figure 3, the PCG signal is much cleaner after being denoised; the first heart sound, the systolic period, the second heart sound and the diastolic period can be clearly identified. The results indicate that the method we have proposed significantly reduces noise and preserves well the characteristics of the PCG signal.
3. Analysis of diastolic murmurs for coronary artery disease based on empirical mode decomposition
Since a novel signal processing algorithm - the Hilbert HuangTransform (HHT) - was proposed by N.E.Huang in 1998 , it has been seen as a data-driven tool for nonlinear and non-stationary signal processing. HHT consists of two parts: the EMD and Hilbert transform. EMD as the important part of the HHT that can adaptively decompose signal into a finite and often a series of small numbers of Intrinsic Mode Functions (IMFs) subjected to the following two conditions:
In the whole dataset, the number of extrema and the number of zero-crossing must either be equal or differ at most by one.
At any time, the mean value of the envelope of the local maxima and the envelope of the local minima must be zero.
These two conditions guarantee the well-behaved Hilbert transform. The IMFs represent the oscillatory modes embedded in the signal. Most signals include more than one oscillatory mode and are not IMFs. EMD is a numerical sifting process to decompose a signal into a finite number of hidden fundamental intrinsic oscillatory modes, i.e., IMFs. Applying the Hilbert transform to each IMF, the instantaneous frequency and amplitude of each IMF can be obtained which constitute the time-frequency-energy distribution of the signal, called the Hilbert spectrum. The Hilbert spectrum provides higher resolution and concentration in the time-frequency plane, and avoids the false high frequency and energy dispersion existent in the Fourier spectrum.
Figure5 shows a classical IMF. The IMFs represent the oscillatory modes embedded in the signal. Each IMF actually is a zero mean monocomponents AM-FM signal with the following form:
with time varying amplitude envelope and phase . The amplitude and phase both have physical and mathematical meaning.
Most signals include more than one oscillatory mode, so they are not IMFs. EMD is a numerical sifting process to disintegrate empirically a signal into a finite number of hidden fundamental intrinsic oscillatory modes, that is, IMFs. The sifting process can be separated into the following steps:
Finding all the local extrema, including maxima and minima; then connecting all the maxima and minima of signal x(t) using smooth cubic splines to get its upper envelope and lower envelope .
Subtracting mean of these two envelopes from the signal to get their difference: .
Regarding the as the new data and repeating steps 1 and 2 until the resulting signal meets the two criteria of an IMF, defined as . The first IMF contains the highest frequency component of the signal. The residual signal is given by .
Regarding as new data and repeating steps (1) (2) (3) until extracting all the IMFs. The sifting procedure is terminated until the m-th residue becomes less than a predetermined small number or becomes monotonic.
The original signal x (t) can thus be expressed as follows:
Heart sounds are recorded from the chest of normal objects and CAD patients using a specially designed high sensitivity cardiac microphone. The ECG signals are also recorded as a time reference to aid in locating the diastolic phase. For each cycle, the central portion of diastole is digitized (sample frequency equals 2.0 kHz).
Figure6 shows the diastolic murmurs of a normal object. Figure7 shows the IMFs of the murmur obtained by EMD. The diastolic murmurs can be decomposed into six IMFs. The Hilbert spectrum is shown in figure 8. The vertical bars on the right of the panel give the relative amplitude scale. Figure6 provides more distinct information on the time-frequency contents of diastolic murmurs, which reveals clearly the dynamic characteristic of murmurs in the time-frequency plane. The Hilbert spectrum contains no energy with frequency above 350Hz. The spectrum appears in the skeleton form and can provide the frequency variations from one instance to the next. Figure 9 shows the marginal spectrum of the diastolic murmurs. It can be clearly seen that the energy mainly concentrates on the lower frequency domain.
Figure 10 shows the diastolic murmurs of the CAD patient, as diagnosed by coronary artery radiography. The left anterior descending artery is stenosed about 60% and the right coronary artery is stenosed about 85%. Figure 11 shows the IMFs of the murmur obtained by EMD. The diastolic cardiac cycle can be decomposed into six IMFs. The Hilbert spectrum is illustrated in figure 12. Figure 13 shows the marginal spectrum of diastolic murmurs. The HHT spectrum has superior temporal and frequency resolutions. The spectrums show precise time-frequency representation of signal. The energies spread over a much wider frequency domain. Much higher spectral energies are concentrated on high frequency compared with those of normal people. More energy distributes in the frequency band over 200Hz and a peak also lies around 350Hz, which often does not appear in diastolic murmurs of normal people. It can be explained as follows: for the CAD patient, the narrowed coronary arteries lead to the blood flow in coronary artery changing from laminar flow to turbulence flow, from simplicity to complexity. Coronary arterial stenosis gives rise to high frequencies of diastolic murmurs. The EMD method makes no assumption about the linearity or stationarity of the signal, and the IMFs are usually easy to interpret and relevant to the underlying dynamic processes being studied.
4. A new method for processing end effect in empirical mode decomposition
In the procedure of EMD, the cubic splines interpolation creates top and bottom envelopes that are implemented in the first step of the above sifting process. It is difficult to interpolate data near the beginnings or ends, where the cubic splines can have swings. The common method to deal with end effect is to consider the end points both as maximum and minimum locations with values unchanged, but this method will give a distorted view of the local mean near the boundaries. We propose a simpler method to restrict the end effect in spline interpolation . The key points are to determine the values and locations of extrema nearby end points. Suppose the length of data x is N, the steps can be implemented as follows:
Finding all the maxima and minima, and considering the end points both as maximum and minimum, that is, maximum= [l maximum N] and minimum= [1 minimum N].
The end points are still considered both as maximum and minimum, whereas their values can be adapted to and . Taking , as the mean of all maximum except for the first and last maximum (the subscript represents location of maximum). Similarly, taking , as the mean of all minimum except for the first and last minimum (the subscript represents location of minimum).
Comparing with x (1), with x (N), with x(1) and with x (N), respectively.
Using cubic splines interpolation to get top and bottom envelopes, and repeating the second step of above sifting process to extract IMF.
The performance of the proposed method is compared with the traditional method where the endpoints are considered both as maximum and minimum with values unchanged. As an example, we decompose a sinusoid signal by the sifting process. Figure 14 shows the signal.
Firstly, we consider the endpoints both as maximum and minimum with value unchanged. Figure 15 shows the top and bottom envelopes calculated by cubic splines interpolation in the sifting process. Top and bottom red dash dot line represent the envelopes. The sinusoid signal is decomposed into six IMFs and one residue by sifting process as depicted in figure 16.
Secondly, applying the proposed method above to restrict the end effect, figure 17 shows the top and bottom envelopes calculated by cubic splines interpolation in the sifting process. Red circles represent the end values predicted. The sinusoid signal is decomposed into one IMF and a residue by the sifting process as depicted in figure 18. The IMF is just the sinusoid and the value of the residue is smaller than 10-6. From figure 18, it can easily be seen that the swings appear near both ends and propagate inwards and produce superfluous IMFs. Actually, the sinusoid signal is an IMF itself in nature because it satisfies the IMF conditions which has the same numbers of zero-crossing and extrema, and can also be local symmetric. Therefore, the sifting process as represented by figure 18 should extract only one IMF. The results indicate that the method we proposed is effective.
5. Instantaneous frequency estimation of diastolic murmurs based on EMD and TEO
Diastolic murmurs can provide clinicians with valuable diagnostic and prognostic information about the function of heart valves. Quantitative analysis of instantaneous frequency (IF) of the murmurs can aid diagnosis [1, 13].
Instantaneous Frequency (IF) is an important signal characteristic, which characterizes the transients and fast changes in frequency as time progresses. The IF of diastolic murmur is used to describe the time-varying spectral contents of the characteristic frequency bands that are of interest for cardiovascular research. The IF of a signal is traditionally obtained by taking the first derivative of the phase of the signal with respect to time using the Hilbert transform. However, this definition is questionable and will mislead interpretation of instantaneous frequency, such as negative frequency. Instantaneous frequency can also be obtained from a time–frequency distribution (TFD) as the first conditional moment in the frequency, suggesting that the instantaneous frequency is the average frequency at each time, whereas the cross terms existing in TFD will lead to a very rapid degradation of performance and severely pollute the instantaneous frequency estimation .
TEO is a powerful nonlinear operator and has been successfully used in a number of applications including speech signal processing, image processing, etc. . TEO can track the modulation energy and estimate the instantaneous amplitude and frequency of AM-FM signals with the form
and are the instantaneous amplitude and frequency respectively.
In continuous time domain, TEO is defined by
For example, for a sinusoid signal , the TEO gives
For a monochromatic signal, the output by TEO is proportional to the squared product of frequency and amplitude.The TEO of the first order derivative of produce the output:
The two results above can be combined to estimate the frequency and amplitude of the signal as follows :
The estimate of instantaneous frequency and amplitude above are also suitable for AM, FM and AM-FM signals.
The discrete-time counterpart of TEO can be defined as:
A discrete-time real value AM-FM signal that is usually used to model time-varying amplitude and frequency patterns can be expressed as:
Where is the time-varying amplitude modulation, is the carrier frequency, is the maximum frequency deviation from the carrier frequency and , is the frequency deviation function and is the initial phase shift. The derivative of the phase , that is, the FM part of the signal is called as instantaneous frequency:
The instantaneous frequency and amplitude of the AM-FM modulated signal at any time instant can be respectively demodulated by applying the TEO to and its difference, which is called the Discrete Energy Separation Algorithm (DESA):
The estimates above are valid under the assumptions that the signal does not vary too fast nor too much compared to the carrier frequency. In general, the first demodulation algorithm (26) ~ (28) is called DESA-1 where ‘1’ implies the approximation of derivatives with a single sample difference. That is, the signal derivative is approximated by the average of forward and backward 1-point differences. The second demodulation algorithm (29) ~ (30) is called DESA-2 where ‘2’ implies a difference between samples whose time indices differ by 2. Both DESA-1 and DESA-2 algorithms yield very small errors and can give the accurate estimate of instantaneous frequency. The DESA-2 algorithm is less computationally complex and has an excellent, almost instantaneous, time resolution which can also lead to a simpler mathematical analysis. In this paper, we focus on the instantaneous frequency rather than the instantaneous amplitude by DESA-2.
Figure 19 shows an AM-FM signal where
The theoretic instantaneous frequency is shown in figure 20. The estimated instantaneous frequency by DESA-2 is shown in figure 21. The estimated amplitude envelope is also illustrated in figure 22. Note that there are no apparent discrepancies between the real values and the DESA-2 calculations. The errors are very slow but less smooth. The results indicate that DESA-2 can be used to track the instantaneous frequency and amplitude accurately.
Another mixture signal is composed of two linear swept-frequency signals shown in figure 23. The frequency of one chirp signal varies from 1Hz to 0.1 Hz and the other varies from 2 Hz to 0.1 Hz. The estimated IF is shown in figure 24. The two chirp signals are also better identified and localized except for near boundaries.
In this paper, we present a novel method to estimate the IF of diastolic murmurs using Empirical Mode Decomposition (EMD) and nonlinear the Teager Energy Operator (TEO). EMD has been analysed as in section 3 and can decompose diastolic murmurs into a series of Intrinsic Mode Functions (IMFs), then accurate IF estimation can be acquired by TEO.
The block diagram of the instantaneous frequency estimate based on EMD-TEO is shown in figure 25 (IF refers to the instantaneous frequency in the block diagram).
The instantaneous frequency of the original signal can be obtained in the following steps:
Decompose the original single into IMFs: j=1…M.
Calculate the instantaneous frequency of each IMF by DESA-2.
Calculate the average instantaneous frequency of the original signal:
It is the average frequency of mainly IMFs at each instant time.
Next we estimate the IF of diastolic murmurs from clinical coronary artery disease (CAD) patient based on the EMD-Teager method. The left anterior descending artery is stenosed about 40% and the right coronary artery is stenosed about 55%, which has already been diagnosed by catherization. Figure 26 shows the diastolic murmurs. Figure 27 shows the IMFs obtained by EMD. The diastolic murmurs can be decomposed into six IMFs and one residue. The amplitudes of IMF5 and IMF6 are smaller compared with the original signal. So IMF5 and IMF6 are abandoned. Figure 28 shows the IF of each effective IMF by DESA-2. Figure 29 shows the average IF of diastolic murmurs. Then some features such as mean, standard deviation, etc., can be extracted from the average IF. For the normal subject, figure 30 shows the IF of each effective IMF and figure 31 shows average IF of diastolic murmurs.
For the CAD object, we can see that both the IF of each IMF and average IF are higher than those for normal subject. The diastolic murmurs contain rich higher frequencies. The mean of average IF is 185Hz and the standard deviation is 40Hz. For the normal subject, the mean of average IF is 140Hz and the standard deviation is 26Hz. These can be explained as follows: for the CAD subject, the narrowed coronary arteries lead to the blood flow in coronary artery changing from laminar flow to turbulence flow, from simplicity to complexity. Coronary arterial stenosis gives rise to high frequencies of diastolic murmurs. The instantaneous frequency features effectively reveal the information as to whether the arteries are blocked and denote the frequency change of diastolic murmurs when the coronary arteries are occluded.
Diastolic murmurs contain the information of coronary artery occlusions which give the basis of CAD diagnosis. The Hilbert Huang Transform is an adaptive powerful method to analyse nonlinear and non-stationary time series. The important part of HHT is the Empirical Mode composition (EMD). In this paper, we firstly studied wavelet shrinkage denoising using the generalized threshold function and the data-driven SURE threshold value, which successfully removed noise from the PCG signal. Secondly, we obtained the Hilbert spectrum and marginal spectrum of diastolic murmurs for normal subjects and CAD patients after EMD. They provide higher resolution and energy concentration in the time-frequency plane. The Hilbert spectrum and marginal spectrum effectively reveal the information as to whether the arteries are blocked and provide a reliable indicator of CAD. For restricting the end effect of EMD, a simple, powerful and effective method is presented. The IF estimation algorithm is studied based on EMD-TEO. The results indicate that the IF of diastolic murmurs effectively reveal the information on whether the arteries are blocked and provide a reliable indicator of CAD and provides a reliable indicator of coronary artery disease.
This paper is partly supported by the National Natural Science Foundation of China (grant no. 61102133) and supported by the Key Project of the Science Technology Department of Zhejiang Province (grant no.2010C11065) and the project of Hangzhou Science and Technology Committee (grant no. 20110833B31).