Моделирование электрокардиосигнала и определение ЧСС

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

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

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

engee.clear()

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

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

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

#gr()
plotlyjs()

Plots.PlotlyJSBackend()

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

В примере используется библиотека Statistics, необходимая для расчёта среднего значения и стандартного отклонения RR-интервалов. Перед выполнением статистического анализа проверяется наличие библиотеки в рабочем окружении.

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

using Statistics

Загрузка и запуск модели

Представленная блочная модель реализует полный тракт формирования и обработки электрокардиосигнала. На первом этапе формируются три стандартных отведения I, II и III. Далее к сигналу добавляются моделируемые помехи измерительного тракта, после чего выполняется цифровая фильтрация. Для II отведения после предварительной обработки выполняется детекция R-пиков, расчёт RR-интервалов и определение частоты сердечных сокращений.

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

Формирование ЭКС – формирование модельного электрокардиосигнала и трёх стандартных отведений;
Моделирование помех – добавление сетевых составляющих и широкополосного шума;
Режекторный фильтр – подавление сетевой составляющей 50 Гц;
ФНЧ – ограничение высокочастотной составляющей сигнала;
Определение ЧСС – детекция R-пиков, вычисление RR-интервалов и частоты сердечных сокращений.

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

modelName = "ECG_HR_Model"
modelPath = joinpath(@__DIR__, modelName * ".engee")

model = engee.load(modelPath; force=true);

Формирование электрокардиосигнала

Подсистема Формирование ЭКС формирует один модельный сердечный цикл как сумму отдельных компонент P, Q, R1, R2, S, ST и T. Такой подход основан на представлении морфологии электрокардиосигнала набором параметризованных несимметричных гауссовых функций и соответствует общей идее примера сообщества Engee «Моделирование электрокардиосигнала».

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

где определяет амплитуду компоненты, – положение её экстремума внутри сердечного цикла, а определяет ширину.

Для получения несимметричной формы ширина задаётся отдельно слева и справа от экстремума:

Суммарная форма одного сердечного цикла имеет вид

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

Период сердечного цикла определяется через частоту сердечных сокращений:

где задаётся в уд/мин, а получается в секундах.

После формирования базовой формы строятся условные потенциалы электродов RA, LA и LL. Стандартные отведения рассчитываются как

Из этих выражений непосредственно следует соотношение Эйнтховена

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

а затем рассчитывается соответствующий интервал . Параметр задаётся в уд/мин. Значение по умолчанию равно 3 уд/мин.

В первом запуске используется постоянная частота сердечных сокращений 75 уд/мин.

Запуск модели с постоянной частотой сердечных сокращений

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

results = engee.run(model; verbose=true);

Building...
Progress 0%
Progress 6%
Progress 20%
Progress 26%
Progress 31%
Progress 37%
Progress 43%
Progress 49%
Progress 55%
Progress 60%
Progress 79%
Progress 84%
Progress 90%
Progress 95%
Progress 100%
Progress 100%

Сформированные стандартные отведения

Отобразим отведения I, II и III на фрагменте модельного сигнала.

src = results["ЭКС"]

t = Float64.(src.time)
lead_I   = Float64.(getindex.(src.value, 1))
lead_II  = Float64.(getindex.(src.value, 2))
lead_III = Float64.(getindex.(src.value, 3))

idx = (t .>= 2.0) .& (t .<= 4.0)

p = plot(t[idx], lead_I[idx],
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    label="I",
    title="Стандартные отведения",
    legend_position=:outertopright)

plot!(p, t[idx], lead_II[idx], label="II")
plot!(p, t[idx], lead_III[idx], label="III")

Полученные кривые содержат одинаковую последовательность P–QRS–T, но отличаются амплитудой вследствие различного способа формирования отведений. Наиболее выраженный R-зубец наблюдается во II отведении, которое далее используется для определения R-пиков.

Моделирование помех

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

Сетевая помеха связана с электромагнитными наводками от сети питания на электроды, соединительные провода и входные каскады измерительной аппаратуры. В модели учитывается основная составляющая 50 Гц и две гармоники 100 и 150 Гц:

В модели заданы

Дополнительно используется ограниченный по полосе белый шум, моделирующий совокупную широкополосную составляющую измерительного тракта. Для блока Band-Limited White Noise задан параметр Cov = 0.2·10⁻⁶.

Суммарный сигнал на выходе подсистемы можно представить в виде

где – исходный электрокардиосигнал, а – шумовая составляющая.

Результат моделирования помех

Для оценки влияния помех сравним исходное и искажённое II отведение.

noisy = results["ЭКС с помехами"]

t_noisy = Float64.(noisy.time)
lead_II_noisy = Float64.(getindex.(noisy.value, 2))

idx_n = (t_noisy .>= 2.0) .& (t_noisy .<= 4.0)

p = plot(t[idx], lead_II[idx],
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    label="Исходное II отведение",
    title="Влияние моделируемых помех",
    legend_position=:outertopright)

plot!(p, t_noisy[idx_n], lead_II_noisy[idx_n],label="II отведение с помехами")

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

Предварительная обработка электрокардиосигнала

Цель предварительной обработки – ослабить сетевые гармоники и широкополосный шум, сохранив форму комплекса QRS, необходимую для последующей детекции R-пиков.

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

Режекторный фильтр

Для подавления основной сетевой составляющей 50 Гц используется БИХ-режекторный фильтр Чебышёва I типа. В текущей модели результирующий порядок фильтра равен 10, граничные частоты составляют 46 и 54 Гц, а допустимая пульсация в полосах пропускания – 0.5 дБ. Частота дискретизации модели равна 1000 Гц.

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

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

частота дискретизации: 1000 Гц;
тип амплитудно-частотной характеристики: режекторный фильтр;
тип фильтра: БИХ;
метод синтеза: Чебышёв I типа;
результирующий порядок: 10;
граничные частоты: 46 и 54 Гц;
пульсация в полосах пропускания: 0.5 дБ.

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

Результат режекторной фильтрации

notch = results["Режекторный фильтр.Y"]

t_notch = Float64.(notch.time)
lead_II_notch = Float64.(getindex.(notch.value, 2))

idx_f = (t_notch .>= 2.0) .& (t_notch .<= 4.0)

p = plot(t_noisy[idx_n], lead_II_noisy[idx_n],
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    label="До фильтра",
    title="Режекторная фильтрация",
    legend_position=:outertopright)

plot!(p, t_notch[idx_f], lead_II_notch[idx_f],label="После режекторного фильтра")

Режекторный фильтр ослабляет узкополосную составляющую в области 50 Гц. При этом гармоники 100 и 150 Гц и широкополосный шум требуют дополнительного ограничения полосы сигнала.

Фильтр нижних частот

По аналогии синтезирован БИХ-фильтр нижних частот Баттерворта. В модели используется порядок 10, частота среза 38.13 Гц и частота дискретизации 1000 Гц.

Фильтр Баттерворта выбран из-за гладкой амплитудно-частотной характеристики в полосе пропускания. Он дополнительно подавляет гармоники 100 и 150 Гц и уменьшает высокочастотную составляющую белого шума. Частота среза выбрана выше основной информативной низкочастотной области модельного электрокардиосигнала и позволяет сохранить выраженный комплекс QRS.

Результат предварительной обработки

filt = results["ФНЧ.Y"]

t_filt = Float64.(filt.time)
lead_II_filt = Float64.(getindex.(filt.value, 2))

idx_lpf = (t_filt .>= 2.0) .& (t_filt .<= 4.0)

p = plot(t_noisy[idx_n], lead_II_noisy[idx_n],
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    label="С помехами",
    title="Результат предварительной обработки",
    legend_position=:outertopright)

plot!(p, t_filt[idx_lpf], lead_II_filt[idx_lpf],label="После фильтрации")

После последовательной фильтрации высокочастотные колебания уменьшаются, а форма P–QRS–T остаётся различимой. Полученный сигнал используется как вход алгоритма определения R-пиков.

Детекция R-пиков и расчёт частоты сердечных сокращений

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

Сначала рассчитывается разность соседних отсчётов:

Разностная операция выделяет участки с высокой скоростью изменения сигнала. Комплекс QRS имеет более крутые фронты по сравнению с зубцами P и T, поэтому его вклад после этой операции выражен сильнее.

Затем полученный сигнал возводится в квадрат:

Возведение в квадрат устраняет знак разности и увеличивает вклад отсчётов с большой скоростью изменения.

Для сглаживания рассчитывается скользящее среднее по 10 отсчётам. При частоте дискретизации 1000 Гц длина окна соответствует 10 мс:

Полученный сигнал используется для выделения области QRS-комплекса.

Порог детекции формируется адаптивно. На окне из 1000 отсчётов, соответствующем 1 с, рассчитывается локальное среднее:

а также среднее квадрата:

Стандартное отклонение рассчитывается по текущим статистическим оценкам:

После этого определяется адаптивный порог:

Область возможного QRS-комплекса определяется условием:

После выделения области QRS положение R-пика уточняется по II отведению после фильтрации. Для этого определяется локальный максимум. Для трёх последовательных отсчётов проверяются условия:

R-кандидат принимается только при одновременном выполнении двух условий: сигнал детекции превышает адаптивный порог, а в II отведении после фильтрации обнаружен локальный максимум.

После принятого R-пика включается рефрактерный интервал длительностью 0.30 с. В течение этого интервала новые R-кандидаты не принимаются, что предотвращает повторную регистрацию одного комплекса QRS.

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

Для двух последовательных детектированных R-пиков рассчитывается RR-интервал:

после чего определяется частота сердечных сокращений:

До появления второго корректно детектированного R-пика RR-интервал и рассчитанная частота сердечных сокращений остаются равными нулю.

Общая последовательность обработки соответствует подходу, рассмотренному в примере сообщества Engee «Детекция R-пиков на ЭКГ сигнале», при этом в рассматриваемой модели детекция выполняется потоково с использованием блоков среды моделирования.

Найденные R-пики

r_sig = results["R-пики"]

t_r = Float64.(r_sig.time)
r_event = Float64.(r_sig.value)

r_idx = findall(r_event .> 0.5)

idx_r = (t_filt .>= 2.0) .& (t_filt .<= 6.0)
r_idx_plot = [i for i in r_idx if 2.0 <= t_r[i] <= 6.0]

p = plot(t_filt[idx_r], lead_II_filt[idx_r],
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    label="II отведение",
    title="Детекция R-пиков",
    legend_position=:outertopright)

scatter!(p,t_r[r_idx_plot],lead_II_filt[r_idx_plot],label="R-пики",markersize=4)

Маркеры располагаются в области максимумов R-зубцов. Для каждого сердечного цикла формируется один детектированный R-пик, а рефрактерный интервал исключает повторную регистрацию одного комплекса QRS.

RR-интервалы и частота сердечных сокращений

hr_model_sig = results["ЧСС модельная"]
hr_measured_sig = results["ЧСС измеренная"]

t_hr_model = Float64.(hr_model_sig.time)
hr_model = Float64.(hr_model_sig.value)

t_hr_measured = Float64.(hr_measured_sig.time)
hr_measured = Float64.(hr_measured_sig.value)

p = plot(t_hr_model, hr_model,
    xlabel="Время, с",
    ylabel="ЧСС, уд/мин",
    label="Заданная ЧСС",
    title="Заданная и измеренная ЧСС",
    legend_position=:outertopright)

plot!(p, t_hr_measured, hr_measured,label="Измеренная ЧСС")

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

Дополнительное временное смещение связано с задержкой, вносимой последовательной цифровой фильтрацией и операциями детекции. После завершения начального переходного участка измеренное значение устанавливается около заданного уровня 75 уд/мин.

Моделирование вариабельности частоты сердечных сокращений

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

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

engee.set_param!(modelName * "/Формирование ЭКС",
    "HR" => 75,
    "HR_VAR" => true,
    "HR_SIGMA" => 3.0)

engee.set_param!(modelName, "StopTime" => 180.0)

results_var = engee.run(model; verbose=true);

Building...
Progress 0%
Progress 5%
Progress 10%
Progress 15%
Progress 20%
Progress 25%
Progress 30%
Progress 35%
Progress 40%
Progress 45%
Progress 50%
Progress 55%
Progress 60%
Progress 65%
Progress 70%
Progress 75%
Progress 80%
Progress 86%
Progress 91%
Progress 96%
Progress 100%
Progress 100%

Модельная и измеренная частота сердечных сокращений

hr_model_var_sig = results_var["ЧСС модельная"]
hr_measured_var_sig = results_var["ЧСС измеренная"]

t_hm = Float64.(hr_model_var_sig.time)
hr_m = Float64.(hr_model_var_sig.value)

t_hd = Float64.(hr_measured_var_sig.time)
hr_d = Float64.(hr_measured_var_sig.value)

p = plot(t_hm, hr_m,
    xlabel="Время, с",
    ylabel="ЧСС, уд/мин",
    label="Модельная ЧСС",
    title="Модельная и измеренная ЧСС",
    legend_position=:outertopright)

plot!(p, t_hd, hr_d,label="Измеренная ЧСС")

Измеренная частота сердечных сокращений изменяется вслед за модельным значением, однако моменты обновления двух зависимостей не совпадают во времени. Модельное значение задаётся в начале сердечного цикла, тогда как измеренная частота сердечных сокращений может быть определена только после детекции следующего R-пика и расчёта завершившегося RR-интервала.

Дополнительное временное смещение вносит тракт предварительной цифровой обработки и последовательность операций детекции R-пика. Поэтому наблюдаемое смещение является следствием структуры измерительного алгоритма и не свидетельствует об ошибке расчёта RR-интервала.

Анализ вариабельности RR-интервалов

На основе последовательности детектированных R-пиков сформируем набор RR-интервалов и выполним их анализ. Используем три представления:

RR-тахограмму – для отображения изменения длительности сердечных циклов во времени;
гистограмму RR-интервалов – для оценки распределения их длительности;
скаттерограмму – для анализа взаимосвязи соседних интервалов и .

Аналогичный набор представлений используется в примере сообщества Engee «Анализ вариабельности сердечного ритма по ЭКГ сигналу».

r_var_sig = results_var["R-пики"]
rr_var_sig = results_var["RR-интервал"]

r_event_var = Float64.(r_var_sig.value)
t_r_var = Float64.(r_var_sig.time)
rr_hold = Float64.(rr_var_sig.value)

r_idx_var = findall(r_event_var .> 0.5)

rr_beats = rr_hold[r_idx_var]
rr_time = t_r_var[r_idx_var]

valid_rr = (rr_beats .> 0.3) .& (rr_beats .< 2.0)

rr_beats = rr_beats[valid_rr]
rr_time = rr_time[valid_rr]
rr_ms = 1000 .* rr_beats
heart_rate_rr = 60.0 ./ rr_beats

mean_rr = mean(rr_ms)
std_rr = std(rr_ms)
mean_hr = mean(heart_rate_rr)

println("Средний RR-интервал: ", round(mean_rr, digits=2), " мс")
println("Стандартное отклонение RR-интервалов: ", round(std_rr, digits=2), " мс")
println("Средняя ЧСС: ", round(mean_hr, digits=2), " уд/мин")

Средний RR-интервал: 802.87 мс
Стандартное отклонение RR-интервалов: 33.02 мс
Средняя ЧСС: 74.86 уд/мин

RR-тахограмма

RR-тахограмма показывает изменение длительности последовательных сердечных циклов во времени. По этой зависимости можно оценить величину и характер моделируемых изменений RR-интервалов.

plot(rr_time, rr_ms,
    xlabel="Время, с",
    ylabel="RR-интервал, мс",
    title="RR-тахограмма",
    legend=false)

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

Гистограмма RR-интервалов

Гистограмма отображает распределение длительности сердечных циклов. По оси X откладываются RR-интервалы в миллисекундах, а по оси Y – количество интервалов в соответствующем диапазоне.

histogram(rr_ms,
    bins=30,
    xlabel="RR-интервал, мс",
    ylabel="Количество интервалов",
    title="Гистограмма RR-интервалов",
    legend=false)

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

Скаттерограмма RR-интервалов

Скаттерограмма показывает взаимосвязь между соседними сердечными циклами. По оси X откладывается текущий интервал , а по оси Y – следующий интервал . Диагональная линия соответствует равенству соседних RR-интервалов.

rr_n = rr_ms[1:end-1]
rr_next = rr_ms[2:end]

p = scatter(rr_n, rr_next,
    xlabel="RR(n), мс",
    ylabel="RR(n+1), мс",
    title="Скаттерограмма RR-интервалов",
    legend=false,
    markersize=3,
    xlim=(0, 1200),
    ylim=(0, 1200),
    aspect_ratio=:equal)

plot!(p, [0, 1200], [0, 1200],linestyle=:dash,linewidth=2)

При небольшой моделируемой вариабельности точки располагаются компактно вблизи диагональной линии. Увеличение разброса по обеим осям соответствует более выраженным различиям между соседними RR-интервалами.

Заключение

В примере рассмотрена единая блочная цепочка формирования и обработки модельного электрокардиосигнала. Параметрическая модель формирует согласованные отведения I, II и III, после чего к сигналу добавляются сетевые гармоники 50, 100 и 150 Гц и широкополосный шум. Режекторный фильтр и фильтр нижних частот выполняют предварительную обработку сигнала перед детекцией R-пиков.

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

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

Моделирование электрокардиосигнала и определение ЧСС