下图是典型的心电图。
运行此演示需要连接互联网。
然后将 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 以下，这里发生了什么？
- 尝试读取其他信号，并应用这个简单的步骤进行分析。
- 你还能想到其他可视化心电图和其他生物医学信号的方法吗？