心电图信号可视化与分析
导言
这个循序渐进的练习介绍了心电图信号的读取、显示和分析过程。 心电图，也称为心电图，是心脏收缩时产生的电信号。 心电图被广泛使用，因为它可以快速识别心脏健康状态，以及各种异常，如心律失常，心脏病发作或心动过速。
下图显示了典型的心电图。
运行此演示需要internet连接
心电图分析
从计算机上的文件下载并显示心电图
第一步是从心电图导入数据。dat文件。
然后我们从ecgdata导入数据。dat文件到ecg变量作为矩阵（1x50000）。
Pkg.add(["CSV", "Findpeaks"])
using CSV, DataFrames
ecg = Matrix( CSV.read("$(@__DIR__)/electrocardiogram.dat", DataFrame, header = 0, delim=';') );
Typeof()命令允许您获取变量的特征。
typeof(ecg)
该文件包含250秒的心电图记录，由于采样率为200Hz，因此该文件包含50,000个样本。 现在我们可以使用plot命令显示心电图。:
using Plots
gr(size=(1700, 600), legend=false)
plot(ecg,label=false)
title!("Электрокардиограмма")
这个信号是从网站上读到的http://people.ucalgary.ca /~ranga/。 在部分讲义，实验室练习和信号数据文件您将发现来自各种数据收集方法（ECG，EEG，EMF等）的大量信号。). 您可以将文件从浏览器保存到硬盘驱动器，然后使用如上所述的下载命令打开它们。
独立工作
尝试从网站或包含此实验室工作材料的文件夹中读取和显示其他文件。
方法
正如您从图表中看到的那样，样本数量会使观察每个心动周期的细节变得困难。 我们可以通过在plot命令中添加ylimits和xlimits来设置图形近似参数。
plot(ecg, label=false, ylimits=(1500,2800), xlimits=(3000, 3800))
独立工作
检查心电图。 您可以更改ylimits和xlimits的值以考虑信号的不同区域。 心电图到处都是一样的吗? 心电图有什么不同的原因？
使用阈值处理检测ECG峰值
这个信号的一个有趣和简单的练习是检测每个周期的峰值-信号超过一定值的区域。 例如，我们可以在2400处任意设置阈值并生成新信号以将ECG与阈值进行比较。:
peaks_1 = ecg .> 2400;
与此同时，我们可以创建一个适合时间的信号，这将在未来有用。:
time_1 = (1:length(ecg))/200;
现在我们可以在时域中绘制原始信号。:
plot(time_1, ecg)
title!("Электрокардиограмма")
我们可以通过将两个信号（原始信号和阈值）绘制在一起来比较它们。 我们可以在命令后使用感叹号将一个图形复盖在另一个图形上。
由于阈值信号的值从0到1，我们需要对其进行缩放，使两个信号处于相同的范围内。 我们可以在cutter命令中执行此操作，方法是将逻辑向量乘以我们感兴趣的峰值。
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变量类型更改为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])
这是一些处理方法的具体细节的一个很好的例子。 峰值检测是对任何高于其邻居的点的检测，这与我们试图检测的不同。 这就引出了一个问题，我们需要什么样的"高峰"？
峰值检测的细化
我们需要相对较大的峰值（不是最有可能与噪声相关的峰值），以及应该与心率密切相关的尖锐和高峰。 这可以导致以下规则：:
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附近有一个突然的峰值，这很可能是由于处理错误。
独立工作:
1）尝试分析代码，以了解为什么会有大约100次/分钟的峰值。
2）尝试改进代码以避免此类问题。
3)除了上述峰值之外，您可能会注意到心率的显着波动。 你认为这是由于病人的心跳不规则吗？ 还是信号处理？ 或者在接收信号的阶段有问题？ 在接近尾声时，有一个值低于20的区域，那里会发生什么？
4）尝试读取其他信号并应用此简单的分步分析。
5）你能想到其他方式来可视化ECG和其他生物医学信号吗？