Detection of systolic FPG peaks using two moving averages and a dynamic threshold
The demonstration example discusses the algorithm for detecting systolic peaks of a photoplethysmographic signal (PPG), proposed by M. Elgendi and co-authors. The method is based on preliminary bandpass filtering, nonlinear signal transformation, calculation of two moving averages with different window lengths and formation of a dynamic threshold.
Unlike the direct search for local maxima, the algorithm first determines the time regions corresponding to potential pulse waves, and then finds the systolic maximum within each accepted region.
A recording fragment from the open set MIMIC PERform AF Dataset is used for the demonstration. The recording of a patient without atrial fibrillation from a file is being considered mimic_perform_non_af_data.mat.
Display Settings
Function engee.clear() cleans up the workspace:
engee.clear()
Two graph rendering modes are available to visualize the results.:
- for static fast rendering, it is used
gr(); - for interactive visualization with the ability to scale and view values , it is used
plotlyjs().
Select one of the modes by commenting out the second one. This example uses interactive mode.
# gr()
plotlyjs()
Installing the necessary libraries
The library is used to download the source data in MAT format. MAT. Library Statistics It is used in the formation of the adaptive threshold of the auxiliary algorithm for detecting R-peaks of the ECG.
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["MAT", "Statistics"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
using MAT, Statistics
Loading the source data
File mimic_perform_non_af_data.mat it must be located in the same folder as the interactive script. After downloading, the structure of the MAT file is displayed. Array data It contains separate recordings, each of which includes PHG, ECG, and additional physiological signals.
data_path = "$(@__DIR__)/mimic_perform_non_af_data.mat"
mat_data = matread(data_path)
The 7th entry from the non-AF group is used for further processing. The FPG and ECG channels are extracted from it. The sampling rate is 125 Hz.
A 15-second fragment of the signal is being considered. This duration allows you to simultaneously observe several consecutive pulse waves and track the change in their amplitude over time.
data = mat_data["data"]
subject_no = 7
subject = data[subject_no]
ppg_all = vec(Float64.(subject["ppg"]["v"]))
ecg_all = vec(Float64.(subject["ekg"]["v"]))
fs = 125.0
start_sample = 60000
duration = 15.0
N = Int(fs * duration)
end_sample = start_sample + N - 1
ppg = ppg_all[start_sample:end_sample]
ecg = ecg_all[start_sample:end_sample]
t = (0:N-1) ./ fs
println("Sampling rate: ", fs, " Hz")
println("Duration of the analyzed fragment: ", duration, " with")
The initial photoplethysmographic signal
First, let's build a temporary representation of the selected fragment of the FPG.
plot(t, ppg,
xlabel="Time, from",
ylabel="The amplitude",
label="FPG",
title="The initial photoplethysmographic signal",
legend_position=:outertopright)
The pulse waves remain distinguishable in the selected area, but their amplitude and shape change. Therefore, for automatic detection, it is convenient to use an algorithm that generates a local changing threshold and does not rely on a single fixed amplitude value.
FPG pre-treatment
The first stage of the algorithm is bandpass filtering of the initial FPG. Elgendi uses a second-order Butterworth filter with a bandwidth from 0.5 to 8 Hz. The lower limit reduces the influence of a slowly changing component, while the upper limit limits high-frequency fluctuations.
Filtering is performed by the function filtfilt() in the forward and reverse directions. This treatment method avoids phase shift, which is especially important in the subsequent determination of the temporal position of systolic maxima.
f_low = 0.5
f_high = 8.0
order = 2
Wn = [f_low, f_high] ./ (fs / 2)
b, a = EngeeDSP.Functions.butter(order, Wn, "bandpass")
ppg_filt = EngeeDSP.Functions.filtfilt(b, a, ppg)
p = plot(t, ppg,
xlabel="Time, from",
ylabel="The amplitude",
label="The original FPG",
title="FPG bandpass filtering",
legend_position=:outertopright)
plot!(p, t, ppg_filt, label="After filtering")
After bandpass filtering, slow changes in the baseline and high-frequency fluctuations decrease. At the same time, the main pulse waves persist and remain localized in time.
After filtering, the negative values of the signal are reset to zero, and the positive part is squared.:
Zeroing out the negative part leaves positive areas of the filtered FPG, and squaring increases the difference between pronounced pulse waves and small fluctuations.
ppg_pos = max.(ppg_filt, 0.0)
z = ppg_pos .^ 2
plot(t, z,
xlabel="Time, from",
ylabel="The amplitude",
label="The converted signal",
title="Nonlinear FPG transformation",
legend_position=:outertopright)
After the nonlinear transformation, the main positive sections of the FPG become pronounced pulse regions. Further, their time structure is analyzed using two moving averages.
Two moving averages
The key stage of the algorithm is the calculation of two moving averages with different window lengths.
The first moving average is calculated in a short window:
A short window highlights local areas corresponding to the systolic part of the pulse wave.
The second moving average is calculated over a longer window.:
The long window changes more slowly and characterizes the local level of the transformed signal on the scale of the cardiac cycle. It is used as the basis for the formation of a dynamic threshold.
In the initial operation, the window durations are about 0.11 s and 0.67 s – the exact optimized values are 111 and 667 ms. At a sampling frequency of 125 Hz, these intervals cannot be accurately represented by an integer number of samples, so the nearest odd window sizes are used.
T_peak = 0.111
T_beat = 0.667
w1 = 2 * floor(Int, T_peak * fs / 2) + 1
w2 = 2 * floor(Int, T_beat * fs / 2) + 1
ma_peak = EngeeDSP.Functions.movmean(z, w1)
ma_beat = EngeeDSP.Functions.movmean(z, w2)
println("The MA_peak window: ", w1, " counts, ", round(1000*w1/fs, digits=0), " ms")
println("The MA_beat window: ", w2, " counts, ", round(1000*w2/fs, digits=0), " ms")
p = plot(t, z,
xlabel="Time, from",
ylabel="The amplitude",
label="The converted signal",
title="Moving averages of the transformed FPG signal",
legend_position=:outertopright)
plot!(p, t, ma_peak, label="MA_peak")
plot!(p, t, ma_beat, label="MA_beat")
On the chart MA_peak It forms local rises near individual pulse waves. MA_beat it changes much more slowly and sets a local level relative to which the presence of a systolic region is further assessed.
Thus, the algorithm uses two time scales: a short one to identify a possible pulse wave, and a long one to adapt the threshold to the current signal level.
Dynamic threshold and areas of interest
The dynamic threshold is formed based on a long moving average:
where
Supplement sets a small offset relative to MA_beat, and the threshold itself changes over time along with the long moving average.
The area of interest is formed when the condition is met
The processing sequence at this stage is as follows:
- compared
MA_peakand the dynamic threshold; - The beginning and the end of the area of interest are determined by the points of intersection.;
- Areas that are too short are discarded.;
- Within each accepted area, the maximum modulo value of the filtered FPG is determined.;
- The position of the found maximum is taken as the position of the systolic peak.
Thus, the maximum modulo value is not searched for the entire record, but only within a pre-allocated time domain.
beta = 0.02
alpha = beta * mean(z)
threshold = ma_beat .+ alpha
boi = ma_peak .> threshold
p = plot(t, ma_peak,
xlabel="Time, from",
ylabel="The amplitude",
label="MA_peak",
title="Dynamic threshold",
legend_position=:outertopright)
plot!(p, t, threshold, label="Dynamic threshold")
In areas where MA_peak If it exceeds the dynamic threshold, candidates for pulse waves are formed. Because the threshold follows slow changes MA_beat, the detection criterion adapts to changes in the local signal level.
edges = diff(Int.(vcat(false, boi, false))) # Determining the moments of entry into and exit from the area of interest
block_start = findall(edges .== 1) # Indexes of the beginning of the areas of interest
block_end = findall(edges .== -1) .- 1 # Indexes of the end of the areas of interest
peak_idx = Int[] # Indices of found systolic peaks
for (i_start, i_end) in zip(block_start, block_end)
if i_end - i_start + 1 >= w1 # Discarding areas that are too short
local_idx = argmax(abs.(ppg_filt[i_start:i_end])) # Searching for the maximum modulo inside the area
push!(peak_idx, i_start + local_idx - 1) # Converting a local index to an index of the entire fragment
end
end
Detection result
Let's display the found systolic peaks on the initial PHG.
p = plot(t, ppg,
xlabel="Time, from",
ylabel="The amplitude",
label="FPG",
title="Detection of systolic peaks",
legend_position=:outertopright)
scatter!(p, t[peak_idx], ppg[peak_idx], label="Systolic peaks", markersize=4)
The markers are located at the maxima of the selected pulse waves. In this case, the position of each systolic peak is determined only after the formation of the area of interest using two moving averages and a dynamic threshold.
Comparison of ECG and AFG
For a visual comparison of cardiac events, consider an ECG from the same time fragment. The R-peaks of the ECG are determined by an auxiliary algorithm and are used only as visual landmarks.
First, the ECG is filtered in the 5-20 Hz band. Next, the difference between neighboring samples is calculated, the result is squared and smoothed by a moving average. This signal highlights the areas of the QRS complexes. Candidates are determined by the adaptive threshold, after which the position of the R-peak is refined by searching for the most pronounced ECG extremum in a small neighborhood of each candidate.
Wn_ecg = [5.0, 20.0] ./ (fs / 2)
b_ecg, a_ecg = EngeeDSP.Functions.butter(2, Wn_ecg, "bandpass")
ecg_filt = EngeeDSP.Functions.filtfilt(b_ecg, a_ecg, ecg)
d_ecg = vcat(0.0, diff(ecg_filt))
sq_ecg = d_ecg .^ 2
w_ecg = 2 * floor(Int, 0.080 * fs / 2) + 1
ecg_det = EngeeDSP.Functions.movmean(sq_ecg, w_ecg)
thr_ecg = Statistics.mean(ecg_det) + 1.5 * Statistics.std(ecg_det)
min_r_distance = Int(round(0.30 * fs))
_, qrs_idx = EngeeDSP.Functions.findpeaks(ecg_det, out=:data, MinPeakHeight=thr_ecg, MinPeakDistance=min_r_distance)
search_radius = Int(round(0.080 * fs))
r_idx = Int[]
for idx in qrs_idx
i1 = max(1, idx - search_radius)
i2 = min(N, idx + search_radius)
local_idx = argmax(abs.(ecg_filt[i1:i2]))
push!(r_idx, i1 + local_idx - 1)
end
r_idx = unique(sort(r_idx));
ECG and detected R-peaks
First, we will display the result of the auxiliary detection of R-peaks on the ECG.
p = plot(t, ecg,
xlabel="Time, from",
ylabel="The amplitude",
label="ECG",
title="Detection of ECG R-peaks",
legend_position=:outertopright)
scatter!(p, t[r_idx], ecg[r_idx], label="R-peaks", markersize=4)
The markers found are located in the area of pronounced QRS complexes and are then used only for visual comparison with the pulse waves of the FPG.
Joint display of ECG and FPG
For comparison, we will plot the ECG and FPG on a common time axis. We note the R peaks on the upper graph, and the systolic peaks of PHG on the lower graph.
p_ecg = plot(t, ecg, ylabel="ECG amplitude", label="ECG", title="ECG", legend_position=:outertopright)
scatter!(p_ecg, t[r_idx], ecg[r_idx], label="R-peaks", markersize=4)
p_ppg = plot(t, ppg, xlabel="Time, from", ylabel="FPG amplitude", label="FPG", title="FPG", legend_position=:outertopright)
scatter!(p_ppg, t[peak_idx], ppg[peak_idx], label="Systolic peaks", markersize=4)
plot(p_ecg, p_ppg, layout=(2, 1), link=:x, size=(1000, 650))
The combined representation allows you to visually compare the sequence of electrical cardiac events on the ECG and the detected pulse waves of the PHG. Each channel uses its own algorithm for identifying characteristic points, so the graph is designed specifically for qualitative comparison of sequences of events.
The absolute time distance between the R-peak of the ECG and the systolic peak of the PHG in this example is not used to calculate the propagation time of the pulse wave.
Conclusion
The demonstration example implements an algorithm for detecting systolic peaks of FPG based on two moving averages and a dynamic threshold.
At the preprocessing stage, band-pass zero-phase filtering and nonlinear signal conversion are performed. Next, the function movmean() Two moving averages are calculated: the short one is used to highlight local areas of pulse waves, and the long one is used to form a time-varying threshold. By intersections MA_peak Areas of interest are determined within which the positions of systolic peaks are determined based on the maximum modulo value of the filtered FPG.
Additionally, auxiliary detection of ECG R-peaks was performed for the same time fragment. The combined display of ECG and PHG allows you to visually compare the sequence of cardiac events and detected pulse waves without calculating heart rate and other derived indicators.
Materials used
-
Elgendi M., Norton I., Brearley M., Abbott D., Schuurmans D. Systolic Peak Detection in Acceleration Photoplethysmograms Measured from Emergency Responders in Tropical Conditions // PLOS ONE. – 2013. – Vol. 8, No. 10. DOI: 10.1371/journal.pone.0076585
-
Charlton P. H. MIMIC PERform Datasets // Zenodo. DOI: 10.5281/zenodo.6807403