残差函数对自适应ECG滤波的影响
差异函数确定自适应滤波器如何调整网络干扰估计。 示例考虑当前残基之间的差异和相邻残基之间的差异。 模型心电图(ECS)显示了这些方法如何影响波形的恢复以及滤波器对等值线位移的灵敏度。
显示设置
功能 engee.clear() 清理工作区:
engee.clear()
两种图形渲染模式可用于可视化结果。:
-对于静态快速渲染，它被使用 gr().
-对于具有缩放和查看值的能力的交互式可视化，它被使用 plotlyjs().
通过注释掉第二种模式来选择其中一种模式。 此示例使用交互模式。
# gr()
plotlyjs()
安装必要的库
库用于生成心率变异性。 Random. 在连接之前检查它的存在。 如果缺少库，则执行 Pkg.add(). 功能 fft 它从AnyMathDSP库中用于计算ECS的幅度谱。
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["Random"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
using Random
import EngeeDSP.Functions: fft
心电信号模型的形成
正在使用源文件 ecg_generator.jl，位于脚本旁边。 该函数将ECS形成为心循环的非对称高斯分量之和。 心率因周期而略有不同。 固定价值 seed 确保再现性。
信号在整个60秒间隔内生成一次。 参数 view_start 和 view_duration 仅确定显示的区域并且不改变计算持续时间。
include(joinpath(@__DIR__, "ecg_generator.jl"))
fs = 1000.0 # 采样率，Hz
duration = 60.0 # 模拟的持续时间，s
hr = 75.0 # 平均心率，节拍/分钟
hr_sigma = 3.0 # 心率标准偏差，节拍/分钟
phase0 = 0.16 # 心周期内的初始位置，具有
t, ecg = generate_ecg(
fs=fs, # 采样率，Hz
duration=duration, # 录音时间，s
hr=hr, # 平均心率，节拍/分钟
amplitude=1.0, # 振幅比例因子
variable_hr=true, # 启用心率变异性
hr_sigma=hr_sigma, # 心率标准偏差，节拍/分钟
phase0=phase0, # 心周期内的初始位置，具有
seed=1234 # 用于再现性的发生器晶粒
)
view_start = 0.0 # 显示的部分的开始，从
view_duration = 5.0 # 显示的部分的持续时间，与
idx = (t .>= view_start) .& (t .< view_start + view_duration)
plot(t[idx], ecg[idx],
xlabel="时间，从",
ylabel="振幅，mV",
label="前",
title="模型心电信号",
ylims=(-0.5, 1.25),
legend_position=:outertopright
)
示波图显示重复的心脏周期。 记录在R波附近开始，因此有用信号的第一个采样明显影响滤波器的初始瞬态。
网络干扰仿真
对于有用的EX，我们添加谐波干扰。:
这里 mV和 赫兹。 产生的输入数组将是两个滤波器共同的。
f0 = 50.0 # 网络干扰的额定频率，Hz
A_net = 0.5 # 网络干扰幅度，mV
d = 0.002 # 自适应校正步骤，mV
net = A_net .* sin.(2π .* f0 .* t)
ecg_noisy = ecg .+ net
p = plot(t[idx], ecg[idx],
xlabel="时间，从",
ylabel="振幅，mV",
label="原前",
title="网络干扰问题",
ylims=(-1.0, 1.5),
legend_position=:outertopright
)
plot!(p, t[idx], ecg_noisy[idx], label="干扰问题")
网络波动叠加在ECS的牙齿上，在心周期之间尤其明显。 自适应滤波器应评估该谐波并将其从输入中减去。
自适应滤波原理
可以从前两个计数预测谐波信号。 因此，下一个干扰值的初步估计计算为
预测后，确定差异。 让我们比较这两个选项:
- 活期存款: .
- 余额差额: .
在第一种情况下，适配对余数的绝对值作出反应，包括有用的EX。 在第二种情况下，余数的恒定分量被减小，并且缓慢的变化被减弱。 然而，QRS复合体的快速变化会影响这两种选择。 此外，恒定的偏移改变了初始条件，因此差异变体并不完全独立于它。
在这两种情况下，预测都会根据差异标志进行调整。:
步骤 mV确定调整的速度。 一个大的步骤加速了校正，但增加了评估和对有用信号的反应的波动。 初始等级设置为原始示例: , . 因此，前两个输出计数为零。
经常性预测和差异差异对应于Frolov等人的公式（4.3）和（4.4）手册。，第49页。 估计校正和干扰减法在那里由公式(4.5)和(4.6)描述。 代码每个计数保留一个校正步骤，如原始演示示例AnyMath.
在当前余额上有差异的过滤器
功能 adaptive_filter_direct() 使用当前余额:
它返回清除的EX、干扰分数和差异。 对于前两个样本，差异没有定义，用 NaN.
# f0为调谐频率，Hz;d为校正步长，mV。
function adaptive_filter_direct(x, fs; f0=50.0, d=0.002)
N = length(x)
e = zeros(Float64, N)
y = zeros(Float64, N)
mismatch = fill(NaN, N)
for k in 1:min(2, N)
e[k] = x[k]
y[k] = x[k] - e[k]
end
c = cos(2π * f0 / fs)
for n in 2:N-1
e_pred = 2 * c * e[n] - e[n-1]
mismatch[n+1] = x[n+1] - e_pred
if mismatch[n+1] > 0
e[n+1] = e_pred + d
elseif mismatch[n+1] < 0
e[n+1] = e_pred - d
else
e[n+1] = e_pred
end
y[n+1] = x[n+1] - e[n+1]
end
return (y=y, e=e, mismatch=mismatch)
end
具有残差的滤波器
功能 adaptive_filter_difference() 比较相邻样本的剩余部分:
预测、修正步骤、初始条件和返回值与第一个选项相同。
# f0为调谐频率，Hz;d为校正步长，mV。
function adaptive_filter_difference(x, fs; f0=50.0, d=0.002)
N = length(x)
e = zeros(Float64, N)
y = zeros(Float64, N)
mismatch = fill(NaN, N)
for k in 1:min(2, N)
e[k] = x[k]
y[k] = x[k] - e[k]
end
c = cos(2π * f0 / fs)
for n in 2:N-1
e_pred = 2 * c * e[n] - e[n-1]
mismatch[n+1] = (x[n+1] - e_pred) - (x[n] - e[n])
if mismatch[n+1] > 0
e[n+1] = e_pred + d
elseif mismatch[n+1] < 0
e[n+1] = e_pred - d
else
e[n+1] = e_pred
end
y[n+1] = x[n+1] - e[n+1]
end
return (y=y, e=e, mismatch=mismatch)
end
过滤结果
让我们将这两个选项应用于同一个输入数组。 调谐频率、音高和初始条件相同。 图表显示了前5秒，包括初始过渡。
y_direct, e_direct, mismatch_direct = adaptive_filter_direct(ecg_noisy, fs; f0=f0, d=d)
y_difference, e_difference, mismatch_difference = adaptive_filter_difference(ecg_noisy, fs; f0=f0, d=d)
p = plot(t[idx], ecg[idx],
xlabel="时间，从",
ylabel="振幅，mV",
label="原前",
title="回收EX的比较",
legend_position=:outertopright
)
plot!(p, t[idx], y_direct[idx], label="活期存款")
plot!(p, t[idx], y_difference[idx], label="余额差额")
在记录开始时，干扰估计包含第一样本的有用分量。 随着您的适应，此错误会减少。 要查看另一个部分，请更改 view_start 并重新执行细胞，从EX的形成开始。 模拟持续时间将保持60秒。
p = plot(t[idx], net[idx],
xlabel="时间，从",
ylabel="振幅，mV",
label="真正的障碍",
title="网络干扰评估",
legend_position=:outertopright
)
plot!(p, t[idx], e_direct[idx], label="活期存款")
plot!(p, t[idx], e_difference[idx], label="余额差额")
与已知干扰的比较显示了初始过渡过程和估计值的偏差。 要详细查看各个周期，您可以使用交互式图表放大所需的区域。
差异函数的比较
让我们构造两个滤波器残差的时间依赖关系。 它们允许您查看当前余额和余额差异如何对缓慢和快速变化的ECS部分做出反应。 该比较解释了校正如何工作，并且根据EX的误差单独评估恢复的准确性。
p = plot(t[idx], mismatch_direct[idx],
xlabel="时间，从",
ylabel="差异，mV",
label="活期存款",
title="形成差异的两种方法",
legend_position=:outertopright
)
plot!(p, t[idx], mismatch_difference[idx], label="余额差额")
差异差异对慢分量的贡献较小。 它的幅度不能与作为质量指标的直接差异进行比较：这些值以不同的方式定义。 质量由EX恢复的误差来评估。
error_direct = y_direct .- ecg
error_difference = y_difference .- ecg
p = plot(t[idx], error_direct[idx],
xlabel="时间，从",
ylabel="错误，mV",
label="活期存款",
title="EX恢复错误",
legend_position=:outertopright
)
plot!(p, t[idx], error_difference[idx], label="余额差额")
# 初始瞬态后的最大模误差。
idx_after_transient = t .>= 2.0 # 误差估计间隔，2-60秒
println("没有异亮氨酸的位移。 2-60s范围内的最大绝对误差，mV")
println("活期存款: ", round(maximum(abs.(error_direct[idx_after_transient])), digits=5))
println("余额差额: ", round(maximum(abs.(error_difference[idx_after_transient])), digits=5))
在不移动等值线的情况下，两个滤波器均可抑制网络波动。 最大恢复误差估计在2-60秒的范围内，没有过渡过程的初始阶段。 QRS复合体附近的误差尖峰反映了对EC快速变化的适应响应。
振幅谱
该函数用于 spect_amp() 从原来的例子。 计算单侧幅度谱，由样本数归一化。 对5秒的显示限制仅适用于时间图。
function spect_amp(x, fs)
N = length(x)
X = fft(x)
K = div(N, 2) + 1
A = abs.(X[1:K]) ./ N
iseven(N) ? (A[2:end-1] .*= 2) : (A[2:end] .*= 2)
f = collect(0:K-1) .* fs ./ N
return f, A
end
f, A_clean = spect_amp(ecg, fs)
_, A_noisy = spect_amp(ecg_noisy, fs)
_, A_direct = spect_amp(y_direct, fs)
_, A_difference = spect_amp(y_difference, fs)
idx_f = f .<= 100.0 # 显示频谱的上限，Hz
p = plot(f[idx_f], A_clean[idx_f],
xlabel="频率，赫兹",
ylabel="振幅，mV",
label="原前",
title="过滤前后的光谱",
legend_position=:outertopright
)
plot!(p, f[idx_f], A_noisy[idx_f], label="干扰问题")
plot!(p, f[idx_f], A_direct[idx_f], label="活期存款")
plot!(p, f[idx_f], A_difference[idx_f], label="余额差额")
降低50Hz线可以抑制干扰。 然而，有用的EX也包含该区域中的分量，因此频谱图应该与相对于已知纯信号的误差一起考虑。
用isoline的偏移过滤EX
让我们为相同的模型值添加一个常量偏移量:
两个滤波器的输入等于 . 参数 offset 在mV中设置一个恒定的偏移量. 的网络干扰保持不变。 滤波器应该抑制它并保持有用信号的偏移。
offset = 0.3 # 等值线的恒定位移EX，mV
ecg_offset = ecg .+ offset
ecg_offset_noisy = ecg_offset .+ net
offset_direct = adaptive_filter_direct(ecg_offset_noisy, fs; f0=f0, d=d)
offset_difference = adaptive_filter_difference(ecg_offset_noisy, fs; f0=f0, d=d)
p = plot(t[idx], ecg_offset[idx],
xlabel="时间，从",
ylabel="振幅，mV",
label="参考EX",
title="等值线的位移为+0.3mV",
legend_position=:outertopright
)
plot!(p, t[idx], offset_direct.y[idx], label="活期存款")
plot!(p, t[idx], offset_difference.y[idx], label="余额差额")
当EX被移位时，滤波器中仍然有明显的波动和当前的剩余部分。 差值版本相对于偏移基准更好地保留波形。 这是因为在计算余数差时减小了常数分量。
要检查结果，错误如下所示 ecg_offset. 所添加的偏移被包括在参考信号中并且不被认为是干扰。
error_offset_direct = offset_direct.y .- ecg_offset
error_offset_difference = offset_difference.y .- ecg_offset
p = plot(t[idx], error_offset_direct[idx],
xlabel="时间，从",
ylabel="错误，mV",
label="活期存款",
title="错误恢复具有偏移量的EX",
legend_position=:outertopright
)
plot!(p, t[idx], error_offset_difference[idx], label="余额差额")
println("与等值线+0.3mV的偏移。 2-60s范围内的最大绝对误差，mV")
println("活期存款: ", round(maximum(abs.(error_offset_direct[idx_after_transient])), digits=5))
println("余额差额: ", round(maximum(abs.(error_offset_difference[idx_after_transient])), digits=5))
当等值线移位时，直接版本保留明显的误差波动。 差变型对余数的恒定分量反应较不强烈。 对于2-60s的相同间隔给出了最大误差。
残差函数选择的影响
这些图表表明，差异的形成方式改变了残余误差的性质。 当前余数考虑了有用信号的绝对值并且即使具有干扰的精确预测也可以指导校正。 残基的差异降低了恒定成分，但QRS复合物的快速变化继续影响适应。
当异亮氨酸移位时，变体之间的差异变得更加明显。 在2-60s范围内的最大绝对误差表征在暂态过程的初始阶段排除之后恢复的EC与参考的最大偏差。 该指标没有描述平均误差和保留波形的所有特征。
获得的结果显示了具有EX和过滤器的指定参数的两个标准的影响。 要在另一项任务中选择准则，必须考虑有用信号的形状，干扰的性质以及对可接受失真的要求。
使用过的来源
- Frolov A.V.，Dik S.K.，Likhachevsky D.V.，Davydov M.V.，Chechik N.M.生物医学信号和图像的数字处理。 明斯克：BGUIR，2016。
- 心电信号小区网络干扰自适应滤波AnyMath.