Детектирование систолических пиков ФПГ с помощью двух скользящих средних и динамического порога

В демонстрационном примере рассматривается алгоритм детектирования систолических пиков фотоплетизмографического сигнала (ФПГ), предложенный M. Elgendi и соавторами. Метод основан на предварительной полосовой фильтрации, нелинейном преобразовании сигнала, расчёте двух скользящих средних с различной длительностью окна и формировании динамического порога.

В отличие от непосредственного поиска локальных максимумов, алгоритм сначала определяет временные области, соответствующие потенциальным пульсовым волнам, а затем находит систолический максимум внутри каждой принятой области.

Для демонстрации используется фрагмент записи из открытого набора MIMIC PERform AF Dataset. Рассматривается запись пациента без фибрилляции предсердий из файла mimic_perform_non_af_data.mat.

Настройка отображения

Функция engee.clear() выполняет очистку рабочего пространства:

engee.clear()

Для визуализации результатов доступны два режима отрисовки графиков:

для статической быстрой отрисовки используется gr();
для интерактивной визуализации с возможностью масштабирования и просмотра значений используется plotlyjs().

Выберите один из режимов, закомментировав второй. В данном примере используется интерактивный режим.

#gr()
plotlyjs()

Plots.PlotlyJSBackend()

Установка необходимых библиотек

Для загрузки исходных данных в формате 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)

Dict{String, Any} with 3 entries:
  "source"  => Dict{String, Any}("matlab_conversion_script"=>"/Users/petercharl…
  "data"    => MatlabStructArray{2}(["ppg", "ekg", "imp", "abp", "fix"], Matrix…
  "license" => Dict{String, Any}("details"=>"This dataset is licensed under the…

Для дальнейшей обработки используется 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, " с")

Частота дискретизации: 125.0 Гц
Длительность анализируемого фрагмента: 15.0 с

Исходный фотоплетизмографический сигнал

Сначала построим временное представление выбранного фрагмента ФПГ.

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), " мс")

Окно MA_peak: 13 отсчётов, 104.0 мс
Окно MA_beat: 83 отсчётов, 664.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