A system for biomedical audio signal processing based on high performance computing techniques

2.1Signal model

The problem considered in this study implies the separation of the lung and heart sounds from the auscultation signal recorded with a microphone placed on the chest wall of the subject. Thus, the captured signal y⁢(m) can be formulated as follows:

where the time-domain sample index is denoted by m, xL⁢(m) and xH⁢(m) are the lung and heart signals, respectively.

Focusing on Eq. (1), the short-time Fourier transform (STFT) of y⁢(m) can be expressed as

where Y⁢(f,t), XH⁢(f,t) and XL⁢(f,t) represent the magnitude spectrograms of y⁢(m), xH⁢(m), and xL⁢(m), respectively. Here, t∈[1,T] denotes the time frame indices and f∈[1,F] represents the frequency bin.Collecting T time frames and F frequency bins, the magnitude spectrogram matrices 𝐘∈ℝ+F×T, 𝐗H∈ℝ+F×T and 𝐗L∈ℝ+F×T are defined, where𝐘=[𝐲⁢(1),…,𝐲⁢(t),…,𝐲⁢(T)] and 𝐲(t)=[Y(1,t),…,Y(f,t),…,Y(F,t)]T, 𝐱H⁢(t) and 𝐱L⁢(t) are defined similarly to 𝐲⁢(t). Finally, 𝐗H and 𝐗L were derived in a manner similar to that of 𝐘.

2.2Baseline method for heart and lung sound separation by NMF

To process heart and lung signals, it is necessary to acoustically cancel the interference between the lungs and heart. Note that both sources are simultaneously active in both time and frequency domains. In particular, the main problem generated by this interference is that heart sounds mask lung sounds and vice versa.

In [28], the separation of sounds generated by the lungs from sounds generated by the heart was performed using an NMF approach. Thus, the separation problem relies on the following factorization model:

where 𝐘^∈ℝ+F×T, 𝐗^H∈ℝ+F×T and 𝐗^L∈ℝ+F×T are the estimated magnitude spectrogram matrices of the mixture, the heart and the lung signals. 𝐇H∈ℝ+KH×T and 𝐖H∈ℝ+F×KH are the activations and bases matrices associated with the heart sounds, respectively, and 𝐇L∈ℝ+KL×T and 𝐖L∈ℝ+F×KL are the activations and bases matrices associated with the lung sounds. Note that K represents the number of NMF bases with K=KH+KL.

The authors proposed dealing with the factorization in Eq. (3) by minimizing the Kullback-Leibler divergence [37] and the sparsity constraint [38]. Therefore, the update rules were computed as follows:

where 𝟏∈ℕF×T represents an all-ones matrix composed of T columns and F rows, 𝐖 and 𝐇 are initialized as random positive matrices and ⊙ denotes the element-wise product.

Canadas-Quesada et al. [28] proposed a clustering process to classify the activations 𝐇 and bases 𝐖 calculated using the NMF procedure. This clustering takes advantage of both the temporal and spectral characteristics of respiratory and cardiac sound sources to distinguish between them. Three methods were used for this purpose. First, the spectral similarity between the bases 𝐖 obtained in Eq. (5) and a dictionary of pre-learned bases obtained from isolated heart-sound signals belonging to a training database is measured. Second, an analysis of the energy distribution of each basis obtained in Eq. (5) was performed to discriminate between heart and lung sounds. This analysis assumes that heart sounds condense their energy below the maximum frequency of 260 Hz, whereas most lung sounds condense their energy above that frequency. Note that this threshold is in line with other works in the state-of-the-art [39, 29, 30, 40]. Then, the authors proposed a third clustering method based on the temporal periodicity of target sounds.

Finally, reconstruction of the target signals is performed using a generalized Wiener filtering strategy. This Wiener filtering strategy minimizes the mean-square error between the estimated and desired signals. The previously predicted parameters were used to estimate the magnitude spectrograms of the lung and heart signals by

Then, a generalized time-frequency mask over the STFT domain is used to obtain the source signals xH⁢(m) and xL⁢(m) from the mixture y⁢(m). In this regard, a Wiener filtering approach was used to ensure that the reconstruction process was conservative [41]. Note that this approach works as a soft mask for the observed mixed signal, which scales the magnitude of the mixed signal at every frequency component with values between 0 and 1 to determine their corresponding frequency component values in the estimated signals.

3.1Hardware prototype

One of the main goals of this study is to design a biomedical audio signal processing system that can be run on portable low-cost systems. Therefore, the main constraint of the development is that low-power microcontrollers do not provide sufficient computational resources, both memory and processing power [42, 43]. Therefore, additional hardware must be considered when computing the proposed method.

To address the target problem, the following requirements for the proposed prototype have been set: 1) low power consumption (suitable to be powered by a battery), 2) low cost (the total price of all components must be less than a few tens of euros), and 3) wireless connectivity, that is, Internet of Things enabled (embedded connectivity is essential to run part of the algorithm on the computing nodes). Based on these requirements, the following components were selected to develop the prototype.

One of the most popular devices for these types of applications is the well-known ESP32 [44] developed by Espressif Systems. This microcontroller offers a cost-effective solution with several GPIO (general-purpose input/output) pins and exposes different communication interfaces to attach peripherals: I2C, SPI, I2S, and UART, in addition to wireless connectivity (Bluetooth and Wi-Fi).

Figure 1.

Block diagram of the proposed approach.

Figure 2.

INPM441 MEMS microphone.

Figure 3.

Low-cost prototype for monitoring and estimating the HR and its main component: (A) INPM441 MEMS microphone, (B) ESP32-based microcontroller and (C) power battery.

Regarding the microphone, the choice was made based on a compromise between the suitability of the device for this application and its price. Although some microphones are specifically designed to capture biomedical signals [45], the main drawback of these devices is their high price, which makes them less attractive for use in the proposed prototype. For this reason, this research focused on finding a general-purpose and low-cost microphone that can be directly attached to the ESP32 using any of the available communication interfaces. A cost-effective option is the INMP441 omnidirectional microphone [46], as displayed in Fig. 2. The most relevant features of this microphone are summarized in Table 1. This microphone was also selected because of its bandwidth, from 60 Hz to 15 kHz. This frequency range guarantees proper performance of the proposed system. To capture the heart and lung signals, an adapter made of cone-shaped plastic was used along with a plastic film to place the microphone in contact with the skin. Note that this is how regular stethoscopes work to match the impedance of the skin to that of air.

As previously mentioned, the ESP32-based microcontroller does not have sufficient computational resources to process the audio stream and to separate the target sounds. Therefore, it serves as an acquisition device that sends captured audio to an external agent. Thus, data from the microphone are digitally sent to ESP32 using an I2S connection. The audio stream is converted from the original 32-bit format into 16-bit for the processing algorithm that runs on computing devices. This is accomplished in a separate task of the FreeRTOS operating system, and is forced to run in a different core from that used in the communication task. Then, ESP32 streams the audio to the network using a multicast connection. This vastly simplifies the entire process, as no IP information is required to receive the data on the different devices, thus enabling the possibility of simultaneous access and recording of the data stream on several devices. Finally, the audio stream is processed on an Android phone (running the algorithm in native C code) or a general-purpose computer.

The prototype is shown in Fig. 3. The total weight of the prototype was 110 g with a battery and 65 g without a battery. Finally, the total cost of the hardware involved (phone and computer devices excluded) is less than 15 euros, which results in a low-cost solution and places it well below current solutions on the market.

3.2Proposed method

The system takes as input the cardiac and respiratory sounds captured by a microphone attached to the surface of the chest wall using previously described hardware. The output of the system is the separation of the lung and heart sounds. As a case study, a final stage for the heart rate prediction has been also designed.

A block diagram of the proposed method is shown in Fig. 1. The framework consists of three main stages: preprocessing, separation, and HR estimation. In the following subsections, the main functions of each stage are described in detail.

3.2.3HR estimation stage

Herein, a simple and efficient method designed for HR estimation is presented. Note that this method takes as input the magnitude spectrogram of the heart signal X^H⁢(f,t) estimated by the separation stage (see Section 2.2). Therefore, substituting this last stage, the developed prototype could also be used to monitor other health parameters based on these audio biosignals. This work is focused on HR estimation for the validation of the prototype, because it is one of the most important parameters for predicting cardiac diseases.

To estimate the HR, the first step consists of transforming this input spectrogram into a more meaningful compact representation from which to infer the heartbeats. The new representation is computed by the Euclidean norm of the estimated spectrum of the heart signal for all temporal frames, as follows:

(20)

γ⁢(t)=∑f=1F|X^H⁢(f,t)|2

Then, a linear interpolation process is performed over γ⁢(t) to use this function in the time domain to obtain γ~⁢(m).

To estimate the HR, the autocorrelation function of the current representation was used to determine the periodicity of the signal obscured by noise. Thus, to achieve an autocorrelation function with clear and well-defined peaks, the values of γ~⁢(m) that are below a specific threshold α must be set to zero to emphasize the strongest and discard the least significant peaks. Thus, the autocorrelation function can be computed as

(21)

Rx⁢x⁢(k)=1N⁢∑0N-1-kγ~⁢(m)⁢γ~⁢(n+k)

k=0,…,N-1

where N is the sample length.

Finally, the strongest peak within the range of possible heart rates for a healthy heart must be selected, that is, [40,190] beats per minute. Therefore, HR can be estimated as

(22)

HR^=60kmax⁢Ts

where kmax is the k index of the strongest peak and Ts is the the sampling period. Figure 4 illustrates the proposed HR estimation model.

Figure 4.

Proposed HR estimation model. a) Input estimated heart signal. b) Magnitude spectrogram of the estimated heart signal. c) Euclidean norm of magnitude spectrogram. d) Autocorrelation function with the beat period.

The proposed method for HR estimation is summarized in Algorithm 4.

: HR estimation algorithm[1] Inputs: Mixture audio signal y⁢(m). Compute 𝐘 using the STFT to obtain the signal representation in the frequency domain. Normalize 𝐘 following the Eq. (7). Compute the SVD of 𝐘 to obtain the orthogonal matrices 𝐔 and 𝐕, respectively. Initialize NMF matrices 𝐖 and 𝐇 using the orthogonal matrices 𝐔 and 𝐕, respectively. 𝑖𝑡𝑒𝑟=1 to N𝑖𝑡𝑒𝑟Update 𝐖 using the Eq. (5). Update 𝐇 using the Eq. (4). i=1 to kDetermine if 𝐰(k)∈𝐖H (see Secion 3.2.2). Compute 𝐗^H using the Eq. (6). Obtain γ⁢(t) from 𝐗^H using Eq. (20). Obtain γ~⁢(m) from the linear interpolation of γ⁢(t). Set the values of γ~⁢(m) below α to zero. Compute the autocorrelation function Rx⁢x⁢(K) using Eq. (21). Find the strongest peak of the autocorrelation function. Estimate the heart rate HR^ using Eq. (22). Outputs: The estimated heart rate HR^.

Concerning the computational complexity of the HR estimation stage, note that the implementation was performed using BLAS routines and OpenMP directives. Thus, the theoretical computational complexity of the parallel version can be computed by considering the property of the sum of the Big O notation as:

(23)

O⁢(T⁢log2⁡(T)p)

4.1Datasets

The performance of the proposed HR estimation approach was assessed using three different databases. The database developed by Yaseen et al. [57] was used first. This database consists of a normal and abnormal set of heart sounds, classified into five categories: normal, aortic stenosis, mitral stenosis, mitral regurgitation, and mitral valve prolapse. For the evaluation of the proposal, a subset of normal heart sounds was used, consisting of 200 audio samples with a duration between 2 and 3 s with an 8 kHz sampling frequency. Because this dataset does not include lung sounds, it was only used for validating the proposed HR estimation method.

Second, a public dataset of heart sounds released in the Classifying Heart Sounds Pascal Challenge competition [58] was used. This competition provides two datasets, Dataset A and Dataset B. Here, Dataset A was used because it was generated in conditions similar to those pursued by the proposal, unlike Dataset B, which is a set collected in a hospital environment. Dataset A was obtained from volunteers and recorded with the iStethoscope (i.e., a mobile application) in real-world conditions. The dataset contains 176 audio samples with a duration between 5 and 10 s in wav format with 44.1 kHz sampling frequency, and is organized as four categories: normal, murmur, extra heart sound, and artifact. Once again, 21 audio samples categorized as normal sounds were selected for the validation of the proposal.

Finally, the database collected by Canadas-Quesada et al. [28] (APAC) is used. This database consists of a selection of 72 audio samples with a duration of 7 s in wav format with 8000 Hz sampling frequency. Audio files were generated using real-world heart and lung sound signals. Therefore, this database is suitable for demonstrating the need to first separate heart and lung sounds before estimating the HR.

These previously described datasets were created for several purposes, including the classification of heart-sound signals and separation of lungs and heart-sounds. However, in this study, the aim was to evaluate the proposed HR estimation task. Hence, cardiologists are required to label all audio files with the ground-truth heart rate. It is important to note that some samples were discarded by specialists because the heartbeat was not distinguishable. However, it was decided to test the proposed algorithm using all audio files.

4.2Experimental setup

In this study, the time-frequency representation is obtained using a 512-point STFT and half overlap between adjacent frames. The sampling rate is 8 kHz. With regard to the signal model, convergence of the NMF decomposition was empirically observed after 100 iterations. Thus, this value was chosen as the number of iterations required for the factorization process. Moreover, it was found that the optimal value of α was 0.2 (see Section 3.2.3).

For experimentation, the NVIDIA Jetson AGX Xavier development kit was used as testbed. This kit consists of a system-on-chip (SoC) with an eight-core ARM v8.2 CPU that operates at 2.26 GHz. This testbed supports several operating modes configured using the NVPModel command tool. This enables simulation of a wide range of mobile devices and embedded systems. The details of all operating modes are listed in Table 2. For each mode, the first equivalent ARM architecture with its market release year is also included. Seven operation modes and four power envelopes were defined. The possible numbers of running cores are eight, six, four, and two with different CPU frequencies.

4.3Evaluation metric

The source-separation performance of the proposal was objectively evaluated using the objective measures provided by the BSS_Eval toolbox [59]. These metrics are commonly accepted and represent a standard approach in the specialized scientific community for testing the quality of separated signals, allowing fair comparison with other state-of-the-art methods. BSS_EVAL provides the following metrics based on the energy ratios for each separated signal: the source to distortion ratio (SDR), the source to interference ratio (SIR) and the source to artifacts ratio (SAR) [59]. SDR reports on the overall quality of the separation process, SIR provides a measure of the presence of lungs in heart sources and vice versa, and SAR reports on the artifacts in the separated signal due to separation and/or resynthesis.

For quantitative evaluation of the reliability of the proposed HR estimator, the relative error metric was chosen. Thus, the estimated heart rate from the audio signal is denoted as HR𝑒𝑠𝑡, the ground truth is denoted as HR𝑔𝑡 and the relative error rate between the two is denoted as HR𝑒𝑟𝑟 and can be expressed as

4.4Algorithms for comparison

To demonstrate the benefits of the proposed method, its separation performance is compared with that of other state-of-the-art algorithms. The different approaches compared here are as follows:

Figure 5.

Heart sound separation results using the BBS_Eval metrics averaged over the APAC dataset [28]. The error bars represent 95% confidence intervals.

Figure 6.

Lung sound separation results using the BBS_Eval metrics averaged over the APAC dataset [28]. The error bars represent 95% confidence intervals.

4.5Results and discussion

This section presents the results obtainedfrom the evaluation of the proposed method in terms of source separation. The proposed HR estimation algorithm was then evaluated using the datasets presented in Section 4.1. Finally, the results in terms of execution times are presented to demonstrate the feasibility of the proposed algorithm in real-time.

4.5.2HR estimation results

The results obtained for each dataset presented in Section 4.1 are shown in Fig. 7. First, the performance of the proposed HR estimation method (see Section 3.2.3) is analyzed. As can be observed, the results obtained for the Yanseen dataset reveal a reliable performance of the HR estimation stage, as the maximum error in the 200 samples that compose the dataset is less than 5%. Note that the separation stage was not required because lung sounds were not included in the dataset. Therefore, this measure informs us about the best estimation that can be achieved in the absence of noise-masking heart sounds, that is, using clean heart sound signals.

Figure 7.

Averaged relative errors in HR estimation results over the datasets described in Section 4.1.

The results for the Pascal and APAC datasets showed the important role of the separation stage in the HR estimation task. As can be seen, the results significantly improve when the separation process is included in the algorithm for both datasets. In the case of the Pascal dataset, 9 of the 21 audio samples obtained a relative error greater than 5% when the separation was not performed. However, when the separation was applied, only 4 audio samples obtained a relative error greater than 5%, and the improvement was not more significant because the audio samples were recorded under highly noisy conditions. Note that the demanding 5% threshold was chosen as a benchmark to validate the reliability of the estimator. In the APAC case, 32 audio samples obtained a relative error greater than 5% when the separation was not performed. On the other hand, the relative error for the 72 audio samples is less than 3% when the separation of lung and heart sounds is performed.

4.5.3Execution time

To assess the possibility of estimating the HR in real time using the developed prototype (i.e., low-cost devices), the execution time of the proposed algorithm was measured over NVIDIA AGX Xavier. As the objective is to measure a subject’s heart rate periodically in real time, it is assumed that the algorithm should provide an estimate at least every second. Moreover, for proper estimation of the HR, an audio sample of approximately ten seconds long was considered in this experiment. Therefore, the HR was estimated every second using the previous ten seconds. Thus, the algorithm will be able to respond in real time if it is able to process ten seconds of audio in less than a second. The choice of this window is justified by the criterion of lower estimation error (the longer the window, the lower the error).

Figure 8.

Experimental results in terms of execution times measured in seconds as a function of the operating mode of the NVIDIA AGX Xavier.

The results obtained in terms of the execution times are summarized in Fig. 8. These results are presented as functions of the operation mode of the NVIDIA AGX Xavier (see Table 2). As expected, the execution time increases as the CPU frequency and the number of cores decrease. Note that real-time is reached for CPU frequencies higher than 1200 MHz, regardless of the number of cores and power budget.

For low-cost devices with four available cores (i.e., Mode 5), the proposed parallel approach allows the execution of the application in real time. It is important to note this configuration because most current mobile devices satisfy this requirement.

5.1Experimental setup

To evaluate the prototype, a real-world scenario was designed in which a subject performed a stress test on a bicycle while his or her HR was estimated. Thus, a setup was defined to adjust the physical strain and consequently modulate the cardiac response. This setup consists of a smart home trainer (Elite Direto) and a conventional bicycle. The main advantage of this trainer is that it is fully software controlled, which means that the resistance to a pedal stroke can be tuned and adjusted to achieve the desired power development and keep it steady over a period of time. The experiment was carried out with four healthy male subjects (authors of this work), aged between 25 and 48 years. It consisted of cycling for four minutes under different levels of pedal resistance (i.e., power zones). Four different power zones were defined and set using trainer software by means of a feature known as the ERG mode [63]. The first power zone (idle or 0 W) was introduced to measure the heart rate of subjects in the absence of physical activity. The other three zones were established at 50 W, 100 W, and 200 W. The duration of each zone was set as one minute.

As a reference, the HR was measured by a standard bluetooth low energy (BLE) HR monitor based on ECG and attached to the chest. These data were recorded, along with the estimation provided by the proposed prototype. Table 3 shows the average HR measured for each power zone and subject in terms of BPM.

The overall setup is illustrated in Fig. 9 and explained as follows: a cycling simulator (Bkool) [64] was installed on a tablet and wirelessly connected to the smart trainer and the BLE HR monitor. The simulator was in charge of displaying the reference HR and controlling the resistance of the smart trainer in response to the provided target power. On the other hand, the microphone of the proposed prototype was attached to the subject’s chest and the ESP32 microcontroller (connected to the network by Wi-Fi) sent the audio stream to the chosen multicast address. Then, an Android phone and/or a computer were also connected to the same network to receive the data and estimate the HR using Algorithm 4.

5.2Results

Figures 10 and 11 summarize the results obtained from the experiment described in Section 5.1. As shown in Fig. 10, the relative error compared with the ECG-based estimator depends on two main factors. First, worse performance was obtained for the first power zone (or idle) because the heart beats with less energy and are masked by external noise. However, for zones of physical activity (i.e., power zones 2, 3, and 4), the deviation from the ECG was well below 5%. Second, it was found that the physical shape of the subjects affected the HR estimation. This can be observed in the estimations of subjects 3 and 4 for zone 1, where the error is higher. Both subjects had a higher body mass, which made it difficult to determine the optimal point for microphone location. Even so, the median error for both subjects is close to 5%, which is acceptable for low-cost devices. Regarding gender bias, previous studies on auscultation signals showed no statistically significant differences between the male and female test groups [65]; therefore, similar results are expected to be obtained with a broader and more diverse sample.

Figure 11 shows the HR estimation of the proposed prototype compared with that of the ECG-based estimator. In general, the results obtained using the proposed method are in line with the ECG-based estimation, and discrepancies are only observed for a few instants of time. It should be noted that this is a strong point of the proposal, as it achieves similar results to commercial proposals, but at a much lower cost.

Finally, to measure the power consumption of the prototype, a 18650 battery (3.7 V) was used (see Fig. 3). The system was then operated uninterrupted at 2200 mAh for 16 h and 36 min until the battery charge was completely depleted.