Pulse rate analysis

This tutorial shows how to extract pulse rate estimates using photoplethysmography (PPG) data and accelerometer data. The pipeline consists of a stepwise approach to determine signal quality, assessing both PPG morphology and accounting for periodic artifacts using the accelerometer. The usage of accelerometer is optional but is recommended to specifically account for periodic motion artifacts. Based on the signal quality, we extract high-quality segments and estimate the pulse rate for every 2 s using the smoothed pseudo Wigner-Ville Distribution.

In this tutorial, we use two days of data from a participant of the Personalized Parkinson Project to demonstrate the functionalities. Since ParaDigMa expects contiguous time series, the collected data was stored in two segments each with contiguous timestamps. Per segment, we load the data and perform the following steps:

1. Preprocess the time series data

2. Extract signal quality features

3. Signal quality classification

4. Pulse rate estimation

We then combine the output of the different segments for the final step:

5. Pulse rate aggregation

Import required modules

import json
from importlib.resources import files
from pathlib import Path

import pandas as pd
import tsdf

from paradigma.classification import ClassifierPackage
from paradigma.config import IMUConfig, PPGConfig, PulseRateConfig
from paradigma.constants import DataColumns
from paradigma.pipelines.pulse_rate_pipeline import (
    aggregate_pulse_rate,
    estimate_pulse_rate,
    extract_signal_quality_features,
    signal_quality_classification,
)
from paradigma.preprocessing import preprocess_ppg_data
from paradigma.util import load_tsdf_dataframe, write_df_data

Load data

This pipeline requires PPG data and can be enhanced with accelerometer data (optional). Here, we start by loading a single contiguous time series (segment), for which we continue running steps 1-4. Below we show how to run these steps for multiple segments. The channel green represents the values obtained with PPG using green light.

In this example we use the internally developed TSDF (documentation) to load and store data [1]. However, we are aware that there are other common data formats. For example, the following functions can be used depending on the file extension of the data:

• .csv: pandas.read_csv() (documentation)

• .json: json.load() (documentation)

# Set the path to where the prepared data is saved and load the data.
# Note: the test data is stored in TSDF, but you can load your data in your own way
path_to_prepared_data =  Path('../../example_data/verily')

ppg_prefix = 'ppg'
imu_prefix = 'imu'

segment_nr = '0001'

df_ppg, metadata_time_ppg, metadata_values_ppg = load_tsdf_dataframe(
    path_to_data=path_to_prepared_data / ppg_prefix,
    prefix=f'PPG_segment{segment_nr}'
)

# Only relevant if you have IMU data available
df_imu, metadata_time_imu, metadata_values_imu = load_tsdf_dataframe(
    path_to_data=path_to_prepared_data / imu_prefix,
    prefix=f'IMU_segment{segment_nr}'
)

time_col = ['time']
acc_cols = ['accelerometer_x', 'accelerometer_y', 'accelerometer_z']
df_acc = df_imu[time_col + acc_cols]

display(df_ppg, df_acc)

Step 1: Preprocess data

The first step after loading the data is preprocessing using the preprocess_ppg_data. This begins by isolating segments containing both PPG and IMU data, discarding portions where one modality (e.g., PPG) extends beyond the other, such as when the PPG recording is longer than the accelerometer data. This functionality requires the starting times (metadata_time_ppg.start_iso8601 and metadata_time_imu.start_iso8601) in iso8601 format as inputs. After this step, the preprocess_ppg_data function resamples the PPG and accelerometer data to uniformly distributed timestamps, addressing the fixed but non-uniform sampling rates of the sensors. If the difference between timestamps is larger than a specified tolerance (config.tolerance, in seconds), it will return an error that the timestamps are not contiguous. If you still want to process the data in this case, you can create segments from discontiguous samples using the function create_segments and analyze these segments consecutively as shown in here. After resampling, a bandpass Butterworth filter (4th-order, bandpass frequencies: 0.4–3.5 Hz) is applied to the PPG signal, while a high-pass Butterworth filter (4th-order, cut-off frequency: 0.2 Hz) is applied to the accelerometer data.

Note: the printed shapes are (rows, columns) with each row corresponding to a single data point and each column representing a data column (e.g.time). The number of rows of the overlapping segments of PPG and accelerometer are not the same due to sampling differences (other sensors and possibly other sampling frequencies).

# Set column names: replace DataColumn.* with your actual column names.
# It is only necessary to set the columns that are present in your data, and
# only if they differ from the default names defined in DataColumns.
imu_column_mapping = {
    'TIME': DataColumns.TIME,
    'ACCELEROMETER_X': DataColumns.ACCELEROMETER_X,
    'ACCELEROMETER_Y': DataColumns.ACCELEROMETER_Y,
    'ACCELEROMETER_Z': DataColumns.ACCELEROMETER_Z,
}

ppg_column_mapping = {
    'TIME': DataColumns.TIME,
    'PPG': DataColumns.PPG,
}

ppg_config = PPGConfig(column_mapping=ppg_column_mapping)
imu_config = IMUConfig(column_mapping=imu_column_mapping)

print(
    f"The tolerance for checking contiguous timestamps is set to "
    f"{ppg_config.tolerance:.3f} seconds for PPG data and "
    f"{imu_config.tolerance:.3f} seconds for accelerometer data."
)
print(
    f"Original data shapes:\n- PPG data: {df_ppg.shape}\n"
    f"- Accelerometer data: {df_imu.shape}"
)
df_ppg_proc, df_acc_proc = preprocess_ppg_data(
    df_ppg=df_ppg,
    ppg_config=ppg_config,
    start_time_ppg=metadata_time_ppg.start_iso8601,  # Optional
    df_acc=df_acc,  # Optional
    imu_config=imu_config,  # Optional
    start_time_imu=metadata_time_imu.start_iso8601  # Optional
)

print(f"Overlapping preprocessed data shapes:\n- PPG data: {df_ppg_proc.shape}\n"
      f"- Accelerometer data: {df_acc_proc.shape}")
display(df_ppg_proc, df_acc_proc)

Step 2: Extract signal quality features

The preprocessed data is windowed into overlapping windows of length ppg_config.window_length_s with a window step of ppg_config.window_step_length_s. From the PPG windows, 10 time- and frequency domain features are extracted to assess PPG morphology. In case of using the accelerometer data (optional), one relative power feature is calculated per window to assess periodic motion artifacts.

The detailed steps are encapsulated in extract_signal_quality_features.

pulse_rate_ppg_config = PulseRateConfig(
    sensor='ppg',
    ppg_sampling_frequency=ppg_config.sampling_frequency,
)

pulse_rate_acc_config = PulseRateConfig(
    sensor='imu',
    imu_sampling_frequency=imu_config.sampling_frequency,
    accelerometer_colnames=imu_config.accelerometer_colnames,
)

print("The default window length for the signal quality feature extraction is "
      f"set to {pulse_rate_ppg_config.window_length_s} seconds.")
print("The default step size for the signal quality feature extraction is "
      f"set to {pulse_rate_ppg_config.window_step_length_s} seconds.")

# Remove optional arguments if you don't have accelerometer data
# (set to None or remove arguments)
df_features = extract_signal_quality_features(
    df_ppg=df_ppg_proc,
    df_acc=df_acc_proc,
    ppg_config=pulse_rate_ppg_config,
    acc_config=pulse_rate_acc_config,
)

display(df_features)

Step 3: Signal quality classification

A trained logistic classifier is used to predict PPG signal quality and returns the pred_sqa_proba, which is the posterior probability of a PPG window to look like the typical PPG morphology (higher probability indicates toward the typical PPG morphology).
If accelerometer is used, the relative power feature from the accelerometer is compared to a threshold for periodic artifacts and therefore pred_sqa_acc_label is used to return a label indicating predicted periodic motion artifacts (label 0) or no periodic motion artifacts (label 1).

The classification step is implemented in signal_quality_classification.

Note on scale sensitivity
The PPG sensor used for developing this pipeline records in arbitrary units. Some features are scale sensitive and require rescaling when applying the pipeline to other datasets or PPG sensors.
In this pipeline, the logistic classifier for PPG morphology was trained on z-scored features, using the means (μ) and standard deviations (σ) from the Personalized Parkinson Project training set (documentation of the training set can be found in the Signal Quality Assessment section here). These μ and σ values are stored in the ppg_quality_classifier_package. When applying the code to another dataset, users are advised to recalculate μ and σ for each feature on their (training) data and update the classifier package accordingly. Example code to recalculate the μ and σ can be found in the code cell below and in section 7.3.1. here.

config = PulseRateConfig()

ppg_quality_classifier_package_filename = 'ppg_quality_clf_package.pkl'
full_path_to_classifier_package = (
    files('paradigma')
    / 'assets'
    / ppg_quality_classifier_package_filename
)

# Load the classifier package
clf_package = ClassifierPackage.load(full_path_to_classifier_package)

# If you use a different sensor or dataset, you have to update the classifier
# package with the mu and sigma calculated on your training data.
# Example code to recalculate the mean and std to update the clf_package:
# import numpy as np

# # create random x_train for demonstration purposes, 100 samples and 10 features
# x_train = np.random.rand(100, 10)
# clf_package.update_scaler(x_train)

df_sqa = signal_quality_classification(
    df=df_features,
    config=config,
    clf_package=clf_package
)

df_sqa

Store as TSDF

The predicted probabilities (and optionally other features) can be stored and loaded in TSDF as demonstrated below.

# Set 'path_to_data' to the directory where you want to save the data
metadata_time_store = tsdf.TSDFMetadata(
    metadata_time_ppg.get_plain_tsdf_dict_copy(),
    path_to_prepared_data
)
metadata_values_store = tsdf.TSDFMetadata(
    metadata_values_ppg.get_plain_tsdf_dict_copy(),
    path_to_prepared_data
)

# Select the columns to be saved
metadata_time_store.channels = ['time']
metadata_values_store.channels = ['pred_sqa_proba', 'pred_sqa_acc_label']

# Set the units
metadata_time_store.units = ['Relative seconds']
metadata_values_store.units = ['Unitless', 'Unitless']
metadata_time_store.data_type = float
metadata_values_store.data_type = float

# Set the filenames
meta_store_filename = f'segment{segment_nr}_meta.json'
values_store_filename = meta_store_filename.replace('_meta.json', '_values.bin')
time_store_filename = meta_store_filename.replace('_meta.json', '_time.bin')

metadata_values_store.file_name = values_store_filename
metadata_time_store.file_name = time_store_filename

write_df_data(
    metadata_time_store,
    metadata_values_store,
    path_to_prepared_data,
    meta_store_filename,
    df_sqa
)

df_sqa, _, _ = load_tsdf_dataframe(
    path_to_prepared_data,
    prefix=f'segment{segment_nr}'
)
df_sqa.head()

Step 4: Pulse rate estimation

For pulse rate estimation, we extract segments of config.tfd_length using estimate_pulse_rate. We calculate the smoothed-pseudo Wigner-Ville Distribution (SPWVD) to obtain the frequency content of the PPG signal over time. We extract for every timestamp in the SPWVD the frequency with the highest power. For every non-overlapping 2 s window we average the corresponding frequencies to obtain a pulse rate per window.

Note: for the test data we set the tfd_length to 10 s instead of the default of 30 s, because the small PPP test data doesn’t have 30 s of consecutive high-quality PPG data.

print(
    "The standard default minimal window length for the pulse rate "
    "extraction is set to", pulse_rate_ppg_config.tfd_length, "seconds."
)

df_pr = estimate_pulse_rate(
    df_sqa=df_sqa,
    df_ppg_preprocessed=df_ppg_proc,
    config=pulse_rate_ppg_config
)

df_pr

Run steps 1 - 4 for multiple segments

If your data is also stored in multiple segments, you can modify segments in the cell below to a list of the filenames of your respective segmented data.

# Set the path to where the prepared data is saved
path_to_prepared_data =  Path('../../example_data/verily')

ppg_prefix = 'ppg'
imu_prefix = 'imu'

# Set the path to the classifier package
ppg_quality_classifier_package_filename = 'ppg_quality_clf_package.pkl'
full_path_to_classifier_package = (
    files('paradigma')
    / 'assets'
    / ppg_quality_classifier_package_filename
)

# Create a list of dataframes to store the estimated pulse rates of all segments
list_df_pr = []

segments = ['0001', '0002'] # list with all available segments

for segment_nr in segments:

    # Load the data
    df_ppg, metadata_time_ppg, _ = load_tsdf_dataframe(
        path_to_data=path_to_prepared_data / ppg_prefix,
        prefix=f'PPG_segment{segment_nr}'
    )
    df_imu, metadata_time_imu, _ = load_tsdf_dataframe(
        path_to_data=path_to_prepared_data / imu_prefix,
        prefix=f'IMU_segment{segment_nr}'
    )

    # Drop the gyroscope columns from the IMU data
    cols_to_drop = df_imu.filter(regex='^gyroscope_').columns
    df_acc = df_imu.drop(cols_to_drop, axis=1)

    # 1: Preprocess the data

    ppg_config = PPGConfig()
    imu_config = IMUConfig()

    df_ppg_proc, df_acc_proc = preprocess_ppg_data(
        df_ppg=df_ppg,
        df_acc=df_acc,
        ppg_config=ppg_config,
        imu_config=imu_config,
        start_time_ppg=metadata_time_ppg.start_iso8601,
        start_time_imu=metadata_time_imu.start_iso8601
    )

    # 2: Extract signal quality features
    ppg_config = PulseRateConfig(
        sensor='ppg',
        ppg_sampling_frequency=ppg_config.sampling_frequency
    )
    acc_config = PulseRateConfig(
        sensor='imu',
        imu_sampling_frequency=imu_config.sampling_frequency,
        accelerometer_colnames=imu_config.accelerometer_colnames,
    )

    df_features = extract_signal_quality_features(
        df_ppg=df_ppg_proc,
        df_acc=df_acc_proc,
        ppg_config=ppg_config,
        acc_config=acc_config,
    )

    # 3: Signal quality classification
    df_sqa = signal_quality_classification(
        df=df_features,
        config=ppg_config,
        clf_package=clf_package
    )

    # 4: Estimate pulse rate
    df_pr = estimate_pulse_rate(
        df_sqa=df_sqa,
        df_ppg_preprocessed=df_ppg_proc,
        config=ppg_config
    )

    # Add the hr estimations of the current segment to the list
    df_pr['segment_nr'] = segment_nr
    list_df_pr.append(df_pr)

df_pr = pd.concat(list_df_pr, ignore_index=True)

Step 5: Pulse rate aggregation

The final step is to aggregate all 2 s pulse rate estimates using aggregate_pulse_rate. In the current example, the mode and 99th percentile are calculated. We hypothesize that the mode gives representation of the resting pulse rate while the 99th percentile indicates the maximum pulse rate. In Parkinson’s disease, we expect that these two measures could reflect autonomic (dys)functioning. The nr_pr_est in the metadata indicates based on how many 2 s windows these aggregates are determined.

pr_values = df_pr['pulse_rate'].values
df_pr_agg = aggregate_pulse_rate(
    pr_values=pr_values,
    aggregates = ['mode', '99p']
)

print(json.dumps(df_pr_agg, indent=2))