Analysis of heart rate variability by ECG signal

The project examines the analysis of heart rate variability based on the ECG signal. Based on the found R-peaks, RR intervals and heart rate are calculated, after which the results are visualized using an RR tachogram, a histogram of RR intervals and a scatterogram. Additionally, a comparison of normal rhythm, bradycardia and tachycardia is performed.

Display Settings

Function engee.clear() cleans up workspaces:

engee.clear()

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

For static, fast rendering, use gr().
For interactive, dynamic visualization with the ability to scale and view values, use plotlyjs().

Select one of the commands by commenting out the second one. Both backends are mutually exclusive, and the last one called will be activated.

# gr()
plotlyjs()

Plots.PlotlyJSBackend()

We will connect using the function include the "read" file.jl" for reading files:

include("$(@__DIR__)/read.jl")

safe_read_ann (generic function with 1 method)

Installing the necessary libraries

At this stage, the necessary Julia libraries are checked and automatically installed if they are not available in the work environment. The example in question uses the library Statistics necessary for performing basic statistical calculations in the analysis of heart rate variability, including processing RR intervals, calculating averages, and estimating data distribution.

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

The library is being activated EngeeDSP and the functions necessary for digital ECG signal processing are imported, including filtering, spectral analysis, peak search, and convolution.

import EngeeDSP.Functions: fft, butter, filter, findpeaks, conv

ECG signal processing and analysis

This section discusses the detection of R-peaks and the analysis of the RR intervals of the ECG signal. Based on the found R-peaks, RR intervals are calculated, after which the basic methods of analyzing heart rate variability are applied.:

RR-tachogram — shows the change in the duration of cardiac cycles over time;
histogram of RR intervals — reflects the distribution of RR intervals on the selected recording fragment;
Scatterogram — allows you to evaluate the relationship between neighboring cardiac cycles.

The combined use of these methods makes it possible to visually analyze the nature of the heart rate and assess its variability.

ECG signal processing function

Let's imagine the function detect_r designed for ECG signal processing and R-peak detection. The main steps described in detail in the example "[Detection of R-peaks on an ECG signal] are performed inside the function (https://engee.com/community/ru/catalogs/projects/detektsiia-r-pikov-na-ekg-signale )": filtering the signal, preparing data for detection, searching for R-peaks, clarifying their position, as well as calculating RR intervals and heart rate.

function detect_r(ecg, fs;
    f_center = 60.0,
    bw = 6.0,
    filter_order = 6,
    smooth_window = 0.020,
    min_peak_distance = 0.30,
    search_window = 0.05
)

    N = length(ecg)

    # Suppression of network interference
    f1 = f_center - bw / 2
    f2 = f_center + bw / 2

    wn1 = f1 / (fs / 2)
    wn2 = f2 / (fs / 2)

    b, a = butter(filter_order, [wn1, wn2], "stop")
    ecg_filt = filter(b, a, ecg)

    # Signal preparation for R-peak detection
    d_ecg = vcat(0.0, diff(ecg_filt))
    sq_ecg = d_ecg .^ 2

    win = max(3, Int(round(smooth_window * fs)))
    kernel = ones(win) ./ win

    det_sig = conv(sq_ecg, kernel)
    det_sig = det_sig[1:N]

    # Search for R-peaks
    thr = mean(det_sig) + 1.5*std(det_sig)
    min_dist = Int(round(min_peak_distance * fs))

    _, locs = findpeaks(
        det_sig,
        out = :data,
        MinPeakDistance = min_dist,
        MinPeakHeight = thr
    )

    # Clarifying the position of the R-peaks
    search_radius = Int(round(search_window * fs))
    r_locs = Int[]

    for idx in locs
        l = max(1, idx - search_radius)
        r = min(N, idx + search_radius)

        seg = ecg_filt[l:r]
        _, imax = findmax(seg)

        push!(r_locs, l + imax - 1)
    end

    r_locs = unique(sort(r_locs))

    # Time and amplitudes of R-peaks
    r_times = (r_locs .- 1) ./ fs
    r_vals = ecg_filt[r_locs]

    # RR-intervals and heart rate
    rr = diff(r_times)
    rr_ms = rr .* 1000
    rr_time = r_times[2:end]
    heart_rate = 60.0 ./ rr

    return ecg_filt, r_times, r_vals, rr_ms, rr_time, heart_rate
end

detect_r (generic function with 1 method)

Signal loading and processing

An ECG fragment is being downloaded from the recording 122 MIT-BIH databases. The first recording channel is selected as the analyzed signal, after which it is transmitted to the function detect_r.

data = "$(@__DIR__)/mitdb/122"

hdr = read_header(data)
fs = Float64(hdr.fs)

# Selecting a signal fragment for analysis
t_start = 10.0 * 60.0       # the beginning of the fragment, from
t_dur   = 10.0 * 60.0      # duration of the fragment, s

# Converting time to countdown numbers
sampfrom = Int(round(t_start * fs))
sampto   = sampfrom + Int(round(t_dur * fs))

# Reading the selected fragment of the ECG signal
_, raw, phys = read_dat_212(data; sampfrom=sampfrom, sampto=sampto)

# Selecting the first ECG signal channel
ecg = phys === nothing ? Float64.(raw[:, 1]) : Float64.(phys[:, 1])

# Formation of the time axis
N = length(ecg)
t = collect(0:N-1) ./ fs

# Signal processing and R-peak detection
ecg_filt, r_times, r_vals, rr_ms, rr_time, heart_rate = detect_r(ecg, fs)

([-0.8085840295479593, -0.6407114219986917, -0.8295744292817472, -1.1115044953253097, -1.218735198190335, -1.0831915944284476, -0.8753605696075456, -0.7877670680815773, -0.8717543420373358, -0.9955104486088345  …  -0.9720761258502096, -0.981535922982182, -0.9813596429430959, -0.974918498829814, -0.9792503169729933, -0.9797968784102892, -0.9649996169615994, -0.9639460503431428, -0.9465707864313421, -0.9284212640303497], [0.225, 0.9722222222222222, 1.7305555555555556, 2.4944444444444445, 3.238888888888889, 3.9583333333333335, 4.736111111111111, 5.4944444444444445, 6.25, 6.963888888888889  …  592.8388888888888, 593.6277777777777, 594.4138888888889, 595.2083333333334, 596.0027777777777, 596.7583333333333, 597.5472222222222, 598.325, 599.1194444444444, 599.9], [0.8073956333793104, 0.8065250858452058, 0.7379285779211336, 0.7061365786747847, 0.7272537751280435, 0.7440087195968599, 0.6796706508104955, 0.7280769766471427, 0.7334193595798135, 0.7459423536834737  …  0.8541559353482163, 0.9772985211844168, 0.8888330583314511, 0.8316965334704701, 0.9346614768511732, 0.9433789791995124, 0.8477715785448744, 0.8200259572797408, 0.9161066785107799, 0.9254050725231524], [747.2222222222223, 758.3333333333334, 763.8888888888888, 744.4444444444445, 719.4444444444446, 777.7777777777773, 758.3333333333337, 755.5555555555555, 713.8888888888886, 750.0  …  780.5555555555657, 788.8888888888914, 786.1111111111541, 794.4444444444798, 794.4444444443661, 755.5555555555884, 788.8888888888914, 777.7777777778283, 794.4444444443661, 780.5555555555657], [0.9722222222222222, 1.7305555555555556, 2.4944444444444445, 3.238888888888889, 3.9583333333333335, 4.736111111111111, 5.4944444444444445, 6.25, 6.963888888888889, 7.713888888888889  …  592.8388888888888, 593.6277777777777, 594.4138888888889, 595.2083333333334, 596.0027777777777, 596.7583333333333, 597.5472222222222, 598.325, 599.1194444444444, 599.9], [80.29739776951672, 79.12087912087911, 78.54545454545455, 80.59701492537313, 83.39768339768338, 77.1428571428572, 79.12087912087908, 79.41176470588235, 84.04669260700393, 80.0  …  76.86832740213424, 76.05633802816877, 76.32508833921844, 75.52447552447217, 75.52447552448297, 79.4117647058789, 76.05633802816877, 77.14285714285214, 75.52447552448297, 76.86832740213424])

Let's construct a filtered ECG signal with the display of the found R-peaks.

plot(t, ecg_filt,
    xlabel="Time, from",
    ylabel="Amplitude, mV",
    label="ECG signal"
)
scatter!(r_times, r_vals, markersize=3, label="R-peaks")

Analysis of RR intervals

Let's calculate the main heart rate indicators: the average RR interval, the standard deviation of the RR intervals, and the average heart rate.

mean_rr = mean(rr_ms)       # average RR interval, ms
std_rr  = std(rr_ms)        # standard deviation of RR intervals, ms
mean_hr = mean(heart_rate)  # average heart rate, beats/min

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), " bpm")

Average RR interval: 742.17 ms
Standard deviation of RR intervals: 37.34 ms
Average heart rate: 81.05 beats/min

RR-tachogram

The RR tachogram reflects the change in the duration of the RR intervals over time, that is, it shows how the duration of consecutive cardiac cycles changes on a selected fragment of the ECG signal. Based on this dependence, it is possible to assess the stability of the heart rate, identify smooth fluctuations in the RR intervals, and also notice sudden changes that may be associated with artifacts, detection errors, or rhythm features.

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

Histogram of RR intervals

The histogram of RR intervals shows the distribution of the duration of cardiac cycles. The values of 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. This visualization allows you to determine near which value the main part of the RR intervals is concentrated, how pronounced their spread is, and whether there are individual intervals that differ markedly from the main group.

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

Scatterogram of RR intervals

The scatterogram of the RR intervals shows the relationship between adjacent cardiac cycles. The current RR interval is plotted on the X-axis. , and on the Y–axis is the following RR interval . This graph allows you to estimate how much the duration of the next cardiac cycle differs from the previous one. If the points are located compactly and close to the diagonal neighboring RR intervals have similar values, and the rhythm is relatively stable. The increased dot spread reflects a more pronounced heart rate variability or the possible presence of irregular contractions.

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

# Axis boundaries for correct display
rr_min = minimum(vcat(rr_n, rr_next))
rr_max = maximum(vcat(rr_n, rr_next))

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

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

Analysis of normal rhythm, bradycardia and tachycardia

This section discusses the comparison of normal rhythm, bradycardia, and tachycardia by RR intervals and heart rate. In general, the normal resting heart rate of an adult is in the range of about 60-100 beats per minute. Bradycardia corresponds to a slow rhythm, in which the heart rate is usually below 60 beats/min, and tachycardia corresponds to a rapid rhythm, in which the heart rate exceeds 100 beats/min.

As part of the example, these types of rhythm are considered as demonstration cases for analyzing ECG signals. For each signal, R-peaks will be detected, RR intervals and average heart rate will be calculated, and an RR tachogram, a histogram of RR intervals, and a scatterogram will be constructed. This will allow us to compare how the duration of cardiac cycles and the distribution of RR intervals change with normal rhythm, bradycardia and tachycardia.

Loading, processing, and displaying ECG signals

For ease of operation, we will create a function load_mitdb, designed to download a selected fragment of an ECG recording from the MIT-BIH database.

function load_mitdb(record_name;
    t_start_min = 5.0,
    t_dur_min = 10.0,
    channel = 1
)
    data = "$(@__DIR__)/mitdb/$record_name"
    # Reading information about a record
    hdr = read_header(data)
    fs = Float64(hdr.fs)
    # Converting time to countdown numbers
    t_start = t_start_min * 60.0
    t_dur = t_dur_min * 60.0
    sampfrom = Int(round(t_start * fs))
    sampto = sampfrom + Int(round(t_dur * fs))
    # Reading the selected fragment of the ECG signal
    _, raw, phys = read_dat_212(data; sampfrom=sampfrom, sampto=sampto)
    # Selecting the ECG signal channel
    ecg = phys === nothing ? Float64.(raw[:, channel]) : Float64.(phys[:, channel])
    # Formation of the time axis
    N = length(ecg)
    t = collect(0:N-1) ./ fs

    return ecg, fs, t
end

load_mitdb (generic function with 1 method)

To demonstrate a normal heart rate, we will upload the recording data. 100 and we will display the result of the processing: a filtered ECG signal with marked R-peaks.

ecg_norm, fs_norm, t_norm = load_mitdb(
    "100";
    t_start_min = 5.0,
    t_dur_min = 5.0,
    channel = 1
)
ecg_filt_norm, r_times_norm, r_vals_norm, rr_ms_norm, rr_time_norm, heart_rate_norm = detect_r(ecg_norm, fs_norm)

p_norm = plot(t_norm, ecg_filt_norm,
    xlabel = "Time, from",
    ylabel = "Amplitude, mV",
    label = "The normal rhythm"
)
scatter!(p_norm, r_times_norm, r_vals_norm,markersize = 2,label = "R-peaks")

We will process the record in the same way 123 selected to demonstrate bradycardia. The graph will show the filtered ECG signal and the found R-peaks, from which the RR intervals are calculated.

ecg_brady, fs_brady, t_brady = load_mitdb(
    "123";
    t_start_min = 5.0,
    t_dur_min = 5.0,
    channel = 1
)
ecg_filt_brady, r_times_brady, r_vals_brady, rr_ms_brady, rr_time_brady, heart_rate_brady = detect_r(ecg_brady, fs_brady)

p_brady = plot(t_brady, ecg_filt_brady,
    xlabel = "Time, from",
    ylabel = "Amplitude, mV",
    label = "Bradycardia"
)
scatter!(p_brady, r_times_brady, r_vals_brady,markersize = 2,label = "R-peaks")

To analyze the rapid heart rate, we use the recording 209. After filtering and detecting the R-peaks, we will display the processed ECG signal, which will visually assess the decrease in the distance between adjacent cardiac cycles.

ecg_tachy, fs_tachy, t_tachy = load_mitdb(
    "209";
    t_start_min = 6.0,
    t_dur_min = 5.0,
    channel = 1
)
ecg_filt_tachy, r_times_tachy, r_vals_tachy, rr_ms_tachy, rr_time_tachy, heart_rate_tachy = detect_r(ecg_tachy, fs_tachy)

p_tachy = plot(t_tachy, ecg_filt_tachy,
    xlabel = "Time, from",
    ylabel = "Amplitude, mV",
    label = "Tachycardia"
)
scatter!(p_tachy, r_times_tachy, r_vals_tachy, markersize = 2,label = "R-peaks")

We will display the processed ECG signals for different types of heart rate.

plot(p_norm, p_brady, p_tachy,
    layout = (3, 1),
    size = (1000, 850),
    margin = 40*Plots.px
)

According to the obtained oscillograms of ECG signals with different types of heart rate, it is possible to notice differences in the distance between neighboring R-peaks. With bradycardia, the intervals between heart contractions increase compared to the normal rhythm, which corresponds to a decrease in heart rate. In tachycardia, on the contrary, the R-peaks are located closer to each other, which indicates an increased heart rate and a decrease in the duration of the RR intervals.

Analysis of RR intervals

For the convenience of displaying numerical indicators, we will form a function print_hr, which calculates and displays the main characteristics of RR intervals and heart rate.

function print_hr(name, rr_ms, heart_rate)
    mean_rr = mean(rr_ms)          # average RR interval, ms
    std_rr  = std(rr_ms)           # standard deviation of RR intervals, ms
    mean_hr = mean(heart_rate)     # average heart rate, beats/min

    println(name)
    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), " bpm")
    println()
end

print_hr("The normal rhythm", rr_ms_norm, heart_rate_norm)
print_hr("Bradycardia", rr_ms_brady, heart_rate_brady)
print_hr("Tachycardia", rr_ms_tachy, heart_rate_tachy)

The normal rhythm
Average RR interval: 771.8 ms
Standard deviation of RR intervals: 43.22 ms
Average heart rate: 77.99 beats/min

Bradycardia
Average RR interval: 1172.63 ms
Standard deviation of RR intervals: 128.66 ms
Average heart rate: 51.79 beats/min

Tachycardia
Average RR interval: 549.12 ms
Standard deviation of RR intervals: 116.93 ms
Average heart rate: 115.26 beats/min

According to the results of the statistical analysis, it can be seen that the values obtained correspond to the expected differences between the types of heart rate. For a normal rhythm, the average RR interval was 771.8 ms, and the average heart rate was 77.99 beats per minute, which corresponds to the range of a normal heart rate.

With bradycardia, the average RR interval increased to 1172.63 ms, and the average heart rate decreased to 51.79 beats/min. This confirms a slowing of the heart rate compared to the normal state. In tachycardia, on the contrary, the average RR interval decreased to 549.12 ms, and the average heart rate increased to 115.26 beats per minute, which corresponds to a rapid heart rate.

Thus, statistical indicators confirm the results of visual analysis of ECG signals.

RR-tachogram

We will construct RR-tachograms for the types of heart rhythm under consideration: normal rhythm, bradycardia and tachycardia.

p_rr_norm = plot(rr_time_norm, rr_ms_norm,
    xlabel = "Time, from",
    ylabel = "RR-interval, ms",
    title = "The normal rhythm",
    legend = false
)

p_rr_brady = plot(rr_time_brady, rr_ms_brady,
    xlabel = "Time, from",
    ylabel = "RR-interval, ms",
    title = "Bradycardia",
    legend = false
)

p_rr_tachy = plot(rr_time_tachy, rr_ms_tachy,
    xlabel = "Time, from",
    ylabel = "RR-interval, ms",
    title = "Tachycardia",
    legend = false
)

plot(p_rr_norm, p_rr_brady, p_rr_tachy,
    layout = (3, 1),
    size = (1000, 850),
    margin = 40*Plots.px
)

Histogram of RR intervals

Let's consider the distribution of RR intervals for the analyzed types of heart rate. Histograms based on a single graph allow you to compare ranges of values and determine where the main part of cardiac cycles is concentrated for each case.

histogram(rr_ms_norm,
    bins = 30,
    alpha = 0.5,
    xlabel = "RR-interval, ms",
    ylabel = "Number of intervals",
    title = "Histograms of RR intervals",
    label = "The normal rhythm"
)

histogram!(rr_ms_brady,
    bins = 30,
    alpha = 0.5,
    label = "Bradycardia"
)

histogram!(rr_ms_tachy,
    bins = 30,
    alpha = 0.5,
    label = "Tachycardia"
)

The histograms show that the distributions of the RR intervals for the analyzed rhythm types occupy different ranges. At a normal rhythm, the main part of the values is concentrated in the region of about 750-800 ms, which corresponds to an average heart rate of about 78 beats/min.

For bradycardia, the distribution is shifted to the region of large RR intervals, approximately 1000-1300 ms. This confirms an increase in the duration of cardiac cycles and a decrease in heart rate. In tachycardia, the distribution, on the contrary, shifts to the region of smaller RR intervals, which corresponds to a rapid heart rate.

Since the analysis is performed on real ECG recordings, the distributions may have an imperfect shape. In particular, two characteristic maxima are noticeable for tachycardia, which indicates the presence of two predominant ranges of RR intervals in the selected recording fragment. This may be due to a change in rhythm within the analyzed area or the characteristics of the actual clinical signal.

Thus, histograms make it possible to visually compare the duration of cardiac cycles and show the differences between normal rhythm, bradycardia and tachycardia according to the distribution of RR intervals.

Scatterogram of RR intervals

Let's consider the scatterograms of RR intervals for the analyzed types of heart rate.

rr_n_norm = rr_ms_norm[1:end-1]
rr_next_norm = rr_ms_norm[2:end]

rr_n_brady = rr_ms_brady[1:end-1]
rr_next_brady = rr_ms_brady[2:end]

rr_n_tachy = rr_ms_tachy[1:end-1]
rr_next_tachy = rr_ms_tachy[2:end]

rr_min = 0
rr_max = 1800

# The normal rhythm
p_scatter_norm = scatter(rr_n_norm, rr_next_norm,
    xlabel = "RR(n), ms",
    ylabel = "RR(n+1), ms",
    title = "The normal rhythm",
    legend = false,
    markersize = 3,
    xlim = (rr_min, rr_max),
    ylim = (rr_min, rr_max)
)

plot!(p_scatter_norm, [rr_min, rr_max], [rr_min, rr_max],
    linestyle = :dash,
    linewidth = 2
)

# Bradycardia
p_scatter_brady = scatter(rr_n_brady, rr_next_brady,
    xlabel = "RR(n), ms",
    ylabel = "RR(n+1), ms",
    title = "Bradycardia",
    legend = false,
    markersize = 3,
    xlim = (rr_min, rr_max),
    ylim = (rr_min, rr_max)
)

plot!(p_scatter_brady, [rr_min, rr_max], [rr_min, rr_max],
    linestyle = :dash,
    linewidth = 2
)

# Tachycardia
p_scatter_tachy = scatter(rr_n_tachy, rr_next_tachy,
    xlabel = "RR(n), ms",
    ylabel = "RR(n+1), ms",
    title = "Tachycardia",
    legend = false,
    markersize = 3,
    xlim = (rr_min, rr_max),
    ylim = (rr_min, rr_max)
)

plot!(p_scatter_tachy, [rr_min, rr_max], [rr_min, rr_max],
    linestyle = :dash,
    linewidth = 2
)

# Displaying scatterograms
plot(p_scatter_norm, p_scatter_brady, p_scatter_tachy,
    layout = (3, 1),
    size = (800, 1400),
    margin = 40*Plots.px
)

For a normal rhythm, the points are located compactly enough and mostly concentrated near the diagonal line. This shows that the neighboring RR intervals have similar values, and the change in the duration of cardiac cycles is moderate.

For bradycardia, the points are shifted to the region of large RR intervals and distributed more widely. There is no clearly defined compact area, which indicates a large variability in the duration of cardiac cycles. This pattern is also consistent with the histogram of RR intervals, where the distribution has a wider range of values.

For tachycardia, two distinct areas of the dots can be seen: one in the area of short RR intervals, the other in the area of longer intervals. This shows that there are two characteristic ranges of cardiac cycle duration in the selected fragment. This result is consistent with the histogram, which also shows two noticeable distribution maxima.

For a more visual comparison, all three tablerograms are additionally displayed on the same graph.

p_scatter_all = scatter(rr_n_norm, rr_next_norm,
    xlabel = "RR(n), ms",
    ylabel = "RR(n+1), ms",
    label = "The normal rhythm",
    markersize = 3,
    xlim = (rr_min, rr_max),
    ylim = (rr_min, rr_max),
    legend_position = :outertopright
)

scatter!(p_scatter_all, rr_n_brady, rr_next_brady,
    label = "Bradycardia",
    markersize = 3
)

scatter!(p_scatter_all, rr_n_tachy, rr_next_tachy,
    label = "Tachycardia",
    markersize = 3
)

plot!(p_scatter_all, [rr_min, rr_max], [rr_min, rr_max],
    linestyle = :dash,
    linewidth = 2,
    label = false
)

The general tablerogram of the RR intervals shows that the three types of rhythm occupy different ranges of values.

Materials used

The work used real experimental electrocardiogram recordings obtained from the open source PhysioNet. The [MIT-BIH Arrhythmia Database] was used as the source data (https://physionet.org/content/mitdb/1.0.0 /). A recording was used for the initial analysis of heart rate variability 122. Additional recordings were used to compare different types of heart rate. 100, 123 and 209 corresponding to demonstration examples of normal rhythm, bradycardia, and tachycardia.

Conclusion

In this example, the analysis of heart rate variability based on the ECG signal was considered. Based on the found R-peaks, RR intervals and heart rate were calculated, after which the main dependencies were constructed to assess the nature of the heart rate.

The work of the example includes the following main steps:

Data download. The initial ECG data was obtained from the open database PhysioNet MIT-BIH Arrhythmia Database. Real experimental electrocardiogram recordings were used for the analysis.

Detection of R-peaks. The function was used to process the ECG signal detect_r, which performs signal filtering, preparation for detection, search for R-peaks and refinement of their position. This function uses the functions of the EngeeDSP library, including butter, filter, conv and findpeaks.

Calculation of RR intervals and heart rate. Based on the time coordinates of the R peaks, RR intervals and heart rate were calculated. Additionally, the average RR interval, the standard deviation of the RR intervals, and the average heart rate were determined.

Analysis of heart rate variability. To assess the nature of the change in heart rate, an RR tachogram, a histogram of RR intervals, and a scatterogram were constructed. These dependences make it possible to estimate the duration of cardiac cycles, the distribution of RR intervals, and the relationship between adjacent intervals.

Comparison of heart rate types. Normal rhythm, bradycardia, and tachycardia were considered to demonstrate the differences. According to the results of the analysis, it can be seen that with bradycardia, the RR intervals increase, and the average heart rate decreases. With tachycardia, on the contrary, the RR intervals become shorter, and the average heart rate increases.

As a result of the example, the sequence of analysis of heart rate variability based on the ECG signal was shown: from the detection of R peaks to the construction of tachograms, histograms and scatterograms of RR intervals.