心电信号建模与心率测定

该示例考虑了形成心电图信号的三个标准导联的块模型，测量路径的干扰建模，数字滤波和r峰的流式检测。基础上发现的，计划的时间和中心。在最短的时间内，可以改变整个结构。

显示设置

功能 engee.clear() 清理工作区:

engee.clear()

两种图形渲染模式可用于可视化结果。:

-对于静态快速渲染，它被使用 gr();
-对于具有缩放和查看值的能力的交互式可视化，它被使用 plotlyjs().

通过注释掉第二种模式来选择其中一种模式。此示例使用交互模式。

# gr()
plotlyjs()

Plots.PlotlyJSBackend()

安装必要的库

示例使用库 Statistics 计划的时间。在执行统计分析之前，会检查库在工作环境中的可用性。

let
    installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
    list_packages = ["Statistics"]
    for pack in list_packages
        pack in installed_packages || Pkg.add(pack)
    end
end

using Statistics

加载和启动模型

所提出的块模型实现了心电信号形成和处理的完整路径。在一个，形成三个标志线，II和III。接下来，将来自测量路径的模拟干扰添加到信号中，之后执行数字滤波。对于先导物II，在预处理后，检测R-的，计算RR间隔，并确定心率。

模型结构包括以下主要子系统:

Формирование ЭКС -形成一个模型心电图和三个标准引线;
Моделирование помех -添加网络组件和宽带噪声;
Режекторный фильтр -管理50hz;
ФНЧ -信号的高频分量的限制;
Определение ЧСС -预测-1，计划和时间。

模型可以直接从文件运行。 .engee 使用模拟环境，或从交互式脚本运行它。在第二种情况下，模型以编程方式加载，之后函数 engee.run 开始计算并返回记录的信号。

modelName = "ECG_HR_Model"
modelPath = joinpath(@__DIR__, modelName * ".engee")

model = engee.load(modelPath; force=true);

心电信号的形成

子系统 Формирование ЭКС 生成一个模型心动周期，作为单个分量p、Q、R1、R2、S、ST和T的总和。这种方法是基于一组参数化的非对称高斯函数表示心电信号的形态，并符合AnyMath社区示例的一般思想"建模электрокардиосигнала".

表达式用于一个组件

哪里确定分量的幅度, -其极值在心动周期内的位置，以及定义宽度。

为了获得不对称的形状，宽度分别设置在极值的左侧和右侧。:

一个心动周期的总形状具有以下形式

该模型使用固定的时间位置和主要组件的宽度。在改变心率主要改变相邻心动周期之间间隔的持续时间，而不是拉伸整个P–QRS–T复合物。

心动周期的周期是通过心率确定的:

哪里它设置在节拍/分钟，和它是在几秒钟内获得的。

在手术过程中,手术过程中,手术过程中. 标准引线计算为

伊恩托文？

子系统掩码设置平均心率。此外，还提供了可变节奏模式。当它打开时，为每个新的心动周期生成一个值。

然后计算相应的间隔。 . 参数它以节拍/分钟为单位设置。默认值为3节拍/分钟。

第一次跑步使用每分钟75次的恒定心率。

运行具有恒定心率的模型

在第一阶段，我们将考虑块模型在每分钟75次心跳的恒定心率下的操作。可变性模式被禁用，这使得可以单独分析心电信号的形成，干扰的影响，预处理和R峰检测算法的操作。模拟的持续时间为10秒。

results = engee.run(model; verbose=true);

Building...
Progress 0%
Progress 5%
Progress 10%
Progress 15%
Progress 22%
Progress 27%
Progress 33%
Progress 38%
Progress 43%
Progress 49%
Progress 55%
Progress 63%
Progress 68%
Progress 73%
Progress 78%
Progress 84%
Progress 89%
Progress 95%
Progress 100%

形成的标准引线

我们八仙在换热的上面显示。

src = results["ЭКС"]

t = Float64.(src.time)
lead_I   = Float64.(getindex.(src.value, 1))
lead_II  = Float64.(getindex.(src.value, 2))
lead_III = Float64.(getindex.(src.value, 3))

idx = (t .>= 2.0) .& (t .<= 4.0)

p = plot(t[idx], lead_I[idx],
    xlabel="时间，从",
    ylabel="振幅，mV",
    label="I",
    title="标准引线",
    legend_position=:outertopright)

plot!(p, t[idx], lead_II[idx], label="II")
plot!(p, t[idx], lead_III[idx], label="III")

但是,由于没有任何问题,没有任何问题. 在第二线路中到达最快的R，然后使用确定的R。

干扰模拟

子系统 Моделирование помех 生成心电图的测量路径的失真特征的模型。

网络干扰与从电源到测量设备的电极、连接导线和输入级的电磁干扰相关联。达到50hz的主观要求是100和150hz。:

模型中设置了以下参数

此外，还使用带限白噪声来模拟测量路径的累积宽带分量。对于块 Band-Limited White Noise 参数已设置 Cov = 0.2·10⁻⁶.

子系统的总输出信号可以表示为

哪里 -初始心电信号，以及 -噪声分量。

干扰仿真结果

为了评估干扰的影响，我们比较了初始和失真的导联II。

noisy = results["ЭКС с помехами"]

t_noisy = Float64.(noisy.time)
lead_II_noisy = Float64.(getindex.(noisy.value, 2))

idx_n = (t_noisy .>= 2.0) .& (t_noisy .<= 4.0)

p = plot(t[idx], lead_II[idx],
    xlabel="时间，从",
    ylabel="振幅，mV",
    label="初始ii铅",
    title="模拟干扰的影响",
    legend_position=:outertopright)

plot!(p, t_noisy[idx_n], lead_II_noisy[idx_n],label="II干扰引线")

失真的信号保留了心动周期的基本形态，但出现高频波动。这种表示允许我们进一步评估数字滤波器对先前已知模型信号的操作。

心电信号的预处理

预处理的目的是衰减网络谐波和宽带噪声，同时保留后续检测R峰所必需的QRS复合物的形状。

该模型使用两个连续级：陷波滤波器和低通滤波器。

陷波滤波器

为50hz的网络流量，使用了iiiiirr。在当前模型中，得到的滤波器阶数为10，边界频率为46和54Hz，带宽中允许的纹波为0.5分贝。模块的频率为1000hz。

要知道这种方法的数量，可以使用这些方法。 Проектирование цифровых фильтров. 该应用程序选择幅度-频率响应的类型，合成方法，阶数和频率参数。在[Anymaths]中进行了描述。https://engee.com/helpcenter/stable/ru-en/tutorial/filter-designer.html )。

以下参数用于有问题的陷波滤波器:

-采样率:1000赫兹;
-幅度-频率响应类型：陷波滤波器;
-延年益寿：波黑;
不合时宜的方法。;
-结果顺序:10;
-边界频率：46和54赫兹;
-分贝：0.5分贝。

合成后，计算的系数用于在子系统中自动生成滤波器。 Режекторный фильтр，它允许您将接收到的参数直接传输到块模型。

陷波滤波的结果

notch = results["Режекторный фильтр.Y"]

t_notch = Float64.(notch.time)
lead_II_notch = Float64.(getindex.(notch.value, 2))

idx_f = (t_notch .>= 2.0) .& (t_notch .<= 4.0)

p = plot(t_noisy[idx_n], lead_II_noisy[idx_n],
    xlabel="时间，从",
    ylabel="振幅，mV",
    label="过滤器之前",
    title="陷波滤波",
    legend_position=:outertopright)

plot!(p, t_notch[idx_f], lead_II_notch[idx_f],label="陷波滤波器后")

2000年9月20日。与此相同，100和150hz的频率。

低通滤波器

新的文艺类,每个都有. 为10hz，为38.13hz，为1000hz。

选择最简单的方法. 它抑制了100和150hz的频率，并没有达到正常的频率。截止频率选择在模型心电图信号的主要信息性低频区域上方，并允许您保持明显的QRS复合体。

预处理的结果

filt = results["ФНЧ.Y"]

t_filt = Float64.(filt.time)
lead_II_filt = Float64.(getindex.(filt.value, 2))

idx_lpf = (t_filt .>= 2.0) .& (t_filt .<= 4.0)

p = plot(t_noisy[idx_n], lead_II_noisy[idx_n],
    xlabel="时间，从",
    ylabel="振幅，mV",
    label="与干扰",
    title="预处理的结果",
    legend_position=:outertopright)

plot!(p, t_filt[idx_lpf], lead_II_filt[idx_lpf],label="过滤后")

在意识中，提倡使用...... 输入信号.

对......

为获得成功。该算法以流式方式工作，因此在仿真期间对每个新样本按顺序执行处理。

首先，计算相邻样本的差异:

差操作突出显示信号变化率高的区域。与P和T齿相比，QRS复合体具有更陡峭的锋面，因此其在此操作后的贡献更为明显。

然后对接收到的信号进行平方:

平方消除了差异的标志，并增加了具有高变化率的计数的贡献。

对于平滑，计算10个计数的移动平均值。在1000hz的频率范围内，有10ms的频率范围内的频率范围内的频率范围内的频率范围内的频率范围内的频率范围内的频率范围内的频率范围内的频率范围:

接收到的信号它可用于LED显示屏。

检测阈值自适应地形成。在对应1s的1000个内部的客户端上，计算的平台。:

还有正方形的平均值:

标准偏差是根据当前统计估计值计算的:

后，确定自适应阈值。:

东风......:

在选择QRS区域后，通过过滤后的第二引线澄清R-峰的位置。为此，确定局部最大值。对于连续三次计数，将检查条件:

仅当同时满足两个条件时才接受R-是：检测信号超过自适应阈值，并且在滤波后在导联II中检测到局部最大值。

灭活0.30s的时间。在此间隔期间，不接受新的r候选，这阻止了单个QRS复合体的重新注册。

在模拟开始时，自适应阈值的统计估计尚未完全形成。因此，在0.5秒的初始间隔之后允许检测，这降低了启动模型时发生误报的可能性。

对于哪个预测到的R-HEYRR时间。:

之后，确定心率。:

直到第二个正确检测到的R峰值出现，RR间隔和计算的心率保持等于零。

整体处理序列对应于AnyMath社区["检测的R-峰上心电信号"]的示例中考虑的方法（https://engee.com/community/ru/catalogs/projects/detektsiia-r-pikov-na-ekg-signale 中而在所考虑的模型中，使用建模环境的块以流式方式执行检测。

[美][美][美]

r_sig = results["R-пики"]

t_r = Float64.(r_sig.time)
r_event = Float64.(r_sig.value)

r_idx = findall(r_event .> 0.5)

idx_r = (t_filt .>= 2.0) .& (t_filt .<= 6.0)
r_idx_plot = [i for i in r_idx if 2.0 <= t_r[i] <= 6.0]

p = plot(t_filt[idx_r], lead_II_filt[idx_r],
    xlabel="时间，从",
    ylabel="振幅，mV",
    label="II铅",
    title="R峰的检测",
    legend_position=:outertopright)

scatter!(p,t_r[r_idx_plot],lead_II_filt[r_idx_plot],label="R峰",markersize=4)

小数位数。每个心动周期形成一个检测到的R-的，难治性间隔排除了一个QRS复合物的重复配准。

管理时间

hr_model_sig = results["ЧСС модельная"]
hr_measured_sig = results["ЧСС измеренная"]

t_hr_model = Float64.(hr_model_sig.time)
hr_model = Float64.(hr_model_sig.value)

t_hr_measured = Float64.(hr_measured_sig.time)
hr_measured = Float64.(hr_measured_sig.value)

p = plot(t_hr_model, hr_model,
    xlabel="时间，从",
    ylabel="ЧСС, уд/мин",
    label="预设心率",
    title="预设和测量的心率",
    legend_position=:outertopright)

plot!(p, t_hr_measured, hr_measured,label="测量心率")

在模拟的初始阶段，测量的心率为零。第一个检测到的R-峰值设置第一个RR间隔的开始，并且其持续时间可以仅在第二个R-峰值出现之后计算。因此，每个人都有自己的需求。

额外的时间偏移是由于串行数字滤波和检测操作引入的延迟。在完成初始过渡部分之后，将测量值设置在75节拍/分钟的预设水平附近。

心率变异性建模

考虑一个制度，其中连续心动周期的持续时间随时间而变化。让我们将平均心率保持在每分钟75次，打开可变性，并将标准偏差设置为每分钟3次。

在每个心动周期开始时形成一次心率的新模型值并存储直到其完成。对于统计分析，我们将模拟持续时间增加到180秒，以获得足够数量的RR间隔。

engee.set_param!(modelName * "/Формирование ЭКС",
    "HR" => 75,
    "HR_VAR" => true,
    "HR_SIGMA" => 3.0)

engee.set_param!(modelName, "StopTime" => 180.0)

results_var = engee.run(model; verbose=true);

Building...
Progress 0%
Progress 5%
Progress 10%
Progress 15%
Progress 15%
Progress 20%
Progress 25%
Progress 30%
Progress 32%
Progress 37%
Progress 42%
Progress 47%
Progress 52%
Progress 57%
Progress 62%
Progress 67%
Progress 72%
Progress 77%
Progress 82%
Progress 87%
Progress 92%
Progress 97%
Progress 100%
Progress 100%

模型和测量的心率

hr_model_var_sig = results_var["ЧСС модельная"]
hr_measured_var_sig = results_var["ЧСС измеренная"]

t_hm = Float64.(hr_model_var_sig.time)
hr_m = Float64.(hr_model_var_sig.value)

t_hd = Float64.(hr_measured_var_sig.time)
hr_d = Float64.(hr_measured_var_sig.value)

p = plot(t_hm, hr_m,
    xlabel="时间，从",
    ylabel="ЧСС, уд/мин",
    label="模型心率",
    title="模型和测量的心率",
    legend_position=:outertopright)

plot!(p, t_hd, hr_d,label="测量心率")

测量的心率跟随模型值变化，但两个依赖关系的更新时间在时间上不重合。模型值在心动周期开始时设置，而测量的心率仅在检测下一个R峰值并计算完成的RR间隔之后才能确定。

通过数字预处理路径和R峰值检测操作的序列引入额外的时间偏移。在观察到的位移是测量算法的结构的结果，并且不指示计算RR间隔的误差。

[医]时间延迟

我们将完成所有的任务。我们使用三种表示:

-Rr-转速表-用于解码的时间;
-Rr的目标--第一时间;
-散射图-用于分析相邻区间的关系和 .

在AnyMath社区的示例中使用了类似的一组表示"ECG心率变异性分析сигналу".

r_var_sig = results_var["R-пики"]
rr_var_sig = results_var["RR-интервал"]

r_event_var = Float64.(r_var_sig.value)
t_r_var = Float64.(r_var_sig.time)
rr_hold = Float64.(rr_var_sig.value)

r_idx_var = findall(r_event_var .> 0.5)

rr_beats = rr_hold[r_idx_var]
rr_time = t_r_var[r_idx_var]

valid_rr = (rr_beats .> 0.3) .& (rr_beats .< 2.0)

rr_beats = rr_beats[valid_rr]
rr_time = rr_time[valid_rr]
rr_ms = 1000 .* rr_beats
heart_rate_rr = 60.0 ./ rr_beats

mean_rr = Statistics.mean(rr_ms)
std_rr = Statistics.std(rr_ms)
mean_hr = Statistics.mean(heart_rate_rr)

println("平均RR间隔: ", round(mean_rr, digits=2), " 女士")
println("RR间隔的标准偏差: ", round(std_rr, digits=2), " 女士")
println("平均心率: ", round(mean_hr, digits=2), " уд/мин")

Средний RR-интервал: 802.98 мс
Стандартное отклонение RR-интервалов: 32.65 мс
Средняя ЧСС: 74.84 уд/мин

Rr-秒速赛车

Rr-表示射手的速度。基本上可以是这样的。

plot(rr_time, rr_ms,
    xlabel="时间，从",
    ylabel="RR-间隔，ms",
    title="RR-转速表",
    legend=false)

RR时对应确定的平台性。值的扩展由心电信号生成单元的掩模中设置的可变性参数确定。

RR影视版

直方图显示心动周期持续时间的分布。以毫秒为单位的RR间隔绘制在X轴上，相应范围内的间隔数绘制在Y轴上。

histogram(rr_ms,
    bins=30,
    xlabel="RR-间隔，ms",
    ylabel="间隔数",
    title="RR间隔直方图",
    legend=false)

Rr区的主体部分在平台上。模型心率的非零标准偏差导致形成连续心动周期持续时间的分布。

RR影视版

散射图显示相邻心动周期之间的关系。当前时间制版。，而且是下载时间。 . 对应关系。

rr_n = rr_ms[1:end-1]
rr_next = rr_ms[2:end]

p = scatter(rr_n, rr_next,
    xlabel="RR(n),ms",
    ylabel="RR（n+1），ms",
    title="RR间隔的散射图",
    legend=false,
    markersize=3,
    xlim=(0, 1200),
    ylim=(0, 1200),
    aspect_ratio=:equal)

plot!(p, [0, 1200], [0, 1200],linestyle=:dash,linewidth=2)

通过较小的模拟可变性，点紧凑地位于对角线附近。沿两个轴传播的增加对应于相邻RR间隔之间更明显的差异。

结论

在示例中，考虑了形成和处理模型心电信号的单个块链。八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八、八陷波滤波器和低通滤波器在检测R峰之前对信号进行预处理。

R峰的检测基于差分处理的平方，滑动平均，自适应阈值，局部最大细化和耐火保护。基础上输出的R-1，计划RR时间和中心。

心电图形成单元的面罩允许您使用恒定或变化的心率。对于具有给定可变性的模式，计算的RR间隔表示为相邻间隔的RR转速表、直方图和散射图。