- Email: [email protected]
An advanced method in fetal phonocardiography
Computer Methods and Programs in Biomedicine 71 (2003) 283 /296 www.elsevier.com/locate/cmpb
An advanced method in fetal phonocardiography Pe´ter Va´rady a, Ludwig Wildt b, Zolta´n Benyo´ c,*, Achim Hein d a
Budapest University of Technology and Economics, Department of Control Engineering and Information Technology, BME-IIT, Pa´zma´ny s. 1/d., Room: B.311, H-1117 Budapest, Hungary b Universita¨tsfrauenklinik */Abt. Gyn. Endokrinologie und Reprodunktionsmedizin, Clinics of Obstetrics and Gynaecology at the University of Erlangen-Nu¨rnberg, Universita¨tsstr, 21 /23. D-91045, Erlangen, Germany c Budapest University of Technology and Economics, Department of Control Engineering and Information Technology, BME-IIT, Pa´zma´ny s. 1/d., Room: B.325, H-1117 Budapest, Hungary d Dr. Hein GmbH, Fu¨rther Str. 212, D-90429 Nurnberg, Germany Received 15 March 2002; received in revised form 5 June 2002; accepted 10 July 2002
Abstract The long-term variability of the fetal heart rate (FHR) provides valuable information on the fetal health status. The routine clinical FHR measurements are usually carried out by the means of ultrasound cardiography. Although the frequent FHR monitoring is recommendable, the high quality ultrasound devices are so expensive that they are not available for home care use. The passive and fully non-invasive acoustic recording called phonocardiography, provides an alternative low-cost measurement method. Unfortunately, the acoustic signal recorded on the maternal abdominal surface is heavily loaded by noise, thus the determination of the FHR raises serious signal processing issues. The development of an accurate and robust fetal phonocardiograph has been since long researched. This paper presents a novel two-channel phonocardiographic device and an advanced signal processing method for determination of the FHR. The developed system provided 83% accuracy compared to the simultaneously recorded reference ultrasound measurements. # 2002 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Phonocardiography; Fetal heart rate; Noise cancellation; Acoustic analysis
1. Introduction
1.1. Background
* Corresponding author. Tel.: /36-1-4631410; fax: /36-14632204 E-mail addresses: [email protected] (P. Va´rady), [email protected] (L. Wildt), [email protected] (Z. Benyo´), [email protected] (A. Hein).
The long-term recording of the fetal heart rate (FHR) is the most frequently used diagnostic measurement to determine the fetal health status [1]. FHR monitoring is usually possible after the 24th gestational week. The routine medical analysis of the recorded FHR diagram is of great
0169-2607/02/$ - see front matter # 2002 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 0 2 ) 0 0 1 1 1 - 6
284
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
clinical importance and since its introduction it has led to the drastic reduction of prenatal and postnatal child mortality [2]. The ultrasound cardiography is the most commonly used non-invasive clinical diagnostic tool which provides the accurate and robust long-term determination of the FHR. During the measurements an ultrasound transducer pulses :/2 MHz ultrasound at 10 mW/cm2 intensity from the maternal abdominal surface towards the fetus. The FHR can be easily calculated from the reflected and received ultrasound, since it is modulated by the fetal heart movements. In the past years combined ultrasound transducers have appeared which are capable of detecting the fetal heart activity and the uterine contractions as well (cardiotocography, CTG). Frequent measurement of FHR is recommendable, although high quality CTG devices are so expensive that this technology is not available for home care use. At the time of writing this paper, a hand-held ultrasound device has just appeared (BFD-100 by BiOSYS Co., Ltd, Korea). Unfortunately, this device displays only the actual heart rate and does not provide an FHR diagram, so it has a limited diagnostic value. Another well-known, but less frequently used non-invasive FHR measurement tool is the fetal electrocardiography (FECG) which uses electrodes on the maternal abdominal surface. This method requires more complex signal processing than CTG, since the FECG signal includes the maternal ECG as well. Although there are numerous methods proposed for the rejection of the maternal signal [3 /5]; the automated evaluation of FECGs is less accurate than CTG. Another problem is that the signal quality highly depends on the proper electrode placement. Therefore long-term recording can be inconvenient for the mother and may require electrode adjustment due to fetal movements. The main benefits of FECG compared to CTG are the lower cost and its passive measurement technology [6].
based on the passive acoustic recording of the fetal heart sounds. Auscultation is one of the oldest medical tool in history which has also been applied in fetal diagnostics using specially-formed (Pinnard-type) wooden stethoscopes. The modern form of auscultation called phonocardiography, provides the non-invasive electronic recording and computer aided analysis of the acoustic cardiac signals. Nonfetal phonocardiography was born in the early 70s and since then numerous useful applications have been introduced, such as electronic stethoscopes and heart noise analysis [7 /10]. Unfortunately, the fetal heart activity produces much less acoustic energy and in addition it is surrounded by a highly noisy environment. Therefore the detection of fetal heart sounds raises serious signal processing issues. Although several methods have been introduced [11 /13]; to the best of our knowledge, currently there are no fetal phonocardiographic devices available yet. By the means of short-time Fourier analysis we found that the dominant part of the acoustic energy produced by the fetal heart activity lies in the band 20 Hz B/f B/200 Hz. Each heart beat can be separated into two subbeats: an S1 (systolic) beat, followed by a softer S2 (diastolic) beat (Fig. 1a). Unfortunately, the analysis of heart sounds measured on the maternal abdominal surface is complicated by the following main factors: . acoustic damping through the fruit water and tissues,
1.2. Preliminaries This paper presents an alternative and fully noninvasive FHR measurement method which is
Fig. 1. Excerpts of a noise-free (a) and a noisy (b) fetal phonocardiogram (length: approximately 0.7 s, bandpass filter used: 35 HzB/f B/200 Hz).
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
285
. acoustic noise produced by the fetal movements, . maternal digestive sounds, . sounds of maternal heart activity, . movements of the measurement head during recording (‘shear noise’), . external noise originating from the environment. An excerpt of a noisy fetal phonocardiogram can be seen in Fig. 1b where the first S2 beat is almost missing, and the signal noise (N ) was caused by fetal movements. In spite of the difficulties, the aim of our work was to develop an advanced phonocardiographic FHR monitoring system which provides accurate and robust real-time FHR measurement (compared to CTG) and can be implemented as a lowcost device for home care use.
2. Signal recording
2.1. Recording device A two-channel recording device with two separated electret microphones (MCE-2000 of Panasonic) was developed for the fetal phonocardiographic measurements. The first channel measured the acoustic signals on the maternal abdominal surface (cA), while the second one recorded the external noise (cE). In the case of the abdominal channel the microphone was inserted in the focus of an aluminum measurement cone with a diameter of 70 mm. The external microphone was fixed in a hole on the top of the device case. Both channels were amplified resulting a stereo line level signal in order to allow low-noise signal transmission. The low-noise precision microphone preamplifier stages were built up using rail-to-rail instrumentation amplifiers (INA-122 of BurrBrown), allowing a 9 V single battery operation. The line level stereo signal was digitized by a notebook computer (fs /11025 Hz sampling rate, 16 bits resolution, see Fig. 2a).
Fig. 2. Block diagram (a) and photos (b) of the phonocardiograph.
2.2. Signal records More than 20 records were taken from different women who were between the 28th and 40th week of singleton pregnancy (the mean value was: 36th week). Unfortunately, only 16 records could be used because of some initial problems (the amplifier stages were overdriven). Nine records were taken by simultaneous CTG monitoring in order to allow a performance comparison. The average duration of the records was about 8 min. Fig. 3 shows a section of a two-channel phonocardiographic record. Please note that the periodic
Fig. 3. Excerpt of a two-channel phonocardiographic record. The scaling of the y axis is in normalized sample units (nu).
286
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
peaks in the abdominal channel are related to the maternal heart sounds.
3. Signal processing The block diagram of the presented signal processing method is depicted in Fig. 4. As realtime signal processing was targeted, a time window was introduced. In each processing step only a section of both input channels is processed. The length of the signal sections is equal to the time window size (Twproc). The processing control algorithm used adaptive window overlapping in order to avoid the loss of S1 or S2 beats near the window borders. Please note that the signals of both channels cA and cE were normalized in the range (/0.5, 0.5) before the entire signal processing. 3.1. External noise The open-air microphone of the external channel recorded only the noise originating from the environment. The abdominal microphone recorded primarily the sounds originating from the maternal abdomen, but these were mixed with a damped version of the external noise. According to our previous studies, the transfer characteristic between the external noise and the abdominal microphone can be modeled as a 2nd order lowpass filter with an upper band limit of approximately 400 Hz [14]. We examined and tested the following methods in order to achieve an optimum external noise cancellation performance: (1) Signal difference: the external channel was low-pass filtered using the 2nd order model
transfer characteristic, and the filtered signal was subtracted from the abdominal channel. This evident method brought only moderate success, since the applied stationary model was only the approximation of a non-stationary and patientspecific one. However this method required the least computing. (2) Adaptive linear filtering with one-step-ahead prediction: this method yielded better results, although the adaptive updating of the filter coefficients with the least squares method required significantly more computation time. (3) Signal separation using higher order statistics: the cancellation of external noise can be stated as a two-source, two-drain problem for blind source separation or independent component analysis. The separation algorithms detailed in Refs. [15] and [16] produced reasonable performance, although the non-stationary external noise required a frequent recalculation of the separation matrices, entailing a lot of computation. (4) Wavelet-based noise cancellation: the 6th order decomposition of both signal channels based on the Coiflet-2 wavelet [17] and subsequent application of adaptive inter-channel coefficient thresholding brought the best trade-off between cancellation performance and speed. This method is detailed in the following subsection.
3.2. Wavelet-based adaptive noise cancellation Wavelet transform is a frequently used tool for denoising signals. The cancellation of the external noise in our case is a two-channel problem. A new method was developed which uses adaptive interchannel thresholding of the wavelet transformed filter coefficients.
Fig. 4. Block diagram of the signal processing.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
The dyadic wavelet transformed of a signal f(t) on the i-th dyadic scale can be written as: 1 W ff (u; 2 )g pffiffiffiffi 2i i
g f (t)c
tu 2i
dtf + c¯ 2i (u)
where c is the orthogonal mother wavelet, and: 1 t ¯ ffiffiffiffi p c2i (t) c2i (t) c 2i 2i Under certain conditions, the discrete implementations of the dyadic wavelet transform can be based on digital filter banks [17]. Let us define: di [n]W ff (n; 2i )g hf (t); c2i (tn)i ai [n]hf (t); c¯ 2i (tn)i where i E/0. For a given filter x [n ], xi [n ] denotes the filter obtained by inserting 2i /1 zeros between every filter coefficients (‘holes algorithm’), a0[n] is the original signal, and its wavelet transformed (decomposition) at the i/1-th scale is: ai1[n ]/ aj + hi [/n], and di1[n ]/dj + gj [/n ] where a defines the approximated signal, and d the signal details at the given scale; h and g are biorthogonal filters corresponding to the selected mother wavelet (Table 1). The reconstruction of the signal is done by the inverse wavelet transform, which has the following form at the i- th scale: 1 ai [n] (a 2
i1 +hi [n]di1
ˆ
+ gˆi [n])
where hˆ and gˆ are the dual filters of h and g, respectively. In case of biorthogonal filters: gˆ/ /h and hˆ/ /g. The Coiflet wavelet with four moments equal to zero (Coiflet-2) was selected as the mother wavelet.
According to the filter bank based realization, the length of the processing time window should be adjusted as Twproc /n + 8192/fs, where n is a positive integer. We selected n /4, resulting a Twproc /3 s. Thus one processing time window consisted of 32768 sample points of both input channels. The two channels were wavelet decomposed up to the 6th order (i/0, 1,. . ., 5), resulting two sets of coefficients: a6A, d6A, d5A, d4A, d3A, d2A, d1A and similarly: a6E, d6E,. . ., d1E, where the A and E letters in the indices refer to the abdominal and the external channel, respectively. The aim of the external noise cancellation is to construct a new set of coefficients: a6F, d6F,. . ., d1F, which can be reconstructed (inverse transformed) as the abdominal signal free of external noise. The algorithm of Scheme 1 was used to obtain this set of coefficients. Please note that only the detail coefficients di A of the abdominal channel are thresholded, whereby the details di E of the external channel are adoptively subtracted from di A.
3.3. Bandpass filtering The adaptive wavelet denoising operates across the whole recorded signal band (0 B/f B/fs/2), however it removes only the external noise. A bandpass digital filter was designed to limit the signal into the band: 35 Hz B/f B/200 Hz. The section of the lower band limit was based on the experimental result that a significant part of the disturbing maternal heart and digestive sounds lies
Table 1 Base coefficients of the equivalent filter bank h0
g0
/0.0007, /0.0018, 0.0056, 0.0237, /0.0594, /0.0765, 0.4170, 0.8127, 0.3861, /0.0674, /0.0415, 0.0164
/0.0164, /0.0415, 0.0674, 0.3861, /0.8127, 0.4170, 0.0765, /0.0594, /0.0237, 0.0056, 0.0018, /0.0007
287
Scheme 1. Wavelet-based inter-channel denoising.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
288
Table 2 Coefficients a and b of the highpass and lowpass filter a /b
Highpass (f /35 Hz)
Lowpass (f B/200 Hz)
a
1.0000, /7.8832, 27.1886, /53.5868, 66.011, /52.044, 25.6459, /7.2217, 0.8897 0.9432, /7.5460, 26.4109, /52.822, 66.027, /52.8219, 26.4109, /7.5460, 0.9432
1.0000, /7.4158, 24.0798, /44.7159, 51.939, /38.642, 17.979, /4.7840, 0.5573 /0, /0, /0, /0, /0, /0, /0, /0, /0
b
below this border, while the fetal heart sounds are not dominantly present there anymore. Due to stability reasons we designed two separate 8th order Butterworth-type IIR digital filters: one low-pass and one high-pass filter. The signal f[n ] is filtered as s[n ]:
s[n]
nX b 1 k1
b[k]f [nk 1]
nX a 1
a[k]s[nk 1]
k2
where a and b are the filter coefficients listed in Table 2, and na /nb /8. The effect of the external noise cancellation and bandpass filtering is illustrated in Fig. 5 using a section of the signal in Fig. 3.
3.4. Locating the amplitude bursts The Wavelet-based external noise cancellation and the bandpass filtering together produce a kind of preprocessed signal s, in which the S1 and S2 beats correspond to amplitude bursts (or simply: bursts, see Fig. 5d). Obviously the band-pass filtering cannot fully remove the internal noise. Disturbances still remain in the filtered signal, mainly caused by fetal movements, fluctuations in the fruit water and ‘shear noise’ of the measurement head. Unfortunately, this noise can cause similar bursts to that of S1 and S2 beats have (Fig. 1b). For the location of bursts several methods were proposed in the field of non-fetal phonocardiography. These methods are usually based on the
Fig. 5. The effect of external noise cancellation and bandpass filtering shown on a section of Fig. 3: (a) abdominal channel, (b) abdominal channel after simple bandpass filtering, (c) abdominal channel after external noise cancellation, (d) abdominal channel after external noise cancellation and bandpass filtering.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
thresholding of some derived signal: either the signal’s absolute value, its squared amplitude, or its Shannon-energy. A fair comparison of the various burst detection methods is provided by Liang et al. [18]. Unfortunately, these methods are only concerned with the magnitude of the derived signal and do not examine the shape and timing of the bursts. As a consequence these methods were insufficient in the case of fetal phonocardiograms where the signal-to-noise ratio is far less than in non-fetal ones. To improve burst detection, the cross-correlation with a reference burst was calculated. This method can measure all three burst parameters at once: magnitude, shape and timing. Since the exact waveform of the bursts varies considerably even during one recording period, not the preprocessed
289
signal, but its envelope was cross-correlated with the reference burst. The envelope e of the preprocessed signal s can be derived from measuring the amplitude difference between the pmin local minimums and the successive pmax local maximums of the signal. The value of each calculated amplitude difference is inserted into e until the next local minimum / maximum pair of s is reached (time spanning). Thus, the resulted envelope signal reflects the amplitude dynamic of the original signal and the same timing is kept. The envelope algorithm is shown in the flow chart in Fig. 6. The section of a preprocessed signal and its envelope is illustrated in Fig. 8a and b, respectively. In the next step the envelope signal was downsampled with the factor 1:10 which allowed faster processing time in the further steps. The downsampling was realized by using a simple moving average filter. Finally, the cross-correlation sum r measured the similarity of the downsampled envelope signal e? to a reference burst envelope y of length k and l, respectively: r[n]
l X
e?[mn]y[m];
n 0; 1; . . . ; k
m0
The cross-correlation of the downsampled envelope signal with the selected reference burst envelope (Fig. 7) produced the signal in Fig. 8c. The peak values in r show the locations of the acoustic amplitude bursts and hence the possible location of the S1 and S2 beats. The peaks were detected with an algorithm that is an extended and modified version of the method used in Ref. [13]. The main steps are:
Fig. 6. Flow chart of the envelope algorithm.
Fig. 7. The reference burst envelope.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
290
t[i]
t[i]m[i] t[i 1]m[i 1] m[i] m[i 1]
. the i/1-th element is deleted and n is adjusted. . eight peaks with the largest values were selected now again from a, and their mean value l is calculated. . a three-state possibility vector p is obtained: if a [i]B/l/16 then p [i]/0, else if a[i] /l /2 then p [i]/1, else p [i]/0.5, where i/1, 2,..., n . return as result those elements of m , t and p , where p[i] /0, i /1, 2,. . ., n . 3.5. Locating the S1 and S2 beats Fig. 8. The internal processing steps: (a) preprocessed signal, (b) envelope signal, (c) cross-correlated signal.
. all local maximums (peaks) of r are located, resulting in three n -dimensional vectors: a , m and t, where a contains the height of the peaks relative to the preceding minimums, m contains the magnitudes of the peaks, and t their locations relative to the start of the recording, respectively. . those jth elements of each a, m and t are deleted where: t[i/1]/t [i ]B/10 ms, i/1, 2,. . ., n/ 1, and if a [i] B/a [i/1] then j/i, else j /i/1. After this n is adjusted. . eight peaks with the largest values were selected from a , and their mean value l is calculated. . those elements of a, m and t are deleted where: a [i] B/l /8, i/1, 2,. . ., n/1. After this n is adjusted. . if t[i/1]/t [i ]B/40 ms, where i/1, 2,. . ., n/ 1, then the i-th element is merged with the i/1-th element according to: a[i]
The resulted m , l and p vectors describe the magnitude, time location and probability of the detected bursts, respectively. Since several bursts can be caused by noise, the next task is to select only the bursts of S1 and S2 beats. The proposed algorithm is depicted in the state machine of Fig. 9. The individual S1 and S2 bursts can be distinguished from the noise bursts by taking various physiological timings into a consideration (Table 3). The result of the detection is a beat information matrix M, where cash row corresponds to an S1 / S2 beat pair: TS1, TS2, PS1, PS2, where TS1, TS2 are the time locations of the S1 and S2 beats relative to the start of recording and PS1, PS2, are their three-state probabilities respectively. The algorithm was constructed with the assumption that the deletion of the acoustically stronger S1 beats is more probable. Therefore the method is
a[i] a[i 1] 2
m[i]
m[i] m[i 1] 2
Fig. 9. State machine of the S1 /S2 detection.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
Table 3 Meaning of the timing values Value
Meaning
TS1S1 TS1S2 TS1S1min TS1S1max TS1S2min TS1S2max TS2S1min TS1S1dev TS1S2dev
Time elapsed between two consecutive S1 beats Time elapsed between an S1 /S2 beat pair Minimal allowed time between two consecutive S1s Maximal allowed time between two consecutive S1s Minimal allowed time between an S1 /S2 pair Maximal allowed time between an S1 /S2 pair Minimal time to be elapsed between an S2 and an S1 Allowed beat-to-beat deviation of TS1S1 Allowed beat-to-beat deviation of TS1S2
291
First S1 found . Appending beat information: h /i (position of first S1 found) if S2 could also be located (p[g ]/0), then append row: t[h]t[g]p [h ]p[g ], and i/g/1, TS1S2 /t[g ]/ t [h ] else append row: t[h] 0 p[h ] 0, and i/i/1, TS1S2 /0 do TR3 Next search
based on the search for S1 of beats. The detection of an S2 beat only upgrade, the reliability of the S1 detection, however, it is possible that no S2 after a given S1 can be detected (i.e. TS2 /PS2 /0). All timing values are adoptively tuned during the detection. If within the given timing no subsequent S1 beat can be found, the algorithm reinitializes itself by setting the default timing values. The following part of this section details the states and the state transitions of the detection algorithm: Start . General initialization: TS1S2dev /40 ms, TS1S2dev /25 ms, TS2S1min / 150 ms, i/1, do TR1 Initial search . Search initialization: TS1S1min /300 ms, TS1S1max /750 ms, TS1S2min /100 ms, TS1S2max /210 ms, s /0 . Checking if the current burst belongs to an S1 beat : start j from i/1, and search for an S2 beat in a range, where TS1S2min B/t [j ]/t [i ]B/TS1S2max and let g the location in this range where p [j] /0 and m [j]B/m [i] if p [i ]/p [g] E/1 then do TR2 else i/i/1 and goto start of this state
. Checking the index: if i /n do TR9 . Checking TS1S1 : if s /1 do TR7 else: TS1S1 /t[i ]/t [h ] (time elasped since last S1 found) if TS1S1 B/TS1S1min or TS1S1 B/TS1S2/ TS2S1min then i/i/1, s/s/1, goto start of this state else if TS1S1 /TS1S1max do TR7, else: . Update TS1S2: if TS1S2 /0 then TS1S2min /TS1S2/TS1S2dev, TS1S2max /TS1S2/ TS1S2dev . Search for the next possible S1 : start j from i/1 and search for an S2 beat in a range, where TS1S2min B/t[j]/t[i]B/TS1S2max and let g the location in this range where p [j]/0 and m [j ]B/m [i ] if p[i]/p [g ]E/1 then do TR4 else i/i/1 and goto start of this state Next S1 found . Checking the index: if i /n do TR8 . Checking TS1S1 : TS1S1 /t[i ]/t [h ] (time elapsed since last S1 found) if TS1S1 B/TS1S1min or TS1S1 /TS1S1max do TR6, else: . Update TS1S1: TS1S1min /TS1S1/TS1S1dev, TS1S1max /TS1S1/ TS1S1dev . Appending beat information : h /i (position of S1 found) if S2 could also be located (p[g ]/0), then
292
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
append row: t[h ]t [g ]p[h ]p [g], and i/g/1, TS1S2 /t[g ]/ t[h ] else append row: t[h ] 0 p [h] 0, and i/i/1, TS1S2 /0 do TR5 In order to keep the notations simple and easyto-understand, the following two features of the real algorithm were not described above: . The search for S1/S2 pairs (in states ‘Initial search’ and ‘Next search’) is continued until t[j]/t[i]B/TS1S1max. If more than one possible S1 locations were found, then the one which has the largest magnitude (m ) is selected. . The TS1S1 and TS1S2 values are calculated using the average of the four subsequently preceding TS1S1 and TS1S2 values, respectively. If no four subsequent values are available (eg. after re-
Fig. 10. Result of the S1 /S2 detection and the calculated FHR values using the signal excerpt in Fig. 3.
initalization), the last single value is used, like in this description above. Fig. 10 illustrates the result of the S1 /S2 pair detection. The top part of this figure shows the cross-correlated envelope signal of the excerpt in Fig. 3. The small dot, on the signal peaks displaying the detected amplitude bursts. The dots in the small squares and circles correspond to the detected S1 and S2 beats, respectively. The left table in the bottom part of the Fig. 10 displays the resulted beat information matrix M . 3.6. Calculation of FHR The last element of the signal processing is the FHR determination. A given FHR value is calculated based on the TS1S1 values inside a time window with a given length (Twfhr /3000 ms). Unfortunately, since some of the S1 beats cannot be obtained, the simple averaging of the TS1S1 times in the selected window is not sufficient. However, we can assume that all beats included in the beat information matrix are valid beats. The proposed FHR calculation method uses only the TS1S1 values and assigns one FHR value to each of them. The method is based on the following main steps: . i indexes the rows of the beat info matrix (i/1, 2,. . ., n ). . k subsequent beats are selected using the time window: TS1S1[i], TS1S1[i/1], . . ., TS1S1[i/k ], where: TS1S1[k ]B/TS1S1[i]/Twfhr . Using the selected k beats an d distance vector is calculated (time elapsed between the beats): d [j]/TS1S1[j/1]/TS1S1[j], where j/1, 2,. . ., k/1 . a validity factor is calculated; v /v1/ v2 where v1 is the maximum number of elements in d , having no more deviation than TS1S1dev and v2 is the number of the neighboring (subsequent) beats amongst them. . if v /3 then the FHR[i] is calculated by averaging the elements of d determined by v1, else no valid FHR could be detected inside Twfhr and FHR[i] is assigned to the last valid FHR value.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
. i/i/1, and goto to the second step. If no valid FHR can be detected over a period of Twfhr/2, then FHR[i ]/0. In the case of the signal in Fig. 10, the resulted FHR values and the pulse values (Pulse[i]/ 60 000/FHR[i]) are displayed in the bottom right table, in ms and bpm units, respectively. Since one FHR value is assigned to each of the S1 beats, the time proportional graphical representations of the FHR values in x, y coordinates are: TS1S1[i], FHR[i].
4. Results The phonocardiographic device was designed and a prototype device was used for the recording. The signal processing was developed in Matlab and after the successful evaluation it was implemented in C//, under Windows NT (Fig. 11). The developed software has the following key features: . real-time signal processing, . input signals either from a wave file or directly from the line level audio input,
293
. numeric display of the actual FHR (in bpm units), . storing and printing of beat information and of the FHR. The evaluation of the developed signal processing method was based on the reference CTG paper records. The paper records were scanned and divided by a grid with a resolution of 3 bpm along the FHR axis, and 5 s along at the time axis. The FHR value of each grid area was compared to the average phonocardiographic FHR value in the corresponding interval. Those positions were not evaluated where the reference CTG record did not include any estimable FHR value. Two of the records yielded more than 90% accuracy, three records were between 80/90%, another two records were within the 70 /80% limit, while the remaining two (rather poor quality) records were below the 70% accuracy. The overall accuracy of the method was: 83%. Fig. 12 shows two FHR sections, both obtained by CTG and by the developed phonocardiographic software. The upper part of the figure shows a steady (possibly sleeping) fetus, while the bottom part illustrates the case of an active and strongly moving fetus. Please note that the CTG has also failed to track the FHR during the strong fetal movements.
5. Discussion and conclusions
Fig. 11. The real-time phonocardiographic software (1 */the calculated FHR diagram, 2/3 */buttons for start/stop processing, 4 */button to open the settings dialog to alter parameters, 5 */exit button, 6 */the actual FHR value, 7 */signal source selection (audio line in or file).
The non-invasive routine measurement of FHR is of great clinical relevance. This paper presented an advanced phonocardiographic method for the determination of the FHR. The phonocardiographic recording technology offers a passive and therefore absolutely non-invasive diagnostic measurement. Fetal phonocardiography has been since long researched, however there is no fetal phonocardiographic device available yet. Existing methods use single channel recording with specialized pressure [11,12] or acoustic sensors [13]. Unfortunately, only the usage of a proper measurement technology itself is not sufficient since the various noise factors make the exact detection of the acoustic
294
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
Fig. 12. FHR diagrams of a steady and a moving fetus obtained by CTG and the phonocardiographic method (vertical scaling is in bpm, horizontal unit is 1 min).
fetal heart beats very difficult. This problem can be handled only by using sophisticated signal processing methods. This paper presented an advanced phonocardiographic system which put the main emphasis not on the acoustic sensors, but on the signal processing issues. The presented phonocardiographic device is built up from commercially available components (high sensitivity electret microphones, low-noise precision instrumentation amplifiers). An innovative solution is the two-channel recording which provided the background for the efficient external noise cancellation.
The determination of the FHR is provided by an advanced real-time signal processing method which consists of several new solutions. The wavelet-based denoising uses an adaptive crosschannel coefficient adjustment in order to cancel external noise more efficiently than any other examined methods. The cross-correlation based amplitude burst detection allowed the precise localization of those bursts, which may correspond with the highest probability to an S1 or S2 beat. The advanced rule-based selection of the S1 and S2 beats by the finite state machine and physiological timing rules is a key point of the entire system. The
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
determination of the periodicity based on the confidence analysis provided an accurate and robust FHR measurement which is comparable to ultrasound CTG. The presented system offers an optimal trade-off between complexity and feasibility and is viable in low-cost, standalone, battery-powered FHR monitoring devices which could be available also for home care use. Although the introduced method produced encouraging results, it must be underlined that a comprehensive evaluation base on a huge number of clinical validation measurements is still needed
6. Future work Our future work will focus on directions of further development The measurement device will be rebuilt including automatic gain, control unit in order in avoid manual gain adjustment and to allow optimal signal quality. Also the signal processing methods can be further enhanced. The selection of the reference burst envelope used for the cross-correlation (Fig. 7) needs a comprehensive validation which should be based on a broad clinical study. The key of the entire signal processing is the proper selection of the S1 and S2 beats. Thus, instead of the finite state machine we intend to use fuzzy or neuro-fuzzy inference system for the selection of the S1 and S2 bursts. The current system presents the FHR diagram without any diagnostic assessment. Although the analysis of the FHR is a complex problem, we would like to develop an automated assessment method which could be useful especially in home care applications. One possible approach is the modelling of the underlying physiological phenomena based on the FHR data [19,20], since the diagnostic information could be obtained from the estimated parameters.
Acknowledgements The phonocardiographic records were taken partly at the Universitats frauenklinik Erlangen,
295
Germany and at the Hospital Universitaire Taher Sfar, Mahdia, Tunesia. The author of this paper would like to thank the supporting help of Andreas Hammel, MD in Erlangen, Mohamed Bannour, MD in Mahdia and Lanouar M. Chouk in Erlangen at taking the records.
References [1] G.S. Dawes, C.R. Houghton, C.W. Redman, G.H. Visser, Pattern of the normal human fetal heart rate, Br. J. Obstet. Gynacol. 89 (1982) (1982) 276 /284. [2] G.S. Dawes, M. Moulden, C.W. Redman, Improvements in computerized fetal heart rate analysis antepartum, J. Perinatal Med. 24 (1996) (1996) 25 /36. [3] Y. Tal, S. Akselrod, Feral heart rate detection by a special transformation method, Proc. IFEE Comp. Soc. Conf. Computers in Cardiology (1990) 275 /278. [4] K.A. Azad, Z.M. Darus, M.A. Ali, Development of a fuzzy rule-based QRS detection algorithm for fetal and maternal heart rate monitoring, Proc. Int. Conf. IEEFE, EMBS 20 (1998) 170-171. [5] J.C. Echeverria, R. Ortiz, N. Ramirez, V. Medina, R. Gonzalez, A reliable method for abdominal ECG signal processing, Computers in Cardiology 25 (1998) 529 /532. [6] C. Lin, H. Wu, T. Liu, N. Lee, T. Kou, S. Young, A portable monitor for fetal heart rate and uterine contraction, IEEE EMBS Magazine (1997 Nov/Dec) 80 /84. [7] J.M. Bentley, P.M. Grant, J.T.E. McDonnel, Time-frequency and time-scale techniques for the classification of native and bioprosthetic heart valve sounds, IEEE Trans. Biomed. Eng. 45 (1998) 125 /128. [8] I. Cathers, Neural network assisted auscultation, Artificial Intelligence in Medicine 7 (1995) 53 /66. [9] J. Ritola, S. Lukkarinen, Comparison of time-frequency distributions in the heart sounds analysis, Medical and Biological Engineering and Computing 34 (1996) 89 /90. [10] X. Zhang, L.G. Durand, L. Senhadji, H.C. Lee, J.L. Coatrieux, Time-frequency scaling transformation of the phonocardiogram based on the matching pursuit method, IEEE Trans. Biomed. Eng. 45 (1998) 972 /979. [11] A.J. Zuckerwar, R.A. Pretlow, J.W. Stoughton, D.A. Baker, Development of a piezopolimer pressure sensor for fetal heart rate monitor, IEEE Trans. Biomed. Eng. 40 (1993) 963 /969. [12] D.G. Talbert, W.L. Davies, F. Johnson, N. Abraham, N. Colley, D.P. Southall, Wide bandwidth fetal phonocardiography using a sensor matched to the compliance of the mother’s abdominal wall, IEEE Trans. Biomed. Eng. 33 (1986) 175 /181. [13] F. Kova´cs, M. To¨ro¨k, I. Habermajer, A rule-based phonocardiographic method for long-term fetal heart rate monitoring, IEEE Trans. Biomed. Eng. 47 (2000) 124 /130.
296
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
[14] P. Va´rady, I. Gross, M.L. Chouk, A. Hein, Analysis of the fetal heart activity by the means of phonocardiography. Proc. IFAC Conference Telematics and Application, Weingarten, Germany, (2001) 41 /46. [15] J.F. Cardoso, A. Soulouminac, Blind beamforming for non-Gaussian signals, IEE-Proc. F 140 (1993) 362 /370. [16] C. Hutten, J. Herault, Blind separation of source, part I: an adaptive algorithm based on neuromimetic architecture, Signal Processing 24 (1991) 1 /10.
[17] G. Strang, T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, 1996. [18] H. Liang, S. Lukkarinen, I. Hartimo, Heart sound segmentation algorithm based on heart sound envelogram, Computers in Cardiology 24 (1997). [19] E.R. Carson, C. Cobelli, Modelling Methodology for Physiology and Medicine, Academic Press, San Diego, 2000. [20] M. Di Rienzo, G. Mancia, Computer Analysis of Cardiovascular Signals, IOS Press, Amsterdam, 1995.
An advanced method in fetal phonocardiography Pe´ter Va´rady a, Ludwig Wildt b, Zolta´n Benyo´ c,*, Achim Hein d a
Budapest University of Technology and Economics, Department of Control Engineering and Information Technology, BME-IIT, Pa´zma´ny s. 1/d., Room: B.311, H-1117 Budapest, Hungary b Universita¨tsfrauenklinik */Abt. Gyn. Endokrinologie und Reprodunktionsmedizin, Clinics of Obstetrics and Gynaecology at the University of Erlangen-Nu¨rnberg, Universita¨tsstr, 21 /23. D-91045, Erlangen, Germany c Budapest University of Technology and Economics, Department of Control Engineering and Information Technology, BME-IIT, Pa´zma´ny s. 1/d., Room: B.325, H-1117 Budapest, Hungary d Dr. Hein GmbH, Fu¨rther Str. 212, D-90429 Nurnberg, Germany Received 15 March 2002; received in revised form 5 June 2002; accepted 10 July 2002
Abstract The long-term variability of the fetal heart rate (FHR) provides valuable information on the fetal health status. The routine clinical FHR measurements are usually carried out by the means of ultrasound cardiography. Although the frequent FHR monitoring is recommendable, the high quality ultrasound devices are so expensive that they are not available for home care use. The passive and fully non-invasive acoustic recording called phonocardiography, provides an alternative low-cost measurement method. Unfortunately, the acoustic signal recorded on the maternal abdominal surface is heavily loaded by noise, thus the determination of the FHR raises serious signal processing issues. The development of an accurate and robust fetal phonocardiograph has been since long researched. This paper presents a novel two-channel phonocardiographic device and an advanced signal processing method for determination of the FHR. The developed system provided 83% accuracy compared to the simultaneously recorded reference ultrasound measurements. # 2002 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Phonocardiography; Fetal heart rate; Noise cancellation; Acoustic analysis
1. Introduction
1.1. Background
* Corresponding author. Tel.: /36-1-4631410; fax: /36-14632204 E-mail addresses: [email protected] (P. Va´rady), [email protected] (L. Wildt), [email protected] (Z. Benyo´), [email protected] (A. Hein).
The long-term recording of the fetal heart rate (FHR) is the most frequently used diagnostic measurement to determine the fetal health status [1]. FHR monitoring is usually possible after the 24th gestational week. The routine medical analysis of the recorded FHR diagram is of great
0169-2607/02/$ - see front matter # 2002 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 0 2 ) 0 0 1 1 1 - 6
284
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
clinical importance and since its introduction it has led to the drastic reduction of prenatal and postnatal child mortality [2]. The ultrasound cardiography is the most commonly used non-invasive clinical diagnostic tool which provides the accurate and robust long-term determination of the FHR. During the measurements an ultrasound transducer pulses :/2 MHz ultrasound at 10 mW/cm2 intensity from the maternal abdominal surface towards the fetus. The FHR can be easily calculated from the reflected and received ultrasound, since it is modulated by the fetal heart movements. In the past years combined ultrasound transducers have appeared which are capable of detecting the fetal heart activity and the uterine contractions as well (cardiotocography, CTG). Frequent measurement of FHR is recommendable, although high quality CTG devices are so expensive that this technology is not available for home care use. At the time of writing this paper, a hand-held ultrasound device has just appeared (BFD-100 by BiOSYS Co., Ltd, Korea). Unfortunately, this device displays only the actual heart rate and does not provide an FHR diagram, so it has a limited diagnostic value. Another well-known, but less frequently used non-invasive FHR measurement tool is the fetal electrocardiography (FECG) which uses electrodes on the maternal abdominal surface. This method requires more complex signal processing than CTG, since the FECG signal includes the maternal ECG as well. Although there are numerous methods proposed for the rejection of the maternal signal [3 /5]; the automated evaluation of FECGs is less accurate than CTG. Another problem is that the signal quality highly depends on the proper electrode placement. Therefore long-term recording can be inconvenient for the mother and may require electrode adjustment due to fetal movements. The main benefits of FECG compared to CTG are the lower cost and its passive measurement technology [6].
based on the passive acoustic recording of the fetal heart sounds. Auscultation is one of the oldest medical tool in history which has also been applied in fetal diagnostics using specially-formed (Pinnard-type) wooden stethoscopes. The modern form of auscultation called phonocardiography, provides the non-invasive electronic recording and computer aided analysis of the acoustic cardiac signals. Nonfetal phonocardiography was born in the early 70s and since then numerous useful applications have been introduced, such as electronic stethoscopes and heart noise analysis [7 /10]. Unfortunately, the fetal heart activity produces much less acoustic energy and in addition it is surrounded by a highly noisy environment. Therefore the detection of fetal heart sounds raises serious signal processing issues. Although several methods have been introduced [11 /13]; to the best of our knowledge, currently there are no fetal phonocardiographic devices available yet. By the means of short-time Fourier analysis we found that the dominant part of the acoustic energy produced by the fetal heart activity lies in the band 20 Hz B/f B/200 Hz. Each heart beat can be separated into two subbeats: an S1 (systolic) beat, followed by a softer S2 (diastolic) beat (Fig. 1a). Unfortunately, the analysis of heart sounds measured on the maternal abdominal surface is complicated by the following main factors: . acoustic damping through the fruit water and tissues,
1.2. Preliminaries This paper presents an alternative and fully noninvasive FHR measurement method which is
Fig. 1. Excerpts of a noise-free (a) and a noisy (b) fetal phonocardiogram (length: approximately 0.7 s, bandpass filter used: 35 HzB/f B/200 Hz).
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
285
. acoustic noise produced by the fetal movements, . maternal digestive sounds, . sounds of maternal heart activity, . movements of the measurement head during recording (‘shear noise’), . external noise originating from the environment. An excerpt of a noisy fetal phonocardiogram can be seen in Fig. 1b where the first S2 beat is almost missing, and the signal noise (N ) was caused by fetal movements. In spite of the difficulties, the aim of our work was to develop an advanced phonocardiographic FHR monitoring system which provides accurate and robust real-time FHR measurement (compared to CTG) and can be implemented as a lowcost device for home care use.
2. Signal recording
2.1. Recording device A two-channel recording device with two separated electret microphones (MCE-2000 of Panasonic) was developed for the fetal phonocardiographic measurements. The first channel measured the acoustic signals on the maternal abdominal surface (cA), while the second one recorded the external noise (cE). In the case of the abdominal channel the microphone was inserted in the focus of an aluminum measurement cone with a diameter of 70 mm. The external microphone was fixed in a hole on the top of the device case. Both channels were amplified resulting a stereo line level signal in order to allow low-noise signal transmission. The low-noise precision microphone preamplifier stages were built up using rail-to-rail instrumentation amplifiers (INA-122 of BurrBrown), allowing a 9 V single battery operation. The line level stereo signal was digitized by a notebook computer (fs /11025 Hz sampling rate, 16 bits resolution, see Fig. 2a).
Fig. 2. Block diagram (a) and photos (b) of the phonocardiograph.
2.2. Signal records More than 20 records were taken from different women who were between the 28th and 40th week of singleton pregnancy (the mean value was: 36th week). Unfortunately, only 16 records could be used because of some initial problems (the amplifier stages were overdriven). Nine records were taken by simultaneous CTG monitoring in order to allow a performance comparison. The average duration of the records was about 8 min. Fig. 3 shows a section of a two-channel phonocardiographic record. Please note that the periodic
Fig. 3. Excerpt of a two-channel phonocardiographic record. The scaling of the y axis is in normalized sample units (nu).
286
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
peaks in the abdominal channel are related to the maternal heart sounds.
3. Signal processing The block diagram of the presented signal processing method is depicted in Fig. 4. As realtime signal processing was targeted, a time window was introduced. In each processing step only a section of both input channels is processed. The length of the signal sections is equal to the time window size (Twproc). The processing control algorithm used adaptive window overlapping in order to avoid the loss of S1 or S2 beats near the window borders. Please note that the signals of both channels cA and cE were normalized in the range (/0.5, 0.5) before the entire signal processing. 3.1. External noise The open-air microphone of the external channel recorded only the noise originating from the environment. The abdominal microphone recorded primarily the sounds originating from the maternal abdomen, but these were mixed with a damped version of the external noise. According to our previous studies, the transfer characteristic between the external noise and the abdominal microphone can be modeled as a 2nd order lowpass filter with an upper band limit of approximately 400 Hz [14]. We examined and tested the following methods in order to achieve an optimum external noise cancellation performance: (1) Signal difference: the external channel was low-pass filtered using the 2nd order model
transfer characteristic, and the filtered signal was subtracted from the abdominal channel. This evident method brought only moderate success, since the applied stationary model was only the approximation of a non-stationary and patientspecific one. However this method required the least computing. (2) Adaptive linear filtering with one-step-ahead prediction: this method yielded better results, although the adaptive updating of the filter coefficients with the least squares method required significantly more computation time. (3) Signal separation using higher order statistics: the cancellation of external noise can be stated as a two-source, two-drain problem for blind source separation or independent component analysis. The separation algorithms detailed in Refs. [15] and [16] produced reasonable performance, although the non-stationary external noise required a frequent recalculation of the separation matrices, entailing a lot of computation. (4) Wavelet-based noise cancellation: the 6th order decomposition of both signal channels based on the Coiflet-2 wavelet [17] and subsequent application of adaptive inter-channel coefficient thresholding brought the best trade-off between cancellation performance and speed. This method is detailed in the following subsection.
3.2. Wavelet-based adaptive noise cancellation Wavelet transform is a frequently used tool for denoising signals. The cancellation of the external noise in our case is a two-channel problem. A new method was developed which uses adaptive interchannel thresholding of the wavelet transformed filter coefficients.
Fig. 4. Block diagram of the signal processing.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
The dyadic wavelet transformed of a signal f(t) on the i-th dyadic scale can be written as: 1 W ff (u; 2 )g pffiffiffiffi 2i i
g f (t)c
tu 2i
dtf + c¯ 2i (u)
where c is the orthogonal mother wavelet, and: 1 t ¯ ffiffiffiffi p c2i (t) c2i (t) c 2i 2i Under certain conditions, the discrete implementations of the dyadic wavelet transform can be based on digital filter banks [17]. Let us define: di [n]W ff (n; 2i )g hf (t); c2i (tn)i ai [n]hf (t); c¯ 2i (tn)i where i E/0. For a given filter x [n ], xi [n ] denotes the filter obtained by inserting 2i /1 zeros between every filter coefficients (‘holes algorithm’), a0[n] is the original signal, and its wavelet transformed (decomposition) at the i/1-th scale is: ai1[n ]/ aj + hi [/n], and di1[n ]/dj + gj [/n ] where a defines the approximated signal, and d the signal details at the given scale; h and g are biorthogonal filters corresponding to the selected mother wavelet (Table 1). The reconstruction of the signal is done by the inverse wavelet transform, which has the following form at the i- th scale: 1 ai [n] (a 2
i1 +hi [n]di1
ˆ
+ gˆi [n])
where hˆ and gˆ are the dual filters of h and g, respectively. In case of biorthogonal filters: gˆ/ /h and hˆ/ /g. The Coiflet wavelet with four moments equal to zero (Coiflet-2) was selected as the mother wavelet.
According to the filter bank based realization, the length of the processing time window should be adjusted as Twproc /n + 8192/fs, where n is a positive integer. We selected n /4, resulting a Twproc /3 s. Thus one processing time window consisted of 32768 sample points of both input channels. The two channels were wavelet decomposed up to the 6th order (i/0, 1,. . ., 5), resulting two sets of coefficients: a6A, d6A, d5A, d4A, d3A, d2A, d1A and similarly: a6E, d6E,. . ., d1E, where the A and E letters in the indices refer to the abdominal and the external channel, respectively. The aim of the external noise cancellation is to construct a new set of coefficients: a6F, d6F,. . ., d1F, which can be reconstructed (inverse transformed) as the abdominal signal free of external noise. The algorithm of Scheme 1 was used to obtain this set of coefficients. Please note that only the detail coefficients di A of the abdominal channel are thresholded, whereby the details di E of the external channel are adoptively subtracted from di A.
3.3. Bandpass filtering The adaptive wavelet denoising operates across the whole recorded signal band (0 B/f B/fs/2), however it removes only the external noise. A bandpass digital filter was designed to limit the signal into the band: 35 Hz B/f B/200 Hz. The section of the lower band limit was based on the experimental result that a significant part of the disturbing maternal heart and digestive sounds lies
Table 1 Base coefficients of the equivalent filter bank h0
g0
/0.0007, /0.0018, 0.0056, 0.0237, /0.0594, /0.0765, 0.4170, 0.8127, 0.3861, /0.0674, /0.0415, 0.0164
/0.0164, /0.0415, 0.0674, 0.3861, /0.8127, 0.4170, 0.0765, /0.0594, /0.0237, 0.0056, 0.0018, /0.0007
287
Scheme 1. Wavelet-based inter-channel denoising.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
288
Table 2 Coefficients a and b of the highpass and lowpass filter a /b
Highpass (f /35 Hz)
Lowpass (f B/200 Hz)
a
1.0000, /7.8832, 27.1886, /53.5868, 66.011, /52.044, 25.6459, /7.2217, 0.8897 0.9432, /7.5460, 26.4109, /52.822, 66.027, /52.8219, 26.4109, /7.5460, 0.9432
1.0000, /7.4158, 24.0798, /44.7159, 51.939, /38.642, 17.979, /4.7840, 0.5573 /0, /0, /0, /0, /0, /0, /0, /0, /0
b
below this border, while the fetal heart sounds are not dominantly present there anymore. Due to stability reasons we designed two separate 8th order Butterworth-type IIR digital filters: one low-pass and one high-pass filter. The signal f[n ] is filtered as s[n ]:
s[n]
nX b 1 k1
b[k]f [nk 1]
nX a 1
a[k]s[nk 1]
k2
where a and b are the filter coefficients listed in Table 2, and na /nb /8. The effect of the external noise cancellation and bandpass filtering is illustrated in Fig. 5 using a section of the signal in Fig. 3.
3.4. Locating the amplitude bursts The Wavelet-based external noise cancellation and the bandpass filtering together produce a kind of preprocessed signal s, in which the S1 and S2 beats correspond to amplitude bursts (or simply: bursts, see Fig. 5d). Obviously the band-pass filtering cannot fully remove the internal noise. Disturbances still remain in the filtered signal, mainly caused by fetal movements, fluctuations in the fruit water and ‘shear noise’ of the measurement head. Unfortunately, this noise can cause similar bursts to that of S1 and S2 beats have (Fig. 1b). For the location of bursts several methods were proposed in the field of non-fetal phonocardiography. These methods are usually based on the
Fig. 5. The effect of external noise cancellation and bandpass filtering shown on a section of Fig. 3: (a) abdominal channel, (b) abdominal channel after simple bandpass filtering, (c) abdominal channel after external noise cancellation, (d) abdominal channel after external noise cancellation and bandpass filtering.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
thresholding of some derived signal: either the signal’s absolute value, its squared amplitude, or its Shannon-energy. A fair comparison of the various burst detection methods is provided by Liang et al. [18]. Unfortunately, these methods are only concerned with the magnitude of the derived signal and do not examine the shape and timing of the bursts. As a consequence these methods were insufficient in the case of fetal phonocardiograms where the signal-to-noise ratio is far less than in non-fetal ones. To improve burst detection, the cross-correlation with a reference burst was calculated. This method can measure all three burst parameters at once: magnitude, shape and timing. Since the exact waveform of the bursts varies considerably even during one recording period, not the preprocessed
289
signal, but its envelope was cross-correlated with the reference burst. The envelope e of the preprocessed signal s can be derived from measuring the amplitude difference between the pmin local minimums and the successive pmax local maximums of the signal. The value of each calculated amplitude difference is inserted into e until the next local minimum / maximum pair of s is reached (time spanning). Thus, the resulted envelope signal reflects the amplitude dynamic of the original signal and the same timing is kept. The envelope algorithm is shown in the flow chart in Fig. 6. The section of a preprocessed signal and its envelope is illustrated in Fig. 8a and b, respectively. In the next step the envelope signal was downsampled with the factor 1:10 which allowed faster processing time in the further steps. The downsampling was realized by using a simple moving average filter. Finally, the cross-correlation sum r measured the similarity of the downsampled envelope signal e? to a reference burst envelope y of length k and l, respectively: r[n]
l X
e?[mn]y[m];
n 0; 1; . . . ; k
m0
The cross-correlation of the downsampled envelope signal with the selected reference burst envelope (Fig. 7) produced the signal in Fig. 8c. The peak values in r show the locations of the acoustic amplitude bursts and hence the possible location of the S1 and S2 beats. The peaks were detected with an algorithm that is an extended and modified version of the method used in Ref. [13]. The main steps are:
Fig. 6. Flow chart of the envelope algorithm.
Fig. 7. The reference burst envelope.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
290
t[i]
t[i]m[i] t[i 1]m[i 1] m[i] m[i 1]
. the i/1-th element is deleted and n is adjusted. . eight peaks with the largest values were selected now again from a, and their mean value l is calculated. . a three-state possibility vector p is obtained: if a [i]B/l/16 then p [i]/0, else if a[i] /l /2 then p [i]/1, else p [i]/0.5, where i/1, 2,..., n . return as result those elements of m , t and p , where p[i] /0, i /1, 2,. . ., n . 3.5. Locating the S1 and S2 beats Fig. 8. The internal processing steps: (a) preprocessed signal, (b) envelope signal, (c) cross-correlated signal.
. all local maximums (peaks) of r are located, resulting in three n -dimensional vectors: a , m and t, where a contains the height of the peaks relative to the preceding minimums, m contains the magnitudes of the peaks, and t their locations relative to the start of the recording, respectively. . those jth elements of each a, m and t are deleted where: t[i/1]/t [i ]B/10 ms, i/1, 2,. . ., n/ 1, and if a [i] B/a [i/1] then j/i, else j /i/1. After this n is adjusted. . eight peaks with the largest values were selected from a , and their mean value l is calculated. . those elements of a, m and t are deleted where: a [i] B/l /8, i/1, 2,. . ., n/1. After this n is adjusted. . if t[i/1]/t [i ]B/40 ms, where i/1, 2,. . ., n/ 1, then the i-th element is merged with the i/1-th element according to: a[i]
The resulted m , l and p vectors describe the magnitude, time location and probability of the detected bursts, respectively. Since several bursts can be caused by noise, the next task is to select only the bursts of S1 and S2 beats. The proposed algorithm is depicted in the state machine of Fig. 9. The individual S1 and S2 bursts can be distinguished from the noise bursts by taking various physiological timings into a consideration (Table 3). The result of the detection is a beat information matrix M, where cash row corresponds to an S1 / S2 beat pair: TS1, TS2, PS1, PS2, where TS1, TS2 are the time locations of the S1 and S2 beats relative to the start of recording and PS1, PS2, are their three-state probabilities respectively. The algorithm was constructed with the assumption that the deletion of the acoustically stronger S1 beats is more probable. Therefore the method is
a[i] a[i 1] 2
m[i]
m[i] m[i 1] 2
Fig. 9. State machine of the S1 /S2 detection.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
Table 3 Meaning of the timing values Value
Meaning
TS1S1 TS1S2 TS1S1min TS1S1max TS1S2min TS1S2max TS2S1min TS1S1dev TS1S2dev
Time elapsed between two consecutive S1 beats Time elapsed between an S1 /S2 beat pair Minimal allowed time between two consecutive S1s Maximal allowed time between two consecutive S1s Minimal allowed time between an S1 /S2 pair Maximal allowed time between an S1 /S2 pair Minimal time to be elapsed between an S2 and an S1 Allowed beat-to-beat deviation of TS1S1 Allowed beat-to-beat deviation of TS1S2
291
First S1 found . Appending beat information: h /i (position of first S1 found) if S2 could also be located (p[g ]/0), then append row: t[h]t[g]p [h ]p[g ], and i/g/1, TS1S2 /t[g ]/ t [h ] else append row: t[h] 0 p[h ] 0, and i/i/1, TS1S2 /0 do TR3 Next search
based on the search for S1 of beats. The detection of an S2 beat only upgrade, the reliability of the S1 detection, however, it is possible that no S2 after a given S1 can be detected (i.e. TS2 /PS2 /0). All timing values are adoptively tuned during the detection. If within the given timing no subsequent S1 beat can be found, the algorithm reinitializes itself by setting the default timing values. The following part of this section details the states and the state transitions of the detection algorithm: Start . General initialization: TS1S2dev /40 ms, TS1S2dev /25 ms, TS2S1min / 150 ms, i/1, do TR1 Initial search . Search initialization: TS1S1min /300 ms, TS1S1max /750 ms, TS1S2min /100 ms, TS1S2max /210 ms, s /0 . Checking if the current burst belongs to an S1 beat : start j from i/1, and search for an S2 beat in a range, where TS1S2min B/t [j ]/t [i ]B/TS1S2max and let g the location in this range where p [j] /0 and m [j]B/m [i] if p [i ]/p [g] E/1 then do TR2 else i/i/1 and goto start of this state
. Checking the index: if i /n do TR9 . Checking TS1S1 : if s /1 do TR7 else: TS1S1 /t[i ]/t [h ] (time elasped since last S1 found) if TS1S1 B/TS1S1min or TS1S1 B/TS1S2/ TS2S1min then i/i/1, s/s/1, goto start of this state else if TS1S1 /TS1S1max do TR7, else: . Update TS1S2: if TS1S2 /0 then TS1S2min /TS1S2/TS1S2dev, TS1S2max /TS1S2/ TS1S2dev . Search for the next possible S1 : start j from i/1 and search for an S2 beat in a range, where TS1S2min B/t[j]/t[i]B/TS1S2max and let g the location in this range where p [j]/0 and m [j ]B/m [i ] if p[i]/p [g ]E/1 then do TR4 else i/i/1 and goto start of this state Next S1 found . Checking the index: if i /n do TR8 . Checking TS1S1 : TS1S1 /t[i ]/t [h ] (time elapsed since last S1 found) if TS1S1 B/TS1S1min or TS1S1 /TS1S1max do TR6, else: . Update TS1S1: TS1S1min /TS1S1/TS1S1dev, TS1S1max /TS1S1/ TS1S1dev . Appending beat information : h /i (position of S1 found) if S2 could also be located (p[g ]/0), then
292
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
append row: t[h ]t [g ]p[h ]p [g], and i/g/1, TS1S2 /t[g ]/ t[h ] else append row: t[h ] 0 p [h] 0, and i/i/1, TS1S2 /0 do TR5 In order to keep the notations simple and easyto-understand, the following two features of the real algorithm were not described above: . The search for S1/S2 pairs (in states ‘Initial search’ and ‘Next search’) is continued until t[j]/t[i]B/TS1S1max. If more than one possible S1 locations were found, then the one which has the largest magnitude (m ) is selected. . The TS1S1 and TS1S2 values are calculated using the average of the four subsequently preceding TS1S1 and TS1S2 values, respectively. If no four subsequent values are available (eg. after re-
Fig. 10. Result of the S1 /S2 detection and the calculated FHR values using the signal excerpt in Fig. 3.
initalization), the last single value is used, like in this description above. Fig. 10 illustrates the result of the S1 /S2 pair detection. The top part of this figure shows the cross-correlated envelope signal of the excerpt in Fig. 3. The small dot, on the signal peaks displaying the detected amplitude bursts. The dots in the small squares and circles correspond to the detected S1 and S2 beats, respectively. The left table in the bottom part of the Fig. 10 displays the resulted beat information matrix M . 3.6. Calculation of FHR The last element of the signal processing is the FHR determination. A given FHR value is calculated based on the TS1S1 values inside a time window with a given length (Twfhr /3000 ms). Unfortunately, since some of the S1 beats cannot be obtained, the simple averaging of the TS1S1 times in the selected window is not sufficient. However, we can assume that all beats included in the beat information matrix are valid beats. The proposed FHR calculation method uses only the TS1S1 values and assigns one FHR value to each of them. The method is based on the following main steps: . i indexes the rows of the beat info matrix (i/1, 2,. . ., n ). . k subsequent beats are selected using the time window: TS1S1[i], TS1S1[i/1], . . ., TS1S1[i/k ], where: TS1S1[k ]B/TS1S1[i]/Twfhr . Using the selected k beats an d distance vector is calculated (time elapsed between the beats): d [j]/TS1S1[j/1]/TS1S1[j], where j/1, 2,. . ., k/1 . a validity factor is calculated; v /v1/ v2 where v1 is the maximum number of elements in d , having no more deviation than TS1S1dev and v2 is the number of the neighboring (subsequent) beats amongst them. . if v /3 then the FHR[i] is calculated by averaging the elements of d determined by v1, else no valid FHR could be detected inside Twfhr and FHR[i] is assigned to the last valid FHR value.
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
. i/i/1, and goto to the second step. If no valid FHR can be detected over a period of Twfhr/2, then FHR[i ]/0. In the case of the signal in Fig. 10, the resulted FHR values and the pulse values (Pulse[i]/ 60 000/FHR[i]) are displayed in the bottom right table, in ms and bpm units, respectively. Since one FHR value is assigned to each of the S1 beats, the time proportional graphical representations of the FHR values in x, y coordinates are: TS1S1[i], FHR[i].
4. Results The phonocardiographic device was designed and a prototype device was used for the recording. The signal processing was developed in Matlab and after the successful evaluation it was implemented in C//, under Windows NT (Fig. 11). The developed software has the following key features: . real-time signal processing, . input signals either from a wave file or directly from the line level audio input,
293
. numeric display of the actual FHR (in bpm units), . storing and printing of beat information and of the FHR. The evaluation of the developed signal processing method was based on the reference CTG paper records. The paper records were scanned and divided by a grid with a resolution of 3 bpm along the FHR axis, and 5 s along at the time axis. The FHR value of each grid area was compared to the average phonocardiographic FHR value in the corresponding interval. Those positions were not evaluated where the reference CTG record did not include any estimable FHR value. Two of the records yielded more than 90% accuracy, three records were between 80/90%, another two records were within the 70 /80% limit, while the remaining two (rather poor quality) records were below the 70% accuracy. The overall accuracy of the method was: 83%. Fig. 12 shows two FHR sections, both obtained by CTG and by the developed phonocardiographic software. The upper part of the figure shows a steady (possibly sleeping) fetus, while the bottom part illustrates the case of an active and strongly moving fetus. Please note that the CTG has also failed to track the FHR during the strong fetal movements.
5. Discussion and conclusions
Fig. 11. The real-time phonocardiographic software (1 */the calculated FHR diagram, 2/3 */buttons for start/stop processing, 4 */button to open the settings dialog to alter parameters, 5 */exit button, 6 */the actual FHR value, 7 */signal source selection (audio line in or file).
The non-invasive routine measurement of FHR is of great clinical relevance. This paper presented an advanced phonocardiographic method for the determination of the FHR. The phonocardiographic recording technology offers a passive and therefore absolutely non-invasive diagnostic measurement. Fetal phonocardiography has been since long researched, however there is no fetal phonocardiographic device available yet. Existing methods use single channel recording with specialized pressure [11,12] or acoustic sensors [13]. Unfortunately, only the usage of a proper measurement technology itself is not sufficient since the various noise factors make the exact detection of the acoustic
294
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
Fig. 12. FHR diagrams of a steady and a moving fetus obtained by CTG and the phonocardiographic method (vertical scaling is in bpm, horizontal unit is 1 min).
fetal heart beats very difficult. This problem can be handled only by using sophisticated signal processing methods. This paper presented an advanced phonocardiographic system which put the main emphasis not on the acoustic sensors, but on the signal processing issues. The presented phonocardiographic device is built up from commercially available components (high sensitivity electret microphones, low-noise precision instrumentation amplifiers). An innovative solution is the two-channel recording which provided the background for the efficient external noise cancellation.
The determination of the FHR is provided by an advanced real-time signal processing method which consists of several new solutions. The wavelet-based denoising uses an adaptive crosschannel coefficient adjustment in order to cancel external noise more efficiently than any other examined methods. The cross-correlation based amplitude burst detection allowed the precise localization of those bursts, which may correspond with the highest probability to an S1 or S2 beat. The advanced rule-based selection of the S1 and S2 beats by the finite state machine and physiological timing rules is a key point of the entire system. The
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
determination of the periodicity based on the confidence analysis provided an accurate and robust FHR measurement which is comparable to ultrasound CTG. The presented system offers an optimal trade-off between complexity and feasibility and is viable in low-cost, standalone, battery-powered FHR monitoring devices which could be available also for home care use. Although the introduced method produced encouraging results, it must be underlined that a comprehensive evaluation base on a huge number of clinical validation measurements is still needed
6. Future work Our future work will focus on directions of further development The measurement device will be rebuilt including automatic gain, control unit in order in avoid manual gain adjustment and to allow optimal signal quality. Also the signal processing methods can be further enhanced. The selection of the reference burst envelope used for the cross-correlation (Fig. 7) needs a comprehensive validation which should be based on a broad clinical study. The key of the entire signal processing is the proper selection of the S1 and S2 beats. Thus, instead of the finite state machine we intend to use fuzzy or neuro-fuzzy inference system for the selection of the S1 and S2 bursts. The current system presents the FHR diagram without any diagnostic assessment. Although the analysis of the FHR is a complex problem, we would like to develop an automated assessment method which could be useful especially in home care applications. One possible approach is the modelling of the underlying physiological phenomena based on the FHR data [19,20], since the diagnostic information could be obtained from the estimated parameters.
Acknowledgements The phonocardiographic records were taken partly at the Universitats frauenklinik Erlangen,
295
Germany and at the Hospital Universitaire Taher Sfar, Mahdia, Tunesia. The author of this paper would like to thank the supporting help of Andreas Hammel, MD in Erlangen, Mohamed Bannour, MD in Mahdia and Lanouar M. Chouk in Erlangen at taking the records.
References [1] G.S. Dawes, C.R. Houghton, C.W. Redman, G.H. Visser, Pattern of the normal human fetal heart rate, Br. J. Obstet. Gynacol. 89 (1982) (1982) 276 /284. [2] G.S. Dawes, M. Moulden, C.W. Redman, Improvements in computerized fetal heart rate analysis antepartum, J. Perinatal Med. 24 (1996) (1996) 25 /36. [3] Y. Tal, S. Akselrod, Feral heart rate detection by a special transformation method, Proc. IFEE Comp. Soc. Conf. Computers in Cardiology (1990) 275 /278. [4] K.A. Azad, Z.M. Darus, M.A. Ali, Development of a fuzzy rule-based QRS detection algorithm for fetal and maternal heart rate monitoring, Proc. Int. Conf. IEEFE, EMBS 20 (1998) 170-171. [5] J.C. Echeverria, R. Ortiz, N. Ramirez, V. Medina, R. Gonzalez, A reliable method for abdominal ECG signal processing, Computers in Cardiology 25 (1998) 529 /532. [6] C. Lin, H. Wu, T. Liu, N. Lee, T. Kou, S. Young, A portable monitor for fetal heart rate and uterine contraction, IEEE EMBS Magazine (1997 Nov/Dec) 80 /84. [7] J.M. Bentley, P.M. Grant, J.T.E. McDonnel, Time-frequency and time-scale techniques for the classification of native and bioprosthetic heart valve sounds, IEEE Trans. Biomed. Eng. 45 (1998) 125 /128. [8] I. Cathers, Neural network assisted auscultation, Artificial Intelligence in Medicine 7 (1995) 53 /66. [9] J. Ritola, S. Lukkarinen, Comparison of time-frequency distributions in the heart sounds analysis, Medical and Biological Engineering and Computing 34 (1996) 89 /90. [10] X. Zhang, L.G. Durand, L. Senhadji, H.C. Lee, J.L. Coatrieux, Time-frequency scaling transformation of the phonocardiogram based on the matching pursuit method, IEEE Trans. Biomed. Eng. 45 (1998) 972 /979. [11] A.J. Zuckerwar, R.A. Pretlow, J.W. Stoughton, D.A. Baker, Development of a piezopolimer pressure sensor for fetal heart rate monitor, IEEE Trans. Biomed. Eng. 40 (1993) 963 /969. [12] D.G. Talbert, W.L. Davies, F. Johnson, N. Abraham, N. Colley, D.P. Southall, Wide bandwidth fetal phonocardiography using a sensor matched to the compliance of the mother’s abdominal wall, IEEE Trans. Biomed. Eng. 33 (1986) 175 /181. [13] F. Kova´cs, M. To¨ro¨k, I. Habermajer, A rule-based phonocardiographic method for long-term fetal heart rate monitoring, IEEE Trans. Biomed. Eng. 47 (2000) 124 /130.
296
P. Va´rady et al. / Computer Methods and Programs in Biomedicine 71 (2003) 283 /296
[14] P. Va´rady, I. Gross, M.L. Chouk, A. Hein, Analysis of the fetal heart activity by the means of phonocardiography. Proc. IFAC Conference Telematics and Application, Weingarten, Germany, (2001) 41 /46. [15] J.F. Cardoso, A. Soulouminac, Blind beamforming for non-Gaussian signals, IEE-Proc. F 140 (1993) 362 /370. [16] C. Hutten, J. Herault, Blind separation of source, part I: an adaptive algorithm based on neuromimetic architecture, Signal Processing 24 (1991) 1 /10.
[17] G. Strang, T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, 1996. [18] H. Liang, S. Lukkarinen, I. Hartimo, Heart sound segmentation algorithm based on heart sound envelogram, Computers in Cardiology 24 (1997). [19] E.R. Carson, C. Cobelli, Modelling Methodology for Physiology and Medicine, Academic Press, San Diego, 2000. [20] M. Di Rienzo, G. Mancia, Computer Analysis of Cardiovascular Signals, IOS Press, Amsterdam, 1995.