Open access

Adaptive Analysis of Diastolic Murmurs for Coronary Artery Disease Based on Empirical Mode Decomposition

Written By

Zhidong Zhao, Yi Luo, Fangqin Ren, Li Zhang and Changchun Shi

Submitted: 30 May 2012 Published: 20 February 2013

DOI: 10.5772/55690

From the Edited Volume

Adaptive Filtering - Theories and Applications

Edited by Lino Garcia Morales

Chapter metrics overview

2,082 Chapter Downloads

View Full Metrics

1. Introduction

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.

Figure 1.

Heart sound signal

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.

Advertisement

2. Wavelet shrinkage method

We consider the following model of a discrete noisy signal:

x=θ+σzE1

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:

  1. Apply discrete wavelet transform to observe noisy signals.

  2. Estimate noise and threshold value, thresholding the wavelet coefficients of the observed signal.

  3. 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’.

δλH(x)={0|x|λx|x|>λE2
δλS(x)={0|x|λxλx>λx+λx<λE3

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:

δλm(x)=xλmxm1m=12E4
λ is threshold value.

When m is an even number:

δλm(x)=xxI(|x|λ)λmxm1I(|x|>λ)E5

When m is odd number:

δλm(x)=xxI(|x|λ)λmxm1I(|x|>λ)sign(x)E6

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.

Figure 2.

Generalized threshold function

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 l2 risk for the generalized threshold function.

Let x~N(θ,1)

Am(θ)=λϕ(xθ)ϕ(x+θ)xmdx Bm(θ)=λϕ(xθ)+ϕ(x+θ)xmdx

ϕ And Φ are density and probability function of standard Gaussian random variable respectively. Then:

Mean:

Mm(λ,θ)=MH(λ,θ)λmAm1(θ)E7

Bias:

SBm(λ,θ)=(Mm(λ,θ)θ)2 E8

Variance:

Vm(λ,θ)=VH(λ,θ)2λmBm2(θ)λ2mAm12(θ)+λ2mB2m2(θ)+2λmMH(λ,θ)Am1(θ)    E9

l2 Risk:

l2E10

Where

ρλm(θ)=E(δλm(x)θ)2=ρλH(θ)2λmBm2(θ)+λ2mB2m2(θ)+2θλmAm1(θ)

MH(λ,θ)=θ+θ[1Φ(λθ)Φ(λ+θ)]+ϕ(λθ)ϕ(λ+θ)VH(λ,θ)=(θ2+1)(2Φ(λθ)Φ(λ+θ)]+(λ+θ)ϕ(λθ)+(λθ)ϕ(λ+θ)MH(λ,θ)2ρλH(θ)=1+(θ21)(Φ(λθ)Φ(λθ))+(λ+θ)ϕ(λ+θ)+(λθ)ϕ(λθ), Mm(λ,θ) ,SBm(λ,θ), Vm(λ,θ) are the mean, bias, variance and risk of generalized threshold function. When m is 1, 2, ρλm(θ), 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) [10] 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, θk=xk+H(xk),k=1...N,H is true risk. Then
ρ(xk,λ)=Eθkθk2E11

is the unbiased risk estimate of ρ(xk,λ)=1+2(ddxkH(xk))+H2(xk):

Proof:

ρ(xk,λ)Eθkθk2=E(xk+H(xk)θk)2=E(zk+H(xk))2=1+2E(zkH(xk))+E(H2(xk))

Where =1+2E(zkH(θk+zk))+E(H2(xk)) and by partial integral

zk=xkθk

Then

E(zkH(θk+zk))=12πzkH(θk+zk)eξ22dξ=12π(ηkθk)H(ηk)exp((ηkθk)22)dηk=12πexp((ηkθk)22)dH(ηk)dηkdηk=E(dH(ηk)dηk|ηk=xk)

So

ρ(xk,λ)=Eθkθk2=1+2E(dH(xk)dxk)+E(H2(xk)) is the unbiased risk estimate of true risk E[ρ(xk,λ)]=1+2(dH(xk)dxk)+H2(xk).

For the generalized threshold function (5) and single coefficient, when m is even,

ρ(xk,λ)

The SURE is

H(xk)=xkI(|xk|λ)λmxkm1I(|xk|>λ)E12

When m is odd,

SURE(xk,λ)=1+(xk22)I(|xk|λ)+(2(m1)λmxkm+λ2mxk2m2)I(|xk|>λ)

The SURE is

H(xk)=xkI(|xk|λ)λmxkm1I(|xk|>λ)sign(xk)E13

Suppose wavelet coefficients are SURE(xk,λ)=1+(xk22)I(|xk|λ)+λ2mxk2m2I(|xk|>λ)+2(m1)λmxkmI(xk>λ)2(m1)λmxkmI(xk<λ), the threshold value λ is set to minimize the estimate of the x1....xN risk for the given data,

l2E14

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.

Figure 3.

a) Noisy PCG signal (b) Phase space diagram of the noisy signal

Figure 4.

a) PCG signal after denoising (b) Phase space diagram of denoised signal

Advertisement

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 [6], 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:

  1. In the whole dataset, the number of extrema and the number of zero-crossing must either be equal or differ at most by one.

  2. 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:

λ=argminλ0k=1NSURE(xk,λ)E15

with time varying amplitude envelope x(t)=a(t)cosϕ(t) and phase a(t). 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:

  1. 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 ϕ(t) and lower envelope xup(t).

  2. Subtracting mean of these two envelopes xlow(t) from the signal to get their difference: m1(t)=(xup(t)+xlow(t))/2.

  3. Regarding the h1(t)=x(t)m1(t) as the new data and repeating steps 1 and 2 until the resulting signal meets the two criteria of an IMF, defined as h1(t). The first IMF c1(t) contains the highest frequency component of the signal. The residual signal c1(t) is given by r1(t).

  4. Regarding r1(t)=x(t)c1(t) as new data and repeating steps (1) (2) (3) until extracting all the IMFs. The sifting procedure is terminated until the m-th residue r1(t) becomes less than a predetermined small number or becomes monotonic.

The original signal x (t) can thus be expressed as follows:

rM(t)E16
x(t)=j=1Mcj(t)+rM(t) is an IMF where j represents the number of corresponding IMFs and cj(t) is residue. The EMD decomposes non-stationary signals into narrow-band components with decreasing frequency. The decomposition is complete, almost orthogonal, local and adaptive. All IMFs form a completely and nearly orthogonal basis for the original signal. The basis comes directly from the signal, which guarantees the inherent characteristic of signal and avoids the diffusion and leakage of signal energy. The sifting process eliminates riding waves, so each IMF is more symmetrical and is actually a zero mean AM-FM component.

Figure 5.

A classical IMF

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 6.

Diastolic murmurs of a normal object

Figure 7.

IMFs of diastolic murmurs from the normal people

Figure 8.

Hilbert spectrum of the diastolic murmurs

Figure 9.

Marginal spectrum of the diastolic murmurs

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.

Figure 10.

Diastolic Murmurs of CAD patient

Figure 11.

Six IMFs of diastolic murmurs from patient

Figure 12.

Hilbert spectrum of the diastolic murmurs from patient

Figure 13.

Marginal spectrum of the diastolic murmurs from patient

Advertisement

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 [11]. 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:

  1. 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].

  2. The end points are still considered both as maximum and minimum, whereas their values can be adapted to rM(t) and δ1,γ1. Taking δN,γN, δ1 as the mean of all maximum except for the first and last maximum (the subscript represents location of maximum). Similarly, taking δN, γ1 as the mean of all minimum except for the first and last minimum (the subscript represents location of minimum).

  3. Comparing γN with x (1), δ1 with x (N), δN with x(1) and γ1 with x (N), respectively.

If γN<x(1) then δ1= x (1);if δ1< x (N) then δN= x (N); if δN> x (1)then γ1=x (1);If γ1>x(N) then γN= x(N).
  1. 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.

Figure 14.

A sinusoid signal

Figure 15.

Cubic splines interpolation in sifting process using the traditional method

Figure 16.

IMFs of the sinusoid 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.

Figure 17.

Cubic splines interpolation in sifting process using the proposed method

Figure 18.

IMF and residue of the sinusoid signal using the proposed method

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.

Advertisement

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 [14].

TEO is a powerful nonlinear operator and has been successfully used in a number of applications including speech signal processing, image processing, etc. [15]. TEO can track the modulation energy and estimate the instantaneous amplitude and frequency of AM-FM signals with the form

γNE17

x(t)=a(t)cos[2π0tω(τ)dτ] and a(t) are the instantaneous amplitude and frequency respectively.

In continuous time domain, TEO is defined by

ω(τ)E18
Ψ(x(t))=[x˙(t)]2x(t)x¨(t) corresponds to continuous signal, x(t) and x˙(t) are the first order and second order time derivatives of x¨(t) respectively.

For example, for a sinusoid signal x(t), the TEO gives

x(t)=Acos(ωt)E19

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 Ψ(x(t))=A2ω2 of x˙(t) produce the output:

x(t)E20

The two results above can be combined to estimate the frequency and amplitude of the signal Ψ(x˙(t))=A2ω4 as follows [14]:

x(t)E21
ω^2(t)=Ψ(x˙(t))Ψ(x(t))E22

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^2(t)|=Ψ2(x(t))Ψ(x˙(t))E23

A discrete-time real value AM-FM signal that is usually used to model time-varying amplitude and frequency patterns can be expressed as:

Ψ(x(n))=x2(n)x(n1)x(n+1)E24

Where x(n)=a(n)cos(ϕ(n))=a(n)cos(ωcn+ωm0nq(k)dk+θ) is the time-varying amplitude modulation, a(n) is the carrier frequency, ωc is the maximum frequency deviation from the carrier frequency and ωm, 0<ωm<ωc is the frequency deviation function and |q(n)|1 is the initial phase shift. The derivative of the phase θ, that is, the FM part of the signal is called as instantaneous frequency:

ϕ(n)E25

The instantaneous frequency ω(n)=dϕ(n)dn=ωc+ωmq(n) and amplitude ω(n) of the AM-FM modulated signal a(n) at any time instant can be respectively demodulated by applying the TEO to x(n) and its difference, which is called the Discrete Energy Separation Algorithm (DESA):

x(n)E26
y(n)=x(n)x(n1)E27
ω(n)=arccos(1Ψ{y(n)}+Ψ{y(n+1)}4Ψ{x(n)})E28

or

|a(n)|=Ψ{x(n)}sin2(ω(n))E29
ω(n)=12arccos(1Ψ{x(n+1)x(n1)}2Ψ{x(n)})E30

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 |a(n)|=2Ψ{x(n)}Ψ{x(n+1)x(n1)} where

x(n)=a(n)cos(ϕ(n))E31

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.

Figure 19.

Original AM-FM signal

Figure 20.

Theoretic instantaneous frequency

Figure 21.

Estimated instantaneous frequency by DESA-2

Figure 22.

Estimated amplitude envelope by DESA-2

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.

Figure 23.

A mixture signal of two chirp signals

Figure 24.

Estimated IF of two IMFs by DESA-2

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.

Figure 25.

Block diagram of Instantaneous Frequency (IF) estimate based on EMD-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:

  1. Decompose the original single into IMFs: a(n)=1+0.6cos(0.01πn)ϕ(n)=π10n+cosπ80n j=1…M.

  2. Calculate the instantaneous frequency cj(t) of each IMF IFj(t) by DESA-2.

  3. Calculate the average instantaneous frequency of the original signal:

cj(t)E32

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.

Figure 26.

Diastolic Murmurs of CAD object

Figure 27.

Six IMFs and one residue by EMD

Figure 28.

Estimated IF of four selective IMFs by DESA-2

Figure 29.

The average instantaneous frequency of diastolic murmurs

Figure 30.

Estimated instantaneous frequency of normal object by DESA-2 algorithm

Figure 31.

Estimated IF of normal object

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.

Advertisement

6. Conclusion

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.

Advertisement

Acknowledgments

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).

References

  1. 1. AkayMet alHarmonic decomposition of diastolic heart sounds associated with coronary artery disease. Signal Processing 19954117990
  2. 2. AkayMet alComparative study of advanced signal processing techniques for detection Coronary Artery Disease. Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society 199121392140
  3. 3. DjebbariABereksi Reguig F. Short-time Fourier transform analysis of the phonocardiogram signal. The 7th IEEE International Conference on Electronics, Circuits and Systems 2002844847
  4. 4. DebbalS. MBereksi Reguig F. Time-frequency analysis of the first and the second heartbeat sounds, Applied Mathematics and Computation 2007184210411052
  5. 5. KhadrLMatalgahMet alThe wavelet transform and its applications to phonocardiogram signal analysis, Medical informatics 1991163221227
  6. 6. HuangN. Eet alThe empirical mode composition and the Hilbert spectrum for non linear and non-stationary time series analysis, Proceedings of the Royal Society of London An1998
  7. 7. ChengJet alResearch on the intrinsic mode function (IMF) criterion in EMD method, Mechanical Systems and Signal Processing 2006204817824
  8. 8. GaoH. YBruceA. GUnderstanding waveshrink: variance and bias estimation, Biometrika 1996834727745
  9. 9. GaoH. YWavelet shrinkage denoising using the non-negative garrote, Journal of Computational and Graphical Statistics 199874469488
  10. 10. SteinCEstimation of the mean of a multivariate normal distribution, Annuals of statics 19819611351151
  11. 11. ZhaoZWangYA new method for processing end effect in Empirical Mode Decomposition, Communications,International Conference on Communications, Circuits, and Systems 2007841845
  12. 12. GauthierDet alSpectral Analysis of Heart Sounds Associated With Coronary Occlusions 6th International Special Topic Conference on Information Technology Applications in Biomedicine 20074952
  13. 13. OliveiraP. MBarrosoVDefinitions of Instantaneous Frequency under physical constraints. Journal of the Franklin Institute 2000
  14. 14. MaragosPet alOn separating amplitude from frequency modulations using energy operators, ICASSP 1992
  15. 15. ZhaoZ. DZhaoZ. Jet alTime-frequency analysis of heart sound based on HHT. International Conference on Communications, Circuits and Systems 2005
  16. 16. ZhaoZ. DPanMInstantaneous Frequency Estimation of Diastolic Murmurs Based on EMD and TEO. 1st International Conference on Bioinformatics and Biomedical Engineering, 2007829832
  17. 17. ZhaoZ. DWavelet shrinkage denoising by generalized threshold function, International Conference on Machine Learning and Cybernetics 200555015506
  18. 18. YoshidaHShinoHYanaKInstantaneous frequency analysis of systolic murmur for phonocardiogram,19th Annual International Conference of the IEEE Engineering in Medicine and Biology Society 1997

Written By

Zhidong Zhao, Yi Luo, Fangqin Ren, Li Zhang and Changchun Shi

Submitted: 30 May 2012 Published: 20 February 2013