数字光体积描记图处理与脉搏波形态分析
该演示演示了光电体积描记图的数字处理过程。 使用实验数据的例子，一致地考虑了脉冲信号分析的主要阶段。 作为示例的一部分，执行信号的预处理和滤波、时间和频谱分析、脉搏波的特征点的检测以及脉搏波的心率和形态参数的计算。
安装必要的库
该演示使用开放的Julia语言库，用于加载、结构化和实验数据的主要处理。
-
图书馆
CSV-用于处理CSV文件的库。 它用于从包含光体积描记录的CSV文件中读取和下载实验数据。 -
图书馆
DataFrames-用于处理表格数据结构的库。 -
图书馆
Statistics-提供基本统计功能的库。 它用于计算数据的统计特性，例如，信号参数分析中使用的平均值。
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["CSV","DataFrames", "Statistics"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
using CSV, DataFrames, Statistics
光体积描记图处理
功能 engee.clear() 清理工作区:
engee.clear()
两种图形渲染模式可用于可视化结果。:
*对于静态，快速渲染，使用 gr().
*对于交互式，具有缩放和查看值的能力的动态可视化，使用 plotlyjs().
通过注释掉第二个命令来选择其中一个命令。 两个后端是互斥的，最后一个被调用的后端将被激活。
# gr()
plotlyjs()
此示例使用CSV文件，其中包含坐姿获得的光电容积图数据和相关信号。 原始记录取自开放数据库[PhysioNet:Pulse Transit Time(PPG)Database](https://physionet.org/content/pulse-transit-time-ppg/1.1.0 /)。
我们用实验光体积描记数据设置文件的路径。 使用建筑 @__DIR__ 允许您创建链接到当前脚本目录的相对文件路径。
data_path = "$(@__DIR__)/s1_sit.csv"
此步骤从CSV文件加载数据。 功能 CSV.read 使用指定的路径读取文件的内容并将其转换为结构 DataFrame.
data = CSV.read(data_path, DataFrame);
在变量中 fs 根据所用数据库的描述，信号的采样频率设置为500赫兹。 接下来，确定光体积描记信号样本的数量并形成时间矢量。
fs = 500.0 # 消化不良的频率，Hz
N = length(data.pleth_1) # 计数数目
t = (0:N-1) ./ fs # 时间向量，用
让我们在时域中绘制光体积描记图信号。
plot(t, data.pleth_1,
xlabel = "时间，从",
ylabel = "的振幅",
label = "pleth_1",
grid = true
)
为了便于分析，我们将在10到70秒的范围内选择记录的时间片段。 根据设定的时间值，计算出相应的采样指标，之后从初始光体积描记信号中选择选定的区域。
t_start = 10.0 # 片段的开头，从
t_end = 70.0 # 片段的结尾，用
n_start = t_start *fs # 片段开头的索引
n_end = t_end * fs # 片段结尾的索引
# 选择FPG信号的片段
x_raw = Float64.(-data.pleth_1[Int64(n_start):Int64(n_end)])
nsampl = length(x_raw) # 片段的计数数
t = range(t_start, step=1/fs, length=nsampl) # 交换向量
从先前获得的分析信号的时间表示，能够注意到信号电平的缓慢变化的存在。 这种行为对于光体积描记录来说是典型的，并且可能与配准条件有关。
因此，使用EngeeDSP库和函数 detrend() 信号中去除所述线性趋势。 接下来，在去除趋势后绘制信号的图形，这简化了对脉搏波形状的后续分析。
x = EngeeDSP.Functions.detrend(x_raw, 1)
plot(t, x,
label="FPG pleth_1 ",
grid=true,
xlabel="时间，从",
ylabel="的振幅",
xlims = (t_start,t_end)
)
EngeeDSP库用于分析信号的频率组成。 使用函数 fft() 计算分析信号的快速傅立叶变换，然后使用函数将结果相对于零频率居中 fftshift(). 结果，形成信号的双向幅度谱。
基于采样频率和信号长度形成频率矢量。 接下来，计算振幅谱，并在有限的频率范围内绘制其曲线图。
fft_1 = EngeeDSP.Functions.fft(x)
fftsh_1 = EngeeDSP.Functions.fftshift(fft_1)
df = fs / nsampl
freq_vec = -fs/2 : df : fs/2 - df
Amp1 = abs.(fftsh_1)./ nsampl
plot(freq_vec, Amp1,
xlabel="频率，赫兹",
ylabel="的振幅",
title="振幅谱",
grid=true,
legend=false,
xlims = (-100, 100)
)
从趋势去除后的信号的曲线图和从先前构造的双向幅度谱中可以看出，记录中包含与脉搏波的形状不相关的高频波动。 对此，我们应用低通滤波器。
为此，在EngeeDSP库中，使用函数 fir1() 计算截止频率fc=30Hz的FIR滤波器的系数。 接下来，由函数执行信号滤波 filter(). 图表显示了趋势去除后的信号与滤波后的信号的比较。
fc = 30.0 # 截止频率，Hz
Wn = fc / (fs/2) # 归一化截止频率
order = 20 # 滤波器FIR顺序
b = EngeeDSP.Functions.fir1(order, Wn)
a = [1.0]
y = EngeeDSP.Functions.filter(b, a, x)
p = plot(t, x,
label="石油气",
xlabel="时间，从",
ylabel="的振幅",
grid=true)
plot!(p, t, y, label="LPF后的FPG fc=△fc Hz")
xlims!(p, 10, 20)
display(p)
在应用低通滤波器后，我们注意到光体积描记信号的高频分量被显着抑制。 结果，降低了原始信号中存在的噪声和微小波动的影响。 滤波后的信号变得更平滑，而脉搏波的基本形状保持不变。 接下来，我们将分析滤波后信号的频率成分如何变化。 为此，使用EngeeDSP库，我们计算滤波信号的频谱，并随后定心。 然后我们将计算双向幅度谱，并将原始信号和滤波后的信号的谱图绘制在同一张图上。
fft_2 = EngeeDSP.Functions.fft(y)
fftsh_2 = EngeeDSP.Functions.fftshift(fft_2)
Amp2 = abs.(fftsh_2)./ nsampl
plot(freq_vec, Amp1,
xlabel="频率，赫兹",
ylabel="的振幅",
title="振幅谱",
label = "原始的FPG",
grid=true,
legend=true,
xlims = (-100, 100)
)
plot!(freq_vec, Amp2,label = "液化石油气后的液化石油气")
EngeeDSP库和函数用于搜索特征脉搏波点。 findpeaks(). 参数 minPeakDistance 设置样本中相邻峰之间的最小允许距离，并防止在同一心动周期内重复检测峰。
首先，对滤波后的光体积描记信号的局部极大值进行搜索。 然后由反转信号确定局部极小值。
minPeakDistance = 300 # 相邻峰之间的最小距离，计数
# 搜索局部脉波最大值
pks_max, locs_max = EngeeDSP.Functions.findpeaks(
y, out=:data,
MinPeakDistance = minPeakDistance
)
# 搜索局部最小值（通过反转信号）
pks_min, locs_min = EngeeDSP.Functions.findpeaks(
-y, out=:data,
MinPeakDistance = minPeakDistance
)
让我们在时域中绘制光电容积图，并标记在其上找到的局部极大值和极小值。
p = plot(t, y,
label="FPG pleth_1 ",
grid=true,
xlabel="时间，从",
ylabel="的振幅"
)
scatter!(p, t[locs_max], y[locs_max], label="高点")
scatter!(p, t[locs_min], y[locs_min], label="最低限额")
xlims!(p, 10, 20)
要计算心率（HR），有必要确定连续脉搏峰值之间的时间间隔。 为此，对找到的最大值的指数进行预排序，之后计算相邻峰值之间的间隔，同时考虑到信号的采样频率。 基于所获得的间隔的平均值，计算平均心率。
# 最大索引数组
locs_max = sort(Int.(locs_max))
# 计算相邻极大值之间的时间间隔，
rr_sec = diff(locs_max) ./ fs
# 平均心率的计算
hr_bpm = 60.0 / mean(rr_sec)
println("平均心率: ", round(hr_bpm, digits=1), " 节拍/分钟")
在演示的第一部分中，考虑了光电体积描记图的数字处理过程，包括信号预处理，时域和频域分析以及心率估计。
接下来，我们将考虑将光电容积图分割为单个心动周期并计算脉搏波的形态参数的方法，例如收缩波的高度，anacrota和catacrota的持续时间以及脉搏周期的总持续时间。
脉搏波形态分析
光体积描记图（PHG）反映了通过光学方法记录的组织血液填充的脉冲变化。 在时域中，AFG是周期性重复波，其形状由心动周期的相位和血管床的状态决定。
一个脉冲周期通常被认为是两个连续最小值之间的信号的一部分。 周期内，分配主收缩压波和下降段，对应于脉压下降的相位。
脉搏波的初始上升段称为anacrota并且对应于活动动脉血流的相位。 Anacrota的特点是信号幅度迅速增加，以达到最大值-收缩波的峰值结束。
波的下降部分，称为大灾变，对应于血液供应减少的阶段。 Catacrota有一个更复杂的形式，可能包括额外的元素，如切口和双质子波，但是，这些元素在本演示中没有分析。
简化形态学参数通常用于光体积描记图的实际分析，可以可靠地自动计算。:
*收缩压波高是最小值和随后的最大值之间的振幅差;
*anacrota的持续时间是从最小到最大的时间间隔;
*大灾变的持续时间是从最大值到下一个最小值的时间间隔;
*脉冲周期持续时间是两个连续最小值之间的时间间隔。
在对光体积描记图进行数字处理后，我们将分析先前考虑的脉搏波的形态特征。 为此，我们将创建阵列，用于存储单个脉冲周期的幅度和时间参数。
然后，在循环中依次考虑各个脉冲周期，其中的每一个被定义为两个相邻最小值之间的信号的一段。 对于每个这样的区域，分配最大脉搏波，之后计算主要形态参数：脉搏波的幅度，anacrota和catacrota的持续时间以及脉搏周期的持续时间。 计算的正确性也被检查。
# 最小索引数组
locs_min = sort(Int.(locs_min))
heights = Float64[] # 收缩压波高
t_ana = Float64[] # 阿纳克罗塔的持续时间,与
t_cat = Float64[] # 大灾难的持续时间，与
t_cyc = Float64[] # 心动周期的持续时间，与
idx_min = Int[] # 周期启动索引
idx_max = Int[] # 周期内的最大索引
idx_min2 = Int[] # 周期结束指数
for k in 1:(length(locs_min) - 1)
i_min = locs_min[k] # 当前最小值的索引
i_min2 = locs_min[k + 1] # 下一个最小值的索引
# 信号在两个极小值之间的部分
sig = y[i_min:i_min2]
# 本网站的最大索引
i_max = i_min - 1 + argmax(sig)
# 脉搏波振幅
h = y[i_max] - y[i_min]
# anacrota的持续时间（从最小到最大）
ta = (i_max - i_min) / fs
# 大灾难的持续时间（从最大到下一个最小值）
tc = (i_min2 - i_max) / fs
# 脉冲周期的持续时间
tp = (i_min2 - i_min) / fs
# 检查值
if h <= 0 || ta <= 0 || tc <= 0 || tp <= 0
continue
end
push!(heights, h)
push!(t_ana, ta)
push!(t_cat, tc)
push!(t_cyc, tp)
push!(idx_min, i_min)
push!(idx_max, i_max)
push!(idx_min2, i_min2)
end
为了便于分析，我们将形成所获得的值的有序表示，并将计算的形态学参数汇总在一个表格中。 为此，为每个脉冲周期分配一个序数，之后创建一个表，其中包含收缩压波的高度，anacrot和catacrot的持续时间以及脉冲周期的持续时间。
n = collect(1:length(heights))
cycles = DataFrame(
周期数=n,
Height_syst_vols=高度,
Anacrota_c=t_ana,
大灾变=t_cat,
周期长度_c=t_cyc
)
使用的材料
这项工作使用了从开源PhysioNet获得的真实实验性光电容积图记录。 将[脉冲传输时间（PPG）数据库]作为源数据（https://physionet.org/content/pulse-transit-time-ppg/1.1.0 /)。
结论
在这个例子中，考虑了使用EngeeDSP库的光电容积图数据的数字处理。
示例的工作包括以下主要步骤:
***数据下载：**初始光体积描记图数据是从包含来自open PhysioNet数据库的真实实验记录的CSV文件中读取的。
***预处理：**使用功能 EngeeDSP.Functions.detrend 信号的线性趋势已被去除。 为了抑制高频振动，使用低通滤波器，使用该功能设计 EngeeDSP.Functions.fir1 并通过实施 EngeeDSP.Functions.filter.
***光谱分析：**使用函数 EngeeDSP.Functions.fft 和 EngeeDSP.Functions.fftshift 对滤波之前和之后的信号的频率组成进行了分析，这使得可以直观地评估滤波器对光体积描记图的光谱的影响。
***特征点的检测：**该函数用于突出脉搏波的局部极大值和极小值 EngeeDSP.Functions.findpeaks. 将发现的特征点在时域中可视化。
***心率的计算：**根据连续脉搏峰值之间的间隔，计算所选记录片段的平均心率。
***脉搏波的形态学分析：**光电容积描记图被分割成单独的脉搏周期。 为每个周期计算形态学参数，包括收缩波的高度，anacrota和catacrota的持续时间以及脉冲周期的总持续时间。
***结果的介绍：**计算的参数使用结构表 DataFrame.