Визуализация и анализ сигнала электрокардиограммы
Вступление
Данное пошаговое упражнение является введением в процессы считывания, отображения и анализа сигнала электрокардиограммы. Электрокардиограмма, также известная как ЭКГ, представляет собой электрический сигнал, который вырабатывается при сокращении сердца. ЭКГ широко используется, поскольку она может быстро выявить состояние здоровья сердца, а также различные аномалии, такие как аритмии, инфаркт или тахикардия.
На изображении ниже показана типичная ЭКГ.
Для запуска этой демонстрации требуется подключения к интернет
Анализ ЭКГ
Загрузка и отображение ЭКГ из файла на вашем компьютере
Первым шагом является импорт данных из файла electrocardiogram.dat.
Затем мы импортируем данные из файла ecgdata.dat в переменную ecg в виде матрицы (1x50000).
using CSV, DataFrames
ecg = Matrix( CSV.read("$(@__DIR__)/electrocardiogram.dat", DataFrame, header = 0, delim=';') );Команда typeof() позволяет получить характеристики переменной.
typeof(ecg)Matrix{Float64} (alias for Array{Float64, 2})Этот файл содержит 250 секунд записи ЭКГ, и поскольку частота дискретизации составляет 200 Гц, файл содержит 50 000 образцов. Теперь мы можем отобразить ЭКГ с помощью команды plot:
using Plots
gr(size=(1700, 600), legend=false)
plot(ecg,label=false)
title!("Электрокардиограмма")
Этот сигнал был прочитан с веб-сайта http://people.ucalgary.ca/~ranga/. В разделе Lecture notes, lab exercises, and signal data files вы найдете большое количество сигналов от различных способов сбора данных (ЭКГ, ЭЭГ, ЭДС и т. д.). Вы можете сохранить файлы из браузера на жесткий диск, а затем открыть с помощью команды загрузки, как указано выше.
Самостоятельная работа
Попробуйте прочитать и отобразить другие файлы с веб-сайта или из папки, содержащей материалы данной лабораторной работы.
Приближение
Как видно из графика, количество образцов может затруднить наблюдение деталей каждого сердечного цикла. Мы можем задать параметры приближения графика, добавив ylimits и xlimits к команде plot.
plot(ecg, label=false, ylimits=(1500,2800), xlimits=(3000, 3800))
Самостоятельная работа
Исследуйте ЭКГ. Вы можете изменить значения ylimits и xlimits, чтобы рассмотреть различные регионы сигнала. ЭКГ везде одинаковая? По каким причинам ЭКГ может выглядеть по-разному?
Обнаружение пиков ЭКГ с помощью пороговой обработки
Одним интересным и простым упражнением с этим сигналом является обнаружение пиков каждого цикла – областей сигнала, которые превышают определенное значение. Например, мы можем произвольно установить порог на уровне сигнала 2400 и сгенерировать новый сигнал на для сравнения ЭКГ с порогом:
peaks_1 = ecg .> 2400;В то же время мы можем создать сигнал, соответствующий времени, что будет полезно в дальнейшем:
time_1 = (1:length(ecg))/200;Теперь мы можем построить график исходного сигнала уже во временном домене:
plot(time_1, ecg)
title!("Электрокардиограмма")
Можем сравнить два сигнала, исходный и пороговый, построив их вместе. Мы можем наложить один график на другой, используя символ восклицательного знака после команды.
Поскольку пороговый сигнал будет иметь значения от 0 до 1, нам нужно масштабировать его, чтобы два сигнала находились в одних и тех же диапазонах. Мы можем сделать это внутри команды catter, умножив логический вектор на интересующее нас значение пиков.
plot(time_1, ecg, label=false, ylimits=(1500,2800), xlimits=(15, 20))
scatter!(time_1, 2600*(peaks_1), mc=:red, ms=2, ma=0.5)
Самостоятельная работа
Измените значение порога (2400) на другие значения. Что произойдет? Вы обнаруживаете больше или меньше пиков?
Размер пика
Обратите внимание, что, поскольку теперь мы используем время, инструкция оси изменилась. Если мы посмотрим внимательно, мы можем обнаружить проблему:
plot(time_1, ecg, label=false, ylimits= (2000, 2800), xlimits = (16.4, 17.2))
scatter!(time_1, 2600*(peaks_1), mc=:red, ms=2, ma=0.5)
Обнаружение пиков с помощью пакета «Findpeaks»
Порог предоставляет значение для каждого отдельного места, где сигнал превышает заданное значение, что не совсем то же самое, что обнаружение пика. К счастью, в Julia есть пакет (всегда проверяйте, есть ли в Julia пакеты для того, что вам нужно. Весьма вероятно, что вы сможете найти пакет, который облегчит вашу работу) для обнаружения пиков, называемый Findpeaks.
Для начала необходимо импортировать пакет Findpeaks.
using Pkg
Pkg.add(url="https://github.com/tungli/Findpeaks.jl") Cloning git-repo `https://github.com/tungli/Findpeaks.jl`
Updating git-repo `https://github.com/tungli/Findpeaks.jl`
Updating registry at `~/.julia/registries/EngeeJuliaRegistry`
Updating git-repo `https://gitlab.kpm-ritm.ru/engee/backend/kernels/simcg/engeejuliaregistry.jl.git`
Updating registry at `~/.julia/registries/ProgrammaticModeling`
Updating git-repo `https://gitlab.kpm-ritm.ru/engee/backend/kernels/ProgrammaticModelingRegistry.jl.git`
Updating registry at `~/.julia/registries/General.toml`
Resolving package versions...
Installed EpollShim_jll ───────── v0.0.20230411+0
Installed DifferentialEquations ─ v7.6.0
Installed PyCall ──────────────── v1.94.1
Installed StaticArrays ────────── v1.5.12
Installed ModelingToolkit ─────── v8.43.0
Updating `~/.julia/environments/v1.8/Project.toml`
[d9b5fb6a] + Findpeaks v0.1.0 `https://github.com/tungli/Findpeaks.jl#master`
Updating `~/.julia/environments/v1.8/Manifest.toml`
[d9b5fb6a] + Findpeaks v0.1.0 `https://github.com/tungli/Findpeaks.jl#master`
[2702e6a9] + EpollShim_jll v0.0.20230411+0
Building PyCall → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/53b8b07b721b77144a0fbbbc2675222ebf40a02d/build.log`
Precompiling project...
✓ EpollShim_jll
✓ Findpeaks
2 dependencies successfully precompiled in 9 seconds. 211 already precompiled.После завершения импорта пакета, нам необходимо поменять тип переменной ecg на AbstractArray и затем пакет Findpeaks можно использовать следующим образом:
using Findpeaks
ecg = vec(ecg::AbstractArray)
peaks_3 = findpeaks(ecg)
plot(time_1, ecg, labels=false)
scatter!(time_1[peaks_3], ecg[peaks_3])
Кажется, что там очень много «пиков», где их быть не должно. Давайте увеличим масштаб:
plot(time_1, ecg, label=false, ylimits=(1500,2800), xlimits=(15, 17))
scatter!(time_1[peaks_3], ecg[peaks_3])
Это хороший пример специфики некоторых методов обработки. Обнаружение пиков — это обнаружение ЛЮБОЙ точки, которая выше своих соседей, что отличается от того, что мы пытаемся обнаружить. Это приводит к вопросу, какие «пики» нам нужны?
Уточнение обнаружения пиков
Нам нужны относительно большие пики (не те, которые, скорее всего, связаны с шумом), а также резкие и высокие пики, которые должны быть тесно связаны с частотой сердечных сокращений. Это можно привести к таким правилам: пики, которые:
а) находятся выше определенного уровня (опять же, 2400 может быть хорошим порогом),
б) не слишком близки друг к другу.
Второе правило может быть не сразу понятно, поэтому давайте рассмотрим их по очереди. Во-первых, давайте отбросим все пики ниже порога и отобразим результат:
peaks_4 = findpeaks(ecg, min_height = 2400.)
plot(time_1, ecg, label=false, ylimits=(1500,3300), xlimits=(20, 23))
scatter!(time_1[peaks_4], ecg[peaks_4])
Теперь мы можем увидеть, как повлияет на результат второе правило, поскольку в одном «пике» обнаружены две позиции. Давайте повторим и найдем только те вершины, которые удалены друг от друга не менее чем на 10 позиций:
peaks_5 = findpeaks(ecg, min_height = 2400., min_dist = 10)
plot(time_1, ecg, label=false, ylimits=(1500,3300), xlimits=(20, 23))
scatter!(time_1[peaks_5], ecg[peaks_5])
В качестве последнего шага мы можем отобразить ЭКГ, а под ней – временные различия между пиками, то есть мы можем увидеть частоту сердечных сокращений, выведенную из этих пиков. Мы можем отобразить два сигнала на одном рисунке с помощью команды subplot следующим образом:
peaks_5 = sort(peaks_5);
s = 60*diff(peaks_5/200);
a = peaks_5[2:end];p1 = plot(time_1, ecg, label=false, title = "ЭКГ");
s1 = scatter!(time_1[peaks_5], ecg[peaks_5]);
p2 = plot(time_1[a], s, legend=false, title = "Сердечный ритм");
plot(p1, p2, layout=(2,1))Наконец, мы отобразили частоту сердечных сокращений, извлеченную из ЭКГ.
Проверка анализа
Обратите внимание, что частота сердечных сокращений относительно стабильна около 40 ударов в минуту, что ниже обычных 60 ударов в минуту. Однако есть один внезапный пик около 100, это, скорее всего, связано с ошибкой в обработке.
Самостоятельная работа
1) Попробуйте проанализировать код, чтобы понять, почему происходит пик около 100 уд/мин. 2) Попробуйте улучшить код, чтобы избежать подобных проблем. 3) Помимо вышеупомянутого пика, вы можете отметить значительные колебания сердечного ритма. Как вы думаете, это связано с нерегулярным сердцебиение пациента? Или с обработкой сигнала? Или с проблемами на этапе получения сигнала? Ближе к концу есть область, где значение падает ниже 20, что происходит там? 4) Попробуйте прочитать другие сигналы и применить этот простой пошаговый анализ. 5) Можете ли вы придумать другие способы визуализации ЭКГ и других биомедицинских сигналов?