AnyMath 文档
Notebook

使用两个移动平均线和动态阈值检测收缩期FPG峰值

演示示例讨论了由M.Elgendi和合作者提出的检测光体积描记信号（PPG）收缩峰的算法。 该方法基于初步的带通滤波，非线性信号变换，计算两个不同窗口长度的移动平均线并形成动态阈值。

与直接搜索局部极大值不同，该算法首先确定与潜在脉搏波对应的时间区域，然后在每个接受区域内找到收缩压最大值。

来自开放集MIMIC Performance AF数据集的记录片段用于演示。 正在考虑从档案中记录没有心房颤动的患者 mimic_perform_non_af_data.mat.

显示设置

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

In [ ]: 
engee.clear()

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

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

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

In [ ]: 
# gr()
plotlyjs()
Out[0]:
Plots.PlotlyJSBackend()

安装必要的库

库用于下载MAT格式的源数据。 MAT. 图书馆 Statistics 它用于形成用于检测ECG的R峰的辅助算法的自适应阈值。

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

using MAT, Statistics

加载源数据

档案 mimic_perform_non_af_data.mat 它必须与交互式脚本位于同一文件夹中。 下载后，将显示MAT文件的结构。 阵列 data 它包含单独的记录，每个记录都包括PHG，ECG和其他生理信号。

In [ ]: 
data_path = "$(@__DIR__)/mimic_perform_non_af_data.mat"
mat_data = matread(data_path)
Out[0]:
Dict{String, Any} with 3 entries:
  "source"  => Dict{String, Any}("matlab_conversion_script"=>"/Users/petercharl…
  "data"    => MatlabStructArray{2}(["ppg", "ekg", "imp", "abp", "fix"], Matrix…
  "license" => Dict{String, Any}("details"=>"This dataset is licensed under the…

来自非AF组的第7个条目用于进一步处理。 从中提取FPG和ECG通道。 采样率为125Hz。

正在考虑信号的15秒片段。 该持续时间允许您同时观察几个连续的脉搏波并跟踪其幅度随时间的变化。

In [ ]: 
data = mat_data["data"]

subject_no = 7
subject = data[subject_no]

ppg_all = vec(Float64.(subject["ppg"]["v"]))
ecg_all = vec(Float64.(subject["ekg"]["v"]))

fs = 125.0

start_sample = 60000
duration = 15.0
N = Int(fs * duration)
end_sample = start_sample + N - 1

ppg = ppg_all[start_sample:end_sample]
ecg = ecg_all[start_sample:end_sample]
t = (0:N-1) ./ fs

println("采样率: ", fs, " 赫兹")
println("分析片段的持续时间: ", duration, " 与")

Частота дискретизации: 125.0 Гц
Длительность анализируемого фрагмента: 15.0 с

初始光体积描记信号

首先，让我们构建fpg的选定片段的临时表示。

In [ ]: 
plot(t, ppg,
    xlabel="时间，从",
    ylabel="的振幅",
    label="石油气",
    title="初始光体积描记信号",
    legend_position=:outertopright)
Out[0]:

脉波在所选区域中保持可区分，但其振幅和形状发生变化。 因此，对于自动检测，使用产生局部变化阈值且不依赖于单个固定幅度值的算法是方便的。

FPG预处理

该算法的第一阶段是初始FPG的带通滤波。 Elgendi使用带宽为0.5至8Hz的二阶巴特沃斯滤波器。 下限减小了缓慢变化的分量的影响，而上限限制了高频波动。

过滤由函数执行 filtfilt() 在正向和反向方向上。 这种治疗方法避免了相移，这在随后确定收缩极大值的时间位置时尤其重要。

In [ ]: 
f_low = 0.5
f_high = 8.0
order = 2

Wn = [f_low, f_high] ./ (fs / 2)

b, a = EngeeDSP.Functions.butter(order, Wn, "bandpass")
ppg_filt = EngeeDSP.Functions.filtfilt(b, a, ppg)

p = plot(t, ppg,
    xlabel="时间，从",
    ylabel="的振幅",
    label="原始的FPG",
    title="FPG带通滤波",
    legend_position=:outertopright)

plot!(p, t, ppg_filt, label="过滤后")
Out[0]:

在带通滤波之后，基线的缓慢变化和高频波动减小。 与此同时，主脉冲波持续存在并在时间上保持局部化。

滤波后，信号的负值被复位为零，正部分被平方。:

将负部分归零会留下滤波后的FPG的正区域，平方会增加明显的脉搏波和微小波动之间的差异。

In [ ]: 
ppg_pos = max.(ppg_filt, 0.0)
z = ppg_pos .^ 2

plot(t, z,
    xlabel="时间，从",
    ylabel="的振幅",
    label="转换后的信号",
    title="非线性FPG变换",
    legend_position=:outertopright)
Out[0]:

非线性变换后，FPG的主要正部分成为明显的脉冲区域。 此外，它们的时间结构使用两个移动平均线进行分析。

两个移动平均线

该算法的关键阶段是计算具有不同窗口长度的两个移动平均线。

第一个移动平均线是在一个短窗口中计算的:

短窗口突出显示对应于脉搏波的收缩部分的局部区域。

第二移动平均线是在较长的窗口上计算的。:

长窗口变化更慢，并表征心动周期尺度上变换信号的局部水平。 它被用作形成动态阈值的基础。

在初始操作中，窗口持续时间约为0.11s和0.67s–精确的优化值为111和667ms。在125hz的采样频率下，这些间隔不能用整数个样本精确表示，因此使用最接近的奇数窗口大小。

In [ ]: 
T_peak = 0.111
T_beat = 0.667

w1 = 2 * floor(Int, T_peak * fs / 2) + 1
w2 = 2 * floor(Int, T_beat * fs / 2) + 1

ma_peak = EngeeDSP.Functions.movmean(z, w1)
ma_beat = EngeeDSP.Functions.movmean(z, w2)

println("MA_peak窗口: ", w1, " 计数, ", round(1000*w1/fs, digits=0), " 女士")
println("MA_beat窗口: ", w2, " 计数, ", round(1000*w2/fs, digits=0), " 女士")

Окно MA_peak: 13 отсчётов, 104.0 мс
Окно MA_beat: 83 отсчётов, 664.0 мс
In [ ]: 
p = plot(t, z,
    xlabel="时间，从",
    ylabel="的振幅",
    label="转换后的信号",
    title="变换后的FPG信号的移动平均线",
    legend_position=:outertopright)

plot!(p, t, ma_peak, label="MA_peak")
plot!(p, t, ma_beat, label="MA_beat")
Out[0]:

在图表上 MA_peak 在个别脉波附近形成局部上升。 MA_beat 它的变化要慢得多，并设置一个相对于收缩区域的存在进一步评估的局部水平。

因此，该算法使用两个时间尺度：一个短的用于识别可能的脉搏波，一个长的用于使阈值适应当前信号电平。

动态阈值和感兴趣的领域

动态阈值是基于长移动平均线形成的:

哪里

补充资料 设置相对于 MA_beat，并且阈值本身随着长移动平均线而随时间变化。

满足条件时形成感兴趣区域

该阶段的处理顺序如下:

  1. 比较 MA_peak 和动态阈值;
  2. 兴趣区域的开始和结束由交叉点确定。;
  3. 太短的区域被丢弃。;
  4. 在每个接受区域内，确定滤波后的FPG的最大模值。;
  5. 找到的最大值的位置作为收缩峰的位置。

因此，不针对整个记录搜索最大模值，而仅在预分配的时域内搜索。

In [ ]: 
beta = 0.02
alpha = beta * mean(z)

threshold = ma_beat .+ alpha
boi = ma_peak .> threshold

p = plot(t, ma_peak,
    xlabel="时间，从",
    ylabel="的振幅",
    label="MA_peak",
    title="动态阈值",
    legend_position=:outertopright)

plot!(p, t, threshold, label="动态阈值")
Out[0]:

在哪里 MA_peak 如果超过动态阈值，则形成脉波候选。 因为阈值跟随缓慢的变化 MA_beat，检测判据适应本地信号电平的变化。

In [ ]: 
edges = diff(Int.(vcat(false, boi, false))) 	# 确定进入和退出感兴趣区域的时刻

block_start = findall(edges .== 1)  # 感兴趣领域开始的索引
block_end = findall(edges .== -1) .- 1  	# 感兴趣的领域结束的索引

peak_idx = Int[]	# 发现的收缩峰指数

for (i_start, i_end) in zip(block_start, block_end)
    if i_end - i_start + 1 >= w1	# 丢弃太短的区域
        local_idx = argmax(abs.(ppg_filt[i_start:i_end]))   # 搜索区域内的最大模
        push!(peak_idx, i_start + local_idx - 1)	# 将本地索引转换为整个片段的索引
    end
end

检测结果

让我们在初始PHG上显示发现的收缩压峰值。

In [ ]: 
p = plot(t, ppg,
    xlabel="时间，从",
    ylabel="的振幅",
    label="石油气",
    title="收缩峰的检测",
    legend_position=:outertopright)

scatter!(p, t[peak_idx], ppg[peak_idx], label="收缩压峰值", markersize=4)
Out[0]:

标记位于所选脉搏波的最大值处。 在这种情况下，每个收缩峰的位置仅在使用两个移动平均线和动态阈值形成感兴趣区域之后确定。

心电图与AFG的比较

对于心脏事件的视觉比较，请考虑来自同一时间片段的心电图。 ECG的R峰由辅助算法确定，并且仅用作视觉地标。

首先，ECG在5-20Hz频段进行滤波。 接下来，计算相邻样本的差异，结果被平方并通过移动平均值平滑。 该信号突出显示了QRS复合体的区域。 通过自适应阈值确定候选者，之后通过在每个候选者的小邻域中搜索最明显的ECG极值来细化R峰的位置。

In [ ]: 
Wn_ecg = [5.0, 20.0] ./ (fs / 2)

b_ecg, a_ecg = EngeeDSP.Functions.butter(2, Wn_ecg, "bandpass")
ecg_filt = EngeeDSP.Functions.filtfilt(b_ecg, a_ecg, ecg)

d_ecg = vcat(0.0, diff(ecg_filt))
sq_ecg = d_ecg .^ 2

w_ecg = 2 * floor(Int, 0.080 * fs / 2) + 1
ecg_det = EngeeDSP.Functions.movmean(sq_ecg, w_ecg)

thr_ecg = Statistics.mean(ecg_det) + 1.5 * Statistics.std(ecg_det)
min_r_distance = Int(round(0.30 * fs))

_, qrs_idx = EngeeDSP.Functions.findpeaks(ecg_det, out=:data, MinPeakHeight=thr_ecg, MinPeakDistance=min_r_distance)

search_radius = Int(round(0.080 * fs))
r_idx = Int[]

for idx in qrs_idx
    i1 = max(1, idx - search_radius)
    i2 = min(N, idx + search_radius)
    local_idx = argmax(abs.(ecg_filt[i1:i2]))
    push!(r_idx, i1 + local_idx - 1)
end

r_idx = unique(sort(r_idx));

ECG和检测到的R峰

首先，我们将在ECG上显示辅助检测R峰的结果。

In [ ]: 
p = plot(t, ecg,
    xlabel="时间，从",
    ylabel="的振幅",
    label="心电图",
    title="心电图R峰的检测",
    legend_position=:outertopright)

scatter!(p, t[r_idx], ecg[r_idx], label="R峰", markersize=4)
Out[0]:

发现的标记位于明显的QRS复合物的区域，然后仅用于与FPG的脉波的视觉比较。

ECG和FPG的联合显示

为了进行比较，我们将在公共时间轴上绘制ECG和FPG。 我们注意到上图上的R峰，下图上的PHG收缩峰。

In [ ]: 
p_ecg = plot(t, ecg, ylabel="心电图振幅", label="心电图", title="心电图", legend_position=:outertopright)
scatter!(p_ecg, t[r_idx], ecg[r_idx], label="R峰", markersize=4)

p_ppg = plot(t, ppg, xlabel="时间，从", ylabel="FPG振幅", label="石油气", title="石油气", legend_position=:outertopright)
scatter!(p_ppg, t[peak_idx], ppg[peak_idx], label="收缩压峰值", markersize=4)

plot(p_ecg, p_ppg, layout=(2, 1), link=:x, size=(1000, 650))
Out[0]:

组合表示允许可视地比较ECG上的电心脏事件的序列和PHG的检测到的脉搏波。 每个通道都使用自己的算法来识别特征点，因此该图是专门为事件序列的定性比较而设计的。

ECG的R峰值与PHG的收缩压峰值之间的绝对时间距离在本例中不用于计算脉搏波的传播时间。

结论

演示示例实现了基于两个移动平均线和动态阈值检测FPG收缩压峰值的算法。

在预处理阶段，进行带通零相位滤波和非线性信号转换。 接下来，功能 movmean() 计算两个移动平均线：短的用于突出脉搏波的局部区域，长的用于形成时变阈值。 通过交叉点 MA_peak 确定感兴趣区域，其中基于滤波后的FPG的最大模值确定收缩峰的位置。

此外，对同一时间片段进行ECG R峰的辅助检测。 ECG和PHG的组合显示使您可以直观地比较心脏事件和检测到的脉搏波的序列，而无需计算心率和其他派生指标。

使用的材料

  1. Elgendi M.，Norton I.，Brearley M.，Abbott D.，Schuurmans D.在热带条件下从紧急响应者测量的加速度光电容积图中的收缩峰检测//PLOS ONE。 – 2013. -卷。 8，10号。 [DOI：10.1371/journal.pone.0076585]（https://doi.org/10.1371/journal.pone.0076585）

  2. Charlton P.H.MIMIC执行数据集//Zenodo。 [DOI：10.5281/zenodo.6807403]（https://doi.org/10.5281/zenodo.6807403）