Novel phonocardiography system for heartbeat detection from various locations

The paper presents evaluation of the proposed phonocardiography (PCG) measurement system designed primarily for heartbeat detection to estimate heart rate (HR). Typically, HR estimation is performed using electrocardiography (ECG) or pulse wave as one of the fundamental diagnostic methodologies for assessing cardiac function. The system includes novel both sensory part and data processing procedure, which is based on signal preprocessing using Wavelet Transform (WT) and Shannon energy computation and heart sounds classification using K-means. Due to the lack of standardization in the placement of PCG sensors, the study focuses on evaluating the signal quality obtained from 7 different sensor locations on the subject’s chest and investigates which locations are most suitable for recording heart sounds. The suitability of sensor localization was examined in 27 subjects by detecting the first two heart sounds (S1, S2). The HR detection sensitivity related to reference ECG from all sensor positions reached values over 88.9 and 77.4% in detection of S1 and S2, respectively. The placement in the middle of sternum showed the higher signal quality with median of the proper S1 and S2 detection sensitivity of 98.5 and 97.5%, respectively.

Wavelet transform.Discrete Wavelet Transform (DWT) is a widely used tool for preprocessing of biological signals.The signal is decomposed using the adequate type of maternal wavelet into detailed and approximate coefficients (wavelet decomposition).The used type of maternal wavelet and level of decomposition differ across studies.Usually, the wavelets Daubechies (db), Discrete Meyer, Coiflet (coif), Morlet, or Symlet (sym) are used.Some studies also verified the performance of the wavelets Harr, Biorthogonal (bior), or Reverse Biorthogonal (rbio), see Table 1.Some of the studies compare the different wavelet types from the perspective of Table 1.Frequently used wavelets for PCG processing..
Signal segmentation.The signal segmentation and associated localization of S1 and S2 are important steps for the PCG signal analysis.The segmentation methods can be characterized as: a. Indirect-PCG is measured together with ECG (or pulse wave) reference, when stamps from ECG signal (moment of occurrence of R and T waves) are used to identify the individual heart cycles and the corresponding heart sounds 2,11,12,22,23 , b. Direct-only PCG is measured, which needs to be transformed into the domain that allows increasing of S1 and S2 amplitude 2 .
The direct segmentation is more complex, see Fig. 1.The segmentation methods in time domain provide direct and the most reliable localization of heart sounds, however, they demand high-quality signal preprocessing due to their high susceptibility to noise 10 .From that reason, using the time-frequency domain is quite popular 2,23,24 , where the increased spectral power on frequencies from 40 Hz to 100 Hz is visible.This approach includes S1 and S2 detection based on signal envelope, which was dealt in 5,7,14,16 .
The direct segmentation is easy to perform in healthy subjects 36 , assuming that systole is shorter than diastole.However, systole and diastole duration is dependent on heart rate 23 , so the detection can be problematic in patients with tachycardia 2 .In some cases, it can be difficult to ensure elimination of extra peaks and preserving of true heart sounds, which is crucial for correct direct segmentation 37 .

Material and methods
Measurement system.The used data acquisition (DAQ) system (see Fig. 2) includes chassis PXIe-1092 38 with a Field-Programmable Gate Array (FPGA) module from National Instruments (NI)-NI PXIe-7862 39 , which allows to measure 16 analogue input signals with 16bit resolution and maximum sampling frequency of 1 MS/s per channel.This module was connected to analogue inputs via terminal block SCB-68A 40 .For data transmission between PC and FPGA module, the direct memory access (DMA) was used.
Seven microphones GRAS 40PP-10 with a frequency range from 10 Hz to 20 kHz and sensitivity of 50 mV/ Pa 41 were used to sense and convert the acoustic signal.A relative connection of microphones was used, which allows the sensing of pressure changes of the chest movement at the location of sensor 5,6 .In order to supply the microphones with the required current of 2 to 20 mA, they were connected to the terminal block via an eightchannel IEPE amplifier M208B, the gain of which was set to 0 dB (gain = 1) and was used only for power supply.
The conduction to the microphone is made through the air-filled tube.An air chamber for conduction of pressure changes to the microphone was created by the 3D-printed sensor of the hemisphere shape with a diameter of 40 mm and PVC tube connecting the sensor with microphone, see Fig. 3.The frequency analysis of the measured signal showed that the proposed system (sensor with a diameter of 40 mm and tube length of 1 m) absorbs signals with frequency higher than 100 Hz.Because the heart sounds S1 and S2 are composed of frequencies lower than 100 Hz, the system is suitable for the HR detection.
The whole sensory part of system that is in contact with patient is made of low-cost and non-electrical materials, which can be advantageous for many applications, such as use within magnetic resonance (MR) environment, when the use of electrical conductive parts could cause an injury of patient 8 .

Dataset.
The measurement was preformed on 27 volunteers (10 women and 17 men) aged from 21 to 29 years.All collected information about subjects, such as age, height, weight, and body mass index (BMI) are summarized in Table 2.One of the probands (#22) suffers from heart murmurs of the intensity of 1/6 with a high blood pressure and tachycardia.All of the remaining signals were without pathology.
All experiments were performed in accordance with the relevant regulations and approved by the Ethics Committee of the Technical University of Ostrava.Experiments were completely safe, using certified HW from National Instruments.An informed consent form was prepared for this study, which all participants signed to   www.nature.com/scientificreports/confirm that they agreed to the publication of results, where all their data will be anonymised.All data generated or analysed during this study are included in this published article (and its supplementary information files).
The study is limited to a small sample of young adult participants, since the aim of this study is to test the usability of the proposed sensors and determine the most suitable position for measuring the PCG signal.It is important to note that having only one pathological sample is not sufficiently indicative, and therefore, a larger sample of participants that encompasses various age groups and includes a wider range of pathological signals will be a goal of future studies to obtain more meaningful results (see Discussion section).
Each volunteer underwent 10-min measurement in the supine position using the proposed measurement system.The sensors were located according to Fig. 4 and secured by three stretch belts.For measuring reference ECG, a clip-on or floating electrodes were used.All the signals were sampled with a frequency of 2.56 kHz.Parameters of the measured signals are summarized in Table 3 and all signal processing steps were realized using LabVIEW software, National Instruments.

PCG signal preprocessing.
For preprocessing, PCG signal was filtered to highlight the desired heart sounds using band-pass filter (BPF) and DWT.Then, calculation of Shannon energy of the signal and its envelope was used for HS detection, see Fig. 1.
Band-pass filter.First, the signal was filtered using Butterworth BPF of the 2nd order with a lower cut-off frequency of 2 Hz and higher cut-off frequency of 100 Hz, see Fig. 5a.Elimination of low frequencies is performed primarily to suppress motion artifacts and compensate isoline of the signal.Frequencies higher than 100 Hz contain undesired noise.www.nature.com/scientificreports/Discrete wavelet transform.As the next step, the signal was processed using DWT.Preliminary experimental measurements were first performed to find out a suitability of different maternal wavelets for the obtained PCG signals.In accordance with a literature (see Table 1), the signal was filtered using available wavelet functions (Advanced Signal Processing Toolkit by NI LabVIEW) from the families of Daubechies, Haar, Biorthogonal, Coiflet, and Symlet.Based on these experimental measurements, Daubechies (db12) wavelet with a decomposition level of 5 was chosen for the final filtering process due to the best results compared to other tested wavelets.The example result of WT filtering is shown in Fig. 5b.
Normalization.The filtered signal was normalized to values from −1 to 1 according to Eq. ( 1): where x n (i) is normalized signal and x(i) is the original one, see Fig. 5c 16, 36 .The normalization was performed within 1 s sections, since its applying on the whole signal can be inaccurate due to presence of undesired high amplitudes, e.g., caused by movement.
Signal energy.From the normalized signal, Shannon energy of the 3rd order was calculated according to Eq. ( 2): where E S (i) is the output energy and x n (i) is the normalized signal 16 .Shannon energy highlights mean values of the signal and suppresses low values together with values close to 1, see Fig. 5d.The benefit of Shannon energy compared to absolute value is the increasing of differences in sounds with middle and high intensity, which simplifies the further localization of heart sounds 5,16,36,42 .
Signal envelope.The envelope was created by filtration of signal energy using Butterworth low-pass filter (LPF) of the 2nd order with zero phase shift and cut-off frequency of 10 Hz, see Fig. 5e.To increase an envelope amplitude, the normalized Shannon energy (see Fig. 5f) was determined using Eq. ( 3): where E N is the resulting energy, µ is mean value of the signal E S and σ is standard deviation of the signal E S .
Heart sound detection.The peaks from this preprocessed signal were detected and further segmentation was applied, as shown in Fig. 1.Further, the intervals between the detected peaks were classified as systole or diastole.The peak before systole was assigned as S1 and the one before diastole as S2.
Peak detection.First, peaks of the envelope located above the defined threshold Th 1 were detected (see Fig. 6): where µ is a mean value and σ is a standard deviation of the signal E S .The value of this threshold is 90 % of the signal shift (see Eq. ( 3)).This threshold allows for the detection of weak heart sounds as well as noise.Peaks detected using this threshold are not primarily labelled as heart sounds but are utilized for further completion based on the criteria described below.Duration of the heart sounds is usually about 150 ms.If the interval between two consecutive peaks was shorter than 150 ms, the peak with lower amplitude was removed.With this process, even low-amplitude false peaks were detected and were confirmed by the further processing.
(1)  23 .The intervals between peaks with higher amplitude than the threshold Th 2 (see Fig. 6) were investigated: where E S is a mean value of the signal E S .Under ideal conditions, the intervals S1-S1 (or S2-S2) take place between even (or odd) peaks, representing heart rate.However, there occurred some extra peaks and some non-detected peaks, see Fig. 7. Therefore, the parameter presumed heart period calculated as modus T S−S of the intervals between even and odd peaks was defined: where T o;n and T e;n are intervals between two odd and even peaks, respectively.This value was used for approxi- mate determination of heart rate and selection of primary centroids for classifying the intervals using K-means method.
For K-means cluster analysis, the number of clusters was selected as k = 3 .Intervals are divided into clus- ters, which can be called "systole", "diastole", and "longer intervals".The primary centroids C were selected as 65 • T S−S , and C 3 = T S−S for systole, diastole, and longer intervals, respectively.Values C 1 and C 2 represent 35% and 65% of T S−S , which is approximate ratio of duration of systole and diastole during rest heart rate 23 .
Peaks located before the interval "systole" are marked as S1.Similarly, the peaks before the interval "diastole" are marked as S2.The peaks located before "longer intervals" or under the threshold Th 2 stayed not classified in this step, see Fig. 8.

Peak confirmation.
Based on two conditions according to 53 , the false peaks were removed and non-detected peaks were assigned, see Fig. 9.Both conditions count with the presumed heart period T S−S : a. Removal of false peaks-if the interval between consecutive peaks S1-S1 (or S2-S2) is shorter than 0.8 • T S−S , the second peak is marked as "non-classified peak".This situation is shown in Fig. 9a.
(5)   9b, the interval is checked whether there is the missed peak.
Around the position, where the peak should occur (position of the assigned peak + T S−S ), window of the length of 300 ms (150 ms in each direction) was located.If there was found a non-classified peak, it was assigned to the same category as the previous one.If there were more peaks than one within the window, the one closer to the center of the window was added.
Processing of reference ECG.The reference ECG was filtered using two Butterworth filters of the 2nd order with zero shift: a. Band-pass filter was used for elimination of artifacts caused by motion and respiratory activity with cut-off frequencies of 1 and 150 Hz. b.Band-stop filter was used to remove power line interference (50 Hz) with cut-off frequencies of 48 and 52 Hz.
The R peaks were detected with using adaptive threshold, which varied for the individual windows of 5 s, when the threshold was set as 70% of the maximum amplitude within the window.Peaks falsely marked as R peak and non-detected peaks were removed or added, respectively, according to the similar conditions like in the case of removal and addition of heart sounds 53 .
After the automatic detection of R peaks, their positions were thoroughly verified by the authors to avoid any incorrect assessments.The chosen detection method demonstrated accurate performance in the examined sample of healthy subjects.For pathological signals, more precise algorithms would be necessary, such as [54][55][56] .
Evaluation parameters.This section describes the parameters used for evaluation of HS detection performance on the individual locations (P1-P7).Also, there are described the methods used for comparison of heart rate values obtained from PCG with the ones from ECG.  www.nature.com/scientificreports/Peak detection.For evaluation of HS detection performance, reference ECG was used.Based on 23 , where the indirect PCG segmentation is described in detail, time windows of the expected occurrence of heart sounds were selected, see Fig. 10.Time duration of these windows depends on heart rate (i.e., it shortens with increasing HR).The S1 should be located in the interval of 0-0.2 • T R−R after the detected R peak, where T R−R is a median of intervals between two consecutive R peaks.Duration of this interval is around 160 ms (median of all measurements is 168 ms).
For S2 detection, the interval was set from 1.2 • T R−T (end of T wave) to 0.6 • T R−R , where T R−T is a median of intervals between R peak and T wave maximum.Duration of this interval was longer than for S1 detection-up to 200 ms (median of all measurements is 189 ms).
Based on the success of peak detection within the selected window, parameters such as True Positive (TP), False Positive (FP), and False Negative (FN) values were calculated: • TP indicates the number of detected peaks within the window, • FP indicates the number of detected peaks out of the window, • FN indicates the number of non-detected peaks within the expected interval.
These parameters were then used for calculation of sensitivity (SE) according to Eq. ( 7), which is one of the most significant statistical parameter stating the probability of the occurrence of the heart sound in the selected interval 2,10 .According to Eq. ( 8), Positive Predictive Value (PPV) is determined, i.e., a probability of the correct peak detection.
Heart rate evaluation.The signal from the position with the best values of SE and PPV was used to determine heart rate HR S1 (or HR S2 ) from the intervals S1-S1 (or S2-S2).Heart rate obtained from R peaks of ECG HR R was used as a reference.Individual HR values were then calculated as follows: where T S1−S1 is median of intervals between two consecutive S1 peaks in the window of 2 s and HR S1 is heart rate with unit BPM (beats per minute, min −1 ).Similarly, HR S2 and HR R can be determined using periods The calculated values of HR S1 and HR S2 were then compared to reference HR R using pair Wilcoxon test, since the normality of data was rejected in some cases (Shapiro-Wilk test, p-value = 0.05) 57 .The p-values less than 0.05 were considered to indicate a statistically significant difference.

Results
Quality of HS detection in the individual locations (P1-P7) was evaluated using SE and PPV.Their results for all measured subjects are presented by range of their minimum and maximum (min;max), median, and interquartile range (IQR), defined as: where Q 1 is the first quartile and Q 3 is the third quartile of the SE and PPV values.All these parameters are summarized in Table 4 and represented by hybrid boxplots in Fig. 11a,c for S1 detection and in Table 5 and Fig. 11b,d for S2 detection.(9) HR S1 = 60

Statistical comparison of HR.
The HR values (see Fig. 13) obtained by PCG ( HR S1 , HR S2 ) from position P3 (with the best performance) were compared to ECG ( HR R ) using Wilcoxon pair test.The statistically signifi- cant difference in medians of pair data was not shown in 24 cases of 27 (88.9%)for HR S1 and in 20 cases of 27 (74.1%)for HR S2 .The results of pair test (p-values) for the individual subjects are presented in Table 6.

Discussion
In this study, the designed low-cost and easy-to-use sensor for PCG measurement was successfully tested.Thanks to the use of non-metallic materials, the sensor is for MRI environments, for which it was primarily designed.The main goal of the study was to determine where on the patient's chest the designed sensor should be placed to obtain a high-quality signal.While this study does not provide a general standardization of sensor   placement for PCG measurement, it identifies the appropriate position for placing the sensors designed in this study.The best signal quality was achieved from position P3 with a median of sensitivity 98.5% for S1 and 97.5% for S2.This position can be chosen for measurements using the designed sensors and also serves as a recommended position for future studies.The feasibility of using these sensors for HR measurement was demonstrated by comparing the HR data obtained from the R-R, S1-S1, and S2-S2 intervals.Wilcoxon's signed-rank test showed that the majority of results are comparable.Cases where the p-value was lower than the selected significance level could arise due to improper sensor attachment, ambient noise, or signal pathologies.Generally, the used system (see Fig. 3) composed from microphone GRAS 40PP-10, PVC tube, and 3D-printed sensor provides a frequency response up to 100 Hz suitable for detection of S1 and S2.However, according to the literature 12 , S2 is composed also with frequencies higher than 100 Hz, whose response could provide higher quality of S2 detection in both healthy and unhealthy patients.

Sensor position SE (%)
In Table 7, the results of the proposed study are compared with the recent studies (from 2017 onwards) that focus on a similar topic (i.e., detection of heart sounds) or provide detailed descriptions of sensor placement.These studies are compared based on used dataset, processing methodology, sensor placement (if specified), and evaluation parameters.The former studies are described and compared in the comprehensive review by Ismail et al. 10 .
Because the mentioned studies use existing databases (PhysioNet 58 , PASCAL 59 , eGeneralMedical 60 , Michigan heart sound and murmur database 61 ) to test their algorithms, they do not focus on the sensor placement itself.It is important to note that the PhysioNet database only describes common auscultation areas, including the aortic area, pulmonic area, tricuspid area, and mitral area.Additionally, in our case, we identified heart sounds based on the R and T wave identification method 23 .Other studies, however, rely on the algorithm 62 , which uses the R peak and the end of the T wave, or employ different methods for sound localization.Furthermore, the evaluation criteria for correct detecting S1 and S2 can vary among studies, with different techniques and metrics being used.In addition to sensitivity (SE) and positive predictive value (PPV), authors may also employ metrics such as accuracy (ACC), specificity (SP), or precision for evaluation.
One of the suitable studies for comparing heart rate (HR) could be study 7 , where the authors use a correlation coefficient for comparing HR S1 and HR R .In conclusion, they also mention that "determining the optimal localization of the PCG sensor is very difficult, and evaluation depends on many parameters".
In this study, DWT was used for PCG signal preprocessing since it is a suitable and very popular tool in the area of non-stationary biological signals.Although many studies endeavour to find an optimal parameter setting for PCG denoising 12,15,[19][20][21] , it is still not clearly defined and types of maternal wavelets differ across studies.In this paper, the results of wavelet db12 are presented as the best achieved from previous testing of wavelets db2-db14 using Advanced Signal Processing Toolkit (NI LabVIEW).In future research, other types of maternal wavelets could be tested, e.g. using software platforms that enable a wider range of wavelets for signal processing  www.nature.com/scientificreports/(e.g., MATLAB).Further research could also explore beyond the discussed DWT, as it is just one of many methods for processing PCG signals.In addition to the mentioned linear filters, PCG signals can be processed using other methods such as Homomorphic filtering 2,42 , Fast Fourier Transform (FFT), STFT 23,24,71,72 , or Empirical Mode Decomposition (EMD) 73,74 .One shortcoming of the proposed method could lie in the differentiation of S1 and S2 based on duration of the intervals between heart sounds (using K-means).This approach is suitable for measuring in rest conditions, when HR is ranged from 60 to 90 bpm.However, during non-resting activities and HR around 130 bpm, duration of systole and diastole is almost the same, see Fig. 14 23 .In this case, heart sounds should be defined according to another criteria, such as amplitude or time of occurrence after R peak using indirect segmentation.Also, when measuring in motion, a closer contact of the sensor with the patient's skin and more sophisticated signal processing would be needed, since the signal would be affected by many motion artifacts.
The motion artifacts occur also due to breathing and can cause a false HS detection in the proposed approach.The sensors were attached to the body via three rubber strips, as shown in Fig. 15.The proper contact of the sensor with the skin is one of primary preconditions for measuring high-quality PCG signal.However, the attachment of many sensors by one strip was very difficult due to different chest size and shape of the individual subjects.Common fastening of more sensors (P4-P7) can cause insufficient sensor-skin contact and movement during breathing.In these cases, another way of the sensor attachment should be chosen to fix this issue.This observation can explain the best performance of the location P3, where the sensor was fastened by the extra belt.
Future research should involve expanding the selected subjects group.This study focused on a relatively small sample, consisting of 27 healthy individuals aged 21-29 years.Within the study, the appropriate sensor placement and usability were verified among the healthy individuals.To further strengthen the study, future research should explore the dependency of HS detection quality on the patients' age.Additionally, the proposed approach should be tested in the context of measuring patients with cardiac conditions, as the effects of heart diseases can significantly affect HS detection quality by reducing the amplitude of the S2 sound (as observed in subject #22).Furthermore, the presence of ambient noise in the signal should be considered.Therefore, further studies should verify the usability of the sensors and the proposed algorithm in noisy environments.The designed algorithm should be optimized, such as by adjusting the thresholds Th 1 and Th 2 , or implementing an adaptive threshold that can dynamically adapt its level based on the amplitudes of the S1 and S2 sounds.For the evaluation of pathological or noisy signals, a commercially available PCG device or an evaluation performed by an expert, serving as a reference standard, could be valuable.

Figure 3 .
Figure 3. 3D printed extension connected to the microphone.
Detail coefficients of 5th level decomposition using db12.(c)Normalized PCG signal (blue).(d)Normalized PCG signal (gray) and its 3rd level Shannon energy (blue).(e)Shannon energy (gray) and envelope of this energy (blue).(f) Original envelope (blue) and its normalized Shannon energy (red).

2 Figure 7 .
Figure 7. Intervals between even ( T e;n ) and odd ( T o;n ) peaks.Correctly detected intervals (which correspond to the HR) are marked in green.Intervals marked in red are detected incorrectly due to non-detected peak (marked with a purple circle); volunteer #4, position P3.
Boxplots for PPV of S2 detection.

Figure 11 .
Figure 11.SE and PPV of HS detection in hybrid boxplots.

Figure 14 .
Figure 14.R-R interval, systole and diastole duration depending on heart rate 23 .
Block diagram of the DAQ system.

Table 2 .
General information about the tested subjects..

Table 3 .
Parameters of measured signals.
(a) Original PCG signal (red) and signal filtered using BPF (blue).

Table 4 .
SE and PPV parameters for S1 detection.Significant are in value[bold].The best performance of HS detection was reached on position P3, both in S1 and S2, as can be seen in Tables4, 5and Fig.11.Contrary, the worst detection quality was obtained from position P7 in both cases.Low values of SE and PPV in positions P1 and P2 are present sporadically, so these locations can be also considered suitable.More often false detection occurred on positions P4, P5, and P6.In one extreme case (subject #25), SE and PPV values of S2 detection approached zero.The S2 was less accurately detected also in subject #22, who suffers from heart disease.The both described cases arised due to decreased S2 amplitude during signal normalization, see Fig.12.

Table 5 .
SE and PPV parameters for S2 detection.Significant are in value [bold].

Table 6 .
Results of a Wilcoxon signed-rank test of HR obtained across subjects from P3; p-values lower than the selected significance level (0.05) are marked in bold.

Table 7 .
Comparison of results with the latest studies focused on PCG signal processing.