An open-access simultaneous electrocardiogram and phonocardiogram database

The 30 s records were recorded during the development phase of the device as the participants were seated in an armchair during their acquisition.

The 30 min recordings (60 recordings) were acquired in an indoor sports center. A structured interview was conducted to confirm that the participants were in good physical condition and none reported symptoms of autonomic or cardiovascular disorder. In preparation for the experiment, volunteers refrained from consuming food, caffeine, alcohol, or smoking for three hours prior to the test. But they were permitted to drink water regularly. Each subject participated in one or a number of the following physical scenarios:

• (i)  
  Scenario A—rest condition: the participant laid horizontally on a bed in a quiet room while the ECG and PCG were recorded for 30 min. A few short samples (30 s) have also been recorded from participants sitting in an armchair.
• (ii)  
  Scenario B—walking condition: the participant walked on a treadmill at a constant speed of 3.7 km h⁻¹ at an incline grade of 1%, where incline grade is:

  $$\begin{eqnarray}&&\mathrm{Grade}\ ( \% )=\left(\displaystyle \frac{\mathrm{vertical}\ \mathrm{rise}}{\mathrm{horizontal}\ \mathrm{run}}\right)\times 100.\end{eqnarray} \tag{ 1 }$$

  The process lasted 30 min at constant speed and incline grade.
• (iii)  
  Scenario C—treadmill stress-test: the modified Bruce protocol was used for the treadmill stress-test (Bruce 1971), as shown in table 1. Our modified stress-test starts at a lower workload as compared to the standard test. Each session lasted 30 min and the increase in speed and treadmill incline grade continued until the subjects reached excessive fatigue, excessive heart rate, or chest pain. Whenever the subject reached this extreme point (subject-dependent), in order to avoid a sudden decrease in the heart rate, the speed was gradually decreased to 3 km h⁻¹, and the incline grade was decreased to a horizontal position. After a few min of walking in this state, the participant stopped and sat in a chair until the end of the test. The signals were acquired up to the end of the acquisition session.
• (iv)  
  Scenario D—bicycle stress-test: the bicycle stress-test protocol is detailed in table 2. Each test lasted 30 min. The participant first rested for 2 min on a stationary exercise bicycle without pedaling, then started pedaling at a workload of 25 Watts min⁻¹. At the beginning of the test, according to table 2, the external load was gradually increased, and the subject was asked to match the power consumption to the power specified in the standard test. It is practically infeasible to pedal at a constant speed (at least for non-athletes). Therefore, the power consumption is only an approximation of the target power consumption. The participant continueed the test until reaching excessive fatigue, excessive heart rate, or chest pain. At this point, the external load level was gradually decreased to 25 Watts min⁻¹ to avoid a sudden decrease in the heart rate. The subject remained in this state for several min, and then stopped and rested on the stationary bicycle until the end of the test session, while the signals were continuously recorded by the EPCG device.

In total, ten subjects participated in the first scenario. Five of the participants of the first scenario did not attend the rest of the stages. With the addition of six new volunteers, a total number of eleven subjects contributed in each of the three other scenarios. Two volunteers failed the test due to physical fatigue in the bicycle exercise stress test (Scenario D).

The EPHNOGRAM database is available online on PhysioNet (Kazemnejad et al 2021). The data files are provided in both MATLAB (ECGPCG00XY.mat where XY = 01, ..., 68) and WFDB (ECGPCG00XY.dat and ECGPCG00XY.hea) formats, with identical base names. The MATLAB files are in double-precision floating point format. Each file was converted into 16 bit WFDB format by using the mat2wfdb.m function from the WFDB Toolbox (Silva and Moody 2014, Moody et al 2022). The accuracy of conversion between MATLAB and WFDB formats was assessed per file and per channel, by comparing the signal-to-noise ratio (SNR) of the original double-precision floating-point MATLAB files versus the WFDB files read by the rdsamp function of WFDB. All 16 bit WFDB files had an SNR of above 60 dB per channel, as compared to the original MATLAB files. Although 60 dB is fully acceptable for most applications, researchers seeking double-precision floating-point accuracy may use the MATLAB files for higher precision.

The description of the corresponding physical activities and the unique IDs of the participants are listed in the spreadsheet ECGPCGSpreadsheet.csv, provided in the database. For basic heart rate extraction and analysis from the ECG and PCG channels a sample MATLAB script TestHeartRateCalculation.m is also provided online. Additional source codes for analyzing this data are available in the OSET (Sameni 2018).

For the ECG channels, a cascade of a moving median filter of length 0.6 s, and a moving average filter of length 0.3 s are used to remove the baseline wander. This combination has been shown to be very effective and robust for ECG baseline wander removal (Sameni et al 2008, Jamshidian-Tehrani and Sameni 2018). As an acoustic signal, the PCG does not physically have a baseline drift, and any non-zero average is associated with electronic circuitry and may be safely removed with a DC blocker, without interfering with the PCG contents. However, to keep the ECG and PCG synchronous, the delay of the two-stage baseline wander filter in the ECG signal path (which is half the total window lengths, or 0.45 s) was added to the PCG signal path, by zero-padding the beginning of the PCG signals during preprocessing ⁵ .

Since the data have been acquired indoors, the power-line noise (50 Hz in this dataset) was inevitable and appeared in variable amplitude throughout the 30 min acquisition sessions, in some cases appearing as short burst noises. In the current study, a fixed second-order IIR notch filter was used to suppress the power-line noise. For more elaborate studies involving ECG morphology analysis, more advanced notch filters such as adaptive Kalman notch filters can be used (Sameni 2012), which adjust the notch filter's quality factor (Q-factor) over time based on the level of power-line interference. This adaptive approach is particularly beneficial for addressing varying levels of power-line interference across different recordings, ensuring more effective noise cancellation. The source code for the Kalman notch filter is online and available in OSET (Sameni 2018).

We used an ECG R-peak detector inspired by the Pan-Tompkins algorithm (Pan and Tompkins 1985). Accordingly, a bandpass FIR filter with a passband between 10 and 40 Hz was applied to the ECG. The filtered signal amplitude was next saturated by a hyperbolic tangent function: $y=\alpha \tanh (x/\alpha )$, with α set to be k times the standard deviation of bandpass filtered signal (k = 10, for the later reported results). The saturation ensures that short-term burst/spike noises do not impact the subsequent peak detection and thresholding steps of the algorithm. Next, the filtered signal's power envelope was calculated over a sliding window of length 75 ms. Finally, the R-peaks were obtained by local peak detection of the power envelopes over windows of 0.5 s. Finally, we built a graphical tool using MATLAB's Signal Labeler app to visualize and tweak the automatically detected R-peaks. A biomedical engineer specializing in ECG analysis used this tool for beat-by-beat review of all the automatically detected R-peaks, making corrections to any missed or inaccurately detected R-peaks. The human overseen and corrected R-peak indexes are also provided the EPHNOGRAM database.

PCG-based beat detection and annotation is more challenging and less-established than ECG annotation. One of the advantages of simultaneous ECG-PCG acquisition and processing is that the ECG-based R-peaks can be used as references for beat detection and segmentation of the PCG components (S1, S2, etc), which are otherwise not trivial to detect. For this study, the ECG R-peaks were first estimated from the ECG and used as initial reference points for estimating the S1 and S2 segments of the PCG, from the PCG power envelope as detailed below.

Typical simultaneous ECG and PCG records from the dataset are shown in figure 5(a). The first and second heart sounds in the PCG, respectively, S1 and S2 components, occur between successive R-peaks. From figure 5(a) that the S1 component is a PCG peak after the R-peak, and the S2 component is a PCG peak after the ECG T-wave, which matches the electrophysiology of the heart. In order to accurately delineate the S1 and S2 waves, we employed wavedet_3D—a wavelet-based algorithm that estimates the ECG fiducial points, like the peaks, onsets and offsets of the key ECG components, including the T-waves (Martínez et al 2004). The function is available online in the open-source ECG-KIT (Llamedo et al 2016). To note, wavedet_3D accepts pre-detected R-peaks as input, or internally detects the R-peaks if they are not provided. For robustness and consistency, we provided the function with the pre-calculated R-peaks that were human overseen by an expert, as detailed in section 4.2.

Zoom In Zoom Out Reset image size

In order to make the S1 and S2 components more prominent and easier to detect, the PCG was filtered with an FIR bandpass with a 20–200 Hz passband, which aligns with the literature that the dominant frequency components of the PCG are below 200 Hz (Mubarak et al 2018). This property was also confirmed through a visual inspection of the spectrograms of our PCG data. After applying this filter, the S1 and S2 components of the PCG manifest as bumpy shapes modulated over narrow-band oscillatory waves. Therefore, unlike the ECG, the local peaks of the PCG are not reliable indicators of the PCG component positions. Instead, we used the Hilbert transform envelope to facilitate the detection of the S1 and S2 components. A typical example of a simultaneous ECG and PCG records and their corresponding spectrograms are shown in figure 5(a). The S1 peak was next detected as the dominant PCG power envelope peak between the corresponding ECG beat's R- and T-wave peaks. Similarly, the S2 peak was identified as the dominant peak of the PCG envelope between the corresponding T-peak and the following R-peak. Similar to the R-peaks, all the detected S1 and S2 components of the PCG were overseen by a human annotator using the Signal Labeler app.

Sample time-series of ECG and PCG-based time intervals are shown in figure 6, for one of the stationary biking tests and one of the treadmill stress tests of the dataset. The left vertical axis in this figure represents the time intervals in milliseconds (ms), and the right vertical axis indicates the equivalent heart rates in beats per minute (bpm). The background shades represent different levels of workload. As expected, the ECG-based, S1-based, and S2-based heart rates have similar trends, with their time-series graphs being superimposed. Upon zooming on the time-series of pedaling on a stationary bicycle in figure 7, it becomes apparent that S2-based heart rate has stronger high-frequency components and randomness compared to ECG-based and S1-based heart rates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 6, we also see that during low-intensity workloads and rest periods, the heart rate exhibits more variability compared to the high-intensity exercise period. Another interesting observation is that heart rate variability (HRV) is more prominent in the S2–S1 (diastolic) interval as compared with the S1–S2, R-S1, or R-S2 intervals. According to figure 6, among the intra-beat intervals S1–S2, S2–S1, R-S1, and R-S2, the R-S1 interval exhibits the least variability during exercise, while the S2–S1 interval demonstrates the most variability. This suggests that the diastolic period, derived from the mechanical behavior of the heart, is the primary parameter that adjusts to the variations in the heart rate at different levels of workload.

There are of course additional observations regarding the heart rate fluctuations in figure 6 that require modeling and statistical analysis. For example, the zoomed-in HR fluctuations in figure 7, implies that S2-based HR has higher frequency fluctuations than ECG-based HR and S1-based HR. To assess this hypothesis, we can calculate the power spectrum and entropy to measure the frequency components and uncertainty in the different HR time-series. The Lomb-Scargle normalized periodogram is a standard method for measuring the power spectrum of HRV time-series, such as R–R intervals (Estévez et al 2016). It is particularly useful for analyzing HRV in cases where the data is irregularly sampled, noisy, or non-stationary (Stewart et al 2020). Figures 8(a) and (b) present the spectral analysis of R–R, S1–S1, and S2–S2 intervals from the sample records from figure 6, under various workload scenarios. Accordingly, the S2-based HR has stronger high-frequency components than ECG-based and S1-based HR, specifically in the higher frequency range (0.2–1 Hz).

Zoom In Zoom Out Reset image size

As a sample nonlinear analysis, the approximate entropy (ApEn) of the HRV time-series was calculated for a range of tolerance parameters (r) to measure uncertainty and randomness (Richman and Moorman 2000). Intuitively, low ApEn values indicate that the system is very persistent, repetitive, and predictive, with apparent patterns that repeat throughout the time-series while higher values reflect independence between the data, fewer repeated patterns, and randomness (Delgado-Bonal and Marshak 2019). The parameter r sets a similarity threshold between data segments, determining whether two sequences are similar based on their differences. Typically, r is a fraction of the data's standard deviation. It highlights pattern regularity, where lower ApEn values suggest predictability and repetition, indicating consistent patterns. In contrast, higher ApEn values imply increased randomness and fewer repeated patterns. A plot of ApEn versus r demonstrates how variations in r affect the system's complexity, specifically showing the point where the signal transitions from predictable to more random patterns, thereby providing an understanding of the signal's structure relative to r.

The ApEn results of the samples of figure 6 are shown in figures 8(c) and 8(d). According to figure 5(a), the amplitude and power of S2 waves are weaker than those of S1. The SNR indeed impacts our precision in detecting these waves, as illustrated in figure 7, where the power spectral density (PSD) for S2 exhibits slightly higher frequency components compared to S1. This is also observed in figures 8(a) and (b). An alternative hypothesis is that there might be more intrinsic randomness in S2, than S1, which is beyond the scope of our current study.

Although the impact of SNR on approximate entropy is less evident, our findings in figures 8(c) and (d) show that ApEn peaks at a slightly higher tolerance parameter r for S2 than for S1, indicating more randomness in the S2 time series. This is consistent with the observed higher entropy in S2-based HRV compared to both ECG-based and S1-based HRVs, showing increased randomness in the S2-based HRV.

We employed the Pearson correlation coefficient, a widely recognized measure of coupling, to investigate the relationship between ECG and PCG-based indices as proxies for the interactions between the electrical and mechanical activities of the heart. The index was calculated for the R–R and R-S1 intervals of all recordings of the dataset, during different physical activities. We compared two categories of workloads: low-intensity activities including lying, sitting, and slow walk (at 3.7 km h⁻¹ speed), and high-intensity activities including compound activity (slow walk, fast walk, sit-to-stand, followed by slow walk), treadmill test, and biking. Figure 9 presents the box plots of the calculated correlation coefficients for these categories. Accordingly, the coupling between R–R and R-S1 is significantly stronger in the high-intensity category compared to the low-intensity category (p-value <10⁻⁸). Specifically, in low-intensity workloads, the R-S1 interval exhibited a very small correlation with the R–R interval, with an average correlation coefficient of 0.04. However, in high-intensity workloads characterized by rapid heart rates, the R-S1 is positively correlated with R–R interval, with an average correlation coefficient of 0.65 across all records.

Zoom In Zoom Out Reset image size

We can also study the relationship between the cardiac electromechanical time intervals and the different physical workloads. As proof of concept, we present the results for subject S009, from whom we have 30 min recordings in three different scenarios: rest (sitting on an armchair), rest (lying on a bed), and exercise (walking at different speeds) conditions. To study the significance of the different ECG-PCG time intervals, we studied the feature importance in classifying the data segments corresponding to sitting, lying, and walking from the combination of the ECG and PCG-based features. The 30 min recordings were divided into 30 s intervals with a sliding window of 5 s, resulting in a total of 1074 segments. Various features were extracted from R–R, S1–S2, S2–S1, R-S1, and R-S2 intervals, including HRV features: normal-to-normal beat intervals mean (NNmean) and standardard deviation (SDNN), and acceleration/deceleration capacity values using phase-rectified signal averaging (Kantelhardt et al 2007), resulting in a total of 20 features for every 30 s segment. The open-source PhysioNet-Cardiovascular-Signal-Toolbox was used to extract these features (Vest et al 2018, 2019).

To evaluate the relative importance of each feature, we used the minimal-redundancy-maximal-relevance (mRMR) criterion to select and rank the features (Peng et al 2005). mRMR is a well-established information-theoretical method for feature analysis, which has also been used for algorithmic heart disease diagnosis (Wang et al 2022). We used the function fscmrmr in MATLAB, setting the parameter corresponding to the prior probabilities for each class to 'empirical', which determines class probabilities from the class frequencies in the class labels. In total, 1074 thirty-second segments were used in the mRMR analysis. Table 3 shows the 10 most significant features out of a total of 20 features that selected and ranked based on mRMR. In this table, RS1-NNmean and RS2-NNmean denote the average of R-peak to S1 and S2 wave intervals of normal beats, respectively. RS2-acceleration and RS2-deceleration denote the acceleration/deceleration capacity values of the R-peak to S2-wave interval.

Accordingly, the most significant feature is RS2-NNmean, which is derived from electro-mechanical intervals. This shows the importance of simultaneous ECG and PCG recording to extract hybrid features. The scatter plot in figure 10 displays the first three features selected by mRMR in three levels of workloads. It shows a clear separation between the physical activities, particularly between sitting and lying.

Zoom In Zoom Out Reset image size

We further used the SHAP (SHapley Additive exPlanations) method to analyze the feature importances. The SHAP method calculates local feature importance for each data point, as opposed to traditional methods which calculate the importance of each feature globally (Lundberg and Lee 2017, Lundberg et al 2020). For feature analysis by the SHAP method, three binary random forest classifiers were used based on the one-vs-rest (OvR) method and applied to the 1074 thirty-second records. The beeswarm plot in figure 11 summarizes all SHAP values for the top 10 features based on their SHAP values. Accordingly, the electromechanical intervals rank higher than the R–R intervals in terms of feature importance. Specifically, figures 11(a) and (b) confirm that RS2-NNmean is the most significant feature for classifying sitting vs other classes and lying vs other classes, in the OvR method. This is consistent with the results of mRMR in figure 10, which identified RS2-NNmean as the most important feature.

Zoom In Zoom Out Reset image size

In summary, both the mRMR and SHAP methods confirm the significance of inter-modal ECG and PCG features R-S1 and R-S2 intervals during different physical activities.

While our study significantly contributes to the understanding of heart rate dynamics and electromechanical coupling, future research should extend these insights to more diverse populations, including both normal and pathological conditions of different sex, age and race groups.

Advancements in EPCG hardware could facilitate the integration of this technology into wearable cardiac monitoring devices. Such development could revolutionize real-time, non-invasive cardiac health assessments for outpatients in everyday life.

The methodological aspects of this research offer several avenues for expansion. Investigating the mutual information between ECG and PCG modalities, along with their unique contributions, is one such direction that can lead into multimodal cardiac monitoring (Lahat et al 2015). Utilizing a combination of time and frequency domain features, we anticipate uncovering detailed mutual and exclusive information between these two modalities. This could prove instrumental in examining micro-variations between the heart's electrical and mechanical mechanisms. Such findings are crucial not only for machine learning applications but also for understanding the cardiac function and its variation with different physical activity levels. This underscores the importance of simultaneous and multimodal heart monitoring for more accurate cardiac diagnosis.

Future research should also consider implementing active noise cancellation algorithms to enhance PCG data quality. As indicated in the EPHNOGRAM data description, certain channels (PCG2, AUX1, and AUX2) exhibit very low amplitude, often at the level of quantization noise. Despite their limited direct utility, these channels, which capture electronic device noises and faint environmental sounds, could be valuable for developing adaptive noise cancellers or multichannel blind and semi-blind source separation algorithms.

Additionally, the dataset could be leveraged to create mathematical ECG-PCG models and digital twins using simultaneous ECG and PCG channels. This approach builds on the past developments in electrocardiography, such as synthetic ECG generators (McSharry 2003), which facilitated the creation of algorithms for Kalman filtering and parameter extraction from noisy ECG recordings (Sameni et al 2007, 2008). Similarly, the present dataset could assist in developing algorithms for denoising and automatic parameter extraction from PCG, or simultaneous ECG and PCG recordings.