Детектирование пиков ФПГ методом двух скользящих средних
Детектирование систолических пиков ФПГ с помощью двух скользящих средних и динамического порога
В демонстрационном примере рассматривается алгоритм детектирования систолических пиков фотоплетизмографического сигнала (ФПГ), предложенный M. Elgendi и соавторами. Метод основан на предварительной полосовой фильтрации, нелинейном преобразовании сигнала, расчёте двух скользящих средних с различной длительностью окна и формировании динамического порога.
В отличие от непосредственного поиска локальных максимумов, алгоритм сначала определяет временные области, соответствующие потенциальным пульсовым волнам, а затем находит систолический максимум внутри каждой принятой области.
Для демонстрации используется фрагмент записи из открытого набора MIMIC PERform AF Dataset. Рассматривается запись пациента без фибрилляции предсердий из файла mimic_perform_non_af_data.mat.
Настройка отображения
Функция engee.clear() выполняет очистку рабочего пространства:
engee.clear()
Для визуализации результатов доступны два режима отрисовки графиков:
- для статической быстрой отрисовки используется
gr(); - для интерактивной визуализации с возможностью масштабирования и просмотра значений используется
plotlyjs().
Выберите один из режимов, закомментировав второй. В данном примере используется интерактивный режим.
#gr()
plotlyjs()
Установка необходимых библиотек
Для загрузки исходных данных в формате MAT используется библиотека MAT. Библиотека Statistics применяется при формировании адаптивного порога вспомогательного алгоритма детектирования R-пиков ЭКГ.
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
Загрузка исходных данных
Файл mimic_perform_non_af_data.mat должен располагаться в одной папке с интерактивным скриптом. После загрузки отображается структура MAT-файла. Массив data содержит отдельные записи, каждая из которых включает ФПГ, ЭКГ и дополнительные физиологические сигналы.
data_path = "$(@__DIR__)/mimic_perform_non_af_data.mat"
mat_data = matread(data_path)
Для дальнейшей обработки используется 7-я запись из группы non-AF. Из неё извлекаются каналы ФПГ и ЭКГ. Частота дискретизации составляет 125 Гц.
Рассматривается 15-секундный фрагмент сигнала. Такая длительность позволяет одновременно наблюдать несколько последовательных пульсовых волн и проследить изменение их амплитуды во времени.
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("Частота дискретизации: ", fs, " Гц")
println("Длительность анализируемого фрагмента: ", duration, " с")
Исходный фотоплетизмографический сигнал
Сначала построим временное представление выбранного фрагмента ФПГ.
plot(t, ppg,
xlabel="Время, с",
ylabel="Амплитуда",
label="ФПГ",
title="Исходный фотоплетизмографический сигнал",
legend_position=:outertopright)
На выбранном участке пульсовые волны остаются различимыми, однако их амплитуда и форма изменяются. Поэтому для автоматического детектирования удобно использовать алгоритм, который формирует локальный изменяющийся порог и не опирается на одно фиксированное значение амплитуды.
Предварительная обработка ФПГ
Первый этап алгоритма – полосовая фильтрация исходной ФПГ. В работе Elgendi используется фильтр Баттерворта второго порядка с полосой пропускания от 0,5 до 8 Гц. Нижняя граница уменьшает влияние медленно изменяющейся составляющей, а верхняя ограничивает высокочастотные колебания.
Фильтрация выполняется функцией filtfilt() в прямом и обратном направлениях. Такой способ обработки позволяет избежать фазового смещения, что особенно важно при последующем определении временного положения систолических максимумов.
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="Время, с",
ylabel="Амплитуда",
label="Исходная ФПГ",
title="Полосовая фильтрация ФПГ",
legend_position=:outertopright)
plot!(p, t, ppg_filt, label="После фильтрации")
После полосовой фильтрации уменьшаются медленные изменения базового уровня и высокочастотные колебания. При этом основные пульсовые волны сохраняются и остаются локализованными во времени.
После фильтрации отрицательные значения сигнала обнуляются, а положительная часть возводится в квадрат:
Обнуление отрицательной части оставляет положительные участки отфильтрованной ФПГ, а возведение в квадрат увеличивает различие между выраженными пульсовыми волнами и малыми колебаниями.
ppg_pos = max.(ppg_filt, 0.0)
z = ppg_pos .^ 2
plot(t, z,
xlabel="Время, с",
ylabel="Амплитуда",
label="Преобразованный сигнал",
title="Нелинейное преобразование ФПГ",
legend_position=:outertopright)
После нелинейного преобразования основные положительные участки ФПГ становятся выраженными импульсными областями. Далее их временная структура анализируется с помощью двух скользящих средних.
Две скользящие средние
Ключевой этап алгоритма – расчёт двух скользящих средних с различной длительностью окна.
Первое скользящее среднее рассчитывается в коротком окне:
Короткое окно выделяет локальные области, соответствующие систолической части пульсовой волны.
Второе скользящее среднее рассчитывается в более длинном окне:
Длинное окно изменяется медленнее и характеризует локальный уровень преобразованного сигнала на масштабе сердечного цикла. Оно используется как основа для формирования динамического порога.
В исходной работе длительности окон составляют около 0,11 с и 0,67 с – точные оптимизированные значения равны 111 и 667 мс. При частоте дискретизации 125 Гц эти интервалы невозможно представить целым числом отсчётов точно, поэтому используются ближайшие нечётные размеры окон.
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("Окно MA_peak: ", w1, " отсчётов, ", round(1000*w1/fs, digits=0), " мс")
println("Окно MA_beat: ", w2, " отсчётов, ", round(1000*w2/fs, digits=0), " мс")
p = plot(t, z,
xlabel="Время, с",
ylabel="Амплитуда",
label="Преобразованный сигнал",
title="Скользящие средние преобразованного ФПГ-сигнала",
legend_position=:outertopright)
plot!(p, t, ma_peak, label="MA_peak")
plot!(p, t, ma_beat, label="MA_beat")
На графике MA_peak формирует локальные подъёмы около отдельных пульсовых волн. MA_beat изменяется значительно медленнее и задаёт локальный уровень, относительно которого далее оценивается наличие систолической области.
Таким образом, алгоритм использует два временных масштаба: короткий – для выделения возможной пульсовой волны, и длинный – для адаптации порога к текущему уровню сигнала.
Динамический порог и области интереса
Динамический порог формируется на основе длинного скользящего среднего:
где
Добавка задаёт небольшое смещение относительно MA_beat, а сам порог изменяется во времени вместе с длинным скользящим средним.
Область интереса формируется при выполнении условия
Последовательность обработки на этом этапе следующая:
- сравниваются
MA_peakи динамический порог; - по моментам пересечения определяются начало и конец области интереса;
- слишком короткие области отбрасываются;
- внутри каждой принятой области определяется максимальное по модулю значение отфильтрованной ФПГ;
- положение найденного максимума принимается за положение систолического пика.
Таким образом, максимальное по модулю значение ищется не по всей записи, а только внутри предварительно выделенной временной области.
beta = 0.02
alpha = beta * mean(z)
threshold = ma_beat .+ alpha
boi = ma_peak .> threshold
p = plot(t, ma_peak,
xlabel="Время, с",
ylabel="Амплитуда",
label="MA_peak",
title="Динамический порог",
legend_position=:outertopright)
plot!(p, t, threshold, label="Динамический порог")
На участках, где MA_peak превышает динамический порог, формируются кандидаты на пульсовые волны. Поскольку порог следует за медленными изменениями MA_beat, критерий детектирования адаптируется к изменению локального уровня сигнала.
edges = diff(Int.(vcat(false, boi, false))) # Определение моментов входа в область интереса и выхода из неё
block_start = findall(edges .== 1) # Индексы начала областей интереса
block_end = findall(edges .== -1) .- 1 # Индексы конца областей интереса
peak_idx = Int[] # Индексы найденных систолических пиков
for (i_start, i_end) in zip(block_start, block_end)
if i_end - i_start + 1 >= w1 # Отбрасывание слишком коротких областей
local_idx = argmax(abs.(ppg_filt[i_start:i_end])) # Поиск максимума по модулю внутри области
push!(peak_idx, i_start + local_idx - 1) # Перевод локального индекса в индекс всего фрагмента
end
end
Результат детектирования
Отобразим найденные систолические пики на исходной ФПГ.
p = plot(t, ppg,
xlabel="Время, с",
ylabel="Амплитуда",
label="ФПГ",
title="Детектирование систолических пиков",
legend_position=:outertopright)
scatter!(p, t[peak_idx], ppg[peak_idx], label="Систолические пики", markersize=4)
Маркеры располагаются в максимумах выделенных пульсовых волн. При этом положение каждого систолического пика определяется только после формирования области интереса с помощью двух скользящих средних и динамического порога.
Сопоставление ЭКГ и ФПГ
Для наглядного сопоставления сердечных событий рассмотрим ЭКГ из того же временного фрагмента. R-пики ЭКГ определяются вспомогательным алгоритмом и используются только как визуальные ориентиры.
Сначала ЭКГ фильтруется в полосе 5–20 Гц. Далее рассчитывается разность соседних отсчётов, результат возводится в квадрат и сглаживается скользящим средним. Такой сигнал подчёркивает области QRS-комплексов. По адаптивному порогу определяются кандидаты, после чего положение R-пика уточняется поиском наиболее выраженного экстремума ЭКГ в небольшой окрестности каждого кандидата.
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));
ЭКГ и найденные R-пики
Сначала отобразим результат вспомогательной детекции R-пиков на ЭКГ.
p = plot(t, ecg,
xlabel="Время, с",
ylabel="Амплитуда",
label="ЭКГ",
title="Детектирование R-пиков ЭКГ",
legend_position=:outertopright)
scatter!(p, t[r_idx], ecg[r_idx], label="R-пики", markersize=4)
Найденные маркеры располагаются в области выраженных комплексов QRS и используются далее только для визуального сопоставления с пульсовыми волнами ФПГ.
Совместное отображение ЭКГ и ФПГ
Для сопоставления построим ЭКГ и ФПГ на общей временной оси. На верхнем графике отметим R-пики, на нижнем – систолические пики ФПГ.
p_ecg = plot(t, ecg, ylabel="Амплитуда ЭКГ", label="ЭКГ", title="ЭКГ", legend_position=:outertopright)
scatter!(p_ecg, t[r_idx], ecg[r_idx], label="R-пики", markersize=4)
p_ppg = plot(t, ppg, xlabel="Время, с", ylabel="Амплитуда ФПГ", label="ФПГ", title="ФПГ", legend_position=:outertopright)
scatter!(p_ppg, t[peak_idx], ppg[peak_idx], label="Систолические пики", markersize=4)
plot(p_ecg, p_ppg, layout=(2, 1), link=:x, size=(1000, 650))
Совместное представление позволяет визуально сопоставить последовательность электрических сердечных событий на ЭКГ и найденных пульсовых волн ФПГ. Для каждого канала используется собственный алгоритм выделения характерных точек, поэтому график предназначен именно для качественного сопоставления последовательностей событий.
Абсолютное временное расстояние между R-пиком ЭКГ и систолическим пиком ФПГ в данном примере не используется для расчёта времени распространения пульсовой волны.
Заключение
В демонстрационном примере реализован алгоритм детектирования систолических пиков ФПГ на основе двух скользящих средних и динамического порога.
На этапе предварительной обработки выполнена полосовая нуль-фазовая фильтрация и нелинейное преобразование сигнала. Далее функцией movmean() рассчитаны две скользящие средние: короткая используется для выделения локальных областей пульсовых волн, а длинная для формирования изменяющегося во времени порога. По пересечениям MA_peak и динамического порога определены области интереса, внутри которых по максимальному по модулю значению отфильтрованной ФПГ определены положения систолических пиков.
Дополнительно для того же временного фрагмента выполнена вспомогательная детекция R-пиков ЭКГ. Совместное отображение ЭКГ и ФПГ позволяет визуально сопоставить последовательность сердечных событий и обнаруженных пульсовых волн без расчёта частоты сердечных сокращений и других производных показателей.
Использованные материалы
-
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