% EPHNOGRAM: A Simultaneous Electrocardiogram and Phonocardiogram Database
% A sample MATLAB script for reading the ECG-PCG data files and extracting the
% R-peaks from the ECG and the S1/S2 peaks of the PCG channels
%
% Dependencies: This script uses some functions from the Open-Source
% Electrophysiological Toolbox (OSET): https://github.com/alphanumericslab/OSET.git
%
% NOTE: This code has not been optimized for any specific dataset and is
% only provided as proof of concept and a starting point to work with the
% EPHNOGRAM database
%
% CITE:
% 1- The EPHNOGRAM on PhysioNet
% 2- A. Kazemnejad, P. Gordany, and R. Sameni. An open-access simultaneous electrocardiogram and phonocardiogram database. bioRxiv, 2021. DOI: https://doi.org/10.1101/2021.05.17.444563
% 3- R. Sameni, The Open-Source Electrophysiological Toolbox (OSET), v 3.14, URL: https://github.com/alphanumericslab/OSET
%
% By: Reza Sameni
% Email: reza.sameni@gmail.com
% May 2021
%
clear;
close all;
clc;
in_folder = '../MATfiles'; % .mat files folder
D = dir([in_folder '/*.mat']); % list of all mat files
% pre-processing parameters
fs = 1000.0; % post-decimation sampling frequency
for m = 1 : length(D) % NOTE: limit the index range to only analyze a subset of the records
% LOAD DATA
in_fname = D(m).name;
in_fname_full = [in_folder '/' in_fname];
dat = load(in_fname_full); % load the data file
% DECIMATION
% decimate the ECG and PCG channel to fs
ECG = decimate(dat.ECG, round(dat.fs/fs));
PCG = decimate(dat.PCG, round(dat.fs/fs));
% NOTCH FILTERING THE ECG
fc = 50.0; % powerline frequency
Qfactor = 45; % Q-factor of the notch filter
Wo = fc/(fs/2); BW = Wo/Qfactor; % nothc filter parameters
[b,a] = iirnotch(Wo, BW); % design the notch filter
ECG = filtfilt(b, a, ECG); % zero-phase non-causal filtering
% ECG PRE-PROCESSING
fl_ECG = 10.0; % lower bandpass cut-off frequency for R-peaks
fh_ECG = 40.0; % upper bandpass cut-off frequency for R-peaks
ECG_hp = ECG - LPFilter(ECG, fl_ECG/fs);
ECG_bp = LPFilter(ECG_hp, fh_ECG/fs);
sigma = 5.0 * std(ECG_bp); % saturate peaks above k-sigma of the ECG
ECG_saturated = sigma * tanh(ECG_bp / sigma);
if(skew(ECG_saturated) > 0) % find the ECG polarity based on skewness sign
ecg_polarity = 1;
else
ecg_polarity = 0;
end
% ECG R-PEAK DETECTION
% First round of R-peak detection (fixed window length peak search)
ff0 = 1.4; % approximate heart rate in Hz (just a rough estimate; will be fine-tuned by the algorithm)
peaks0 = PeakDetection(ECG_saturated, ff0/fs, ecg_polarity);
I0 = find(peaks0);
ff1 = median(fs ./ diff(I0));
% Second round of R-peak detection (adaptive window length peak search)
peaks_ECG = PeakDetectionAdaptiveHR(ECG_saturated, ff1, fs, ecg_polarity);
I_ECG_peaks = find(peaks_ECG);
RR_intervals_ecg = diff(I_ECG_peaks); % RR-interval time series in samples
ff_ecg = fs ./ RR_intervals_ecg;
HR_ecg = 60.* ff_ecg; % ECG-based heart rate
% PCG PRE-PROCESSING
fl_PCG = 10.0; % lower bandpass cut-off frequency for PCG peaks
fh_PCG = 100.0; % upper bandpass cut-off frequency for PCG peaks
PCG_hp = PCG - LPFilter(PCG, fl_PCG/fs);
PCG_bp = LPFilter(PCG_hp, fh_PCG/fs);
wlen = round(0.015 * fs); % PCG envelope detector window length
PCG_envelope = sqrt(filtfilt(ones(1, wlen), wlen, PCG_bp.^2));
% PCG PEAK DETECTION
search_wlen = round(0.15 * fs);
% find two dominant local peaks in the PCG envelope between successive
% ECG R-peaks (gives both S1 and S2)
num_peaks = 2;
peaks_S1S2_PCG = IntermediatePeakDetection(PCG_envelope, peaks_ECG, search_wlen, num_peaks, 1, 1);
% find the first dominant local peak in the PCG envelope between
% successive ECG R-peaks (gives S1 only)
num_peaks = 1;
peaks_S1_PCG = IntermediatePeakDetection(PCG_envelope, peaks_ECG, search_wlen, num_peaks, 1, 2);
I_PCG_S1_peaks = find(peaks_S1_PCG);
S1S1_intervals_pcg = diff(I_PCG_S1_peaks); % S1S1-interval time series in samples
ff_pcg = fs ./ S1S1_intervals_pcg;
HR_pcg_S1 = 60.* ff_pcg; % PCG-based heart rate based on S1S1 time intervals
peaks_S2_PCG = abs(peaks_S1S2_PCG - peaks_S1_PCG); % find the S2 peaks (exclude S1 from the two-peak time sequence)
I_PCG_S2_peaks = find(peaks_S2_PCG);
S2S2_intervals_pcg = diff(I_PCG_S2_peaks); % S2S2-interval time series in samples
ff_pcg_S2 = fs ./ S2S2_intervals_pcg;
HR_pcg_S2 = 60.* ff_pcg_S2; % PCG-based heart rate based on S2S2 time intervals
% INTERPOLATE THE HEART RATE SEQUENCES
t = (0 : length(ECG)-1)/fs;
% Add start and end points to the HR time sequences before interpolation
tt_ecg = t;
if(t(I_ECG_peaks(1)) > t(1))
tt_ecg = cat(2, t(1), tt_ecg(I_ECG_peaks));
RR_intervals_ecg = cat(2, RR_intervals_ecg(1), RR_intervals_ecg(1), RR_intervals_ecg);
end
if(t(I_ECG_peaks(end)) < t(end))
tt_ecg = cat(2, tt_ecg, t(end));
RR_intervals_ecg = cat(2, RR_intervals_ecg, RR_intervals_ecg(end));
end
tt_pcg_S1 = t;
if(t(I_PCG_S1_peaks(1)) > t(1))
tt_pcg_S1 = cat(2, t(1), tt_pcg_S1(I_PCG_S1_peaks));
S1S1_intervals_pcg = cat(2, S1S1_intervals_pcg(1), S1S1_intervals_pcg(1), S1S1_intervals_pcg);
end
if(t(I_PCG_S1_peaks(end)) < t(end))
tt_pcg_S1 = cat(2, tt_pcg_S1, t(end));
S1S1_intervals_pcg = cat(2, S1S1_intervals_pcg, S1S1_intervals_pcg(end));
end
tt_pcg_S2 = t;
if(t(I_PCG_S2_peaks(1)) > t(1))
tt_pcg_S2 = cat(2, t(1), tt_pcg_S2(I_PCG_S2_peaks));
S2S2_intervals_pcg = cat(2, S2S2_intervals_pcg(1), S2S2_intervals_pcg(1), S2S2_intervals_pcg);
end
if(t(I_PCG_S2_peaks(end)) < t(end))
tt_pcg_S2 = cat(2, tt_pcg_S2, t(end));
S2S2_intervals_pcg = cat(2, S2S2_intervals_pcg, S2S2_intervals_pcg(end));
end
% interpolate the HRs
RR_interval_ecg_interpolated = interp1(tt_ecg, RR_intervals_ecg, t);
S1S1_interval_pcg_interpolated = interp1(tt_pcg_S1, S1S1_intervals_pcg, t);
S2S2_interval_pcg_interpolated = interp1(tt_pcg_S2, S2S2_intervals_pcg, t);
RR_S1S1_diff = (RR_interval_ecg_interpolated - S1S1_interval_pcg_interpolated) / fs;
RR_S2S2_diff = (RR_interval_ecg_interpolated - S2S2_interval_pcg_interpolated) / fs;
beta = 0.05; % saturate differences above this amound of time (s)
RR_S1S1_diff_sat = beta * tanh(RR_S1S1_diff / beta);
RR_S2S2_diff_sat = beta * tanh(RR_S2S2_diff / beta);
HR_ecg_interpolated = interp1(t(I_ECG_peaks), [HR_ecg(1), HR_ecg], t);
HR_pcg_S1_interpolated = interp1(t(I_PCG_S1_peaks), [HR_pcg_S1(1), HR_pcg_S1], t);
HR_pcg_S2_interpolated = interp1(t(I_PCG_S2_peaks), [HR_pcg_S2(1), HR_pcg_S2], t);
HR_ecg_pcg_S1_diff = HR_ecg_interpolated - HR_pcg_S1_interpolated;
HR_ecg_pcg_S2_diff = HR_ecg_interpolated - HR_pcg_S2_interpolated;
alpha = 10.0; % saturate differences above this number of BPMs
HR_ecg_pcg_S1_diff_sat = alpha * tanh(HR_ecg_pcg_S1_diff / alpha);
HR_ecg_pcg_S2_diff_sat = alpha * tanh(HR_ecg_pcg_S2_diff / alpha);
% PLOT THE RESULTS
figure
lg = {};
subplot(511);
plot(t, ECG); lg = cat(2, lg, {'ECG'});
hold on
plot(t, ECG_saturated); lg = cat(2, lg, {'saturated ECG'});
plot(t, PCG); lg = cat(2, lg, {'PCG'});
plot(t, PCG_envelope); lg = cat(2, lg, {'PCG envelope'});
plot(t(I_ECG_peaks), ECG(I_ECG_peaks), 'ro'); lg = cat(2, lg, {'R-peaks'});
plot(t(I_PCG_S1_peaks), PCG_envelope(I_PCG_S1_peaks), 'gx'); lg = cat(2, lg, {'S1-peaks'});
plot(t(I_PCG_S2_peaks), PCG_envelope(I_PCG_S2_peaks), 'kx'); lg = cat(2, lg, {'S2-peaks'});
legend(lg);
grid
subplot(512);
lg = {};
plot(t(I_ECG_peaks), [HR_ecg(1), HR_ecg]); lg = cat(2, lg, {'ECG-based HR'});
hold on
plot(t(I_PCG_S1_peaks), [HR_pcg_S1(1), HR_pcg_S1]); lg = cat(2, lg, {'S1-based HR'});
plot(t(I_PCG_S2_peaks), [HR_pcg_S2(1), HR_pcg_S2]); lg = cat(2, lg, {'S2-based HR'});
legend(lg);
ylabel('BPM');
grid
subplot(513);
lg = {};
plot(t, HR_ecg_interpolated); lg = cat(2, lg, {'ECG-based HR'});
hold on
plot(t, HR_pcg_S1_interpolated); lg = cat(2, lg, {'S1-based HR'});
plot(t, HR_pcg_S2_interpolated); lg = cat(2, lg, {'S2-based HR'});
legend(lg)
grid
ylabel('Interpolated HRs (BPM)');
subplot(514);
lg = {};
plot(t, HR_ecg_pcg_S1_diff_sat); lg = cat(2, lg, {'RR-S1S1'});
hold on
plot(t, HR_ecg_pcg_S2_diff_sat); lg = cat(2, lg, {'RR-S2S2'});
legend(lg)
grid
ylabel('RR-SS HR diff saturated (BPM)');
subplot(515);
lg = {};
plot(t, 1000 * RR_S1S1_diff_sat); lg = cat(2, lg, {'RR-S1S1'});
hold on
plot(t, 1000 * RR_S2S2_diff_sat); lg = cat(2, lg, {'RR-S1S1'});
legend(lg)
grid
ylabel('RR-SS diff saturated (ms)');
% NOTE: Add a keyboard hit or pause here to check the results!
end