Quasi-periodic modeling for heart sound localization and suppression in lung sounds

Biomedical Signal Processing and Control 8 (2013) 586–595

Contents lists available at SciVerse ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Quasi-periodic modeling for heart sound localization and suppression in lung sounds Miroslav Zivanovic ∗ , Miriam González-Izal Dept. Ingeniería Eléctrica y Electrónica, Universidad Pública de Navarra, Campus Arrosadía, 31006 Pamplona, Spain

a r t i c l e

i n f o

Article history: Received 10 December 2012 Received in revised form 6 May 2013 Accepted 5 June 2013 Keywords: Heart sound localization Lung sound denoising Quasi-periodic modeling

a b s t r a c t We present a novel approach to treating lung sounds contaminated by heart sounds by means of quasiperiodic signal modeling. Heart sounds are described through the long and short-term quasi-periodicity generated by cardiac cycle repetition and oscillations caused by valves openings and closures, respectively. In terms of signals, these quasi-periodicities drive time-variant HS envelope and phase which are modeled by single and piecewise time polynomials, respectively. Single polynomials account for slow and continuous envelope time variations, while piecewise polynomials capture fast and abrupt phase changes in short time intervals. Such a compact signal description provides an efficient way to fundamental heart sound (FHS) components localization and posterior removal from lung sounds. The results show that the proposed method outperforms two reference methods for medium (15 ml/s/kg) and high (22.5 ml/s/kg) air flow rates. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction Lung sound (LS) analysis is a key issue for accurate medical diagnosis as it provides valuable information about lungs psychologies and pathologies. The major general problem, however, remains the presence of heart sounds (HS) which interferes with LS in 20–150 Hz approximately [1–3]. Incorrectly localized and removed fundamental heart sound components (FHS) can usually render LS unfit for diagnosis purposes. Denoising of LS is often viewed as a two-step process, namely, (1) HS localization and (2) suppression. In the literature we find a number of approaches to HS localization based on time–frequency filtering [4], higher-order statistics [5], wavelet decomposition [6,7], sound entropy [8] and singular spectrum analysis [9]. Additional efforts have been made to suppress FHS components from LS with the aim of preserving most of the LS signal content in the HS frequency range. Approaches like independent component analysis [10], adaptive filtering [11–13], linear and non-linear prediction [14,15] and blind source separation [16] were dealt with. A thorough analysis would be needed to evaluate the performance of the aforementioned methods and bring out their virtues and shortcomings. However, a general remark regarding limitations

∗ Corresponding author. Tel.: +34 948169024; fax: +34 948169720. E-mail addresses: [email protected] (M. Zivanovic), [email protected] (M. González-Izal). 1746-8094/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.bspc.2013.06.003

usually stems from either an a priori assumption on HS stationarity or a lack of robustness in conditions of middle and high air flow rates. Herein we present a novel approach to HS characterization by explicit modeling of the long and short-term quasi-periodicities, which arise from cardiac cycle repetition and vibrations caused by openings and closures of the valves, respectively. In the signal space the former account for the instantaneous envelope while the latter address the issue of the time-variant phase. Under the assumption of slow and continuous periodicity variation, the envelope is modeled by a harmonic structure whose coefficients are low degree single time polynomials. Phase changes, on the contrary, are fast and abrupt and occur in short time intervals corresponding to FHS components; accordingly, the phase is modeled as a single sinusoid scaled by uniformly distributed piecewise polynomials. Such a modeling allows for a robust HS characterization in middle and high air flow rates scenarios. Moreover, the method is computationally highly efficient because the model parameters are estimated from an overdetermined linear system of equations. Note that single and piecewise polynomials have already been utilized in biomedical signal processing e.g. [17,18] but to the best of our knowledge never for the present application. The present paper is organized as follows: Section 2 is dedicated to HS modeling in terms of time-variant envelope and phase; in Section 3 we present experimental results in the context of a twofold comparative study; finally, the conclusions appear in Section 4.

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

587

The denominator in (5) stands for the fact that f(t) is the derivative of the instantaneous phase. By inserting (4) and (5) into (3) and after some elementary trigonometric manipulations (omitted for the sake of simplicity) A(t) becomes a linear combination of sine and cosine product terms of the non-linear argument 2nfm tm+1 /(m + 1), m = 1, . . ., M. Assuming slow time variations in an (t) and f(t), the small argument approximation can be applied:



sin

 cos

2nfm t m+1 m+1 2nfm t m+1 m+1



≈

2nfm t m+1 , m+1

(6)

 ≈ 1.

(7)

The combination of (3)–(7) yields the final model for the HS signal envelope: A(t) =

N 

sn (t) sin(2nf0 t) + cn (t) cos(2nf0 t),

(8)

n=1

Fig. 1. Illustration of the quasi-periodicities: (a) a HS signal with the long-term quasi-period 1/f(t) and (b) a FHS component with the short-term quasi-period 1/F(t).



2M+1

sn (t) =

 (n) ai

i−2 

− 2n



2M+1

cn (t) =

2.1. HS signal model

 (n) bi

i−2 

+ 2n

In general terms we can represent a HS as a narrow-band signal h(t) with a time-variant envelope and phase: h (t) = A(t) cos ϕ(t).

(1)

A(t) characterizes the instantaneous RMS amplitude of the h(t) and it exhibits prominent peaks localized at the FHS components. The repetitive nature of the cardiac cycle induces quasi-periodic patterns in A(t) with the quasi-period of the order of 1 s (longterm quasi-periodicity). The phase ϕ(t) describes the time-variant frequency content of the FHS components with the quasi-period typically in the range 5–50 ms (short-term quasi-periodicity). For the sake of clarity both quasi-periodic trends are illustrated in Fig. 1. 2.2. HS envelope modeling Due to its quasi-periodic nature A(t) can be modeled by a set of N time-variant harmonics: A(t) =

N 

(2)

The last expression can be rewritten in the following way: N 

− (an (t) sin n ) sin(2nf (t)t)

n=1

+ (an (t) cos n ) cos(2nf (t)t).

(3)

The time-variant parameters in (3) are modeled by polynomials of order M under the assumption of continuous time variation: −an (t) sin n =

M 

(n) am t m ,

m=0

f (t) =

M  fm m=0

i−1

ti,

(9)

ti.

(10)



sn (t) and cn (t) are polynomials of degree 2M + 1 whose internal structure reveals the way the harmonic amplitudes are related to the frequency. This structure is transparent to the envelope modeling process; nevertheless, it can provide valuable information about cardiac cycle time evolution, which is related to the issue of heart rate variability (out of the scope of the present paper). From (8) to (10) we observe that if A(t) was periodic, then sn (t) = a0 (n) , cn (t) = b0 (n) , f(t) = f0 and (8) would become a set of harmonics of the fundamental cardiac period 1/f0 describing a stationary HS. Otherwise, long-term quasi-periodicity is captured by means of the non-zero polynomial coefficients. Accordingly, (8) correctly characterizes the envelope of a HS signal allowing variations in energy and time distance between the FHS components. Moreover, A(t) conveys necessary information for LS–HS localization, which can be readily obtained e.g. by adaptive thresholding [9]. 2.3. HS phase modeling

an (t) cos[2nf (t) + n ].

n=1

A(t) =

(n) f am i−m−1

m=0

i=0



i−1

m=0

i=0

2. Methods

(n) f bm i−m−1

m+1

an (t) cos n =

M 

The term cos ϕ(t) in (1) characterizes the oscillatory behavior of a HS signal in short time intervals (typically 40–100 ms) centered at the FHS components. The character of these oscillations has been studied in the past [19–21], where it was found that they resemble chirp signals. Accordingly, we model ϕ(t) by means of a linear trend instantaneous frequency F(t): ϕ(t) = 2F(t)t + ˛,

F(t) = F0 +

(4)

m=0

(11)

By inserting (11) into (1) and assuming A(t) is known we obtain:





h(t) = A(t) cos 2 F0 + (n) bm t m ,

F1 t. 2





F1 t t+˛ . 2

(12)

After trigonometrically decomposing (12) and applying the approximations (6) and (7) we can write: h(t) = (− sin ˛ − F1 t 2 cos ˛)A(t) sin 2F0 t

tm.

(5)

+ (cos ˛ + F1 t 2 sin ˛)A(t) cos 2F0 t.

(13)

588

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

The last expression shows that the short-term quasi-periodic behavior of the FHS components can be described by two parabolas in a short time interval. Observe that for F1 = 0 the oscillatory character of the FHS components becomes periodic with the period 1/F0 . It turns out, however, that (13) needs additional information about the FHS components localization for the sake of implementation (recall that (13) is limited to short time intervals between 40 ms and 100 ms approximately). Accordingly, we would first have to localize and delimit the FHS components by examining A(t) and then estimate the model parameters for each FHS in the analysis window (compare this approach to (8)–(10) which model A(t) in the whole analysis window). Moreover, the computational time would strongly depend on the number of FHS components detected in the analysis window, rendering the approach unfit for the present application. An adequate way to overcome the aforementioned shortcomings is by using piecewise polynomial approximation in the whole analysis window. Piecewise polynomials, among which uniform splines are the most popular, have very convenient properties out of which we mention a few: (1) an accurate modeling in a short time interval with a low-degree polynomial, (2) continuity of derivatives of a certain order in the data points, (3) efficient calculation due to uniform spacing [22–24]. For the present application we have chosen quadratic splines, which are represented as a linear combination of uniformly time-distributed third-order B-splines. Such a spline forms an approximation basis which perfectly fits to the short-time quadratic polynomials in (13). Each third-order B-spline is constructed over three adjacent uniformly spaced time knots as follows:

Br,2 (t) =

⎧ ⎪ 0.5(t − tr )2 , ⎪ ⎨ ⎪ ⎪ ⎩

t ∈ [tr+1 , tr+2 ]

0.5(1 − (t − tr+2 ))2 ,

t ∈ [tr+2 , tr+3 ],

(14)

where Br,2 (t) is a third-order B-spline, r = 1, . . ., R is a knot index and R is a number of knots in the analysis window. Accordingly, (13) can be rewritten as: h(t) =

R 

A = (S(1) S(2) . . . S(N) C(1) C(2) . . . C(N) ) × (s1 s2 . . . sN c1 c2 . . . cN )T , (16) (n)

lmod(2M+1)

sin(2nf0 tk ),

(17)

(n)

lmod(2M+1)

cos(2nf0 tk ),

(18)

Sk,l = tk

Ck,l = tk

for l = 0, . . ., 2N(M + 2), n = 1, . . ., N and k = 0, . . ., K − 1 where K is the number of data samples and T is the transpose operator. The vector A contains the output of the 10 Hz cut-off low-pass Butterworth filter while the subvectors s1 , . . ., cN contain the polynomial coefficients. The splines coefficients  r and r are estimated from the following system written in the matrix form: T

HS = (1 2 . . . R 1 2 . . . R ) × (␰) ,

(19)

k,r = Br,2 (tk )A(tk ) sin 2F0 tk

(20)

k,r = Br,2 (tk )A(tk ) cos 2F0 tk

(21)

where r = 1, . . ., R and the vector HS contains the input data. Both overdetermined linear systems (16) and (19) can be represented as a linear least-squares problem which can easily be solved using, for instance, the QR factorization or singular value decomposition. 3. Results

t ∈ [tr , tr+1 ]

−(t − tr+1 )2 + (t − tr+1 ) + 0.5,

cut-off band-pass for F0 estimation. In both cases a third-order Butterworth filter was used. The rest of the parameters are estimated from two linear sets of Eqs. (8) and (15), respectively. The coefficients in sn (t) and cn (t) are estimated from the following matrix system:

In order to best assess the validity of the proposed modeling method under different analysis circumstances, we have carried out a thorough comparative study in terms of HS localization as well as LS denoising by means of both artificial and real data mixtures. In the following sections we describe in detail all the issues related to the experiments and discuss the results. 3.1. Reference methods

r Br,2 (t)A(t) sin 2F0 t +

r Br,2 (t)A(t) cos

2F0 t,

(15)

r=1

where the parameters  r and r are the B-spline coefficients. Regarding the last expression we give a few remarks. First, observe that (15) provides not only the phase model but also the full model for HS signals. Accordingly, HS removal from LS is carried out by a simple h(t) subtraction from the original LS signal. Second, (15) can be useful in applications where HS recordings are corrupted by external noises (skin contact, muscular activity, surrounding speech) e.g. [25]. Third, the choice of the uniform time-distribution of Br,2 (t) has been motivated by its low computational complexity, although some studies show that non-uniform distribution can achieve more optimal signal approximations in terms of the residual norm [26–28]. 2.4. Parameter estimation In order to fully characterize a HS signal the following set of parameters must be estimated: the coefficients of the polynomials sn (t) and cn (t) in (8), the scalars  r and r in (15) and the frequencies f0 and F0 . The frequencies are straightforwardly estimated by the classical auto-correlation based approach [29] after the HS has been adequately filtered: 10 Hz cut-off low-pass for f0 and 20–150 Hz

Regarding HS localization in LS signals we have compared the proposed method to [9]. This method relies on a singular spectrum analysis approach (SSA) which aims at decomposing HS–LS mixtures onto two signal subspaces by means of detecting different trends in the eigenspectra. As shown in [9], the SSA approach seems to be more suitable for HS localization than the multiresolution decomposition and entropy-based approach. As for LS denoising we have chosen [14] as a reference method. The method performs HS localization and LS denoising through multiscale product (MP) and linear prediction (LP), respectively (in the present study we will refer to this method as the MP–LP). First, it completely removes segments which contain HS, thus creating gaps in the wavelet decomposition of the LS signal; next, the gaps are “filled” through LP under the assumption of a non-stationary LS. The results show that the MP–LP yields a fairly good HS suppression in terms of power spectral density for low and medium flow rates. 3.2. Signals For the present study we have used two sets of test signals. The first set comprises real HS and LS recordings from [30,31] which we have artificially mixed in order to simulate different analysis scenarios. By adjusting the Signal-to-noise ratio (SNR) between the LS

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

1 0.8

Normalized amplitude

and HS we could simulate different levels of flow rates. In particular, medium flow rates were simulated by the SNR in 5–7 dB while for high flow rates the SNR was in the range 7–9 dB. Moreover, the knowledge of the underlying HS signal provided a quantitative study in terms of a set of parameters enumerated in the following subsection. An illustration of these test signals is given in Fig. 2 where the following time-waveforms are plotted: (a) heart signal, (b) lung signal and (c) their mixture for a medium flow rate. The second set of test signals included LS recorded from four healthy subjects aged between 30 and 45. The signals were captured by a piezoelectric accelerometer Murata which was placed between the second and third intercostal spaces at the left sternal border. Five recordings (approximately 30 s each) were made for all subjects, which were instructed to breathe at medium and high flow rates. The signals were digitized at 44.1 kHz with 16 bits precision and later downsampled to fs = 4410 Hz.

0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0

1

2

3

3.3. Performance evaluation parameters

0.8 0.6

t =

t h2 (t ) k=q k FHS k , p h2 (t ) k=q FHS k



P

t hˆ 2 (t ) k=Q k FHS k

P

hˆ 2 (t ) k=Q FHS k

6

7

8

9

,

(23)

6

7

8

9

6

7

8

9

1

Normalized amplitude

(22)

p

tˆ =

5

(a)

• Duration error[samples] = fs−1 |(tp − tP ) − (tq − tQ ). • Position error[samples] =

4

Time [s]

3.3.1. HS localization Let us denote an original FHS signal component as hFHS (t), tq ≤ t ≤ tp and the corresponding component detected from a mixˆ FHS (t), tQ ≤ t ≤ tP . The quality of localization for a given FHS ture as h is measured by means of the following parameters:

fs−1 |t − tˆ|,

589

0.4 0.2 0 -0.2 -0.4 -0.6

where the operator “<·>” stands for the mean time or center of gravity of the detected component [32]. The center of gravity is more accurate than a simple geometric center because it points at the maximum energy time-concentration of the underlying FHS component.

-0.8 -1 0

1

2

3

4

5

Time [s]

⎛

q

h2 (t ) k=p FHS k

SDRdB = 10 log10 ⎝



2 q ˆ  (hFHS (tk ) − hFHS (tk )) k=p

(b)

⎞ ⎠

(24) 1

which measures time-waveform similarity between the original and estimated FHS component. • False positive and false negative errors [%] which evaluate the methods’ performance as a FHS component detector. 3.3.2. HS suppression Denoting the PSD of a clean LS as Y(fu ) and the PSD of a denoised ˆ (fu ), we measure the quality of HS removal from a LS signal by LS as Y the Power spectrum density difference (PSDD) in dB evaluated in a set of logarithmically spaced frequency bands. For a band delimited by the frequency indices [u1 , u2 ] the PSDD is calculated as follows:

 u

u=u1

Y (fu )

u=u1

Yˆ (fu )

2

PSDDdB = 10 log10

u2



0.8

Normalized amplitude

• Signal-to-distortion ratio (SDR) [dB]:

0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1

(25)

3.4. HS localization – comparative study The first part of the study was performed on the artificially mixed signals by means of the duration error, position error and

0

1

2

3

4

5

Time [s]

(c) Fig. 2. Signals used for the HS localization comparative study: (a) heart signal, (b) lung signal and (c) their mixture for a medium flow rate.

590

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

MEDIUM FLOW

HIGH FLOW SSA modeling

0.8

0.6

0.4

0.2

0

SSA modeling

1

Norm. histogram

Norm. histogram

1

0.8

0.6

0.4

0.2

0

5

10

15

20

25

30

35

0

40

0

10

20

Position error [samples]

SSA modeling

SSA modeling

1

0.8

Norm. histogram

Norm. histogram

60

HIGH FLOW

1

0.6

0.4

0.2

0.8

0.6

0.4

0.2

0

20

40

60

80

0

100

0

50

Duration error [samples]

100

150

200

Duration error [samples]

(b)

(e)

MEDIUM FLOW

HIGH FLOW SSA modeling

1

SSA modeling

1

0.8

Norm. histogram

Norm. histogram

50

(d)

MEDIUM FLOW

0.6

0.4

0.2

0

40

Position error [samples]

(a)

0

30

0.8

0.6

0.4

0.2

2

4

6

8

10

SDR [dB]

(c)

12

14

16

18

0

0

2

4

6

8

10

12

14

16

18

SDR [dB]

(f)

Fig. 3. Position error, duration error and SDR histograms obtained for the SSA and proposed modeling method: (a–c) correspond to medium flow rates; (d–f) correspond to high flow rates.

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

SIGNAL

591

SIGNAL SEGMENT

2.5

2

2

1.5

1.5 1

Amplitude

Amplitude

1 0.5 0 -0.5

0.5 0 -0.5

-1 -1

-1.5 -1.5

-2 -2.5

0

1

2

3

4

5

6

7

8

-2

9

5

5.5

Time [s]

(a)

SSA 1

2

0.8

1.5

0.6

1

0.4

Amplitude

Amplitude

LS ESTIMATED

0.5 0 -0.5

0.2 0 -0.2

-1

-0.4

-1.5

-0.6

-2

-0.8

0

1

2

3

4

5

6

7

8

-1

9

5

5.5

Time [s]

MODELING 1

0.8

0.8

0.6

0.6

0.4

0.4

Amplitude

Amplitude

HS ESTIMATED

0.2 0 -0.2

0.2 0 -0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

1

2

3

4

5

Time [s]

(c)

6.5

(e)

1

0

6

Time [s]

(b)

-1

6.5

(d)

2.5

-2.5

6

Time [s]

6

7

8

9

-1

5

5.5

6

6.5

Time [s]

(f)

Fig. 4. High flow HS–LS separation illustration: (a) mixture signal, (b and c) estimated LS and HS by the modeling method, (d) original signal segment, (e) the segment processed by the SSA, and (f) the segment processed by the modeling method.

592

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

Table 1 Position error, duration error and SDR for the artificial mixtures. Position error (samples) SSA (low flow) Modeling (low flow) SSA (medium flow) Modeling (medium flow) SSA (high flow) Modeling (high flow)

± ± ± ± ± ±

4.88 2.49 15.07 4.20 23.33 6.98

Duration error (samples)

3.06 1.91 8.72 3.84 13.33 5.57

6.21 5.98 21.11 6.33 68.55 6.66

± ± ± ± ± ±

2.66 2.80 34.44 7.58 32.62 8.73

Table 4 Mean execution time (s) as a function of the number of samples on an Intel 3.19 GHz processor.

SDR (dB) 7.93 13.22 5.49 11.14 5.90 8.34

± ± ± ± ± ±

1.80 2.71 2.13 2.23 1.72 2.22

Table 2 False positives and negatives (%) for the artificial mixtures. False positives (%) SSA (low flow) Modeling (low flow) SSA (medium flow) Modeling (medium flow) SSA (high flow) Modeling (high flow)

0.86 0.91 2.84 1.58 52.53 3.21

± ± ± ± ± ±

0.30 0.53 1.72 1.27 22.28 2.78

False negatives (%) 0.70 0.62 1.84 1.40 15.84 2.67

± ± ± ± ± ±

0.23 0.41 1.05 1.52 8.36 1.90

Number of samples

10,000

20,000

30,000

40,000

SSA Modeling

2.05 0.27

3.54 0.53

4.79 0.65

6.23 0.82

percentage errors. For each subject and flow rate these parameters were calculated and the results are shown in Tables 2 and 3. We can see that in general for medium flow rates both methods exhibit statistically large scores (small false errors percentages). However, for high flow rates the SSA undergoes an important degradation due to a stronger impact of the LS onto HS localization. The reason behind this effect is the principal component (PC) pairs selection criterion which interprets the smallest eigenvalues as noise and largest as signal. It turns out that this criterion does not hold for high flow rates, as the noise starts generating larger eigenvalues

MEDIUM FLOW

Table 3 False positives and negatives (%) for the real mixtures from four subjects.

False positives (%) SSA (low flow) Modeling (low flow) SSA (medium flow) Modeling (medium flow) SSA (high flow) Modeling (high flow) False negatives (%) SSA (low flow) Modeling (low flow) SSA (medium flow) Modeling (medium flow) SSA (high flow) Modeling (high flow)

Subject 1

Subject 2

Subject 3

Subject 4

0.87 0.75 2.16 1.27 62.76 3.82

0.95 0.82 1.86 1.76 58.93 2.20

1.13 1.04 3.06 0.91 46.80 1.68

0.79 0.80 1.88 1.31 55.48 2.09

0.68 0.54 1.98 1.46 19.74 3.05

0.61 0.48 1.86 1.36 14.17 1.28

0.93 0.73 2.88 0.72 13.29 2.81

0.55 0.68 1.41 1.50 22.10 4.72

4

Modeling MP-LP

3.5

Average PSDD [dB]

3 2.5 2 1.5 1 0.5 0

1

2

3

4

Frequency band

(a) HIGH FLOW 6

Modeling MP-LP 5

Average PSDD [dB]

SDR. A total of 100 10-s mixtures were generated representing different medium and high-flow rate scenarios. These signals were processed by the SSA and proposed modeling method and the outcomes are represented as normalized histograms. For the modeling method we used M = 1, R = 20 (expressions (4) and (15)) while the number of harmonics in A(t) was chosen according to the following criterion. The lower limit on FHS duration is assumed 40 ms: hence, the envelope bandwidth spans up to 1/40 ms = 25 Hz approximately. By taking the mean cardiac period as 1 s, we conclude that A(t) can be properly modeled with N = 25 harmonics. In addition, we used the same adaptive thresholding criterion as in [9] for delimiting components in the estimated HS signals. Regarding the SSA method, the same parameters as those in [9] were used. The histograms are shown in Fig. 3 where the left column (a–c) corresponds to the medium and right column (d–f) to high flow rates. Judging by the histograms’ shape and localization on the horizontal axes we can say that the modeling method outperforms the SSA regarding the three performance evaluation parameters. For both medium and high flow rates the modeling method yields position and duration error histograms with a smaller mean and deviation than those of the corresponding SSA histograms. Regarding the SDR both methods give rise to a similar deviation but the modeling method achieves larger mean values. The concrete mean and deviation values for the histograms are given in Table 1. The second part of the present study concerned the artificial and real mixtures processed in terms of false positive and negative

4

3

2

1

0

1

2

3

4

Frequency band

(b) Fig. 5. Average PSDD for the MP–LP and modeling method under (a) medium and (b) high flow rates. The frequency bands 1–4 correspond to: 20–40, 40–70, 70–150 and 150–300 Hz.

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

which strongly overlap with the signal’s eigenvalues. Consequently, the selected eigenvalue pairs contain a significant amount of noise which hampers the S1–S2 component localization by means of adaptive thresholding [9]. According to Table 3 the false positives detection increases significantly over the false negatives detection. This effect is probably due to the envelope distortion through a number of spikes clustered around the local energy maxima of LS. The modeling method, nevertheless, proves to be more robust against the LS, as its statistical measures clearly indicate. Fig. 4 shows a high flow rate time-waveform example in order to illustrate the methods’ performance. In subfigures (a–c) the original respiratory signal, the estimated LS and the estimated HS for the modeling method, respectively are shown. We can see that the estimated HS still contains certain LS residual which cannot be completely removed. Nevertheless, the FHS components can be clearly identified and used for further processing. In subfigures (d–f) we show, for the sake of comparison, a short segment (around 2 cardiac cycles) of the respiratory signal and the output of the SSA and modeling algorithm, respectively. We see that the modeling method is more robust against LS than the SSA regarding both localization and characterization of FHS components. Furthermore, we have compared the SSA and modeling method in terms of mean execution time as a function of the number of

593

samples to be processed. The results are obtained on an Intel 3.19 GHz processor and given in Table 4. Observe that the modeling method outperforms the SSA by almost an order of magnitude independently on the record length. Although this outcome might come as a surprise, it is in fact due to the singular value decomposition (SVD) computation of a large trajectory matrix in the SSA method [9]. 3.5. HS suppression – comparative study This study was performed on the same signals from the previous subsection, but herein the modeling method was compared to the MP–LP in terms of the PSDD parameter. The PSDs were calculated by the Welch method using a 200 ms Hamming window with 50% overlap. The PSDD was evaluated in the frequency subbands 20–40, 40–70, 70–150 and 150–300 Hz and averaged over all trials. Regarding the artificially mixed signals the results are shown in Fig. 5 for (a) medium and (b) high flow rates. Both methods seem to obey a similar trend which exhibits the minimum HS suppression capacity in 40–70 Hz (the maximum PSDD). This is due to the fact that this band contains most of the HS energy; therefore, it will enclose the major part of the HS residual after denoising. The modeling method is superior in the whole analysis range by

SIGNAL

MP-LP

1.5

1 0.8

1

0.6 0.4

Amplitude

Amplitude

0.5

0

-0.5

0.2 0 -0.2 -0.4 -0.6

-1

-0.8 -1.5

0

1

2

3

4

5

6

7

8

-1 0

9

1

2

3

Time [s]

(a)

0.8

0.8

0.6

0.6

0.4

0.4

0.2 0 -0.2

-0.6

-0.8

-0.8

5

Time [s]

(b)

9

6

7

8

9

-0.2 -0.4

4

8

0

-0.6

3

7

0.2

-0.4

2

6

MODELING 1

Amplitude

Amplitude

LS

1

5

(c)

1

-1 0

4

Time [s]

6

7

8

9

-1 0

1

2

3

4

5

Time [s]

(d)

Fig. 6. Illustration of LS estimation: (a) mixture signal, (b) LS, (c) estimated LS by the MP–LP method, (d) estimated LS by the modeling method.

594

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595

MEDIUM FLOW -30

Clean LS MP-LP Modeling

Average PSD [dB]

-35 -40 -45 -50 -55 -60 -65

signal estimation, as most of the energy in the bursts is preserved (compare (c) and (d) to (b)). In order to support the above discussion we have compared average PSDs of the real mixtures, clean LS and estimated LS across the methods. The PSDs of the clean LS were obtained by processing the segments of the LS signals which did not contain HS. The PSDs for medium and high flow rates are shown in Fig. 7(a and b), respectively. We can see that below 200 Hz approximately the MP–LP introduces important differences between the clean and estimated LS for medium flow rates. For high flow rates this effect spreads up to 300 Hz. The LS estimated by the modeling method resembles the clean LS in almost whole analysis range, except for the small differences at low frequencies (below 100 Hz). Therefore, we state that the modeling method is more robust against LS than the MP–LP in terms of HS suppression for medium and high flow rates.

2

10

4. Conclusions

Frequency [Hz]

(a) HIGH FLOW -25 Clean LS MP-LP Modeling

Average PSD [dB]

-30 -35 -40 -45 -50 -55 -60

2

10

Frequency [Hz]

(b) Fig. 7. Average PSD from the subjects for the MP–LP and modeling method under (a) medium and (b) high flow rates.

We have presented a novel approach to HS localization and suppression in LS by means of quasi-periodic modeling. We have shown that inherent quasi-periodicities in HS (envelope and phase) can be compactly modeled by a joint action of single and piecewise polynomials. Such a modeling scheme allows for a correct HS characterization even in conditions of large LS. The model parameters are estimated in an efficient fashion from an overdetermined linear system of equations. The HS model has proven to be capable of fitting well to data. Considering LS localization we have shown that the proposed modeling method outperforms a singular spectrum analysis method in terms of fundamental HS component position, duration and signalto-distortion ratio. The differences become even more prominent for high flow rates, indicating a very good robustness against LS for the proposed method. As for HS suppression in LS, the proposed method has proven to be better than an approach based on multiscale product and linear prediction. It was shown that in terms of average PSD difference between clean and denoised LS for both medium and high flow rates, the proposed method is situated 2–5 dB below the other method. In view of these results, we think that the proposed method is a good candidate for clinical applications which include LS sound analysis. Competing interests None declared.

approximately 3 dB for low and 4 dB for high flow rates. Such a difference in performance is a result of the errors generated by both HS localization and suppression. The MP–LP FHS localization is based on the multiscale product of the wavelet coefficients of the original signal. Although this approach works well for low flow rates, its performance degrades for medium and especially high flow rates. In addition, the HS suppression “gap-filling” based approach is also error prone, because connecting two signal segments to new data introduces errors in the whole frequency range, not only around the HS frequencies. The modeling method, on the contrary, provides a very good localization and characterization of the FHS components (Section 2) under both medium and high flow rate conditions (Section 3). Accordingly, HS suppression will only have a small impact on the LS signal in the HS frequencies. As an illustration, we have shown in Fig. 6 the following waveforms: (a) original mixture, (b) original LS, (c) estimated LS by the MP–LP method, and (d) estimated LS by the modeling method. We see that the synthesized LS signal by the MP–LP method turns to be seriously degraded at the position of the bursts. This effect is due to HS detection mechanism which is very sensible to low signal-tonoise ratios. The modeling method, however, generates a better LS

Funding None declared. References [1] S. Reichert, R. Gass, C. Brandt, E. Andres, Analysis of respiratory sounds: state of the art, clinical medicine insight, Circulatory, Respiratory and Pulmonary Medicine 2 (2008) 45–58. [2] P.J. Arnott, G.W. Pfeiffer, M.E. Tavel, Spectral analysis of heart sounds: relationships between some physical characteristics and frequency spectra of first and second sounds in normal and hypertensives, Journal of Biomedical Engineering 6 (2) (1984) 121–128. [3] A. Sovijarvi, L. Malmberg, G. Charbonneau, F. Vanderschoot, J. Dalmasso, C. Sacco, M. Rossi, J. Earis, Characteristic of breath sounds and adventitious respiratory sounds, Journal of European Respiratory Review 10 (2000) 591–596. [4] M.T. Pourazad, Z. Moussavi, G. Thomas, Heart sound cancellation from lung sound recordings using time–frequency filtering, Medical & Biological Engineering & Computing 44 (2006) 216–225. [5] L.J. Hadjileontiadis, S.M. Panas, Separation of discontinuous adventitious sounds from vesicular sounds using a wavelet-based filter, IEEE Transactions on Biomedical Engineering 44 (12) (1997) 1269–1281. [6] Z. Moussavi, D. Flores, G. Thomas, Heart sound cancellation based on multiscale products and linear prediction, in: Proceedings of IEEE Conference on Engineering in Medicine & Biological Society, 2004, pp. 3840–3843.

M. Zivanovic, M. González-Izal / Biomedical Signal Processing and Control 8 (2013) 586–595 [7] C. Ahlstrom, T. Lanne, P. Ask, A. Johansson, A method for accurate localization of the first heart sound and possible applications, Physiological Measurement 29 (3) (2008) 417–428. [8] A. Yadollahi, Z.M.K. Moussavi, A robust method for heart sounds localization using lung sound entropy, IEEE Transactions on Biomedical Engineering 53 (3) (2006) 497–502. [9] F. Ghaderi, H.R. Mohseni, S. Sanei, Localizing heart sounds in respiratory signals using singular spectrum analysis, IEEE Transactions on Biomedical Engineering 58 (12) (2011) 3360–3367. [10] M.T. Pourazadi, Z. Moussavi, F. Farahmand, R.K. Ward, Heart sounds separation from lung sounds using independent component analysis, in: Proceedings of IEEE Conference on Engineering in Medicine & Biological Society, 2005, pp. 2736–2739. [11] J. Gnitecki, I. Hossain, H. Pasterkamp, Z. Moussavi, Qualitative and quantitative evaluation of heart sound reduction from lung sound recordings, IEEE Transactions on Biomedical Engineering 52 (10) (2005) 1788–1792. [12] J. Gnitecki, Z. Moussavi, Separating heart sounds from lung sounds: accurate diagnosis of respiratory disease depends on understanding noises, IEEE Engineering in Medicine and Biology Magazine 26 (1) (2007) 20–29. [13] K. Basak, S. Mandal, M. Manjunatha, J. Chatterjee, A.K. Ray, Phonocardiogram signal analysis using adaptive line enhancer methods on mixed signal processor, in: Proceedings of International Conference on Signal Processing and Communications, 2010, pp. 1–5. [14] D. Flores-Tapias, Z.M.K. Moussavi, G. Thomas, Heart sound cancellation based on multiscale products and linear prediction, IEEE Transactions on Biomedical Engineering 54 (2) (2007) 234–243. [15] C. Ahlstrom, O. Liljefeldt, P. Hult, P. Ask, Heart sound cancellation from lung sound recordings using recurrence time statistics and nonlinear prediction, IEEE Transactions on Biomedical Engineering 12 (12) (2005) 812–815. [16] T. Tsalaile, S.M. Naqvi, K.S. Nazarpour, J.A. Chambers, Blind source extraction of heart sound signals from lung sound recordings exploiting periodicity of the heart sound, in: Proceedings of International Conference on Acoustics Speech and Signal Processing, 2008, pp. 461–464.

595

[17] M. Zivanovic, M. Gonzalez-Izal, Non-stationary harmonic modelling for ECG removal in surface EMG signals, IEEE Transactions on Biomedical Engineering 59 (6) (2012) 1633–1640. [18] A. Ahmadian, S. Karimifard, H. Sadoughi, M. Abdoli, An efficient piecewise modelling of ECG signals based on hermitian basis functions, in: Proceedings of IEEE Conference on Engineering in Medicine & Biological Society, 2007, pp. 3180–3183. [19] J. Xu, L.G. Durand, P. Pibarot, Nonlinear transient chirp signal modeling of the aortic and pulmonary components of the second heart sound, IEEE Transactions on Biomedical Engineering 47 (10) (2000) 1328–1335. [20] J. Xu, L.G. Durand, P. Pibarot, Extraction of the aortic and pulmonary components of the second heart sound using a nonlinear transient chirp signal model, IEEE Transactions on Biomedical Engineering 48 (3) (2001) 277–283. [21] A. Djebbari, F. Bereksi-Reguig, A new chirp-based wavelet for heart sound time–frequency analysis, Journal on Communications Antenna and Propagation 1 (1) (2011) 92–102. [22] M. Unser, A. Aldroubi, M. Eden, B-spline signal processing: part I. Theory, IEEE Transactions on Signal Processing 41 (2) (1993) 821–833. [23] M. Unser, A. Aldroubi, M. Eden, B-spline signal processing: part II. Efficient design and applications, IEEE Transactions on Signal Processing 41 (2) (1993) 834–848. [24] C. de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978. [25] T. Li, H. Tang, T. Qiu, Y. Park, Best subsequence selection of heart sound recording based on degree of sound periodicity, Electronic Letters 47 (15) (2010) 841–843. [26] L. Chua, A. Deng, Canonical piece-wise representation, IEEE Transactions on Circuits and Systems 35 (1) (1988) 4425–4431. [27] G. Bucci, M. Faccio, C. Landi, New ADC with piecewise characteristic: case study. Implementation of smart humidity sensor, IEEE Transactions on Instrumentation and Measurement 49 (6) (2000) 1154–1166. [28] T. Blu, P. Thevenaz, M. Unser, Linear interpolation revitalized, IEEE Transactions on Image Processing 13 (5) (2004) 471–475. [29] W.J. Hess, Pitch Determination of Speech Signals, Springer, New York, 1993. [30] http://www.med.umich.edu/lrc/psb/heartsounds/index.htm [31] http://darwin.unmc.edu/hnp/m2/breathsounds.htm [32] L. Cohen, Time–Frequency Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1995.