心电信号分析心率变异性
该项目检查了基于ECG信号的心率变异性分析。 基于发现的R峰,计算RR间隔和心率,之后使用RR转速表,rr间隔直方图和散射图可视化结果。 此外,进行正常节律,心动过缓和心动过速的比较。
显示设置
功能 engee.clear() 清理工作区:
engee.clear()
两种图形渲染模式可用于可视化结果。:
*对于静态,快速渲染,使用 gr().
*对于交互式,具有缩放和查看值的能力的动态可视化,使用 plotlyjs().
通过注释掉第二个命令来选择其中一个命令。 两个后端是互斥的,最后一个被调用的后端将被激活。
# gr()
plotlyjs()
我们将使用函数进行连接 include 的"读"文件。jl"用于读取文件:
include("$(@__DIR__)/read.jl")
安装必要的库
在这个阶段,必要的Julia库将被检查并自动安装,如果它们在工作环境中不可用。 有问题的示例使用库 Statistics 在心率变异性分析中执行基本统计计算所必需的,包括处理RR间隔、计算平均值和估计数据分布。
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
图书馆正在激活 EngeeDSP 并导入了数字心电信号处理所需的功能,包括滤波、频谱分析、峰值搜索和卷积。
import EngeeDSP.Functions: fft, butter, filter, findpeaks, conv
心电信号处理与分析
本节讨论R峰的检测和ECG信号的RR间隔的分析。 基于发现的R峰,计算RR间隔,之后应用分析心率变异性的基本方法。:
*RR-tachogram-显示心动周期持续时间随时间的变化;
*rr间隔直方图-反映rr间隔在所选记录片段上的分布;
*散射图-允许您评估相邻心动周期之间的关系。
这些方法的联合使用使得可以直观地分析心率的性质并评估其变异性。
心电信号处理功能
让我们想象一下功能 detect_r 专为心电信号处理和R峰值检测而设计。 例"[在ECG信号上检测R峰]中详细描述的主要步骤在函数内部执行(https://engee.com/community/ru/catalogs/projects/detektsiia-r-pikov-na-ekg-signale )":过滤信号,准备用于检测的数据,搜索R峰,澄清它们的位置,以及计算RR间隔和心率。
function detect_r(ecg, fs;
f_center = 60.0,
bw = 6.0,
filter_order = 6,
smooth_window = 0.020,
min_peak_distance = 0.30,
search_window = 0.05
)
N = length(ecg)
# 抑制网络干扰
f1 = f_center - bw / 2
f2 = f_center + bw / 2
wn1 = f1 / (fs / 2)
wn2 = f2 / (fs / 2)
b, a = butter(filter_order, [wn1, wn2], "stop")
ecg_filt = filter(b, a, ecg)
# 用于R峰检测的信号准备
d_ecg = vcat(0.0, diff(ecg_filt))
sq_ecg = d_ecg .^ 2
win = max(3, Int(round(smooth_window * fs)))
kernel = ones(win) ./ win
det_sig = conv(sq_ecg, kernel)
det_sig = det_sig[1:N]
# 搜索R峰
thr = mean(det_sig) + 1.5*std(det_sig)
min_dist = Int(round(min_peak_distance * fs))
_, locs = findpeaks(
det_sig,
out = :data,
MinPeakDistance = min_dist,
MinPeakHeight = thr
)
# 澄清R峰的位置
search_radius = Int(round(search_window * fs))
r_locs = Int[]
for idx in locs
l = max(1, idx - search_radius)
r = min(N, idx + search_radius)
seg = ecg_filt[l:r]
_, imax = findmax(seg)
push!(r_locs, l + imax - 1)
end
r_locs = unique(sort(r_locs))
# R峰的时间和振幅
r_times = (r_locs .- 1) ./ fs
r_vals = ecg_filt[r_locs]
# RR-间隔和心率
rr = diff(r_times)
rr_ms = rr .* 1000
rr_time = r_times[2:end]
heart_rate = 60.0 ./ rr
return ecg_filt, r_times, r_vals, rr_ms, rr_time, heart_rate
end
信号加载和处理
正在从记录中下载ECG片段 122 麻省理工学院-波黑数据库。 第一记录通道被选择为被分析的信号,之后它被传送到功能 detect_r.
data = "$(@__DIR__)/mitdb/122"
hdr = read_header(data)
fs = Float64(hdr.fs)
# 选择用于分析的信号片段
t_start = 10.0 * 60.0 # 片段的开头,从
t_dur = 10.0 * 60.0 # 片段的持续时间,s
# 将时间转换为倒计时数字
sampfrom = Int(round(t_start * fs))
sampto = sampfrom + Int(round(t_dur * fs))
# 读取ECG信号的选定片段
_, raw, phys = read_dat_212(data; sampfrom=sampfrom, sampto=sampto)
# 选择第一心电信号通道
ecg = phys === nothing ? Float64.(raw[:, 1]) : Float64.(phys[:, 1])
# 时间轴的形成
N = length(ecg)
t = collect(0:N-1) ./ fs
# 信号处理和R峰值检测
ecg_filt, r_times, r_vals, rr_ms, rr_time, heart_rate = detect_r(ecg, fs)
让我们用发现的R峰的显示来构建一个滤波的ECG信号。
plot(t, ecg_filt,
xlabel="时间,从",
ylabel="振幅,mV",
label="心电信号"
)
scatter!(r_times, r_vals, markersize=3, label="R峰")
RR区间分析
让我们计算主要的心率指标:平均RR间隔,rr间隔的标准偏差和平均心率。
mean_rr = mean(rr_ms) # 平均RR间隔,ms
std_rr = std(rr_ms) # RR间隔的标准偏差,ms
mean_hr = mean(heart_rate) # 平均心率,节拍/分钟
println("平均RR间隔: ", round(mean_rr, digits=2), " 女士")
println("RR间隔的标准偏差: ", round(std_rr, digits=2), " 女士")
println("平均心率: ", round(mean_hr, digits=2), " bpm")
RR-转速表
RR转速表反映了RR间隔的持续时间随时间的变化,即它显示了连续心动周期的持续时间如何在ECG信号的选定片段上变化。 基于这种依赖性,可以评估心率的稳定性,识别rr间隔中的平滑波动,并且还注意到可能与伪像,检测错误或节奏特征相关联的突然变化。
plot(rr_time, rr_ms,
xlabel="时间,从",
ylabel="RR-间隔,ms",
title="RR-转速表",
legend=false
)
RR间隔直方图
RR间隔的直方图显示心动周期持续时间的分布。 以毫秒为单位的RR间隔的值绘制在X轴上,相应范围内的间隔数绘制在Y轴上。 此可视化允许您确定RR间隔的主要部分集中在哪个值附近,它们的传播有多明显,以及是否存在与主要组明显不同的单个间隔。
histogram(rr_ms,
bins=30,
xlabel="RR-间隔,ms",
ylabel="间隔数",
title="RR间隔直方图",
legend=false
)
RR间隔的散射图
RR间隔的散射图显示了相邻心动周期之间的关系。 当前RR间隔绘制在X轴上。 ,并且在Y轴上是以下RR间隔 . 此图表允许您估计下一个心动周期的持续时间与前一个心动周期的持续时间有多大差异。 如果点位于紧凑和接近对角线 相邻的RR区间具有相似的值,并且节奏相对稳定。 点扩散的增加反映了更明显的心率变异性或不规则收缩的可能存在。
rr_n = rr_ms[1:end-1]
rr_next = rr_ms[2:end]
# 正确显示的轴边界
rr_min = minimum(vcat(rr_n, rr_next))
rr_max = maximum(vcat(rr_n, rr_next))
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)
)
plot!([0, 1200], [0, 1200],
linestyle=:dash,
linewidth=2
)
正常节律,心动过缓和心动过速的分析
本节讨论了rr间隔和心率对正常节律,心动过缓和心动过速的比较。 一般来说,成人的正常静息心率在每分钟约60-100次的范围内。 心动过缓对应于慢节奏,其中心率通常低于60次/分钟,心动过速对应于快速节奏,其中心率超过100次/分钟。
作为示例的一部分,这些类型的节律被视为用于分析ECG信号的示范案例。 对于每个信号,将检测R峰值,计算RR间隔和平均心率,并构建RR转速表,rr间隔直方图和散射图。 这将使我们能够比较心动周期的持续时间和RR间隔的分布如何随着正常节律,心动过缓和心动过速而变化。
加载、处理和显示ECG信号
为了便于操作,我们将创建一个函数 load_mitdb,旨在从MIT-BIH数据库下载ECG记录的选定片段。
function load_mitdb(record_name;
t_start_min = 5.0,
t_dur_min = 10.0,
channel = 1
)
data = "$(@__DIR__)/mitdb/$record_name"
# 阅读有关记录的信息
hdr = read_header(data)
fs = Float64(hdr.fs)
# 将时间转换为倒计时数字
t_start = t_start_min * 60.0
t_dur = t_dur_min * 60.0
sampfrom = Int(round(t_start * fs))
sampto = sampfrom + Int(round(t_dur * fs))
# 读取ECG信号的选定片段
_, raw, phys = read_dat_212(data; sampfrom=sampfrom, sampto=sampto)
# 选择心电信号通道
ecg = phys === nothing ? Float64.(raw[:, channel]) : Float64.(phys[:, channel])
# 时间轴的形成
N = length(ecg)
t = collect(0:N-1) ./ fs
return ecg, fs, t
end
为了演示正常的心率,我们将上传记录数据. 100 我们将显示处理的结果:具有标记R峰的滤波ECG信号。
ecg_norm, fs_norm, t_norm = load_mitdb(
"100";
t_start_min = 5.0,
t_dur_min = 5.0,
channel = 1
)
ecg_filt_norm, r_times_norm, r_vals_norm, rr_ms_norm, rr_time_norm, heart_rate_norm = detect_r(ecg_norm, fs_norm)
p_norm = plot(t_norm, ecg_filt_norm,
xlabel = "时间,从",
ylabel = "振幅,mV",
label = "正常的节奏"
)
scatter!(p_norm, r_times_norm, r_vals_norm,markersize = 2,label = "R峰")
我们将以同样的方式处理记录 123 选择来证明心动过缓。 该图表将显示过滤后的ECG信号和找到的R峰,从中计算RR间隔。
ecg_brady, fs_brady, t_brady = load_mitdb(
"123";
t_start_min = 5.0,
t_dur_min = 5.0,
channel = 1
)
ecg_filt_brady, r_times_brady, r_vals_brady, rr_ms_brady, rr_time_brady, heart_rate_brady = detect_r(ecg_brady, fs_brady)
p_brady = plot(t_brady, ecg_filt_brady,
xlabel = "时间,从",
ylabel = "振幅,mV",
label = "心动过缓"
)
scatter!(p_brady, r_times_brady, r_vals_brady,markersize = 2,label = "R峰")
为了分析快速心率,我们使用录音 209. 在过滤和检测R峰后,我们将显示处理后的ECG信号,这将直观地评估相邻心动周期之间距离的减少。
ecg_tachy, fs_tachy, t_tachy = load_mitdb(
"209";
t_start_min = 6.0,
t_dur_min = 5.0,
channel = 1
)
ecg_filt_tachy, r_times_tachy, r_vals_tachy, rr_ms_tachy, rr_time_tachy, heart_rate_tachy = detect_r(ecg_tachy, fs_tachy)
p_tachy = plot(t_tachy, ecg_filt_tachy,
xlabel = "时间,从",
ylabel = "振幅,mV",
label = "心动过速"
)
scatter!(p_tachy, r_times_tachy, r_vals_tachy, markersize = 2,label = "R峰")
我们将针对不同类型的心率显示处理后的ECG信号。
plot(p_norm, p_brady, p_tachy,
layout = (3, 1),
size = (1000, 850),
margin = 40*Plots.px
)
根据获得的具有不同类型心率的ECG信号的示波图,能够注意到相邻R峰之间的距离的差异。 对于心动过缓,与正常节律相比,心脏收缩之间的间隔增加,这对应于心率的降低。 在心动过速中,相反,R-峰位于彼此更靠近的位置,这表明心率增加和rr间隔持续时间的减少。
RR区间分析
为了方便显示数字指标,我们将形成一个函数 print_hr,其中计算并显示RR间隔和心率的主要特征。
function print_hr(name, rr_ms, heart_rate)
mean_rr = mean(rr_ms) # 平均RR间隔,ms
std_rr = std(rr_ms) # RR间隔的标准偏差,ms
mean_hr = mean(heart_rate) # 平均心率,节拍/分钟
println(name)
println("平均RR间隔: ", round(mean_rr, digits=2), " 女士")
println("RR间隔的标准偏差: ", round(std_rr, digits=2), " 女士")
println("平均心率: ", round(mean_hr, digits=2), " bpm")
println()
end
print_hr("正常的节奏", rr_ms_norm, heart_rate_norm)
print_hr("心动过缓", rr_ms_brady, heart_rate_brady)
print_hr("心动过速", rr_ms_tachy, heart_rate_tachy)
根据统计分析的结果,可以看出所获得的值对应于心率类型之间的预期差异。 对于正常节律,平均RR间隔为771.8ms,平均心率为每分钟77.99次,这对应于正常心率的范围。
随着心动过缓,平均RR间隔增加至1172.63ms,平均心率下降至51.79次/分钟。 这证实了与正常状态相比心率的减慢。 在心动过速中,相反,平均RR间隔减少到549.12ms,平均心率增加到每分钟115.26次,这对应于快速心率。
因此,统计指标确认ECG信号的视觉分析的结果。
RR-转速表
我们将针对所考虑的心律类型构建RR-tachogram:正常节律,心动过缓和心动过速。
p_rr_norm = plot(rr_time_norm, rr_ms_norm,
xlabel = "时间,从",
ylabel = "RR-间隔,ms",
title = "正常的节奏",
legend = false
)
p_rr_brady = plot(rr_time_brady, rr_ms_brady,
xlabel = "时间,从",
ylabel = "RR-间隔,ms",
title = "心动过缓",
legend = false
)
p_rr_tachy = plot(rr_time_tachy, rr_ms_tachy,
xlabel = "时间,从",
ylabel = "RR-间隔,ms",
title = "心动过速",
legend = false
)
plot(p_rr_norm, p_rr_brady, p_rr_tachy,
layout = (3, 1),
size = (1000, 850),
margin = 40*Plots.px
)
RR间隔直方图
让我们考虑分析的心率类型的RR间隔分布。 基于单个图表的直方图允许您比较值的范围,并确定每种情况下心动周期的主要部分集中在哪里。
histogram(rr_ms_norm,
bins = 30,
alpha = 0.5,
xlabel = "RR-间隔,ms",
ylabel = "间隔数",
title = "RR间隔直方图",
label = "正常的节奏"
)
histogram!(rr_ms_brady,
bins = 30,
alpha = 0.5,
label = "心动过缓"
)
histogram!(rr_ms_tachy,
bins = 30,
alpha = 0.5,
label = "心动过速"
)
直方图显示分析的节奏类型的RR间隔分布占据不同的范围。 在正常节律下,值的主要部分集中在约750-800ms的区域,其对应于约78次/分钟的平均心率。
对于心动过缓,分布转移到大RR间隔的区域,大约为1000-1300ms。这证实了心动周期持续时间的增加和心率的降低。 相反,在心动过速中,分布转移到较小RR间隔的区域,这对应于快速心率。
由于分析是在真实的ECG记录上进行的,因此分布可能具有不完美的形状。 特别是,对于心动过速,两个特征极大值是明显的,这表明在所选择的记录片段中存在rr间隔的两个主要范围。 这可能是由于所分析区域内的节律变化或真实临床信号的特征。
因此,直方图使得可以直观地比较心动周期的持续时间,并根据RR间隔的分布显示正常节律,心动过缓和心动过速之间的差异。
RR间隔的散射图
让我们考虑分析的心率类型的RR间隔的散射图。
rr_n_norm = rr_ms_norm[1:end-1]
rr_next_norm = rr_ms_norm[2:end]
rr_n_brady = rr_ms_brady[1:end-1]
rr_next_brady = rr_ms_brady[2:end]
rr_n_tachy = rr_ms_tachy[1:end-1]
rr_next_tachy = rr_ms_tachy[2:end]
rr_min = 0
rr_max = 1800
# 正常的节奏
p_scatter_norm = scatter(rr_n_norm, rr_next_norm,
xlabel = "RR(n),ms",
ylabel = "RR(n+1),ms",
title = "正常的节奏",
legend = false,
markersize = 3,
xlim = (rr_min, rr_max),
ylim = (rr_min, rr_max)
)
plot!(p_scatter_norm, [rr_min, rr_max], [rr_min, rr_max],
linestyle = :dash,
linewidth = 2
)
# 心动过缓
p_scatter_brady = scatter(rr_n_brady, rr_next_brady,
xlabel = "RR(n),ms",
ylabel = "RR(n+1),ms",
title = "心动过缓",
legend = false,
markersize = 3,
xlim = (rr_min, rr_max),
ylim = (rr_min, rr_max)
)
plot!(p_scatter_brady, [rr_min, rr_max], [rr_min, rr_max],
linestyle = :dash,
linewidth = 2
)
# 心动过速
p_scatter_tachy = scatter(rr_n_tachy, rr_next_tachy,
xlabel = "RR(n),ms",
ylabel = "RR(n+1),ms",
title = "心动过速",
legend = false,
markersize = 3,
xlim = (rr_min, rr_max),
ylim = (rr_min, rr_max)
)
plot!(p_scatter_tachy, [rr_min, rr_max], [rr_min, rr_max],
linestyle = :dash,
linewidth = 2
)
# 显示散射图
plot(p_scatter_norm, p_scatter_brady, p_scatter_tachy,
layout = (3, 1),
size = (800, 1400),
margin = 40*Plots.px
)
对于一个正常的节奏,点位于足够紧凑和主要集中在对角线附近。 这表明邻近的RR间隔具有相似的值,并且心动周期持续时间的变化是适度的。
对于心动过缓,点转移到大RR间隔的区域并且分布更广泛。 没有明确定义的紧凑区域,这表明心动周期持续时间的大可变性。 这种模式也与RR间隔的直方图一致,其中分布具有更宽的值范围。
对于心动过速,可以看到两个不同的点区域:一个在短RR间隔区域,另一个在较长间隔区域。 这表明在所选片段中存在心动周期持续时间的两个特征范围。 该结果与直方图一致,该直方图还显示出两个明显的分布极大值。
为了进行更直观的比较,所有三个表格图都另外显示在同一个图形上。
p_scatter_all = scatter(rr_n_norm, rr_next_norm,
xlabel = "RR(n),ms",
ylabel = "RR(n+1),ms",
label = "正常的节奏",
markersize = 3,
xlim = (rr_min, rr_max),
ylim = (rr_min, rr_max),
legend_position = :outertopright
)
scatter!(p_scatter_all, rr_n_brady, rr_next_brady,
label = "心动过缓",
markersize = 3
)
scatter!(p_scatter_all, rr_n_tachy, rr_next_tachy,
label = "心动过速",
markersize = 3
)
plot!(p_scatter_all, [rr_min, rr_max], [rr_min, rr_max],
linestyle = :dash,
linewidth = 2,
label = false
)
RR间隔的一般表格图显示三种类型的节奏占据不同的值范围。
使用的材料
这项工作使用了从开源PhysioNet获得的真实实验心电图记录。 将[MIT-BIH心律失常数据库]作为源数据(https://physionet.org/content/mitdb/1.0.0 /)。 记录用于心率变异性的初始分析 122. 另外的记录被用来比较不同类型的心率. 100, 123 和 209 对应于正常节律、心动过缓和心动过速的示范实例。
结论
在这个例子中,考虑了基于ECG信号的心率变异性的分析。 基于发现的R峰,计算RR间隔和心率,之后构建主要依赖关系以评估心率的性质。
示例的工作包括以下主要步骤:
**数据下载。**从开放数据库PhysioNet MIT-BIH心律失常数据库获得初始ECG数据。 采用真实的实验心电图记录进行分析.
R峰的检测。该功能用于处理心电信号 detect_r,其中执行信号滤波、准备检测、搜索R峰和细化它们的位置。 此函数使用EngeeDSP库的功能,包括 butter, filter, conv 和 findpeaks.
**计算RR间隔和心率。**基于R峰值的时间坐标,计算RR间隔和心率。 此外,还确定了平均RR间隔,RR间隔的标准偏差和平均心率。
**心率变异性分析。**为了评估心率变化的性质,构建了RR转速表,rr间隔直方图和散射图。 这些依赖性使得能够估计心动周期的持续时间、RR间隔的分布以及相邻间隔之间的关系。
**心率类型的比较。**正常节律,心动过缓和心动过速被认为证明了差异。 根据分析结果可以看出,随着心动过缓,RR间隔增加,平均心率降低。 随着心动过速,相反,RR间隔变短,平均心率增加。
作为示例的结果,示出了基于ECG信号的心率变异性的分析序列:从r峰值的检测到rr间隔的测速图、直方图和散射图的构建。