Electrocardiosignal modeling and heart rate determination

The example considers a block model of the formation of three standard leads of an electrocardiosignal, modeling of interference of the measuring path, digital filtering and streaming detection of R-peaks. Based on the found R-peaks, RR intervals and heart rate are calculated. In the final part, the mode of variable heart rate is considered and the distribution of RR intervals is analyzed.

Display Settings

Function engee.clear() cleans up the workspace:

engee.clear()

Two graph rendering modes are available to visualize the results.:

for static fast rendering, it is used gr();
for interactive visualization with the ability to scale and view values , it is used plotlyjs().

Select one of the modes by commenting out the second one. This example uses interactive mode.

# gr()
plotlyjs()

Plots.PlotlyJSBackend()

Installing the necessary libraries

The example uses the library Statistics necessary for calculating the average value and standard deviation of the RR intervals. Before performing statistical analysis, the library is checked for availability in the working environment.

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

using Statistics

Loading and launching the model

The presented block model implements the complete path of electrocardiosignal formation and processing. At the first stage, three standard leads I, II and III are formed. Next, simulated interference from the measuring path is added to the signal, after which digital filtering is performed. For lead II, after pretreatment, R-peaks are detected, RR intervals are calculated, and heart rate is determined.

The model structure includes the following main subsystems:

Формирование ЭКС – formation of a model electrocardiosignal and three standard leads;
Моделирование помех – adding network components and broadband noise;
Режекторный фильтр – suppression of the 50 Hz network component;
ФНЧ – limitation of the high-frequency component of the signal;
Определение ЧСС – detection of R-peaks, calculation of RR intervals and heart rate.

The model can be run directly from the file. .engee using the simulation environment or run it from an interactive script. In the second case, the model is loaded programmatically, after which the function engee.run starts the calculation and returns the recorded signals.

modelName = "ECG_HR_Model"
modelPath = joinpath(@__DIR__, modelName * ".engee")

model = engee.load(modelPath; force=true);

Formation of an electrocardiosignal

Subsystem Формирование ЭКС generates one model cardiac cycle as the sum of the individual components P, Q, R1, R2, S, ST, and T. This approach is based on the representation of the morphology of the electrocardiosignal by a set of parameterized asymmetric Gaussian functions and corresponds to the general idea of the example of the Engee community "Modeling электрокардиосигнала".

The expression is used for one component

where determines the amplitude of the component, – the position of its extremum inside the cardiac cycle, and defines the width.

To obtain an asymmetrical shape, the width is set separately to the left and right of the extremum.:

The total shape of one cardiac cycle has the form

The model uses fixed time positions and widths of the main components. Therefore, changing the heart rate primarily changes the duration of the interval between adjacent cardiac cycles, rather than stretching the entire P–QRS–T complex.

The period of the cardiac cycle is determined through the heart rate:

where it is set in beats/min, and it is obtained in seconds.

After forming the basic shape, the conditional potentials of the electrodes RA, LA and LL are constructed. Standard leads are calculated as

The Einthoven relation follows directly from these expressions.

The subsystem mask sets the average heart rate. Additionally, a variable rhythm mode is provided. When it is turned on, a value is generated for each new cardiac cycle.

and then the corresponding interval is calculated. . Parameter It is set in beats/min. The default value is 3 beats/min.

The first run uses a constant heart rate of 75 beats per minute.

Running a model with a constant heart rate

At the first stage, we will consider the operation of the block model at a constant heart rate of 75 beats per minute. The variability mode is disabled, which makes it possible to separately analyze the formation of the electrocardiosignal, the effects of interference, preprocessing, and the operation of the R-peak detection algorithm. The duration of the simulation is 10 seconds.

results = engee.run(model; verbose=true);

Building...
Progress 0%
Progress 5%
Progress 10%
Progress 15%
Progress 22%
Progress 27%
Progress 33%
Progress 38%
Progress 43%
Progress 49%
Progress 55%
Progress 63%
Progress 68%
Progress 73%
Progress 78%
Progress 84%
Progress 89%
Progress 95%
Progress 100%

Formed standard leads

We will display leads I, II, and III on a fragment of the model signal.

src = results["ЭКС"]

t = Float64.(src.time)
lead_I   = Float64.(getindex.(src.value, 1))
lead_II  = Float64.(getindex.(src.value, 2))
lead_III = Float64.(getindex.(src.value, 3))

idx = (t .>= 2.0) .& (t .<= 4.0)

p = plot(t[idx], lead_I[idx],
    xlabel="Time, from",
    ylabel="Amplitude, mV",
    label="I",
    title="Standard leads",
    legend_position=:outertopright)

plot!(p, t[idx], lead_II[idx], label="II")
plot!(p, t[idx], lead_III[idx], label="III")

The resulting curves contain the same P–QRS–T sequence, but differ in amplitude due to the different method of lead formation. The most pronounced R-wave is observed in the second lead, which is then used to determine the R-peaks.

Interference simulation

Subsystem Моделирование помех generates a model of distortions characteristic of the measuring path of electrocardiography.

Network interference is associated with electromagnetic interference from the power supply to the electrodes, connecting wires and input stages of measuring equipment. The model takes into account the main component of 50 Hz and two harmonics of 100 and 150 Hz.:

The following parameters are set in the model

Additionally, band-limited white noise is used to simulate the cumulative broadband component of the measuring path. For the block Band-Limited White Noise the parameter is set Cov = 0.2·10⁻⁶.

The total output signal of the subsystem can be represented as

where – the initial electrocardiosignal, and – the noise component.

Interference simulation result

To assess the effect of interference, we compare the initial and distorted lead II.

noisy = results["ЭКС с помехами"]

t_noisy = Float64.(noisy.time)
lead_II_noisy = Float64.(getindex.(noisy.value, 2))

idx_n = (t_noisy .>= 2.0) .& (t_noisy .<= 4.0)

p = plot(t[idx], lead_II[idx],
    xlabel="Time, from",
    ylabel="Amplitude, mV",
    label="Initial II lead",
    title="The effect of simulated interference",
    legend_position=:outertopright)

plot!(p, t_noisy[idx_n], lead_II_noisy[idx_n],label="II Interference lead")

The distorted signal retains the basic morphology of the cardiac cycle, but high-frequency fluctuations appear. This representation allows us to further evaluate the operation of digital filters on a previously known model signal.

Pre-treatment of the electrocardiosignal

The purpose of preprocessing is to attenuate network harmonics and broadband noise, while preserving the shape of the QRS complex necessary for subsequent detection of R-peaks.

The model uses two successive stages: a notch filter and a low-pass filter.

Notch filter

To suppress the main 50 Hz network component, a type I Chebyshev IIR-notch filter is used. In the current model, the resulting filter order is 10, the boundary frequencies are 46 and 54 Hz, and the allowable ripple in the bandwidth is 0.5 dB. The sampling frequency of the model is 1000 Hz.

To quickly obtain the coefficients of such a filter, you can use the built-in Engee application. Проектирование цифровых фильтров. The application selects the type of amplitude-frequency response, synthesis method, order, and frequency parameters. The application is described in the Engee documentation.

The following parameters are used for the notch filter in question:

Sampling rate: 1000 Hz;
type of amplitude-frequency response: notch filter;
filter type: BIH;
synthesis method: Chebyshev type I;
resulting order: 10;
boundary frequencies: 46 and 54 Hz;
Ripple in bandwidth: 0.5 dB.

After synthesis, the calculated coefficients are used to automatically generate the filter in the subsystem. Режекторный фильтр, which allows you to transfer the received parameters directly to the block model.

The result of notch filtering

notch = results["Режекторный фильтр.Y"]

t_notch = Float64.(notch.time)
lead_II_notch = Float64.(getindex.(notch.value, 2))

idx_f = (t_notch .>= 2.0) .& (t_notch .<= 4.0)

p = plot(t_noisy[idx_n], lead_II_noisy[idx_n],
    xlabel="Time, from",
    ylabel="Amplitude, mV",
    label="Before the filter",
    title="Notch filtering",
    legend_position=:outertopright)

plot!(p, t_notch[idx_f], lead_II_notch[idx_f],label="After the notch filter")

The notch filter attenuates the narrowband component in the 50 Hz region. At the same time, harmonics of 100 and 150 Hz and broadband noise require additional bandwidth limitation.

Low-pass filter

By analogy, Butterworth's IIR low-pass filter is synthesized. The model uses an order of 10, a cutoff frequency of 38.13 Hz, and a sampling frequency of 1000 Hz.

The Butterworth filter was chosen because of the smooth amplitude-frequency response in the passband. It additionally suppresses harmonics of 100 and 150 Hz and reduces the high-frequency component of white noise. The cutoff frequency is selected above the main informative low-frequency region of the model electrocardiosignal and allows you to maintain a pronounced QRS complex.

The result of preprocessing

filt = results["ФНЧ.Y"]

t_filt = Float64.(filt.time)
lead_II_filt = Float64.(getindex.(filt.value, 2))

idx_lpf = (t_filt .>= 2.0) .& (t_filt .<= 4.0)

p = plot(t_noisy[idx_n], lead_II_noisy[idx_n],
    xlabel="Time, from",
    ylabel="Amplitude, mV",
    label="With interference",
    title="The result of preprocessing",
    legend_position=:outertopright)

plot!(p, t_filt[idx_lpf], lead_II_filt[idx_lpf],label="After filtering")

After sequential filtering, the high–frequency fluctuations decrease, and the P–QRS-T shape remains distinguishable. The received signal is used as an input to the R-peak detection algorithm.

Detection of R-peaks and calculation of heart rate

To detect R-peaks, lead II is used after pretreatment. The algorithm works in a streaming manner, so processing is performed sequentially for each new sample during the simulation.

First, the difference of neighboring samples is calculated:

The difference operation highlights areas with a high rate of signal change. The QRS complex has steeper fronts compared to the P and T teeth, so its contribution after this operation is more pronounced.

Then the received signal is squared:

Squaring eliminates the sign of the difference and increases the contribution of the counts with a high rate of change.

For smoothing, a moving average of 10 counts is calculated. At a sampling frequency of 1000 Hz, the window length corresponds to 10 ms:

The received signal It is used to highlight the area of the QRS complex.

The detection threshold is formed adaptively. On a window of 1000 samples corresponding to 1 s, the local average is calculated.:

and also the average of the square:

The standard deviation is calculated based on current statistical estimates:

After that, the adaptive threshold is determined.:

The area of a possible QRS complex is determined by the condition:

After selecting the QRS region, the position of the R-peak is clarified by the second lead after filtration. For this, a local maximum is determined. For three consecutive counts, the conditions are checked:

The R-candidate is accepted only if two conditions are met simultaneously: the detection signal exceeds the adaptive threshold, and a local maximum is detected in lead II after filtering.

After the received R-peak, a refractory interval of 0.30 s is activated. During this interval, no new R candidates are accepted, which prevents the re-registration of a single QRS complex.

At the beginning of the simulation, statistical estimates of the adaptive threshold have not yet been fully formed. Therefore, detection is allowed after an initial interval of 0.5 seconds, which reduces the likelihood of false alarms when starting the model.

The RR interval is calculated for two consecutive detected R-peaks.:

after that, the heart rate is determined.:

Until the second correctly detected R-peak appears, the RR interval and the calculated heart rate remain equal to zero.

The overall processing sequence corresponds to the approach considered in the example of the Engee community ["Detection of R-peaks on an ECG signal"] (https://engee.com/community/ru/catalogs/projects/detektsiia-r-pikov-na-ekg-signale ), while in the considered model, detection is performed in a streaming manner using blocks of the modeling environment.

Found R-peaks

r_sig = results["R-пики"]

t_r = Float64.(r_sig.time)
r_event = Float64.(r_sig.value)

r_idx = findall(r_event .> 0.5)

idx_r = (t_filt .>= 2.0) .& (t_filt .<= 6.0)
r_idx_plot = [i for i in r_idx if 2.0 <= t_r[i] <= 6.0]

p = plot(t_filt[idx_r], lead_II_filt[idx_r],
    xlabel="Time, from",
    ylabel="Amplitude, mV",
    label="II lead",
    title="Detection of R-peaks",
    legend_position=:outertopright)

scatter!(p,t_r[r_idx_plot],lead_II_filt[r_idx_plot],label="R-peaks",markersize=4)

The markers are located in the area of the maxima of the R-teeth. One detected R-peak is formed for each cardiac cycle, and the refractory interval excludes the repeated registration of one QRS complex.

RR-intervals and heart rate

hr_model_sig = results["ЧСС модельная"]
hr_measured_sig = results["ЧСС измеренная"]

t_hr_model = Float64.(hr_model_sig.time)
hr_model = Float64.(hr_model_sig.value)

t_hr_measured = Float64.(hr_measured_sig.time)
hr_measured = Float64.(hr_measured_sig.value)

p = plot(t_hr_model, hr_model,
    xlabel="Time, from",
    ylabel="ЧСС, уд/мин",
    label="Preset Heart rate",
    title="Preset and measured heart rate",
    legend_position=:outertopright)

plot!(p, t_hr_measured, hr_measured,label="Measured heart rate")

At the initial stage of the simulation, the measured heart rate is zero. The first detected R-peak sets the beginning of the first RR interval, and its duration can be calculated only after the appearance of the second R-peak. Therefore, the first correct value of the heart rate is formed after the completion of the first complete RR interval.

The additional time offset is due to the delay introduced by serial digital filtering and detection operations. After completing the initial transition section, the measured value is set near the preset level of 75 beats/min.

Modeling of heart rate variability

Consider a regime in which the duration of consecutive cardiac cycles varies over time. Let's keep the average heart rate at 75 beats per minute, turn on the variability, and set the standard deviation to 3 beats per minute.

A new model heart rate value is formed once at the beginning of each cardiac cycle and stored until its completion. For statistical analysis, we will increase the simulation duration to 180 seconds to obtain a sufficient number of RR intervals.

engee.set_param!(modelName * "/Формирование ЭКС",
    "HR" => 75,
    "HR_VAR" => true,
    "HR_SIGMA" => 3.0)

engee.set_param!(modelName, "StopTime" => 180.0)

results_var = engee.run(model; verbose=true);

Building...
Progress 0%
Progress 5%
Progress 10%
Progress 15%
Progress 15%
Progress 20%
Progress 25%
Progress 30%
Progress 32%
Progress 37%
Progress 42%
Progress 47%
Progress 52%
Progress 57%
Progress 62%
Progress 67%
Progress 72%
Progress 77%
Progress 82%
Progress 87%
Progress 92%
Progress 97%
Progress 100%
Progress 100%

Model and measured heart rate

hr_model_var_sig = results_var["ЧСС модельная"]
hr_measured_var_sig = results_var["ЧСС измеренная"]

t_hm = Float64.(hr_model_var_sig.time)
hr_m = Float64.(hr_model_var_sig.value)

t_hd = Float64.(hr_measured_var_sig.time)
hr_d = Float64.(hr_measured_var_sig.value)

p = plot(t_hm, hr_m,
    xlabel="Time, from",
    ylabel="ЧСС, уд/мин",
    label="Model Heart rate",
    title="Model and measured heart rate",
    legend_position=:outertopright)

plot!(p, t_hd, hr_d,label="Measured heart rate")

The measured heart rate changes following the model value, but the update times of the two dependencies do not coincide in time. The model value is set at the beginning of the cardiac cycle, whereas the measured heart rate can be determined only after detecting the next R peak and calculating the completed RR interval.

An additional time offset is introduced by the digital preprocessing path and the sequence of R-peak detection operations. Therefore, the observed displacement is a consequence of the structure of the measuring algorithm and does not indicate an error in calculating the RR interval.

RR interval variability analysis

Based on the sequence of detected R peaks, we will form a set of RR intervals and analyze them. We use three representations:

RR-tachogram – for displaying changes in the duration of cardiac cycles over time;
a histogram of RR intervals – to estimate the distribution of their duration;
a scatterogram – for analyzing the relationship of neighboring intervals and .

A similar set of representations is used in the example of the Engee community "ECG heart Rate Variability Analysis сигналу".

r_var_sig = results_var["R-пики"]
rr_var_sig = results_var["RR-интервал"]

r_event_var = Float64.(r_var_sig.value)
t_r_var = Float64.(r_var_sig.time)
rr_hold = Float64.(rr_var_sig.value)

r_idx_var = findall(r_event_var .> 0.5)

rr_beats = rr_hold[r_idx_var]
rr_time = t_r_var[r_idx_var]

valid_rr = (rr_beats .> 0.3) .& (rr_beats .< 2.0)

rr_beats = rr_beats[valid_rr]
rr_time = rr_time[valid_rr]
rr_ms = 1000 .* rr_beats
heart_rate_rr = 60.0 ./ rr_beats

mean_rr = Statistics.mean(rr_ms)
std_rr = Statistics.std(rr_ms)
mean_hr = Statistics.mean(heart_rate_rr)

println("The average RR interval: ", round(mean_rr, digits=2), " ms")
println("Standard deviation of RR intervals: ", round(std_rr, digits=2), " ms")
println("Average heart rate: ", round(mean_hr, digits=2), " уд/мин")

Средний RR-интервал: 802.98 мс
Стандартное отклонение RR-интервалов: 32.65 мс
Средняя ЧСС: 74.84 уд/мин

RR-tachogram

The RR-tachogram shows the change in the duration of consecutive cardiac cycles over time. Based on this dependence, it is possible to estimate the magnitude and nature of the modeled changes in the RR intervals.

plot(rr_time, rr_ms,
    xlabel="Time, from",
    ylabel="RR-interval, ms",
    title="RR-tachogram",
    legend=false)

The RR intervals vary around the average level corresponding to a given average heart rate. The spread of values is determined by the variability parameter set in the mask of the electrocardiosignal generation unit.

Histogram of RR intervals

The histogram shows the distribution of the duration of cardiac cycles. The RR intervals in milliseconds are plotted on the X-axis, and the number of intervals in the corresponding range is plotted on the Y–axis.

histogram(rr_ms,
    bins=30,
    xlabel="RR-interval, ms",
    ylabel="Number of intervals",
    title="Histogram of RR intervals",
    legend=false)

The main part of the RR intervals is concentrated near the average value. A non-zero standard deviation of the model heart rate leads to the formation of a distribution of the durations of consecutive cardiac cycles.

Scatterogram of RR intervals

The scatterogram shows the relationship between adjacent cardiac cycles. The current interval is plotted on the X-axis. , and on the Y – axis is the following interval . The diagonal line corresponds to the equality of neighboring RR intervals.

rr_n = rr_ms[1:end-1]
rr_next = rr_ms[2:end]

p = scatter(rr_n, rr_next,
    xlabel="RR(n), ms",
    ylabel="RR(n+1), ms",
    title="Scatterogram of RR intervals",
    legend=false,
    markersize=3,
    xlim=(0, 1200),
    ylim=(0, 1200),
    aspect_ratio=:equal)

plot!(p, [0, 1200], [0, 1200],linestyle=:dash,linewidth=2)

With a small simulated variability, the points are compactly located near the diagonal line. An increase in the spread along both axes corresponds to more pronounced differences between neighboring RR intervals.

Conclusion

In the example, a single block chain of forming and processing a model electrocardiosignal is considered. The parametric model generates matched leads I, II, and III, after which network harmonics of 50, 100, and 150 Hz and broadband noise are added to the signal. The notch filter and the low-pass filter pre-process the signal before detecting the R-peaks.

The detection of R-peaks is based on difference processing, squaring, sliding averaging, adaptive threshold, local maximum refinement, and refractory protection. The sequence of the found R-peaks is used to calculate the RR intervals and heart rate.

The mask of the electrocardiosignal formation unit allows you to use a constant or changing heart rate. For a mode with a given variability, the calculated RR intervals are represented as an RR tachogram, histogram, and scatterogram of adjacent intervals.