Engee 文档
Notebook

使用连续辐射雷达和深度学习的心电监测（上）

导言

现代雷达越来越多地用于非接触监测重要功能，其中主要是呼吸和心脏活动。 与接触登记（ECG）不同，雷达不需要与身体直接相互作用，并且在存在障碍物和患者运动的情况下也起作用。

心电图是诊断学的"黄金标准"，但传统的切除对长期注册来说是不方便的。 该演示考察了从雷达信号重建心电图形状的任务，该雷达信号捕获胸部微运动和心脏活动的混合物。 在这种解释中，胸部的运动是一个障碍。 结果，任务减少到在噪音和呼吸的背景下突出微弱的心脏信号。

演示项目由两部分组成：在本演示示例中，分析和比较ECG信号和反射回波雷达的信号（峰值、间隔、心率），在第二个演示示例中，训练一个神经网络，从雷达信号重建ECG的形状，然后将其与参考ECG信号进行比较。

关于数据

该论文以Schellenberger等人提出的开放数据集为例。 （科学数据，2020）。 该数据集包含来自30名健康志愿者的雷达和参考传感器（ECG，阻抗心电图，连续BP）的~24小时同步记录。
从正在调查的数据集中提取了一个变异

导入库

我们从Julia和Python（skimage，scipy，neurokit2）导入必要的库，并导入utils。jl文件到脚本，其中包含所有必要的功能。

In [ ]: 
include("$(@__DIR__)/install_packages.jl");

为了使脚本工作，您需要安装几个Python包，其功能通过PyCall调用。 要安装，请转到命令提示符并单击";"符号。 接下来，shell打开，您需要在其中注册:

pip安装-r要求.txt脧脗脭脴

重要的一个位于工作目录中。

In [ ]: 
# ENV["PYTHON"] = '/user/.user_venvs/py3.11/bin'
using PyCall
using Measures
using Peaks
using Plots
using DSP
using Statistics
using Wavelets
skimage = pyimport("skimage");
nk = pyimport("neurokit2")
scipy = pyimport("scipy")
include("$(@__DIR__)/utils.jl");

上传数据

使用函数从文件下载数据 load_data，在哪里 i -雷达信号的共模分量, q -雷达信号的正交分量, fs_radar -雷达信号的采样频率, ecg_data1, ecg_data2 -使用位于受试者身体上的两个不同电极记录的ECG信号 fs_ecg -心电信号的采样率

In [ ]: 
Dirpath = "$(@__DIR__)"
filepath = raw"GDN0011_1_Resting.mat"
order_filter = 6;
full_path = joinpath(Dirpath, filepath);

i, q, fs_radar, ecg_data1, ecg_data2, fs_ecg = load_data(full_path);
cur_ecg = vec(ecg_data2); # сигнал ЭКГ, который будем анализировать

我们可视化原始雷达和心电信号

In [ ]: 
p = plot(
  layout        = (4,1),
  size          = (800,800),
  margin        = 5mm,
  bottom_margin = 15mm,     
)

plot!(p, i[1:10000],     subplot=1, label="I",   xlabel="Отсчет",        ylabel="Амплитуда, В")
plot!(p, q[1:10000],     subplot=2, label="Q",   xlabel="Отсчет",        ylabel="Амплитуда, В")
plot!(p, ecg_data1[1:10000], subplot=3, label="1 отведение", xlabel="Отсчет", ylabel="U, мВ")
plot!(p, ecg_data2[1:10000], subplot=4, label="2 отведение", xlabel="Отсчет", ylabel="U, мВ")
Out[0]:

让我们看看没有预处理原始信号的雷达数据的相位差和胸部位移是什么样的。

In [ ]: 
complex_radar_Data = Complex.(i, q);
sigma = angle.(complex_radar_Data);
lambda_ = 12.5e-3;
range_ = unwrap(sigma) ./ (4 * pi) * lambda_;
plot(layout = (2, 1), size=(600,600),margin=5mm)

plot!(sigma[10000:20000] .* (180/π), subplot=1, label="Разность фаз", xlabel="Отсчет", ylabel="Разность фаз, градусы", legendfontsize = 10)
plot!(range_[10000:20000] .*1e3, subplot=2, label="Смещение грудной клетки", xlabel="Отсчет", ylabel="Расстояние, мм", legendfontsize = 10)
Out[0]:

雷达信号分析

功能 EllipseReconstruction(i, q) 标准化雷达的IQ信号：用椭圆近似它，居中，旋转它，并将其缩放成一个圆圈，然后提取相位并将其转换为目标的微位移。 这使您可以获得干净的位移信号，然后使用该功能隔离呼吸和脉搏 getVitalSignals.

In [ ]: 
I_comp, Q_comp, sigma, range_ = EllipseReconstruction(i, q);
radar_respiration, radar_pulse, radar_heartsound, radar_heartsound_denoise = getVitalSignals(range_, fs_radar, order_filter);

由于信号足够大，我们在间隔中可视化提取的脉搏，呼吸和愤怒的声音以进行视觉评估。 [1, 10000]

In [ ]: 
plot(layout = (2, 2), size=(800,500), margin=5mm)

plot!(radar_respiration[1:10000], subplot=1, label="Дыхание", xlabel="Отсчет", ylabel="Амплитуда, мВ")
plot!(radar_pulse[1:10000], subplot=2, label="Пульс", xlabel="Отсчет", ylabel="Амплитуда, мВ")
plot!(radar_heartsound[1:10000], subplot=3, label="Звуки сердца", xlabel="Отсчет", ylabel="Амплитуда, мВ")
plot!(radar_heartsound_denoise[1:10000], subplot=4, label="Звуки сердца после фильтра", xlabel="Отсчет", ylabel="Амплитуда, мВ")
Out[0]:

左上图显示了一个缓慢正弦形状的呼吸周期，心脏活动声音的雷达信号（下图）具有明显的高频峰值。 右上图显示了脉搏分量和经过滤波的心率声音。

搜索雷达信号上的R峰值

接下来，我们计算胸部运动的加速度，在带上进行过滤 [4, 20] 赫兹

In [ ]: 
acc = diff(diff(radar_pulse)) * fs_radar^2;
radar_pulse_acc = Filter(acc, 4, 20, fs_radar, 6);
time_acc = range(0, step=1/fs_radar, length=length(radar_pulse_acc));

我们在间隔处可视化胸部加速度信号 [10000, 40000] 赫兹

In [ ]: 
plot(time_acc[10000:40000], radar_pulse_acc[10000:40000], xlabel="Время (с)", ylabel="Ускорение, м/с^2", title="Ускорение грудной клетки", label="")
Out[0]:

接下来，使用希尔伯特变换从雷达脉冲信号计算振幅包络，之后搜索心脏峰值—仅搜索那些足够高（高于中位数的20％）并且间隔至少0.76秒的峰值。

In [ ]: 
env = abs.(hilbert(radar_pulse_acc));
radar_peacks, _ = scipy.signal.find_peaks(env,
                      distance=Int(fs_radar * 0.76), height=median(env) * 0.2);

让我们显示雷达信号的片段：黑线显示包络，橙线显示脉冲分量（胸部加速度），红点标记所选时间窗口内找到的R峰。

In [ ]: 
time_start = 80000; time_end = 100000;
r_subset = filter(x -> time_start ≤ x ≤ time_end, radar_peacks)
radar_part = env[time_start:time_end]
time_part = time_acc[time_start:time_end]

time_peaks = time_acc[r_subset]

plot(time_part, radar_part;
    color = :black,
    lw = 1.5,
    label = "Огибающая",
    xlabel = "Время, с",
    ylabel = "Амплитуда, м/с^2",
    title = "Полный фрагмент"
)

plot!(time_part, radar_pulse_acc[time_start:time_end];
    label = "Ускорение грудной клетки"
)

scatter!(time_peaks, env[r_subset];
    marker = (:circle, 6),
    mc = :red,
    label = "R-пики"
)
Out[0]:

心电图分析

接下来，我们将分析可用的心电信号。

首先从ECG信号中去除噪声 (ecg_clean)，然后找到R峰的位置 (ecg_peaks)，之后执行峰值搜索 P, Q, S 和 T 使用方法 peak

In [ ]: 
ecg_cleaned = nk.ecg_clean(cur_ecg, sampling_rate=fs_ecg);
signals, rpeaks = nk.ecg_peaks(ecg_cleaned, sampling_rate=fs_ecg);
_, waves = nk.ecg_delineate(ecg_cleaned, rpeaks, sampling_rate=fs_ecg, method="peak");

可视化具有峰值的ECG部分 P, Q, R, S, T 使用文件中定义的自定义函数的波 utils.py

In [ ]: 
plot_ecg_with_peaks(ecg_cleaned, waves, rpeaks; range_start=90000, range_end=100000)
Out[0]:
In [ ]: 
rr_intervals_ecg = diff(rpeaks["ECG_R_Peaks"]) / fs_ecg;
times_ecg = cumsum([0.0; rr_intervals_ecg]);

rr_intervals_radar = diff(radar_peacks) / fs_radar;
times_radar = cumsum([0.0; rr_intervals_radar]);

scatter(
    times_radar[2:end], rr_intervals_radar;
    label   = "RR радар",
    xlabel  = "Отсчет",
    ylabel  = "RR",
    title   = "RR радар & ЭКГ",
    color   = :green
)

scatter!(
    times_ecg[2:end], rr_intervals_ecg;
    label = "RR ЭКГ",
    color   = :red
)
Out[0]:
In [ ]: 
print("Кол-во RR интервалов на радарном сигнале: ", length(rr_intervals_radar), "\n")
print("Кол-во RR интервалов на ЭКГ: ", length(rr_intervals_ecg))

Кол-во RR интервалов на радарном сигнале: 621
Кол-во RR интервалов на ЭКГ: 659

在韵律图上，红点（ECG RR）和绿点（雷达RR）大多彼此靠近，这表明间隔的总体良好收敛。 与此同时，雷达具有散射向上发射（高达1.9s），ECG具有向下发射（高达0.6s），这表明在检测过程中罕见的遗漏或误报。 在中间部分（从大约100到400个样品），簇特别密集和一致。 总体而言，RR_radar和RR_ECG之间具有高一致性。

接下来，我们计算雷达和ECG数据的瞬时心率，然后从HR数据计算基本统计数据-平均值，中位数，最小值和最大值，这使我们能够比较平均脉搏和两种方法中

In [ ]: 
HR_radar = 60 .÷ rr_intervals_radar;
HR_ecg = 60 .÷ rr_intervals_ecg;
mean_HRRadar   = mean(HR_ecg)
median_HRRadar = median(HR_ecg)
min_HRRadar    = minimum(HR_ecg)
max_HRRadar    = maximum(HR_ecg)

mean_HRECG   = mean(HR_radar)
median_HRECG = median(HR_radar)
min_HRECG    = minimum(HR_radar)
max_HRECG    = maximum(HR_radar);
In [ ]: 
println("Статистика HR (уд/мин)")
println("──────────────────────────────")
println("             |     ECG     |   Radar")
println("─────────────┼─────────────┼────────────")
println("Mean HR      | ", round(mean_HRECG, digits=2), "       | ", round(mean_HRRadar, digits=2))
println("Median HR    | ", round(median_HRECG, digits=2), "        | ", round(median_HRRadar, digits=2))
println("Min HR       | ", round(min_HRECG, digits=2), "        | ", round(min_HRRadar, digits=2))
println("Max HR       | ", round(max_HRECG, digits=2), "        | ", round(max_HRRadar, digits=2))

Статистика HR (уд/мин)
──────────────────────────────
             |     ECG     |   Radar
─────────────┼─────────────┼────────────
Mean HR      | 58.18       | 62.1
Median HR    | 58.0        | 60.0
Min HR       | 29.0        | 42.0
Max HR       | 77.0        | 112.0

一般来说，一致性还不错-平均和中值心率值接近（≈58vs per62每分钟节拍），但雷达方法具有明显更大的传播：最小值和最大值要高得多，这表明检测中的假峰值或

心率变异性(HRV)分析

您可以获得一整套HRV指标，从简单的统计指标到使用该函数的高级非线性和多分形指标 hrv 从图书馆 neuroit2.

In [ ]: 
hrv_indices = nk.hrv(rpeaks, sampling_rate=fs_ecg, show=false)
Out[0]:
HRV_MeanNN HRV_SDNN HRV_SDANN1 HRV_SDNNI1 HRV_SDANN2 HRV_SDNNI2 HRV_SDANN5 HRV_SDNNI5 HRV_RMSSD HRV_SDSD HRV_CVNN HRV_CVSD HRV_MedianNN HRV_MadNN HRV_MCVNN HRV_IQRNN HRV_SDRMSSD HRV_Prc20NN HRV_Prc80NN HRV_pNN50 HRV_pNN20 HRV_MinNN HRV_MaxNN HRV_HTI HRV_TINN HRV_ULF HRV_VLF HRV_LF HRV_HF HRV_VHF HRV_TP HRV_LFHF HRV_LFn HRV_HFn HRV_LnHF HRV_SD1 HRV_SD2 HRV_SD1SD2 HRV_S HRV_CSI ... HRV_C1a HRV_SD1d HRV_SD1a HRV_C2d HRV_C2a HRV_SD2d HRV_SD2a HRV_Cd HRV_Ca HRV_SDNNd HRV_SDNNa HRV_DFA_alpha1 HRV_MFDFA_alpha1_Width HRV_MFDFA_alpha1_Peak HRV_MFDFA_alpha1_Mean HRV_MFDFA_alpha1_Max HRV_MFDFA_alpha1_Delta HRV_MFDFA_alpha1_Asymmetry HRV_MFDFA_alpha1_Fluctuation HRV_MFDFA_alpha1_Increment HRV_DFA_alpha2 HRV_MFDFA_alpha2_Width HRV_MFDFA_alpha2_Peak HRV_MFDFA_alpha2_Mean HRV_MFDFA_alpha2_Max HRV_MFDFA_alpha2_Delta HRV_MFDFA_alpha2_Asymmetry HRV_MFDFA_alpha2_Fluctuation HRV_MFDFA_alpha2_Increment HRV_ApEn HRV_SampEn HRV_ShanEn HRV_FuzzyEn HRV_MSEn HRV_CMSEn HRV_RCMSEn HRV_CD HRV_HFD HRV_KFD HRV_LZC
0 982.672231 138.257783 17.203166 134.240934 14.416785 136.780236 NaN NaN 174.439493 174.572008 0.140696 0.177515 999.5 105.2646 0.105317 140.25 0.792583 909.5 1080.0 63.429439 84.218513 532.0 1420.0 26.36 445.3125 NaN 0.011477 0.036597 0.070714 0.007943 0.126731 0.517539 0.288778 0.557983 -2.649114 123.441051 151.645602 0.81401 58808.391731 1.228486 ... 0.402532 95.415128 78.31781 0.29099 0.70901 81.802871 127.689776 0.413133 0.586867 88.870007 105.920627 0.731215 2.380941 1.043536 1.665046 -1.553554 -2.035018 -0.238965 0.003628 0.44613 0.526557 0.592637 0.639657 0.637235 0.591616 0.056814 -0.504086 0.000063 0.020011 1.361697 1.581122 8.666462 1.230945 0.872857 1.439033 1.934975 1.740408 1.928947 3.763396 0.952044

1 rows × 91 columns

让我们可视化获得的几个统计数据

In [ ]: 
time_keys  = ["HRV_SDNN", "HRV_RMSSD", "HRV_SDSD", "HRV_pNN50", "HRV_pNN20"]
freq_keys  = ["HRV_VLF", "HRV_LF", "HRV_HF", "HRV_LFHF"]
poinc_keys = ["HRV_SD1", "HRV_SD2"]
geom_keys  = ["HRV_TINN", "HRV_HTI"]
ent_keys   = ["HRV_SampEn", "HRV_ApEn"]


getval(key) = Float64(PyCall.getindex(hrv_indices[key], 1))


time_vals  = [getval(k) for k in time_keys]
freq_vals  = [getval(k) for k in freq_keys]
poinc_vals = [getval(k) for k in poinc_keys]
geom_vals  = [getval(k) for k in geom_keys]
ent_vals   = [getval(k) for k in ent_keys]


p = plot(layout = (1, 2), size=(500, 500), theme=:default, framestyle=:box, grid=true)


bar!(
    p[1, 1],
    1:length(time_keys), time_vals;
    xticks=(1:length(time_keys), time_keys),
    rotation=45,
    xlabel="Временной домен",
    title="Временной HRV",
    legend=false
)


bar!(
    p[1, 2],
    1:length(freq_keys), freq_vals;
    xticks=(1:length(freq_keys), freq_keys),
    rotation=45,
    xlabel="Частотный домен", 
    title="Частотный HRV",
    legend=false)
Out[0]:

结论

在本演示示例中，分析了分析心电信号和预处理雷达数据的方法。