Engee 文档
Notebook

心电图信号的可视化和分析

导言

本逐步练习介绍如何读取、显示和分析心电图信号。心电图又称心电图,是心脏收缩时产生的电信号。心电图能快速检测心脏健康状况以及各种异常情况,如心律失常、心脏病发作或心动过速,因此被广泛使用。

下图是典型的心电图。

ecg_001.jpg

运行此演示需要连接互联网

心电图分析

从计算机上的文件加载和显示心电图

第一步是从心电图.dat文件导入数据。

然后将 ecgdata.dat 文件中的数据以矩阵(1x50000)形式导入 ecg 变量。

In [ ]:
Pkg.add(["CSV", "Findpeaks"])
In [ ]:
using CSV, DataFrames
ecg = Matrix( CSV.read("$(@__DIR__)/electrocardiogram.dat", DataFrame, header = 0, delim=';') );

通过 typeof() 命令,我们可以获得变量的特征。

In [ ]:
typeof(ecg)
Out[0]:
Matrix{Float64} (alias for Array{Float64, 2})

该文件包含 250 秒的心电图记录,由于采样率为 200 Hz,因此文件包含 50,000 个样本。现在我们可以使用 plot 命令绘制心电图:

In [ ]:
using Plots
gr(size=(1700, 600), legend=false)
plot(ecg,label=false)
title!("Электрокардиограмма")
Out[0]:

该信号从 http://people.ucalgary.ca/~ranga/ 读取。在讲义、实验练习和信号数据文件下,您可以找到大量来自各种数据采集方法(心电图、脑电图、EDS 等)的信号。您可以将这些文件从浏览器保存到硬盘,然后使用上述下载命令打开它们。

独立工作

尝试阅读和显示本网站或包含本实验材料的文件夹中的其他文件。

近似

从图中可以看出,由于样本数量较多,很难观察到每个心动周期的细节。我们可以通过在 plot 命令中添加 ylimitsxlimits 来设置图表近似参数。

In [ ]:
plot(ecg, label=false, ylimits=(1500,2800), xlimits=(3000, 3800))
Out[0]:

独立工作

检查心电图。您可以更改 ylimitsxlimits 来检查信号的不同区域。各处的心电图是否相同?心电图看起来不同的原因是什么?

使用阈值检测心电图峰值

使用该信号的一个有趣而简单的方法是检测每个周期的峰值,即信号中超过某一特定值的区域。例如,我们可以任意将阈值设定为 2400 信号电平,并在此电平处生成一个新信号,将心电图与阈值进行比较:

In [ ]:
peaks_1 = ecg .> 2400;

同时,我们还可以生成一个时间匹配信号,这在以后会很有用:

In [ ]:
time_1 = (1:length(ecg))/200;

现在,我们可以在时域中绘制原始信号:

In [ ]:
plot(time_1, ecg)
title!("Электрокардиограмма")
Out[0]:

我们可以将信号源和阈值这两个信号绘制在一起进行比较。我们可以使用命令后的感叹号符号将一个图形叠加到另一个图形上。

由于阈值信号的值介于 0 和 1 之间,我们需要对其进行缩放,使两个信号处于相同的范围内。我们可以在 catter 命令中将逻辑矢量乘以相关峰值来实现这一点。

In [ ]:
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)
Out[0]:

独立工作

将阈值 (2400) 改为其他值。会发生什么?检测到的峰值是多了还是少了?

峰值大小

请注意,由于我们现在使用的是时间,因此坐标轴指令发生了变化。如果我们仔细观察,就能发现问题所在:

In [ ]:
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)
Out[0]:

使用 Findpeaks 软件包进行峰值检测

阈值为信号超过给定值的每一个地方提供一个值,这与峰值检测不太一样。幸运的是,Julia有一个软件包(请务必检查Julia是否有满足您需要的软件包。你很有可能找到一个能让你的工作更轻松的软件包),用于峰值检测,名为Findpeaks

要开始使用,您需要导入Findpeaks软件包。

In [ ]:
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变量类型更改为抽象数组,然后就可以按如下方式使用Findpeaks软件包了:

In [ ]:
using Findpeaks
ecg = vec(ecg::AbstractArray)
peaks_3 = findpeaks(ecg)
plot(time_1, ecg, labels=false)
scatter!(time_1[peaks_3], ecg[peaks_3])
Out[0]:

在不应该有 "峰值 "的地方似乎出现了很多 "峰值"。让我们放大看看:

In [ ]:
plot(time_1, ecg, label=false, ylimits=(1500,2800), xlimits=(15, 17))
scatter!(time_1[peaks_3], ecg[peaks_3])
Out[0]:

这是一个很好的例子,说明了某些处理方法的特殊性。峰值检测是检测任何一个高于其相邻点的点,这与我们想要检测的内容不同。这就引出了一个问题:我们想要什么样的 "峰值"?

澄清峰值检测

我们需要相对较大的峰值(而不是可能与噪音有关的峰值),我们需要与心率密切相关的尖而高的峰值。这可以概括为以下规则:峰值必须是:

(a) 高于某一水平(同样,2400 可能是一个很好的阈值)、

b) 相互之间不要太接近。

第二条规则可能不是很明显,所以让我们一条一条地看。首先,我们舍弃所有低于阈值的峰值并显示结果:

In [ ]:
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])
Out[0]:

现在我们可以看到第二条规则对结果的影响,因为在一个 "峰 "中发现了两个位置。让我们重复一遍,只找出相距至少 10 个位置的峰值:

In [ ]:
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])
Out[0]:

最后一步,我们可以显示心电图,并在其下方显示峰值之间的时间差,也就是说,我们可以看到从这些峰值推断出的心率。我们可以使用 subplot 命令在一张图中显示两个信号,如下所示:

In [ ]:
peaks_5 = sort(peaks_5);
s = 60*diff(peaks_5/200);
a = peaks_5[2:end];
In [ ]:
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 bpm 附近出现峰值。
  2. 尝试改进代码以避免出现此类问题。
  3. 除了上述峰值外,您可能还会注意到心率的明显波动。您认为这是患者心律不齐造成的吗?还是信号处理问题?还是信号采集阶段出现了问题?在接近终点时,有一个区域的数值降到了 20 以下,这里发生了什么?
  4. 尝试读取其他信号,并应用这个简单的步骤进行分析。
  5. 你还能想到其他可视化心电图和其他生物医学信号的方法吗?