Digital photoplethysmogram processing and pulse wave morphology analysis

The demo shows the process of digital processing of a photoplethysmogram. Using the example of experimental data, the main stages of pulse signal analysis are consistently considered. As part of the example, preprocessing and filtering of the signal, temporal and spectral analysis, detection of characteristic points of the pulse wave, as well as calculation of the heart rate and morphological parameters of the pulse wave are performed.

Installing the necessary libraries

The demo uses open Julia language libraries designed for loading, structuring, and primary processing of experimental data.

Library CSV - a library for working with CSV files. It is used to read and download experimental data from CSV files containing photoplethysmogram recordings.
The Library DataFrames - a library for working with tabular data structures.
The Library Statistics - a library that provides basic statistical functions. It is used to calculate the statistical characteristics of data, for example, the average values used in the analysis of signal parameters.

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

using CSV, DataFrames, Statistics

Photoplethysmogram processing

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

This example uses a CSV file with photoplethysmogram data and associated signals obtained in a sitting position. The original records are taken from the open database [PhysioNet: Pulse Transit Time (PPG) Database](https://physionet.org/content/pulse-transit-time-ppg/1.1.0 /).

We set the path to the file with experimental photoplethysmogram data. Using the construction @__DIR__ allows you to create a relative file path linked to the directory of the current script.

data_path = "$(@__DIR__)/s1_sit.csv"

"/user/DemoPublic/biomedical/ppg_signal_processing/s1_sit.csv"

This step loads the data from the CSV file. Function CSV.read reads the contents of the file using the specified path and converts it into a structure DataFrame.

data = CSV.read(data_path, DataFrame);

In the variable fs The sampling frequency of the signals is set to 500 Hz, in accordance with the description of the database used. Next, the number of photoplethysmogram signal samples is determined and a time vector is formed.

fs = 500.0  # frequency of dysretisation, Hz
N  = length(data.pleth_1)    # number of counts

t = (0:N-1) ./ fs   # time vector, with

0.0:0.002:508.05

Let's plot the photoplethysmogram signal in the time domain.

plot(t, data.pleth_1,
    xlabel = "Time, from",
    ylabel = "The amplitude",
    label = "pleth_1",
    grid = true
)

For the convenience of analysis, we will select a time fragment of the recording in the range from 10 to 70 seconds. According to the set time values, the corresponding sampling indices are calculated, after which the selected area is selected from the initial photoplethysmogram signal.

t_start = 10.0  # the beginning of the fragment, from
t_end = 70.0    # end of the fragment, with

n_start = t_start *fs   # index of the beginning of the fragment
n_end = t_end * fs  # index of the end of the fragment

# selecting a fragment of the FPG signal
x_raw = Float64.(-data.pleth_1[Int64(n_start):Int64(n_end)])

nsampl = length(x_raw)  # number of counts of the fragment
t =  range(t_start, step=1/fs, length=nsampl)   # The exchange vector

10.0:0.002:70.0

From the previously obtained time representation of the analyzed signal, it is possible to notice the presence of a slow change in the signal level. This behavior is typical for photoplethysmographic recordings and may be related to registration conditions.

Therefore, using the EngeeDSP library and the function detrend() the linear trend is removed from the signal. Next, a graph of the signal is plotted after the trend is removed, which simplifies the subsequent analysis of the pulse wave shape.

x = EngeeDSP.Functions.detrend(x_raw, 1)

plot(t, x, 
    label="FPG pleth_1 ",
    grid=true,
    xlabel="Time, from",
    ylabel="The amplitude",
    xlims = (t_start,t_end)
)

The EngeeDSP library is used to analyze the frequency composition of the signal. Using the function fft() a fast Fourier transform of the analyzed signal is calculated, after which the result is centered relative to the zero frequency using the function fftshift(). As a result, a two-way amplitude spectrum of the signal is formed.

A frequency vector is formed based on the sampling frequency and the length of the signal. Next, the amplitude spectrum is calculated and its graph is plotted in a limited frequency range.

fft_1 = EngeeDSP.Functions.fft(x)
fftsh_1 = EngeeDSP.Functions.fftshift(fft_1)

df = fs / nsampl
freq_vec = -fs/2 : df : fs/2 - df

Amp1 = abs.(fftsh_1)./ nsampl

plot(freq_vec, Amp1,
    xlabel="Frequency, Hz",
    ylabel="The amplitude",
    title="The amplitude spectrum",
    grid=true,
    legend=false,
    xlims = (-100, 100)
)

From the graph of the signal after trend removal and from the previously constructed two-way amplitude spectrum, it can be seen that the recording contains high-frequency fluctuations that are not related to the shape of the pulse wave. In this regard, a low-pass filter is used.

To do this, in the EngeeDSP library, use the function fir1() The coefficients of the FIR filter with a cutoff frequency fc = 30 Hz are calculated. Next, the signal filtering is performed by the function filter(). The chart shows a comparison of the signal after trend removal and the signal after filtering.

fc = 30.0   # cutoff frequency, Hz
Wn = fc / (fs/2)    # normalized cutoff frequency

order = 20      # filter FIR order
b = EngeeDSP.Functions.fir1(order, Wn)
a = [1.0]

y = EngeeDSP.Functions.filter(b, a, x)

p = plot(t, x,
        label="FPG",
         xlabel="Time, from",
         ylabel="The amplitude",
         grid=true)
plot!(p, t, y, label="FPG after LPF fc = $fc Hz")

xlims!(p, 10, 20)  
display(p)

After applying a low-pass filter, we note that the high-frequency components of the photoplethysmogram signal are significantly suppressed. As a result, the effect of noise and small fluctuations present in the original signal is reduced. The filtered signal becomes smoother, while the basic shape of the pulse wave remains. Next, we will analyze how the frequency composition of the signal has changed after filtering. To do this, using the EngeeDSP library, we calculate the spectrum of the filtered signal with subsequent centering. Then we will calculate the two-way amplitude spectrum and plot the spectra of the original signal and the signal after filtering on the same graph.

fft_2 = EngeeDSP.Functions.fft(y)
fftsh_2 = EngeeDSP.Functions.fftshift(fft_2)

Amp2 = abs.(fftsh_2)./ nsampl

plot(freq_vec, Amp1,
    xlabel="Frequency, Hz",
    ylabel="The amplitude",
    title="The amplitude spectrum",
    label = "The original FPG",
    grid=true,
    legend=true,
    xlims = (-100, 100)
)
plot!(freq_vec, Amp2,label = "FPG after LPF")

The EngeeDSP library and the function are used to search for characteristic pulse wave points. findpeaks(). Parameter minPeakDistance sets the minimum allowable distance between adjacent peaks in the samples and prevents the repeated detection of peaks within the same cardiac cycle.

First, a search is performed for local maxima of the filtered photoplethysmogram signal. Then the local minima are determined by the inverted signal.

minPeakDistance = 300 # minimum distance between adjacent peaks, counts

# search for local pulse wave maxima
pks_max, locs_max = EngeeDSP.Functions.findpeaks(
    y, out=:data,
    MinPeakDistance = minPeakDistance
)

# search for local minima (by inverted signal)
pks_min, locs_min = EngeeDSP.Functions.findpeaks(
    -y, out=:data,
    MinPeakDistance = minPeakDistance
)

(Ypk = [17.584343668110932, -17.39531901266868, -18.669050763309713, -1.94542846968509, 33.74814704376922, 30.829367013334327, 25.84233903755211, 55.10538690758822, 69.8434389451998, 55.96626074426286  …  83.94463029754985, 84.82209067984081, 74.2801343867382, 88.52240598255462, 76.69576038233366, 84.17858772403335, 86.21498933712168, 97.85916644294964, 89.52541179013434, 88.90142660511627], Xpk = [140, 509, 863, 1242, 1651, 2067, 2464, 2871, 3305, 3721  …  26356, 26756, 27151, 27560, 27960, 28369, 28775, 29189, 29589, 29979], Wpk = [47.32674470024982, 148.10917709249253, 151.45931980024034, 159.00063126558462, 158.54474148630038, 196.83710033698867, 157.79178842043348, 170.28400609777964, 194.41276238628552, 166.93745546182754  …  195.92966275264916, 185.3148943352062, 189.9316956490511, 196.05989918619525, 192.93180929202936, 190.39849655982107, 189.35220067773844, 184.01907248867064, 191.21180624231783, 64.07862430432942], Ppk = [49.7927053863731, 119.12776439623973, 140.21705732007086, 142.2028851101346, 156.92264570424751, 148.04180549297627, 142.9548167946231, 156.68669000047524, 170.38020646509636, 146.51604671055836  …  180.93364223768032, 184.0924373176546, 173.06960737918646, 189.21725033818433, 159.4380899750749, 176.79651341589843, 185.6646479372822, 182.45101470253513, 170.93938020617196, 53.46875268894843])

Let's plot the photoplethysmogram in the time domain and mark the local maxima and minima found on it.

p = plot(t, y,
        label="FPG pleth_1 ",
        grid=true,
        xlabel="Time, from",
        ylabel="The amplitude"
)
scatter!(p, t[locs_max], y[locs_max], label="highs")
scatter!(p, t[locs_min], y[locs_min], label="minimums")

xlims!(p, 10, 20)

To calculate the heart rate (HR), it is necessary to determine the time intervals between consecutive pulse peaks. To do this, the indices of the found maxima are pre-sorted, after which the intervals between neighboring peaks are calculated, taking into account the sampling frequency of the signal. Based on the average value of the obtained intervals, the average heart rate is calculated.

# array of maximum indexes
locs_max = sort(Int.(locs_max))

# calculation of time intervals between neighboring maxima, with
rr_sec = diff(locs_max) ./ fs

# calculation of the average heart rate
hr_bpm = 60.0 / mean(rr_sec)

println("Average heart rate: ", round(hr_bpm, digits=1), " beats/min")

Средняя ЧСС: 74.5 уд/мин

In the first part of the demo, the process of digital processing of a photoplethysmogram was considered, including signal preprocessing, time and frequency domain analysis, and heart rate estimation.

Next, we will consider methods for segmenting photoplethysmograms into individual cardiac cycles and calculating the morphological parameters of the pulse wave, such as the height of the systolic wave, the duration of anacrota and catacrota, as well as the total duration of the pulse cycle.

Analysis of pulse wave morphology

A photoplethysmogram (PHG) reflects pulse changes in blood filling of tissues recorded by the optical method. In the time domain, the AFG is a periodically repeating wave, the shape of which is determined by the phases of the cardiac cycle and the state of the vascular bed.

One pulse cycle is usually considered as a portion of the signal between two consecutive minima. Within the cycle, the main systolic wave and a descending section are allocated, corresponding to the phase of the pulse pressure decrease.

The initial ascending segment of the pulse wave is called anacrota and corresponds to the phase of active arterial blood flow. Anacrota is characterized by a rapid increase in the amplitude of the signal, ending with reaching a maximum - the peak of the systolic wave.

The descending section of the wave, called cataclysm, corresponds to the phase of decrease in blood supply. Catacrota has a more complex form and may include additional elements such as incision and a dicrotic wave, however, these elements are not analyzed in this demo.

Simplified morphological parameters are often used for practical analysis of photoplethysmograms, which can be reliably calculated automatically.:

systolic wave height is the amplitude difference between the minimum and the subsequent maximum;
the duration of anacrota is the time interval from minimum to maximum;
the duration of the cataclysm is the time interval from the maximum to the next minimum;
Pulse cycle duration is the time interval between two consecutive minima.

After digital processing of the photoplethysmogram, we will analyze the previously considered morphological features of the pulse wave. To do this, we will create arrays designed to store the amplitude and time parameters of individual pulse cycles.

Then, individual pulse cycles are sequentially considered in the cycle, each of which is defined as a section of the signal between two adjacent minima. For each such area, the maximum pulse wave is allocated, after which the main morphological parameters are calculated: the amplitude of the pulse wave, the duration of anacrota and catacrota, as well as the duration of the pulse cycle. The correctness of the calculation is also checked.

# array of minimum indexes
locs_min = sort(Int.(locs_min))

heights = Float64[] # systolic wave height
t_ana   = Float64[] # duration of anacrota, with
t_cat  = Float64[]  # duration of cataclysm, with
t_cyc   = Float64[] # duration of the cardiac cycle, with

idx_min  = Int[]    # cycle start index
idx_max  = Int[]    # the maximum index within the cycle
idx_min2 = Int[]    # end-of-cycle index

for k in 1:(length(locs_min) - 1)

    i_min  = locs_min[k]    # the index of the current minimum
    i_min2 = locs_min[k + 1]    # the index of the next minimum

    # the portion of the signal between two minima
    sig = y[i_min:i_min2]

    # the maximum index for this site
    i_max = i_min - 1 + argmax(sig)

    # pulse wave amplitude
    h = y[i_max] - y[i_min]

    # duration of anacrota (from minimum to maximum)
    ta = (i_max - i_min) / fs

    # duration of cataclysm (from maximum to next minimum)
    tc = (i_min2 - i_max) / fs

    # pulse cycle duration
    tp = (i_min2 - i_min) / fs

    # checking the values
    if h <= 0 || ta <= 0 || tc <= 0 || tp <= 0
        continue
    end

    push!(heights, h)
    push!(t_ana, ta)
    push!(t_cat, tc)
    push!(t_cyc, tp)

    push!(idx_min, i_min)
    push!(idx_max, i_max)
    push!(idx_min2, i_min2)
end

For the convenience of analysis, we will form an ordered representation of the obtained values and summarize the calculated morphological parameters in a table. To do this, an ordinal number is assigned to each pulse cycle, after which a table is created containing the height of the systolic wave, the duration of anacrot and catacrot, as well as the duration of the pulse cycle.

n = collect(1:length(heights))

cycles = DataFrame(
    Number of the cycle = n,
    Height_syst_vols = heights,
    Anacrota_c = t_ana,
    Cataclysm_c = t_cat,
    Cycle length_c = t_cyc
)

Materials used

The work used real experimental photoplethysmogram recordings obtained from the open source PhysioNet. The [Pulse Transit Time (PPG) Database] was used as the source data (https://physionet.org/content/pulse-transit-time-ppg/1.1.0 /).

Conclusion

In this example, the digital processing of photoplethysmogram data using the EngeeDSP library was considered.

The work of the example includes the following main steps:

Data Download: The initial photoplethysmogram data was read from a CSV file containing real experimental recordings from the open PhysioNet database.
Pre-processing: Using the function EngeeDSP.Functions.detrend the linear trend of the signal has been removed. To suppress high-frequency vibrations, a low-pass filter is used, designed using the function EngeeDSP.Functions.fir1 and implemented through EngeeDSP.Functions.filter.
Spectral analysis: Using the functions EngeeDSP.Functions.fft and EngeeDSP.Functions.fftshift The analysis of the frequency composition of the signal before and after filtering was performed, which made it possible to visually assess the effect of the filter on the spectrum of the photoplethysmogram.
Detection of characteristic points: The function was used to highlight the local maxima and minima of the pulse wave EngeeDSP.Functions.findpeaks. The found characteristic points were visualized in the time domain.
Calculation of heart rate: Based on the intervals between consecutive pulse peaks, the average heart rate for the selected recording fragment was calculated.
Morphological analysis of the pulse wave: The photoplethysmogram was segmented into separate pulse cycles. Morphological parameters were calculated for each cycle, including the height of the systolic wave, the duration of anacrota and catacrota, as well as the total duration of the pulse cycle.
Presentation of results: The calculated parameters were tabulated using the structure DataFrame.

Row	Номер_цикла	Высота_сист_волы	Анакрота_с	Катакрота_с	Длительность_цикла_с
	Int64	Float64	Float64	Float64	Float64
1	1	154.107	0.124	0.614	0.738
2	2	141.491	0.12	0.588	0.708
3	3	148.613	0.114	0.644	0.758
4	4	142.203	0.128	0.69	0.818
5	5	150.961	0.118	0.714	0.832
6	6	154.004	0.122	0.672	0.794
7	7	142.955	0.12	0.694	0.814
8	8	156.687	0.142	0.726	0.868
9	9	160.393	0.104	0.728	0.832
10	10	156.503	0.118	0.694	0.812
11	11	159.537	0.126	0.7	0.826
12	12	155.849	0.112	0.678	0.79
13	13	151.217	0.13	0.666	0.796
⋮	⋮	⋮	⋮	⋮	⋮
63	63	173.556	0.106	0.74	0.846
64	64	186.218	0.104	0.736	0.84
65	65	187.327	0.136	0.678	0.814
66	66	183.215	0.116	0.684	0.8
67	67	185.517	0.124	0.666	0.79
68	68	173.07	0.128	0.69	0.818
69	69	171.265	0.114	0.686	0.8
70	70	169.314	0.124	0.694	0.818
71	71	187.603	0.116	0.696	0.812
72	72	185.665	0.14	0.688	0.828
73	73	179.273	0.116	0.684	0.8
74	74	174.117	0.118	0.662	0.78