The algorithm to be presented is based on ideas proposed by Johansen for the detection of cavitation in mechanical heart valve patients [18, 21]. In what follows, we set this algorithm on a firm mathematical foundation via a rigorous analysis and we show that theoretically, it should quantify the quantities we seek to capture. The algorithm is then implemented on real heart-beat acoustic signals obtained via a stethoscope or high frequency hydrophone with healthy subjects. We seek to evaluate the algorithm with real-life signals and their inherent difficulties, rather than 'fake' test signals. In particular, slight changes in the length of each heart-beat and typical levels of noise could pose difficulties. Low-pass filtering for noise removal is not performed as it would also likely remove the signal-of-interest (usually a high frequency component). Thus the algorithm will have to be able to deal with higher levels of noise than might be desirable.
By using signals from healthy subjects, no cavitation or abnormalities are expected to occur in the signal. Since no cavitation or other non-deterministic features of interest are expected in the signal then a quantitative measure of the "goodness" of the algorithm is lower calculated levels of "potential features of interest". A change to the algorithm that lowers the levels of non-deterministic markers is considered to be improving the algorithm. This is an important point to make as this allows us to work with real heart signals.
For implementation of the algorithm, both stethoscope and hydrophone signals were used. Some of the stethoscope signals that were used in the algorithms were pre-recorded stethoscope test signals that were obtained from a CD provided with the Littmann stethoscope. Other signals that were used in the algorithms are in-vivo stethoscope signals recorded in our lab on ourselves, representing healthy subjects. Those recordings were made with the same electronic stethoscope, sound card, and laptop as described below. The acquisition of these signals was in compliance with the University of Ottawa Ethics Guidelines.
For further analysis with other types of signals, in-vitro recordings were also obtained using a left-heart simulator in a laboratory at the University of Ottawa Heart Institute [23]. The left-heart simulator is a mechanical system that simulates the activity in the left portion of the heart and is designed for in vitro testing of bioprosthetic and mechanical heart valves. The left heart simulator used in this experiment was purchased from ViVitro Labs Inc. (Victoria, Canada). Designed to be physiologically realistic, it can assess the function of heart valves and other devices under simulated cardiac conditions. The simulator consists of a Superpump system, a viscoelastic impedance adapter, a left heart model, flow and pressure measuring systems, a waveform generator, and a PC data acquisition system. The valves used in this experiment were a bioprosthetic valve of type Magna from Edwards Lifescience (Irvine, California, USA), and a St. Jude Medical bileaflet mechanical heart valve. The bioprosthetic valve was an aortic trileaflet bovine pericardial valve with a commercial denomination of 27 mm. The mechanical heart valve had pyrolytic carbon occluders. Its interior diameter was 18 mm and its exterior diameter including the cuff was 28 mm. Its commercial denomination was 23 mm. Cavitation was not visually observed with the use of either valve. Full details are given in [24].
The high-frequency pressure fluctuations were also measured using a miniature hydrophone. The hydrophone used is the 8103 miniature hydrophone by Brüel & Kjaer, Naerum, Denmark. Its frequency range is 0.1 Hz to 180 kHz. The hydrophone was connected to the Nexus Conditioning Amplifier by Brüel & Kjaer. To help segment the hydrophone signal, an electronic stethoscope was positioned beside the hydrophone, and the stethoscope recording was done simultaneously with the hydrophone recording. The electronic stethoscope used was the Welch Allyn Elite Electronic Stethoscope, New York, USA. Its frequency range is 20 Hz to 20 kHz. The stethoscope and hydrophone were connected to a Y-adapter, which was connected to the sound card. The Y-adapter used was a dual mono jack to stereo plug adapter. The sound card used in this research was the Creative Sound Blaster Audigy 2 ZS Notebook sound card. It was capable of recording with a sampling rate of up to 96 kHz. To maintain a manageable data set, a sampling rate of 44.1 kHz was used for healthy subjects since no cavitation was expected. The sampling rate could be changed to 96 kHz when cavitation was expected. Finally, the sound card was plugged into the laptop (ASUS A3E) where the data were stored. The data was recorded with software called Creative Smart Recorder which came with the sound card. Figure 1 illustrates the equipment used and the experimental setup.
2.1. The Algorithm
2.1.1. Rationale
It has been suggested in [18] and [21] by Johansen et al. that the cavitation near mechanical heart valves can be quantified by separating the acoustic pressure signal into deterministic and non-deterministic components. The deterministic component represents the valve closing sound assumed deterministic since valve closure is cyclic. The non-deterministic component is the information of interest. This was assumed to originate from cavitation but may stem from other 'signals of interest' such as heart-murmurs, thus widening the applicability of this algorithm. The non-deterministic component of the signal also contains unwanted random noise from various sources such as the detection equipment used.
Based on this, the algorithm follows on the assumption that the information of interest is quantified by the level of non-deterministic energy, defined as the difference between the deterministic energy and the total energy:
(1)
Here E
non-det represents the non-deterministic signal energy which is the energy of interest, attributed to the random, non-deterministic events of interest. E
total represents the total signal energy containing all energies and E
det represents the deterministic signal energy of the repeating components of the signal that we seek to eliminate, the heart beat itself. How these energies are defined mathematically will be presented subsequently. Since the quantity of interest is a often a high-frequency component [16, 18, 21], no low-pass filtering is applied to the original signal.
2.1.2. The algorithm
The first step in the algorithm is to segment the total signal into individual heartbeats and then line up all the heart beats in time. Then, each heart beat is truncated so that all beats have the same length since vectors of different lengths cannot be averaged. The end part of the beats is truncated since no important information is located there. The truncated heart beats are then superimposed one over the other and ensemble averaged to obtain an average heart beat signal. This theoretically eliminates unwanted noise as well as reduces the signal parts that do not repeat from beat to beat.
The next step in this algorithm is then to find the deterministic energy. This is defined as
(2)
where N is the number of samples in each heart beat, p
ea
[n] is the ensemble average of the heart beats, and F is the Fourier Transform. The ensemble average is the average heartbeat calculated according to
(3)
where HB is the number of heart beats measured, and p
xy
[n, m] represents the nth sample of the mth heart beat in the total signal.
The third step in this algorithm consists of determining the total energy. The energy is determined for each segmented and truncated heart beat and then these energies are averaged, which represents the total energy. Mathematically,
(4)
where |F(p
xy
[n, i])|2is the amplitude spectrum squared of the ith heart beat.
The final step consists of subtracting the deterministic energy from the total energy to obtain the non-deterministic energy. Figure 2 is a block diagram summarizing the algorithm.
2.1.3. Segmentation algorithm
The first step of the algorithm is to segment the signal into individual heartbeats. The segmentation problem is a topic of research in itself and we do not attempt a literature review of the subject here. Many algorithms have been proposed in this area, for example [25–36] and any suitable segmentation algorithm can be chosen for this step. The segmentation algorithm used in this paper was the method developed in [37]. This algorithm was used instead the cross-correlation method proposed by Johansen since the results obtained by Johansen [18] could not be reproduced.
The Courtemanche segmentation algorithm used herein proposed the use of an adaptive threshold wavelet transform filtering technique used with Shannon energy, physiological factors and heart rate approximation to properly identify the first heart sounds (S1) and segment the heart signal. However, this method can still present some errors when faced with complex signals. Therefore, the addition of a Mel-Scaled Wavelet Transform (MSWT) validation step was proposed. The MSWT is a modified Mel-Frequency Cepstral Coefficient (MFCC) algorithm with the Discrete Wavelet Transform (DWT), and it was used to reduce the impact of noise on the coefficients. The preliminary results obtained in the paper indicated that the MSWT is less prone to noise than the MFCC and can distinguish S1 sounds from others when faced with complex signals [37].
The segmentation points obtained using the Courtemanche segmentation algorithm [37] were compared with those obtained using a different segmentation algorithm as suggested in [38]. The stethoscope signals were independently segmented by the first author of [38], and it was found that the segmentation points obtained were very similar to the ones obtained using the Courtemanche algorithm of [37]. This served to verify the calculation of the segmentation points.
2.2. Mathematical analysis of the algorithm
To ensure that the algorithm should indeed quantify the non-deterministic signals of interest, a rigorous mathematical analysis is presented. To simplify the mathematical analysis, Parseval's relation is used so that the energy of a signal can be equivalently found from the time or frequency domain.
2.2.1. Calculations with a continuous signal x(t)
Deterministic energy
Analytically, the segmentation and truncation of the original signal is accomplished by multiplying the continuous signal x(t) by a rectangular window w
n
(t) at time t
n
and having a length of one truncated heart beat T. The truncated heart beats are then shifted, superimposed and then averaged to obtain an average heart beat signal given by
(5)
where each time-shifted heart beat y
n
(t) = y(t + t
n
) is found by shifting each heart beat back to time zero, after multiplying the original signal with the window y(t) = x(t)·w
n
(t) to obtain each beat. The deterministic energy is the energy of the ensemble-averaged heart beat:
(6)
These steps can be illustrated with n heartbeats but for brevity, the calculation is illustrated with three heart beats. For three heart beats, we have that
(7)
Hence the deterministic energy is given by
(8)
If y
1(t) = y
2(t) = y
3(t), implying that the three heart beats are identical and perfectly superimposed, then
(9)
This pattern is the same for n heartbeats. As observed in (8) compared to (9), if the heart beats are not perfectly superimposed, cross-terms appear in the deterministic energy result, which impacts both the deterministic and thus non-deterministic energy results. This demonstrates the importance of the correct segmentation of the original signal. Imperfect segmentation leads to additional terms in the deterministic energy calculation, and the source of these terms is entirely from the imperfect segmentation and not noise or cavitation. This would, in turn, affect the non-deterministic energy not because of any true additional non-deterministic energy in the signal but rather through imperfect processing of the signal.
Total energy
The total energy is defined as the average of individual heart beat energies:
(10)
To demonstrate for three heart beats, the total energy is given by
(11)
Then if the heartbeats are identical and perfectly superimposed so that y
1(t) = y
2(t) = y
3(t), it follows that
(12)
Similar results are obtained with n beats.
Non-deterministic energy
The non-deterministic energy is obtained by subtracting the deterministic energy from the total energy. Continuing our example for three heart beats, the non-deterministic energy is
(13)
For beats that are identical, then E
non
-det = 0. This clearly demonstrates that the proposed algorithm, so far in theory, works perfectly. Namely, if a signal consists of heart beats that repeat perfectly and can be segmented perfectly, then the non-deterministic energy should be zero.
Suppose that an additional component is added to the signal so that the beats are not identical but are 'similar'. This is done by considering one beat to be fundamentally the same as another but with the addition of an extra signal, which may be cavitation or noise, etc. This is modelled mathematically as
(14)
where η
i
(t)represent cavitation, noise or some other signal of interest. Then the non-deterministic energy calculation gives
(15)
The previous equation again demonstrates analytically that the proposed algorithm is a quantitative measure that contains only the non-repeating elements of a signal since the repeating portions of the signal do not appear in the expression for non-deterministic energy. Therefore, the theory behind the algorithm is sound.
2.2.2. Effect of beat misalignment
In the case when the heart beats are not perfectly superimposed, the non-deterministic energy result is affected. The impact is demonstrated analytically with two heart beats. We start with the non-deterministic energy but now consider the case where ε represents the factor by which the second heart beat is identical to the first but misaligned with respect to it so that y
2(t) = y
1(t + ε). This gives
(16)
A signal of interest (such as cavitation) η(t) is then added to the heart signal so that . Thus equation (16) becomes
(17)
From the non-deterministic energy result in (17), one can observe that if ε = 0, then
(18)
which implies that the non-deterministic energy represents only the non-repeating portions of the signal, as expected and previously shown. However, if ε ≠ 0, then ε appears in the non-deterministic energy result as shown in (17), meaning that the non-deterministic energy is not a measure of η(t)alone. Therefore, if the beats are not properly superimposed, those beats that are not lined up will have an impact on the non-deterministic result, which might be falsely interpreted as a greater level of the non-repeating portion in the signal.
2.2.3. Calculations with a test sine signal
To calculate the relative sizes of contributions to the energies from the heartbeat portion compared with the non-deterministic portion, test sine signals were used since they are periodic and allow for simple calculations in closed form. A test signal was constructed consisting of a simple low-frequency sine signal, which represents the heartbeat, along with a higher frequency sine signal representing an extremely simple cavitation signal.
(19)
Here, A
i
is the amplitude of the i th heart beat signal, C
i
is the amplitude of the i th cavitation signal, Ω is the frequency of the heart signal, ω
1 and ω
2 are the (higher) frequencies of the cavitation signal, and φ
i
is the phase shift for each cavitation signal. A phase shift is required to ensure the cavitation signal will not be in-phase with the heartbeat. Using the parameters, , ω
1 = 35000Ω Hz, ω
2 = 36000Ω Hz so that the cavitation has a much higher frequency than the heartbeat, and also A
1 = A
2 = A, C
1 = C
2 = C, φ
1 = φ
2 = 0, y
1(t) ≠ y 2(t), the deterministic energy is found to be
(20)
The value of ω
1 was chosen as 35 kHz times the value of Ω since the high frequency pressure fluctuations due to cavitation occur in the frequency range of 35 kHz to 350 kHz [39]. The frequency ω
2 was assigned a different value than ω
1 to account for possible variations in frequency between heart beats. The value of Ω is approximately 1 Hz since the duration of a heart beat is approximately 1 second. If the value of ω
1 and ω
2 is changed to any integer between the frequency range mentioned previously, the deterministic energy equation in (20) will remain the same. From equation (20), it is noted that the amplitudes squared of both the heart-beat and cavitation signals appear, with the heart-beat itself having a stronger presence due to the factor 2.For the chosen signals, the total energy can be shown to be given by
(21)
The non-deterministic energy is obtained by subtracting the deterministic energy from the total energy:
(22)
It is observed in (22) that E
non-det does not depend on A, the amplitude of the heart signal, although both the deterministic and total energies do. It depends on C, the amplitude of the cavitation signal. Therefore, for perfectly aligned heart beats, the non-deterministic energy represents the cavitation in the signal confirming the algorithm.
2.2.4. Superposition problem
When heart beats are not perfectly superimposed, there are two different possibilities after truncation:
The first possibility is that y
2(t), the second heart beat, will be zero for the first few points since it is slightly shifted to the right due to its misalignment. This will resemble a misplaced heart beat since the value before the first heart sound in a heartbeat should be near zero. The second possibility is that the beginning of y
2(t) will be the part of the beat that was cut off at the end by the truncation, since it is slightly shifted to the right. Some calculations are demonstrated below to compare the results obtained with the two possibilities mentioned previously. The calculations are done using two heart beats.
1st case
The test signal used is identical to the one used in (19).
Now, the second heart beat is "misaligned" so that y
2(t) = y
1(t - ε)u(t - ε), where u(t - ε)is the unit step function and the heart beats are truncated to have length T.
With the chosen form of the test signal in (19), and the parameters ω
1 = 35000 Hz, Ω = 1, φ
1 = φ
2 = 0, ε = 10ms, this gives
(23)
Note that with misalignment of the beats explicitly modelled as a shift in one of the heart beats, the non-deterministic energy now depends on the magnitude of the cavitation signal as well as the magnitude of the heart beats themselves and the dependence on both is about the same, with both factors having a π/4 dependence. When the signals were perfectly aligned, then using the same parameters, equation (22) indicates that the nondeterministic energy was whereas the nondeterministic energy with misaligned beats is . The contribution of the cavitation portion of the signal has effectively halved and the heart-beat portion now has an equally weighted contribution. This is potentially particularly problematic as in general the heart beat is louder than any 'signal of interest', implying that A will be much larger than C in the previous equation. Thus, in this case the level of non-deterministic energy will rise - and not due to an increase in non-heartbeat signals but due to a misalignment of the beats.
2nd case
If we assume that y
2(t) = y
1(t-ε), then using the same form of signals and same parameters as the previous calculation, the non-deterministic energy becomes
(24)
Once again, with the effect of the misalignment of the beats accounted for, it can be seen that the non-deterministic energy depends on the magnitudes of both the cavitation and heart beat signals. The result obtained with the first possibility in (23) seems better than the result obtained with the second possibility in (24) since A, the amplitude of the heart beat signal, contributes less to the non-deterministic energy (0.7854A2 vs. 2.8888A2). In equation (23), the weighting of the contribution of A2 and C2 to the nondeterministic energy is equal whereas in equation (24), the contribution of A2 is about twice as much as the contribution of C2. Ideally, the component coming from the heart beat should not be present in the non-deterministic energy since the latter is supposed to represent the non-heartbeat portion of the signal only.
Since the quality of the segmentation of the heart signal has a large impact on how well the heart beats are superimposed and thus the calculation of energies of interest, we propose to implement an algorithm that aligns the heart beats in the signal after the initial segmentation has been performed.
To demonstrate the need for further beat alignment, even after the initial segmentation is performed; Figure 3 shows the superimposed heart beats of the pre-recorded stethoscope test signal (Pre-recorded PCG 1) after the initial segmentation was performed. The first large peak of activity in the heart beats of Figure 3 is the first heart sounds (S1) and the second peak of activity is the second heart sounds (S2). The heart beats in Figure 3 are not perfectly lined up, which is not surprising since perfect alignment is practically impossible as the heart does not produce a perfectly periodic heart signal. Therefore, the signals were processed so as to align all the S1's in the heart signal after all the heart beats have been superimposed. To do this, the maximum peak of the first quarter of each heart beat, which is the peak of S1, was determined. The S1 is usually contained in the first quarter of the heart beat. The algorithm searches for the maximum peak in the first quarter of the heart beat instead of in the entire heart beat since the S2 peak can sometimes be higher than the S1 peak, in which case Matlab would detect the S2 peak as the maximum peak instead of S1. After the S1 peaks were found for each individual heart beat, each heart beat was shifted to the left or to the right to line up with the S1 of the first heart beat of the heart signal. Once all the S1's were lined up, the beginning or the end of the heart beats were truncated so that all the beats have the same length. The larger the shift, the more the heart beat is going to be truncated. The beats that require a very large shift compared to the other beats may eventually removed from the heart signal since it is considered a 'bad' heart beat. Too much of that heart beat would be truncated after shifting; hence, too much information would be lost. The more signals are lined up, the higher the ensemble average will be.