Engee documentation
Notebook

# ECG monitoring using continuous radiation radar and deep learning (Part 1)

## Introduction

Modern radars are increasingly being used for non-contact monitoring of vital functions, the main of which are respiration and cardiac activity. Unlike contact registration (ECG), radar does not require direct interaction with the body, and also works in the presence of obstacles and movements of the patient.

An ECG is the "gold standard" in diagnostics, but its traditional removal is inconvenient for long-term registration. This demonstration examines the task of reconstructing the shape of an ECG from a radar signal that captures a mixture of chest micro-movement and heart activity. In this interpretation, the movement of the chest is a hindrance. As a result, the task is reduced to highlighting weak cardiac signals against the background of noise and breathing.

The demo project consists of two parts: in this demo example, the signs of the ECG signal and the reflected echo radar (peaks, intervals, heart rate) are analyzed and compared, and in the second demo example, a neural network is trained that reconstructs the shape of the ECG from the radar signal and then compares it with the reference ECG signal.

## About the data

The paper uses as an example the open data set presented by Schellenberger et al. (Scientific Data, 2020)​. This dataset contains ~24 hours of synchronous radar and reference sensor recordings (ECG, impedance cardiography, continuous BP) from 30 healthy volunteers.
One variation was extracted from the dataset, which was being investigated

## Importing libraries

We import the necessary libraries from Julia and Python (skimage, scipy, neurokit 2), and also import the utils.jl file into the script, which contains all the necessary functions.

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

For the script to work, you need to install several Python packages, the functions of which are called via PyCall. To install, go to the command prompt and click the ";" symbol. Next, the shell opens, where you need to register:

</no-translation>` pip install -r requirements.txt <no-translation>`

the important one is located in the working directory.

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");

## Uploading data

Download the data from the file using the function load_data, where i - common-mode component of the radar signal, q - quadrature component of the radar signal, fs_radar - sampling frequency of the radar signal, ecg_data1, ecg_data2 - ECG signals recorded using two different electrodes located on the body of the subject fs_ecg - sampling rate of the ECG signal

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); # the ECG signal that we will analyze

We visualize raw radar and ECG signals

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

plot!(p, i[1:10000],     subplot=1, label="I",   xlabel="The countdown",        ylabel="The amplitude, In")
plot!(p, q[1:10000],     subplot=2, label="Q",   xlabel="The countdown",        ylabel="The amplitude, In")
plot!(p, ecg_data1[1:10000], subplot=3, label="1 lead", xlabel="The countdown", ylabel="U, mV")
plot!(p, ecg_data2[1:10000], subplot=4, label="2 lead", xlabel="The countdown", ylabel="U, mV")
Out[0]:

Let's see what the phase difference and chest displacement look like for radar data without preprocessing the original signals.

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="Phase difference", xlabel="The countdown", ylabel="Phase difference, degrees", legendfontsize = 10)
plot!(range_[10000:20000] .*1e3, subplot=2, label="Displacement of the chest", xlabel="The countdown", ylabel="Distance, mm", legendfontsize = 10)
Out[0]:

## Radar signal analysis

Function EllipseReconstruction(i, q) Normalizes the radar's IQ signal: approximates it with an ellipse, centers it, rotates it, and scales it into a circle, after which it extracts the phase and translates it into micro-displacements of the target. This allows you to get a clean displacement signal, from which the respiration and pulse are then isolated using the function 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);

Since the signals are large enough, we visualize the extracted pulse, breathing, and angry sounds in the interval for visual evaluation. [1, 10000]

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

plot!(radar_respiration[1:10000], subplot=1, label="Breath", xlabel="The countdown", ylabel="Amplitude, mV")
plot!(radar_pulse[1:10000], subplot=2, label="Pulse", xlabel="The countdown", ylabel="Amplitude, mV")
plot!(radar_heartsound[1:10000], subplot=3, label="Sounds of the heart", xlabel="The countdown", ylabel="Amplitude, mV")
plot!(radar_heartsound_denoise[1:10000], subplot=4, label="Heart sounds after the filter", xlabel="The countdown", ylabel="Amplitude, mV")
Out[0]:

The upper-left graph shows a breathing cycle with a slow sinusoidal shape, a radar signal of the sound of cardiac activity (below) with pronounced high-frequency peaks. The upper-right graph shows the pulse component and filtered heart rate sounds.

## Search for R peaks on the radar signal

Next, we calculate the acceleration of the chest movement, filter it on the band [4, 20] Hz

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));

We visualize the chest acceleration signal at the interval [10000, 40000] Hz

In [ ]: 
plot(time_acc[10000:40000], radar_pulse_acc[10000:40000], xlabel="Time (s)", ylabel="Ускорение, м/с^2", title="Chest acceleration", label="")
Out[0]:

Next, the amplitude envelope is calculated from the radar pulse signal using the Hilbert transform, after which cardiac peaks are searched for — only those that are high enough (above 20% of the median) and separated by at least 0.76 seconds.

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);

Let's display a fragment of the radar signal: the black line shows the envelope, the orange line shows the pulse component (chest acceleration), and the red dots mark the found R-peaks within the selected time window.

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 = "The envelope",
    xlabel = "Time, from",
    ylabel = "Амплитуда, м/с^2",
    title = "The full fragment"
)

plot!(time_part, radar_pulse_acc[time_start:time_end];
    label = "Chest acceleration"
)

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

## ECG analysis

Next, we will analyze the available ECG signal.

Noise is first removed from the ECG signal (ecg_clean), then the positions of the R-peaks are found (ecg_peaks), after which a peak search is performed P, Q, S and T using the method 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");

Visualize an ECG section with peaks P, Q, R, S, T waves using a custom function defined in the file 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 radar",
    xlabel  = "The countdown",
    ylabel  = "RR",
    title   = "RR Radar & ECG",
    color   = :green
)

scatter!(
    times_ecg[2:end], rr_intervals_ecg;
    label = "RR ECG",
    color   = :red
)
Out[0]:
In [ ]: 
print("Number of RR intervals on the radar signal: ", length(rr_intervals_radar), "\n")
print("Number of RR intervals per ECG: ", length(rr_intervals_ecg))

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

On the rhythmogram, the red dots (ECG RR) and green dots (radar RR) mostly lie close to each other around, which indicates an overall good convergence of the intervals. At the same time, the radar has scattered upward emissions (up to ~1.9 s), and the ECG has downward emissions (up to ~0.6 s), which indicates rare omissions or false alarms during detection. In the middle part (from about 100 to 400 samples), the clusters are particularly dense and consistent. Overall, there is high consistency between RR_radar and RR_ECG.

Next, we calculate the instantaneous heart rate for radar and ECG data, and then basic statistics are calculated from the HR data - the average, median, minimum and maximum values, which allows us to compare the average pulse and the heart rate spread in both methods.

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

In general, the consistency is not bad — the average and median heart rate values are close (≈58 vs ≈62 beats per minute), but the radar method has a noticeably greater spread: the minima and maxima are much higher, which indicates false peaks or gaps in detection.

## Heart Rate Variability (HRV) Analysis

You can get a complete set of HRV indicators, from simple statistical metrics to advanced nonlinear and multi—fractal indexes using the function hrv from the library 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

Let's visualize several statistics obtained

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="Temporary Domain",
    title="Temporary HRV",
    legend=false
)


bar!(
    p[1, 2],
    1:length(freq_keys), freq_vals;
    xticks=(1:length(freq_keys), freq_keys),
    rotation=45,
    xlabel="Frequency domain", 
    title="Frequency HRV",
    legend=false)
Out[0]:

# Conclusions

In this demo example, the method of analyzing ECG signals and preprocessing radar data was analyzed.