Цифровая обработка фотоплетизмограммы и анализ морфологии пульсовой волны

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

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

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

Библиотека CSV - библиотека для работы с файлами формата CSV. Применяется для чтения и загрузки экспериментальных данных из CSV-файлов, содержащих запись фотоплетизмограммы.
Библиотека DataFrames - библиотека для работы с табличными структурами данных.
Библиотека Statistics - библиотека, предоставляющая базовые статистические функции. Применяется для вычисления статистических характеристик данных, например средних значений, используемых при анализе параметров сигнала.

let
    installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
    list_packages = ["CSV","DataFrames", "Statistics"]
    for pack in list_packages
        pack in installed_packages || Pkg.add(pack)
    end
end

using CSV, DataFrames, Statistics

Обработка фотоплетизмограммы

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

engee.clear()

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

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

Выберите одну из команд, закомментировав вторую. Оба бэкенда являются взаимоисключающими, и активирован будет последний вызванный.

#gr()
plotlyjs()

Plots.PlotlyJSBackend()

В этом примере используется CSV-файл с данными фотоплетизмограммы и сопутствующих сигналов, полученными в положении сидя. Исходные записи взяты из открытой базы данных PhysioNet: Pulse Transit Time (PPG) Database.

Задаём путь к файлу с экспериментальными данными фотоплетизмограммы. Использование конструкции @__DIR__ позволяет формировать относительный путь к файлу, привязанный к директории текущего скрипта.

data_path = "$(@__DIR__)/s1_sit.csv"

"/user/DemoPublic/biomedical/ppg_signal_processing/s1_sit.csv"

На этом шаге выполняется загрузка данных из CSV-файла . Функция CSV.read считывает содержимое файла по заданному пути и преобразует его в структуру DataFrame.

data = CSV.read(data_path, DataFrame);

В переменной fs задаётся частота дискретизации сигналов, равная 500 Гц, в соответствии с описанием используемой базы данных. Далее определяется количество отсчётов сигнала фотоплетизмограммы и формируется временной вектор.

fs = 500.0  # частота дисретизации, Гц
N  = length(data.pleth_1)    # число отсчётов

t = (0:N-1) ./ fs   # вектор времени, с

0.0:0.002:508.05

Построим график сигнала фотоплетизмограммы во временной области.

plot(t, data.pleth_1,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    label = "pleth_1",
    grid = true
)

Для удобства анализа выберем временной фрагмент записи в интервале от 10 до 70 секунд. По заданным значениям времени рассчитываются соответствующие индексы отсчётов, после чего из исходного сигнала фотоплетизмограммы выделяется выбранный участок.

t_start = 10.0  # начало фрагмента, с
t_end = 70.0    # конец фрагмента, с

n_start = t_start *fs   # индекс начала фрагмента
n_end = t_end * fs  # индекс конца фрагмента

# выделение фрагмента сигнала ФПГ
x_raw = Float64.(-data.pleth_1[Int64(n_start):Int64(n_end)])

nsampl = length(x_raw)  # число отсчетв фрагмента
t =  range(t_start, step=1/fs, length=nsampl)   #   врменной вектор

10.0:0.002:70.0

Из ранее полученного временного представления анализируемого сигнала можно заметить наличие медленного изменения уровня сигнала. Такое поведение характерно для фотоплетизмографических записей и может быть связано с условиями регистрации.

Поэтому с использованием библиотеки EngeeDSP и функции detrend() из сигнала удаляется линейный тренд. Далее строится график сигнала после удаления тренда, что упрощает последующий анализ формы пульсовой волны

x = EngeeDSP.Functions.detrend(x_raw, 1)

plot(t, x, 
    label="ФПГ pleth_1 ",
    grid=true,
    xlabel="Время, с",
    ylabel="Амплитуда",
    xlims = (t_start,t_end)
)

Для анализа частотного состава сигнала используется библиотека EngeeDSP. С помощью функции fft() вычисляется быстрое преобразование Фурье анализируемого сигнала, после чего результат центрируется относительно нулевой частоты с использованием функции fftshift(). В результате формируется двухсторонний амплитудный спектр сигнала.

На основе частоты дискретизации и длины сигнала формируется вектор частот. Далее рассчитывается амплитудный спектр и строится его график в ограниченном частотном диапазоне

fft_1 = EngeeDSP.Functions.fft(x)
fftsh_1 = EngeeDSP.Functions.fftshift(fft_1)

df = fs / nsampl
freq_vec = -fs/2 : df : fs/2 - df

Amp1 = abs.(fftsh_1)./ nsampl

plot(freq_vec, Amp1,
    xlabel="Частота, Гц",
    ylabel="Амплитуда",
    title="Амплитудный спектр",
    grid=true,
    legend=false,
    xlims = (-100, 100)
)

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

Для этого в библиотеке EngeeDSP с помощью функции fir1() рассчитываются коэффициенты КИХ-фильтра с частотой среза fc = 30 Гц. Далее фильтрация сигнала выполняется функцией filter(). На графике приводится сравнение сигнала после удаления тренда и сигнала после фильтрации.

fc = 30.0   # частота среза, Гц
Wn = fc / (fs/2)    # нормированная частота среза

order = 20      # порядок КИХ фильтра
b = EngeeDSP.Functions.fir1(order, Wn)
a = [1.0]

y = EngeeDSP.Functions.filter(b, a, x)

p = plot(t, x,
        label="ФПГ",
         xlabel="Время, с",
         ylabel="Амплитуда",
         grid=true)
plot!(p, t, y, label="ФПГ после ФНЧ fc = $fc Гц")

xlims!(p, 10, 20)  
display(p)

После применения фильтра нижних частот заметим, что высокочастотные составляющие сигнала фотоплетизмограммы существенно подавляются. В результате уменьшается влияние шумов и мелких колебаний, присутствующих в исходном сигнале. Отфильтрованный сигнал становится более сглаженным, при этом основная форма пульсовой волны сохраняется. Далее проанализируем, как изменился частотный состав сигнала после фильтрации. Для этого с использованием библиотеки EngeeDSP вычислим спектр отфильтрованного сигнала с последующим центрированием. Затем рассчитаем двухсторонний амплитудный спектр и построим на одном графике спектры исходного сигнала и сигнала после фильтрации.

fft_2 = EngeeDSP.Functions.fft(y)
fftsh_2 = EngeeDSP.Functions.fftshift(fft_2)

Amp2 = abs.(fftsh_2)./ nsampl

plot(freq_vec, Amp1,
    xlabel="Частота, Гц",
    ylabel="Амплитуда",
    title="Амплитудный спектр",
    label = "Исходная ФПГ",
    grid=true,
    legend=true,
    xlims = (-100, 100)
)
plot!(freq_vec, Amp2,label = "ФПГ после ФНЧ")

Для поиска характерных точек пульсовой волны используется библиотека EngeeDSP и функция findpeaks(). Параметр minPeakDistance задаёт минимальное допустимое расстояние между соседними пиками в отсчётах и позволяет исключить повторное обнаружение пиков внутри одного сердечного цикла.

Сначала выполняется поиск локальных максимумов отфильтрованного сигнала фотоплетизмограммы. Затем локальные минимумы определяются по инвертированному сигналу.

minPeakDistance = 300 # минимальное расстояние между соседними пиками, отсчёты

# поиск локальных максимумов пульсовой волны 
pks_max, locs_max = EngeeDSP.Functions.findpeaks(
    y, out=:data,
    MinPeakDistance = minPeakDistance
)

# поиск локальных минимумов (по инвертированному сигналу)
pks_min, locs_min = EngeeDSP.Functions.findpeaks(
    -y, out=:data,
    MinPeakDistance = minPeakDistance
)

(Ypk = [17.584343668110932, -17.39531901266868, -18.669050763309713, -1.94542846968509, 33.74814704376922, 30.829367013334327, 25.84233903755211, 55.10538690758822, 69.8434389451998, 55.96626074426286  …  83.94463029754985, 84.82209067984081, 74.2801343867382, 88.52240598255462, 76.69576038233366, 84.17858772403335, 86.21498933712168, 97.85916644294964, 89.52541179013434, 88.90142660511627], Xpk = [140, 509, 863, 1242, 1651, 2067, 2464, 2871, 3305, 3721  …  26356, 26756, 27151, 27560, 27960, 28369, 28775, 29189, 29589, 29979], Wpk = [47.32674470024982, 148.10917709249253, 151.45931980024034, 159.00063126558462, 158.54474148630038, 196.83710033698867, 157.79178842043348, 170.28400609777964, 194.41276238628552, 166.93745546182754  …  195.92966275264916, 185.3148943352062, 189.9316956490511, 196.05989918619525, 192.93180929202936, 190.39849655982107, 189.35220067773844, 184.01907248867064, 191.21180624231783, 64.07862430432942], Ppk = [49.7927053863731, 119.12776439623973, 140.21705732007086, 142.2028851101346, 156.92264570424751, 148.04180549297627, 142.9548167946231, 156.68669000047524, 170.38020646509636, 146.51604671055836  …  180.93364223768032, 184.0924373176546, 173.06960737918646, 189.21725033818433, 159.4380899750749, 176.79651341589843, 185.6646479372822, 182.45101470253513, 170.93938020617196, 53.46875268894843])

Построим график фотоплетизмограммы во временной области и отметим на нём найденные локальные максимумы и минимумы.

p = plot(t, y,
        label="ФПГ pleth_1 ",
        grid=true,
        xlabel="Время, с",
        ylabel="Амплитуда"
)
scatter!(p, t[locs_max], y[locs_max], label="максимумы")
scatter!(p, t[locs_min], y[locs_min], label="минимумы")

xlims!(p, 10, 20)

Для расчёта частоты сердечных сокращений (ЧСС) необходимо определить временные интервалы между последовательными пульсовыми пиками. Для этого индексы найденных максимумов предварительно сортируются, после чего рассчитываются интервалы между соседними пиками с учётом частоты дискретизации сигнала. На основе среднего значения полученных интервалов вычисляется средняя частота сердечных сокращений.

# массив индексов максимумов
locs_max = sort(Int.(locs_max))

# вычисление временных интервалов между соседними максимумами, с
rr_sec = diff(locs_max) ./ fs

# расчёт средней частоты сердечных сокращений
hr_bpm = 60.0 / mean(rr_sec)

println("Средняя ЧСС: ", round(hr_bpm, digits=1), " уд/мин")

Средняя ЧСС: 74.5 уд/мин

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

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

Анализ морфологии пульсовой волны

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

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

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

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

Для практического анализа фотоплетизмограммы часто используются упрощённые морфологические параметры, которые могут быть надёжно вычислены автоматически:

высота систолической волны - разность амплитуд между минимумом и последующим максимумом;
длительность анакроты - временной интервал от минимума до максимума;
длительность катакроты - временной интервал от максимума до следующего минимума;
длительность пульсового цикла — интервал времени между двумя последовательными минимумами.

После цифровой обработки фотоплетизмограммы проведём анализ ранее рассмотренных морфологических признаков пульсовой волны. Для этого создадим массивы, предназначенные для хранения амплитудных и временных параметров отдельных пульсовых циклов.

Затем в цикле последовательно рассматриваются отдельные пульсовые циклы, каждый из которых определяется как участок сигнала между двумя соседними минимумами. Для каждого такого участка выделяется максимум пульсовой волны, после чего рассчитываются основные морфологические параметры: амплитуда пульсовой волны, длительности анакроты и катакроты, а также длительность пульсового цикла. Также выполняется проверка корректности расчета.

# массив индексов минимумов
locs_min = sort(Int.(locs_min))

heights = Float64[] # высота систолической волны
t_ana   = Float64[] # длительность анакроты, с
t_cat  = Float64[]  # длительность катакроты, с
t_cyc   = Float64[] # длительность сердечного цикла, с

idx_min  = Int[]    # индекс начала цикла
idx_max  = Int[]    # индекс максимума внутри цикла
idx_min2 = Int[]    # индекс конца цикла

for k in 1:(length(locs_min) - 1)

    i_min  = locs_min[k]    # индекс текущего минимума
    i_min2 = locs_min[k + 1]    # индекс следующего минимума

    # участок сигнала между двумя минимумами
    sig = y[i_min:i_min2]

    # индекс максимума на данном участке
    i_max = i_min - 1 + argmax(sig)

    # амплитуда пульсовой волны
    h = y[i_max] - y[i_min]

    # длительность анакроты (от минимума до максимума)
    ta = (i_max - i_min) / fs

    # длительность катакроты (от максимума до следующего минимума)
    tc = (i_min2 - i_max) / fs

    # длительность пульсового цикла
    tp = (i_min2 - i_min) / fs

    # проверка значений
    if h <= 0 || ta <= 0 || tc <= 0 || tp <= 0
        continue
    end

    push!(heights, h)
    push!(t_ana, ta)
    push!(t_cat, tc)
    push!(t_cyc, tp)

    push!(idx_min, i_min)
    push!(idx_max, i_max)
    push!(idx_min2, i_min2)
end

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

n = collect(1:length(heights))

cycles = DataFrame(
    Номер_цикла = n,
    Высота_сист_волы = heights,
    Анакрота_с = t_ana,
    Катакрота_с = t_cat,
    Длительность_цикла_с = t_cyc
)

Используемые материалы

В работе использовались реальные экспериментальные записи фотоплетизмограммы, полученные из открытого ресурса PhysioNet. В качестве исходных данных применялась база Pulse Transit Time (PPG) Database.

Заключение

В данном примере была рассмотрена цифровая обработка данных фотоплетизмограммы с использованием библиотеки EngeeDSP.

Работа примера включает следующие основные этапы:

Загрузка данных: Исходные данные фотоплетизмограммы были считаны из CSV-файла, содержащего реальные экспериментальные записи из открытой базы данных PhysioNet.
Предварительная обработка: С использованием функции EngeeDSP.Functions.detrend выполнено удаление линейного тренда сигнала. Для подавления высокочастотных колебаний применён фильтр нижних частот, спроектированный с помощью функции EngeeDSP.Functions.fir1 и реализованный через EngeeDSP.Functions.filter.
Спектральный анализ: С помощью функций EngeeDSP.Functions.fft и EngeeDSP.Functions.fftshift выполнен анализ частотного состава сигнала до и после фильтрации, что позволило наглядно оценить влияние фильтра на спектр фотоплетизмограммы.
Детекция характерных точек: Для выделения локальных максимумов и минимумов пульсовой волны использовалась функция EngeeDSP.Functions.findpeaks. Найденные характерные точки были визуализированы во временной области.
Расчёт частоты сердечных сокращений: На основе интервалов между последовательными пульсовыми пиками была рассчитана средняя частота сердечных сокращений для выбранного фрагмента записи.
Морфологический анализ пульсовой волны: Фотоплетизмограмма была сегментирована на отдельные пульсовые циклы. Для каждого цикла рассчитаны морфологические параметры, включая высоту систолической волны, длительности анакроты и катакроты, а также полную длительность пульсового цикла.
Представление результатов: Рассчитанные параметры были сведены в таблицу с использованием структуры DataFrame.

Row	Номер_цикла	Высота_сист_волы	Анакрота_с	Катакрота_с	Длительность_цикла_с
	Int64	Float64	Float64	Float64	Float64
1	1	154.107	0.124	0.614	0.738
2	2	141.491	0.12	0.588	0.708
3	3	148.613	0.114	0.644	0.758
4	4	142.203	0.128	0.69	0.818
5	5	150.961	0.118	0.714	0.832
6	6	154.004	0.122	0.672	0.794
7	7	142.955	0.12	0.694	0.814
8	8	156.687	0.142	0.726	0.868
9	9	160.393	0.104	0.728	0.832
10	10	156.503	0.118	0.694	0.812
11	11	159.537	0.126	0.7	0.826
12	12	155.849	0.112	0.678	0.79
13	13	151.217	0.13	0.666	0.796
⋮	⋮	⋮	⋮	⋮	⋮
63	63	173.556	0.106	0.74	0.846
64	64	186.218	0.104	0.736	0.84
65	65	187.327	0.136	0.678	0.814
66	66	183.215	0.116	0.684	0.8
67	67	185.517	0.124	0.666	0.79
68	68	173.07	0.128	0.69	0.818
69	69	171.265	0.114	0.686	0.8
70	70	169.314	0.124	0.694	0.818
71	71	187.603	0.116	0.696	0.812
72	72	185.665	0.14	0.688	0.828
73	73	179.273	0.116	0.684	0.8
74	74	174.117	0.118	0.662	0.78