Детекция R-пиков на ЭКГ сигнале

В проекте рассматривается обработка ЭКГ сигнала для детекции R-пиков. Выполняется загрузка сигнала, фильтрация сетевой помехи, подготовка сигнала к детекции и отображение найденных R-пиков на ЭКГ сигнале.

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

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

engee.clear()

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

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

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

#gr()
plotlyjs()

Plots.PlotlyJSBackend()

Подключим с помощью функции include файл "read.jl" для чтения файлов:

include("$(@__DIR__)/read.jl")

safe_read_ann (generic function with 1 method)

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

На данном этапе выполняется проверка наличия необходимых библиотек Julia и их автоматическая установка при отсутствии в рабочем окружении. В рассматриваемом примере используется библиотека 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

Подключается библиотека EngeeDSP и импортируются функции, необходимые для цифровой обработки ЭКГ сигнала, включая фильтрацию, спектральный анализ, поиск пиков и выполнение свертки.

import EngeeDSP.Functions: fft, butter, filter, findpeaks, conv, freqz

Обработка и анализ ЭКГ сигнала

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

Функция расчета спектра сигнала

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

function calc_fft(x, fs)
    len = length(x)
    S = fft(x)  # ДПФ
    
    # Односторонний спектр
    K = len ÷ 2
    S = S[1:K]      
    freq = (0:K-1) .* fs ./ len     # частотная ось от 0 до fs/2

    # Односторонний амплитудный спектр
    S_amp = 2 .* abs.(S) ./ len
    S_amp[1] = 0.0

    # Односторонний амплитудный спектр в дБ
    S_dB = 20 .* log10.(S_amp .+ 1e-12)

    return freq, S_amp, S_dB
end

calc_fft (generic function with 1 method)

Загрузка и подготовка данных

Выполняется загрузка выбранного фрагмента ЭКГ сигнала из базы данных MIT-BIH. На основе информации из заголовочного файла определяется частота дискретизации записи, после чего задаются начало и длительность анализируемого участка сигнала. Далее производится чтение соответствующего временного фрагмента ЭКГ записи и выбор первого канала сигнала, соответствующего второму стандартному отведению ЭКГ, которое широко используется при анализе сердечного ритма и детекции R-пиков.

data = "$(@__DIR__)/mitdb/100"

# Считывание заголовочного файла и получение частоты дискретизации
hdr = read_header(data)
fs = Float64(hdr.fs)

# Выбор фрагмента сигнала для анализа
t_start = 5.0 * 60.0        # начало фрагмента - с 5 мин.
t_dur   = 5.0  * 60.0      # длительность фрагмента - 5 мин.

# Перевод времени из секунд в номера отсчетов
sampfrom = Int(round(t_start * fs))
sampto   = sampfrom + Int(round(t_dur * fs))

# Чтение выбранного фрагмента ЭКГ сигнала
_, raw, phys = read_dat_212(data; sampfrom=sampfrom, sampto=sampto)

# Выбор канала отведения для анализа
ecg = phys === nothing ? Float64.(raw[:,1]) : Float64.(phys[:,1])

# Расчет количества отсчетов и формирование временной оси
N = length(ecg)
t = collect(0:N-1) ./ fs

108000-element Vector{Float64}:
   0.0
   0.002777777777777778
   0.005555555555555556
   0.008333333333333333
   0.011111111111111112
   0.013888888888888888
   0.016666666666666666
   0.019444444444444445
   0.022222222222222223
   0.025
   0.027777777777777776
   0.030555555555555555
   0.03333333333333333
   ⋮
 299.96666666666664
 299.96944444444443
 299.97222222222223
 299.975
 299.97777777777776
 299.98055555555555
 299.98333333333335
 299.9861111111111
 299.9888888888889
 299.9916666666667
 299.99444444444447
 299.9972222222222

Отображение исходного ЭКГ сигнала

Отображается исходный фрагмент ЭКГ сигнала во временной области.

p1 = plot(t, ecg,
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    label="Исходный сигнал"
)

Спектральный анализ исходного ЭКГ сигнала

С использованием ранее созданной функции calc_fft выполняется расчет спектра ЭКГ сигнала и построение одностороннего амплитудного спектра. Полученный график позволяет оценить распределение амплитуд сигнала по частотам и определить основные спектральные компоненты ЭКГ сигнала.

freq0, S0_amp, _ = calc_fft(ecg, fs)

p2 = plot(freq0, S0_amp,
    xlabel="Частота, Гц",
    ylabel="Амплитуда, мВ",
    label="Исходный сигнал",
    xlim=(0, fs/2)
)

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

Фильтрация ЭКГ сигнала

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

f_center = 60.0     # резонансная частота, Гц
bw = 6.0        # полоса пропускания фильтра, Гц
f1 = f_center - bw/2      # нижняя граничная частота, Гц
f2 = f_center + bw/2      # верхняя граничная частота, Гц

# Нормировання частота относительно частоты дискретизации
wn1 = f1 / (fs / 2)
wn2 = f2 / (fs / 2)

order = 6      # порядок проектируемого фильтра

b, a = butter(order, [wn1, wn2], "stop")    # синтез режекторного фильтра Баттерворта

2-element Vector{Any}:
 [0.8167515449979387 -4.907234464272076 … -4.907234464272074 0.8167515449979387]
 [1.0 -5.805671517442012 … -4.14311738707587 0.6670830862565154]

Для анализа свойств синтезированного фильтра выполняется расчет и построение его амплитудно-частотной характеристики. С использованием функции freqz определяется комплексный коэффициент передачи фильтра, после чего выполняется перевод угловой частоты в герцы и расчет АЧХ в логарифмическом масштабе.

# Расчет амплитудно-частотной характеристики (АЧХ) фильтра
h, w = freqz(b, a, 2048)
f_resp = w .* fs ./ (2π)    # перевод угловой частоты в Гц

# Расчет АЧХ в дБ
amp_resp = 20 .* log10.(abs.(h) .+ 1e-12)

plot(f_resp, amp_resp,
    xlabel="Частота, Гц",
    ylabel="АЧХ, дБ",
    legend=false,
    xlim=(0, fs/2)
)

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

С использованием ранее синтезированного режекторного фильтра выполняется фильтрация ЭКГ сигнала. В процессе обработки ослабляется спектральная составляющая в области 60 Гц, соответствующая сетевой помехе. Это позволяет уменьшить влияние внешних электрических наводок и повысить качество ЭКГ сигнала перед дальнейшей детекцией R-пиков и расчетом RR-интервалов.

ecg_filt = filter(b, a, ecg)        # фильтрация ЭКГ сигнала

108000-element Vector{Float64}:
 -0.2613604943993404
 -0.21250240167607903
 -0.26573066807757795
 -0.36614738177362005
 -0.4161172656214138
 -0.39324744502909514
 -0.32258636917030514
 -0.2986311375931025
 -0.3313726245576049
 -0.3962893283973373
 -0.40059139122729925
 -0.3890850840151229
 -0.36649244512364276
  ⋮
 -0.34206831520401354
 -0.3327934850479681
 -0.3324339179822744
 -0.347472287278352
 -0.3435299337450428
 -0.3383511370436336
 -0.33223499162723336
 -0.33419189908414015
 -0.31892828917862315
 -0.3344302495989344
 -0.320462905333482
 -0.3240091843406425

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

freq1, S1_amp, _ = calc_fft(ecg_filt, fs)

plot!(p2, freq1, S1_amp, label="После фильтрации")

plot!(p1, t, ecg_filt, label="После фильтрации")

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

Детекция R-пиков

Для повышения выраженности QRS-комплексов и R-пиков выполняется расчет производной ЭКГ сигнала. Производная позволяет выделить резкие изменения сигнала, характерные для области комплекса QRS. После применения функции diff количество отсчетов уменьшается на единицу, поэтому с использованием функции vcat в начало массива добавляется нулевое значение, что позволяет сохранить исходную длину сигнала.

d_ecg = vcat(0.0, diff(ecg_filt))   # расчет производной ЭКГ сигнала

plot(t, d_ecg,
    xlabel="Время, с",
    ylabel="Амплитуда производной, отн. ед.",
    title="Производная сигнала",
    legend=false
)

После расчета производной выполняется возведение сигнала в квадрат. Такая операция позволяет усилить области с большими амплитудами, соответствующие QRS-комплексам и R-пикам, а также дополнительно подавить низкоамплитудные составляющие сигнала. Кроме того, после возведения в квадрат все значения сигнала становятся положительными, что упрощает дальнейшую обработку и детекцию R-пиков.

sq_ecg = d_ecg .^ 2     # возведение производной сигнала в квадрат

plot(t, sq_ecg,
    xlabel="Время, с",
    ylabel="Амплитуда, отн. ед.",
    title="Квадрат производной",
    legend=false
)

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

# Сглаживание сигнала методом скользящего среднего
win = max(3, Int(round(0.010 * fs)))
kernel = ones(win) ./ win

det_sig = conv(sq_ecg, kernel)  # сглаживание сигнала с использованием свертки
det_sig = det_sig[1:N]  # приведение длины сигнала к исходному размеру

plot(t, det_sig,
    xlabel="Время, с",
    ylabel="Амплитуда, отн. ед.",
    title="Сигнал для детекции R-пиков",
    legend=false
)

Для детекции R-пиков используется подготовленный сигнал det_sig, в котором области QRS-комплексов становятся более выраженными. Порог детекции рассчитывается адаптивно как сумма среднего значения сигнала и 1,5 стандартных отклонений. Такой способ лучше фиксированного порога от максимума, так как меньше зависит от одиночных выбросов и учитывает общий уровень сигнала.

Параметр MinPeakDistance задает минимальное расстояние между соседними найденными пиками. Это необходимо для того, чтобы внутри одного QRS-комплекса не было обнаружено несколько максимумов. Для анализа нормального ритма, брадикардии и тахикардии используется расстояние около 0,30 с, так как значение 0,50 с ограничивает обнаружение ритмов с частотой выше 120 уд/мин и может приводить к пропуску пиков при тахикардии.

Функция findpeaks из библиотеки EngeeDSP выполняет поиск локальных максимумов, удовлетворяющих заданным условиям по высоте и расстоянию между пиками. В результате формируются массив амплитуд найденных пиков pks и массив их положений locs , которые далее используются для расчета RR-интервалов.

thr = mean(det_sig) + 1.5 * std(det_sig)    # адаптивный порог детекции
min_dist = Int(round(0.5 * fs))   #  минимальное расстояние между R-пикам

# Поиск R-пиков
pks, locs = findpeaks(
    det_sig, out=:data,
    MinPeakDistance = min_dist,
    MinPeakHeight  = thr
)

(Ypk = [0.07705275825787196, 0.0916773825278909, 0.10834663376541775, 0.09481633031944532, 0.09963496404605718, 0.10460719300050016, 0.08650497487378789, 0.08319442366179106, 0.09247478104048162, 0.11295427417467355  …  0.1004942866529982, 0.08518680051583521, 0.08025327798569036, 0.08981339449670496, 0.0967951803374813, 0.08915943080450855, 0.08684800494412531, 0.08676418960626395, 0.09614281275354113, 0.08352697256693856], Xpk = [52, 349, 650, 932, 1206, 1493, 1780, 2083, 2381, 2680  …  105334, 105614, 105901, 106188, 106472, 106744, 107012, 107288, 107570, 107857], Wpk = [4.3850150608365865, 4.442932153811, 4.291772746874244, 4.367091061929614, 4.328827355369185, 4.38801284823262, 4.733936150335694, 4.452001752022625, 4.403056166576334, 4.284659123849451  …  4.387227236977196, 4.42038689427136, 4.44097986187262, 4.375775133899879, 4.294981876926613, 4.350528891707654, 4.336650176890544, 4.633688054411323, 4.2678951177804265, 4.363150334407692], Ppk = [0.07704879058713884, 0.09167628535186181, 0.10834571777219215, 0.09481181465597753, 0.09963224995512014, 0.10460611620966644, 0.08650225879249397, 0.0831901731489352, 0.09247192121249741, 0.11295348704301397  …  0.10049394311552719, 0.0851814030130436, 0.08024837323306237, 0.08981076469939231, 0.09679261413663849, 0.08915510722090965, 0.08684034685887856, 0.08675754683512554, 0.09613799667239774, 0.08351632852143341])

Отображаются результаты первичной детекции R-пиков

plot(t, det_sig,
    xlabel="Время, с",
    ylabel="Амплитуда, отн. ед.",
    title="Найденные R-пики",
    legend=false
)
r_times = (locs .- 1) ./ fs   # перевод индексов найденных пиков в секунды

scatter!(r_times, pks, markersize=3)  # отображение найденных пиков на сигнале детекции

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

Так как первичный поиск выполняется не по самому ЭКГ сигналу, а по сигналу det_sig, положение найденного максимума может быть немного смещено относительно реального R-пика. Поэтому далее выполняется уточнение координат: около каждого найденного положения анализируется небольшой участок отфильтрованного ЭКГ сигнала, и за R-пик принимается его локальный максимум.

# Уточнение положения R-пиков 
halfwin = Int(round(0.05 * fs))   # окно поиска
r_locs = Int[]

for idx in locs
    # Границы участка вокруг найденного пика
    l = max(1, idx - halfwin)
    r = min(N, idx + halfwin)

     # Поиск максимума внутри выбранного участка
    seg = ecg_filt[l:r]
    _, imax = findmax(seg)

    push!(r_locs, l + imax - 1)
end

r_locs = unique(sort(r_locs))   # удаление повторяющихся индексов

r_times = (r_locs .- 1) ./ fs   # время найденных R-пиков
r_vals = ecg_filt[r_locs]       # амплитуды найденных R-пиков

plot(t, ecg_filt,
    xlabel="Время, с",
    ylabel="Амплитуда, мВ",
    title="Найденные R-пики",
    legend=false
)

scatter!(r_times, r_vals, markersize=3)

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

Расчет RR-интервалов

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

rr = diff(r_times)  # разность между соседними R-пиками, с

rr_ms = rr.* 1e3  # перевод RR-интервалов в мс

rr_time = r_times[2:end]    # временная ось для RR-интервалов

heart_rate = 60.0 ./ rr # Расчет частоты сердечных сокращений, уд/мин

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

На графике отображается изменение RR-интервалов во времени. По нему можно оценить равномерность сердечного ритма и увидеть, насколько длительность сердечных циклов изменяется на выбранном фрагменте ЭКГ сигнала.

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

plot(rr_time, heart_rate,
    xlabel="Время, с",
    ylabel="ЧСС, уд/мин",
    title="Частота сердечных сокращений",
    legend=false
)

Полученный график позволяет оценить динамику сердечного ритма на выбранном фрагменте ЭКГ сигнала. При относительно стабильном ритме значения ЧСС изменяются плавно и находятся в близком диапазоне. Резкие скачки на графике могут быть связаны как с реальными изменениями сердечного ритма, так и с ошибками детекции R-пиков или артефактами в сигнале.

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

В работе использовались реальные экспериментальные записи электрокардиограммы, полученные из открытого ресурса PhysioNet. В качестве исходных данных применялась база MIT-BIH Arrhythmia Database. В рассматриваемом примере используется запись 100, представленная файлами 100.hea и 100.dat.

Заключение

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

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

Загрузка данных. Исходные данные ЭКГ были получены из открытой базы данных PhysioNet MIT-BIH Arrhythmia Database. Для анализа использовался фрагмент записи 100, содержащий реальный экспериментальный ЭКГ сигнал.
Спектральный анализ. Для оценки частотного состава исходного сигнала был построен односторонний амплитудный спектр. По спектру была выявлена выраженная составляющая в области 60 Гц, соответствующая сетевой помехе.
Фильтрация сигнала. Для подавления сетевой помехи был синтезирован режекторный фильтр Баттерворта. С использованием функций butter, freqz, filter из библиотек EngeeDSP.
Подготовка сигнала к детекции. Для выделения областей QRS-комплексов была рассчитана производная отфильтрованного ЭКГ сигнала, выполнено возведение производной в квадрат и сглаживание методом скользящего среднего. Сглаживание реализовано с помощью функции conv из библиотеки EngeeDSP.
Детекция R-пиков. Для поиска локальных максимумов использовалась функция findpeaks из библиотеки EngeeDSP. Детекция выполнялась с учетом адаптивного порога и минимального расстояния между соседними пиками, что позволило избежать повторных срабатываний внутри одного сердечного цикла.
Уточнение положения R-пиков. После первичной детекции положения R-пиков были уточнены на отфильтрованном ЭКГ сигнале. Это позволило связать результаты поиска по подготовленному сигналу с реальной формой ЭКГ сигнала.
Расчет RR-интервалов и ЧСС. На основе временных координат найденных R-пиков были рассчитаны RR-интервалы и мгновенная частота сердечных сокращений. Полученные графики позволяют оценить изменение длительности сердечных циклов и динамику сердечного ритма на выбранном фрагменте записи.

В результате выполнения примера была реализована последовательность обработки ЭКГ сигнала от загрузки исходных данных до детекции R-пиков, расчета RR-интервалов и частоты сердечных сокращений.