BVP Signal Analysis - A Complete Tour |
Physiological signals are informational resources containing a great potential of knowledge to be extracted, however, this knowledge is not directly accessible after signal acquisition.
In order to convert information into knowledge, physiological signals should be pre-processed, ensuring noise attenuation from unwanted sources and removal of unexpected artifacts. After pre-processing the signal, researcher can be totally secure that his valuable informational resource is ready to be processed, through the extraction of parameters.
Each one of the extracted parameters define an objective way to analyse the acquired and pre-processed physiological data. Subsequent interpretation of the results/extracted parameters will be extremely important to understand which type of events/mechanisms are taking place inside the human organism and, at this point, an amazing journey is reaching its fundamental stage with the horizon of achievable discoveries becoming observable by the eyes of knowledge.
☌ The current Jupyter Notebook belongs to a set of Notebooks entitled "Complete Tour", where it is presented, in a sequential way, a guide containing recommended steps to be followed, namely: 1) sensor/electrode placement; 2) pre-acquisition notes; 3) pre-processing methods (filtering and motion artifact reduction); 4) detection of specific events and 5) parameter extraction procedures. |
A - Blood Volume Pulse (BVP) sensor | Finger Probe Placement
Taking into consideration the great density of capillaries and great translucency in the hand fingers, these become the ideal site to place the BVP sensor and our finger probe, as demonstrated in the illustrative example of the following figure:B - Pre-Acquisition Requirements
In research literature it is considered that the typical informational band of photoplethysmographic BVP signal is in the range:
Accordingly to these articles, informational content of BVP signal reach, at the most broader criteria, the 10 Hz frequency components. So, applying the widely known Nyquist Theorem , before starting the data collection, the signal acquisition device must be configured to support sampling rates of at least 20 Hz, avoiding aliasing problems during the analog-to-digital conversion procedure.
C - Pre-Processing Stage
C0 - Importing Packages and load signal
# biosignalsnotebooks own package.
import biosignalsnotebooks as bsnb
# Scientific programming package.
from numpy import average, std, diff, max, pad, log2, ceil, absolute, power
from numbers import Number
# Pacakge dedicated to Wavelet decomposition algorithms.
from pywt import swt, iswt
# Deep copy function to clone lists.
from copy import deepcopy
# Load entire acquisition data.
data, header = bsnb.load("../../signal_samples/bvp_motion_artifact.txt", get_header=True)
# Store the desired channel (CH2) in the "signal" variable
signal = data["CH2"]
# Sampling rate definition.
sr = header["sampling rate"]
from numpy import linspace
time = linspace(0, len(signal) / sr, len(signal))
bsnb.plot([time], [signal])
C1 - Filtering Stage 1
In this filtering stage it will be used a conventional bandpass filter, intended to attenuate all signal components outside the typical BVP frequency bands.C1.1 - Generation of filtered signals
# Conventional 2nd order Butterworth bandpass filtering
# [Option RES]
signal_res = bsnb.bandpass(signal, 0.02, 2.1, order=2, fs=sr, use_filtfilt=True)
# Conventional 2nd order Butterworth bandpass filtering
# [Option EXT]
signal_ext = bsnb.bandpass(signal, 0.6875, 10, order=2, fs=sr, use_filtfilt=True)
C1.2 - Comparison between filtered and unfiltered signals
# Generation of power spectra.
# [Original]
freq_axis, power_spect = bsnb.plotfft(signal - average(signal), sr)
# [Option RES]
freq_axis_res, power_spect_res = bsnb.plotfft(signal_res - average(signal_res), sr)
# [Option EXT]
freq_axis_ext, power_spect_ext = bsnb.plotfft(signal_ext - average(signal_ext), sr)
from bokeh.layouts import gridplot
from bokeh.plotting import show
from numpy import min
fig_list_1 = bsnb.plot([time, time, time], [signal, signal_res + average(signal), signal_ext + average(signal)], y_axis_label="RAW Units", legend_label=["Acquired BVP signal", "Filtered signal [0.02; 2.1] Hz", "Filtered signal [0.6875; 10] Hz"], x_range=(time[0], time[-1]),
y_range=(0.80 * min(signal_ext + average(signal)), 1.20 * max(signal_ext + average(signal))), get_fig_list=True, show_plot=False)
fig_list_3 = bsnb.plot([freq_axis, freq_axis_res, freq_axis_ext], [power_spect, power_spect_res, power_spect_ext], y_axis_label="Relative Intensity (a.u.)", x_axis_label="Frequency (Hz)", legend_label=["Acquired BVP signal", "Filtered signal [0.02; 2.1] Hz", "Filtered signal [0.6875; 10] Hz"], get_fig_list=True, show_plot = False)
# Highlight motion artifacts.
from bokeh.models import BoxAnnotation
fill_color_1 = bsnb.opensignals_color_pallet()
motion_art_1 = BoxAnnotation(left=3.8, right=7, fill_color=fill_color_1, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_1, fill_alpha=0.1, legend_label="1st Motion Artifact")
fig_list_1[0].add_layout(motion_art_1)
fill_color_2 = bsnb.opensignals_color_pallet()
motion_art_2 = BoxAnnotation(left=13.5, right=16, fill_color=fill_color_2, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_2, fill_alpha=0.1, legend_label="2nd Motion Artifact")
fig_list_1[0].add_layout(motion_art_2)
show(fig_list_1[0])
# Zoomed figures.
fig_list_2 = bsnb.plot([time[10400:12000], time[10400:12000], time[10400:12000]], [signal[10400:12000], (signal_res + average(signal))[10400:12000], (signal_ext + average(signal))[10400:12000]], y_axis_label="RAW Units", legend_label=["Acquired BVP signal", "Filtered signal [0.02; 2.1] Hz", "Filtered signal [0.6875; 10] Hz"], get_fig_list=True, grid_plot=True, grid_columns=3, grid_lines=1)
#grid_plot_1 = gridplot([[fig_list_1[0]], [fig_list_2[0], fig_list_2[1], fig_list_2[2]]], **bsnb.opensignals_kwargs("gridplot"))
#show(grid_plot_1)
With this first filtering stage we achieved a reasonable result:
Taking this information into consideration, it will be necessary to include an additional filtering stage. As a consequence of the preservation of all relevant information, in the subsequent steps only the signal_ext will be analysed.
The following section is dedicated to minimize the artifacts influence on our BVP signal.
C2 - Filtering Stage 2
After smoothing the BVP signal, removing high-frequency noise components, in the following steps it will be adapted an algorithm proposed by C. M. Lee et al. ( research article ) to minimize the influence of motion artifacts in the BVP signal.C2.1 - Decomposition of signal_ext using Stationary Wavelet Transform (SWT)
SWT has the advantage, when comparing with Discrete Wavelet Transform (DWT), of not producing additional "ringing effects in the vicinity of a discontinuity", accordingly to Chen et al. ( research article ).# SWT 7th level Decomposition using "Daubechies" mother wavelet [https://ieeexplore.ieee.org/abstract/document/6908199].
signal_ext = pad(signal_ext, (0, 2**int(ceil(log2(len(signal_ext)))) - len(signal_ext)), "constant") # Ensure that the number of signal samples is a power of 2.
swt_orig_coeffs = swt(signal_ext[:32768], "db4", level=7)
C2.2 - Removal of motion artifact wavelet coefficients (SWT scaled/approximation coefficients) between decomposition levels 1 and 7 (list entries 0 to 6)
A reference threshold (for detail coefficients) and a range (for scaling/approximation coefficients) are determined for each decomposition level, consisting in:The removed coefficients are related with motion artifacts, taking into consideration that they are outliers outside the acceptability ranges ($[\mu_j^d; +\infty]$ and $[T_{low}; T_{high}]$).
# "level" iteration for loop.
swt_coeffs_copy = deepcopy(swt_orig_coeffs)
for level in range(0, 7):
thr_avg_dt = average(swt_orig_coeffs[level][1]) # Only "detail" coefficients are being taken into account.
thr_avg_sc_low = average(swt_orig_coeffs[level][0]) - 3 * std(swt_orig_coeffs[level][0])
thr_avg_sc_high = average(swt_orig_coeffs[level][0]) + 3 * std(swt_orig_coeffs[level][0])
# "coefficients" iteration loop.
for coeff_nbr in range(0, len(swt_orig_coeffs[level][1])):
if swt_orig_coeffs[level][1][coeff_nbr] > thr_avg_dt: # Motion artifact coefficients.
swt_orig_coeffs[level][1][coeff_nbr] = 0
if swt_orig_coeffs[level][0][coeff_nbr] < thr_avg_sc_low or swt_orig_coeffs[level][0][coeff_nbr] > thr_avg_sc_high: # Motion artifact coefficients.
swt_orig_coeffs[level][0][coeff_nbr] = 0
else: # Storage of noise artifact coefficients in a separate list.
swt_coeffs_copy[level][0][coeff_nbr] = 0
C2.3 - Signal reconstruction through an inverse SWT
rec_signal = iswt(swt_orig_coeffs, "db4")
from numpy import max
fig_list_1 = bsnb.plot([time], [rec_signal[:len(time)] / max(rec_signal[:len(time)])], y_axis_label="Normalized Electrodermal Response (a.u.)", legend_label=["Reconstructed BVP signal"], show_plot=False, get_fig_list=True, x_range=(time[0], time[-1]),
y_range=(-1, 1))
fig_list_2 = bsnb.plot([time], [signal_ext[:len(time)] / max(signal_ext[:len(time)])], y_axis_label="Normalized Electrodermal Response (a.u.)", legend_label=["Peak Artifact Signal"], show_plot=False, get_fig_list=True, x_range=(time[0], time[-1]),
y_range=(-1, 1))
# Highlight motion artifacts.
fill_color_1 = bsnb.opensignals_color_pallet()
motion_art_1 = BoxAnnotation(left=3.8, right=7, fill_color=fill_color_1, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_1, fill_alpha=0.1, legend_label="1st Motion Artifact")
fig_list_1[0].add_layout(motion_art_1)
fill_color_2 = bsnb.opensignals_color_pallet()
motion_art_2 = BoxAnnotation(left=13.5, right=16, fill_color=fill_color_2, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_2, fill_alpha=0.1, legend_label="2nd Motion Artifact")
fig_list_1[0].add_layout(motion_art_2)
motion_art_1 = BoxAnnotation(left=3.8, right=7, fill_color=fill_color_1, fill_alpha=0.1)
fig_list_2[0].circle([-100], [0], fill_color=fill_color_1, fill_alpha=0.1, legend_label="1st Motion Artifact")
fig_list_2[0].add_layout(motion_art_1)
motion_art_2 = BoxAnnotation(left=13.5, right=16, fill_color=fill_color_2, fill_alpha=0.1)
fig_list_2[0].circle([-100], [0], fill_color=fill_color_2, fill_alpha=0.1, legend_label="2nd Motion Artifact")
fig_list_2[0].add_layout(motion_art_2)
show(fig_list_2[0])
show(fig_list_1[0])
In spite of a considerable improve in the quality of our signal, a parcel of the artifact still exists.
Fortunately, the wavelet coefficients of the noise artifact was stored on step C2.2 in the swt_coeffs_copy variable.
Now we can reconstruct the noise artifact component.
C3 - Noise artifact isolation
C3.1 - Signal reconstruction
art_signal = iswt(swt_coeffs_copy, "db4")
fig_list_1 = bsnb.plot([time], [art_signal[:len(time)] / max(art_signal[:len(time)])], y_axis_label="Normalized Electrodermal Response (a.u.)", legend_label=["Reconstructed Noise Artifact signal"], show_plot=False, get_fig_list=True, x_range=(time[0], time[-1]),
y_range=(-1, 1))
# Highlight motion artifacts.
fill_color_1 = bsnb.opensignals_color_pallet()
motion_art_1 = BoxAnnotation(left=3.8, right=7, fill_color=fill_color_1, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_1, fill_alpha=0.1, legend_label="1st Motion Artifact")
fig_list_1[0].add_layout(motion_art_1)
fill_color_2 = bsnb.opensignals_color_pallet()
motion_art_2 = BoxAnnotation(left=13.5, right=16, fill_color=fill_color_2, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_2, fill_alpha=0.1, legend_label="2nd Motion Artifact")
fig_list_1[0].add_layout(motion_art_2)
show(fig_list_1[0])
C3.2 - Enhancement of noise complex
In order to isolate the noise artifact complex the following steps will be applied:# Signal rectification.
rect_art_signal = absolute(art_signal)
# Signal smoothing stage.
smooth_art_signal = bsnb.smooth(rect_art_signal, 2 * sr)
# Increase separation between variant and quasi-invariant components of the signal.
power_art_signal = power(smooth_art_signal, 2)
from numpy import power
fig_list_1 = bsnb.plot([time], [power_art_signal[:len(time)]], y_axis_label="Relative Electrodermal Response (a.u.)", legend_label=["Smooth Signal"], show_plot=False, get_fig_list=True, x_range=(time[0], time[-1]),
y_range=(-0.20*max(power_art_signal[:len(time)]), 1.2*max(power_art_signal[:len(time)])), hor_lines=[[average(power_art_signal[:len(time)])]])
# Highlight motion artifacts.
fill_color_1 = bsnb.opensignals_color_pallet()
motion_art_1 = BoxAnnotation(left=3.8, right=7, fill_color=fill_color_1, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_1, fill_alpha=0.1, legend_label="1st Motion Artifact")
fig_list_1[0].add_layout(motion_art_1)
fill_color_2 = bsnb.opensignals_color_pallet()
motion_art_2 = BoxAnnotation(left=13.5, right=16, fill_color=fill_color_2, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_2, fill_alpha=0.1, legend_label="2nd Motion Artifact")
fig_list_1[0].add_layout(motion_art_2)
show(fig_list_1[0])
C3.3 - Signal binarisation
All sample values greater than 1 (this value may not be appropriate for all BVP acquisitions) will be equal to 1 while the remaining samples will be updated to 0.noise_binary_complex = []
for i in range(0, len(power_art_signal)):
if power_art_signal[i] > 1:
noise_binary_complex.append(1)
else:
noise_binary_complex.append(0)
fig_list_1 = bsnb.plot([time, time], [art_signal[:len(time)] / max(art_signal[:len(time)]), noise_binary_complex[:len(time)]], y_axis_label="Normalized Electrodermal Response (a.u.)", legend=["Reconstructed Noise Artifact signal", "Binary Signal (Noise Complexes)"], show_plot=False, get_fig_list=True, x_range=(time[0], time[-1]),
y_range=(-1, 1.2))
# Highlight motion artifacts.
fill_color_1 = bsnb.opensignals_color_pallet()
motion_art_1 = BoxAnnotation(left=3.8, right=7, fill_color=fill_color_1, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_1, fill_alpha=0.1, legend_label="1st Motion Artifact")
fig_list_1[0].add_layout(motion_art_1)
fill_color_2 = bsnb.opensignals_color_pallet()
motion_art_2 = BoxAnnotation(left=13.5, right=16, fill_color=fill_color_2, fill_alpha=0.1)
fig_list_1[0].circle([-100], [0], fill_color=fill_color_2, fill_alpha=0.1, legend_label="2nd Motion Artifact")
fig_list_1[0].add_layout(motion_art_2)
show(fig_list_1[0])
C3.4 - Identification of the start and end of each noise artifact complex
Upward transitions define the start of noise complexes while downward transitions are the respective ends.diff_signal = diff(noise_binary_complex)
noise_complexes = [(None, None)]
for i in range(0, len(diff_signal)):
# Search for a starting period.
if diff_signal[i] == 1 and noise_complexes[-1][0] == None:
# Convert tuple to list.
tuple_list = list(noise_complexes[-1])
# Change entry reserved to the "end" value.
tuple_list[0] = i
# Update noise_complexes content.
noise_complexes[-1] = tuple(tuple_list)
# Search for an ending period.
if diff_signal[i] == -1 and noise_complexes[-1][1] == None and noise_complexes[-1][0] != None:
# Convert tuple to list.
tuple_list = list(noise_complexes[-1])
# Change entry reserved to the "end" value.
tuple_list[1] = i
# Update noise_complexes content.
noise_complexes[-1] = tuple(tuple_list)
# Preparing for the next searching iteration.
noise_complexes.append((None, None))
# Remove last tuple if it does not have an end.
if noise_complexes[-1][1] == None:
noise_complexes.pop()
print("Noise Complexes tuples:")
print(noise_complexes)
Noise Complexes tuples: [(3405, 7048), (13383, 16581)]
C4 - Removal of the remaining influence of noise complexes from the BVP signal
for i in range(0, len(noise_complexes)):
noise_period_start = noise_complexes[i][0]
noise_period_end = noise_complexes[i][1]
rec_signal[noise_period_start:noise_period_end] = average(rec_signal)
bsnb.plot([time[:noise_complexes[0][0]], time[noise_complexes[0][1]:noise_complexes[1][0]], time[noise_complexes[1][1]:]], [rec_signal[:noise_complexes[0][0]], rec_signal[noise_complexes[0][1]:noise_complexes[1][0]], rec_signal[noise_complexes[1][1]:noise_complexes[1][1] + len(time[noise_complexes[1][1]:])]], y_axis_label="Electrodermal Response (a.u.)", grid_plot=True, grid_columns=3, grid_lines=1)
D - Parameter Extraction
From the filtered BVP signals it is possible to execute a deep heart rate variability analysis, similar to the one presented in another Jupyter Notebook (ECG Analysis - Heart Rate Variability Parameters" ).A more specific analysis focused on the BVP morphology will be soon available, but, for now, we invite you to explore a very interesting research article from Mohamed Elgendi (source link ).
With the steps described on the current Jupyter Notebook you are able to acquire a broad perspective of the entire methodology involved on BVP signal analysis, starting in the acquisition stage and reaching the parameter extraction.
This is a mere illustrative example aiming to show some possible paths, but, you are free to explore this amazing world by your own.
We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don't stop your journey and learn more with the remaining Notebooks !
Auxiliary Code Segment (should not be replicated by the user)
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()
.................... CSS Style Applied to Jupyter Notebook .........................