下图是典型的心电图。
运行此演示需要连接互联网。
然后将 ecgdata.dat 文件中的数据以矩阵(1x50000)形式导入 ecg 变量。
Pkg.add(["CSV", "Findpeaks"])
using CSV, DataFrames
ecg = Matrix( CSV.read("$(@__DIR__)/electrocardiogram.dat", DataFrame, header = 0, delim=';') );
通过 typeof() 命令,我们可以获得变量的特征。
typeof(ecg)
该文件包含 250 秒的心电图记录,由于采样率为 200 Hz,因此文件包含 50,000 个样本。现在我们可以使用 plot 命令绘制心电图:
using Plots
gr(size=(1700, 600), legend=false)
plot(ecg,label=false)
title!("Электрокардиограмма")
该信号从 http://people.ucalgary.ca/~ranga/ 读取。在讲义、实验练习和信号数据文件下,您可以找到大量来自各种数据采集方法(心电图、脑电图、EDS 等)的信号。您可以将这些文件从浏览器保存到硬盘,然后使用上述下载命令打开它们。
独立工作¶
尝试阅读和显示本网站或包含本实验材料的文件夹中的其他文件。
近似¶
从图中可以看出,由于样本数量较多,很难观察到每个心动周期的细节。我们可以通过在 plot 命令中添加 ylimits 和 xlimits 来设置图表近似参数。
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")
导入软件包完成后,我们需要将ecg变量类型更改为抽象数组,然后就可以按如下方式使用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])
这是一个很好的例子,说明了某些处理方法的特殊性。峰值检测是检测任何一个高于其相邻点的点,这与我们想要检测的内容不同。这就引出了一个问题:我们想要什么样的 "峰值"?
澄清峰值检测¶
我们需要相对较大的峰值(而不是可能与噪音有关的峰值),我们需要与心率密切相关的尖而高的峰值。这可以概括为以下规则:峰值必须是:
(a) 高于某一水平(同样,2400 可能是一个很好的阈值)、
b) 相互之间不要太接近。
第二条规则可能不是很明显,所以让我们一条一条地看。首先,我们舍弃所有低于阈值的峰值并显示结果:
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 左右有一个突然的峰值,这很可能是处理错误造成的。
独立作业:¶
- 尝试分析代码,了解为何在 100 bpm 附近出现峰值。
- 尝试改进代码以避免出现此类问题。
- 除了上述峰值外,您可能还会注意到心率的明显波动。您认为这是患者心律不齐造成的吗?还是信号处理问题?还是信号采集阶段出现了问题?在接近终点时,有一个区域的数值降到了 20 以下,这里发生了什么?
- 尝试读取其他信号,并应用这个简单的步骤进行分析。
- 你还能想到其他可视化心电图和其他生物医学信号的方法吗?