# 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. 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 ## Load data This pipeline requires accelerometer and PPG data to run. Here, we start by loading a single contiguous time series (segment), for which we continue running steps 1-4. [Below](#multiple_segments_cell) 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 interally developed `TSDF` ([documentation](https://biomarkersparkinson.github.io/tsdf/)) to load and store data [[1](https://arxiv.org/abs/2211.11294)]. 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](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html)) - _.json_: `json.load()` ([documentation](https://docs.python.org/3/library/json.html#json.load)) ```python from pathlib import Path from paradigma.util import load_tsdf_dataframe # 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') 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}' ) 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}' ) # 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) display(df_ppg, df_acc) ```
time green
0 0.00000 262316
1 0.03340 262320
2 0.06680 262446
3 0.10020 262770
4 0.13360 262623
... ... ...
1029370 34339.49720 1049632
1029371 34339.53056 1049632
1029372 34339.56392 1049632
1029373 34339.59728 1049632
1029374 34339.63064 1020788

1029375 rows × 2 columns

time accelerometer_x accelerometer_y accelerometer_z
0 0.000000 -0.474641 -0.379426 0.770335
1 0.009933 -0.472727 -0.378947 0.765072
2 0.019867 -0.471770 -0.375598 0.766986
3 0.029800 -0.472727 -0.375598 0.770335
4 0.039733 -0.475120 -0.379426 0.772249
... ... ... ... ...
3455326 34339.561333 -0.257895 -0.319139 -0.761244
3455327 34339.571267 -0.555502 -0.153110 -0.671292
3455328 34339.581200 -0.286124 -0.263636 -0.981340
3455329 34339.591133 -0.232536 -0.161722 -0.832536
3455330 34339.601067 0.180383 -0.368421 -1.525837

3455331 rows × 4 columns

## Step 1: Preprocess data The first step after loading the data is preprocessing using the [preprocess_ppg_data](https://github.com/biomarkersParkinson/paradigma/blob/main/src/paradigma/preprocessing.py#:~:text=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. After this, 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). ```python from paradigma.config import PPGConfig, IMUConfig from paradigma.preprocessing import preprocess_ppg_data ppg_config = PPGConfig() imu_config = IMUConfig() print(f"Original data shapes:\n- PPG data: {df_ppg.shape}\n- Accelerometer data: {df_imu.shape}") 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 ) print(f"Overlapping preprocessed data shapes:\n- PPG data: {df_ppg_proc.shape}\n- Accelerometer data: {df_acc_proc.shape}") display(df_ppg_proc, df_acc_proc) ``` Original data shapes: - PPG data: (1029375, 2) - Accelerometer data: (3455331, 7) Overlapping preprocessed data shapes: - PPG data: (1030188, 2) - Accelerometer data: (3433961, 4)
time green
0 0.000000 -26.315811
1 0.033333 91.335299
2 0.066667 181.603416
3 0.100000 225.760466
4 0.133333 219.937282
... ... ...
1030183 34339.433333 224556.234611
1030184 34339.466667 210075.529517
1030185 34339.500000 163811.629247
1030186 34339.533333 94537.897763
1030187 34339.566667 12915.304284

1030188 rows × 2 columns

time accelerometer_x accelerometer_y accelerometer_z
0 0.00 -0.002324 -0.001442 -0.002116
1 0.01 -0.000390 -0.000914 -0.007396
2 0.02 0.000567 0.002474 -0.005445
3 0.03 -0.000425 0.002414 -0.002099
4 0.04 -0.002807 -0.001408 -0.000218
... ... ... ... ...
3433956 34339.56 -0.402941 0.038710 0.461449
3433957 34339.57 -0.659832 0.098696 0.817136
3433958 34339.58 -0.464138 0.033607 0.471552
3433959 34339.59 -0.389065 0.108485 0.622471
3433960 34339.60 -0.082625 -0.014490 0.119875

3433961 rows × 4 columns

## Step 2: Extract signal quality features The preprocessed data (PPG & accelerometer) 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 and from the accelerometer windows one relative power feature is calculated to assess periodic motion artifacts. The detailed steps are encapsulated in `extract_signal_quality_features` (documentation can be found [here](https://github.com/biomarkersParkinson/paradigma/blob/main/src/paradigma/pipelines/pulse_rate_pipeline.py#:~:text=extract_signal_quality_features)). ```python from paradigma.config import PulseRateConfig from paradigma.pipelines.pulse_rate_pipeline import extract_signal_quality_features ppg_config = PulseRateConfig('ppg') acc_config = PulseRateConfig('imu') print("The default window length for the signal quality feature extraction is set to", ppg_config.window_length_s, "seconds.") print("The default step size for the signal quality feature extraction is set to", ppg_config.window_step_length_s, "seconds.") df_features = extract_signal_quality_features( df_ppg=df_ppg_proc, df_acc=df_acc_proc, ppg_config=ppg_config, acc_config=acc_config, ) df_features ``` The default window length for the signal quality feature extraction is set to 6 seconds. The default step size for the signal quality feature extraction is set to 1 seconds.
time var mean median kurtosis skewness signal_to_noise auto_corr f_dom rel_power spectral_entropy acc_power_ratio
0 0.0 1.145652e+05 282.401234 238.829637 2.170853 0.107401 3.320049 0.544165 0.585938 0.138454 0.516336 0.026409
1 1.0 1.102401e+05 271.582177 236.891936 2.251393 -0.029309 3.041878 0.491829 0.585938 0.160433 0.511626 0.023402
2 2.0 1.061479e+05 262.348604 225.915756 2.415221 0.216631 2.818552 0.469092 0.585938 0.167007 0.525025 0.028592
3 3.0 9.514719e+04 245.089445 203.417715 2.481465 0.110420 2.677071 0.415071 0.585938 0.170626 0.550495 0.019296
4 4.0 7.393010e+04 218.379138 187.583266 2.405921 0.084566 2.796140 0.338369 0.585938 0.121113 0.595214 0.020083
... ... ... ... ... ... ... ... ... ... ... ... ...
34329 34329.0 8.176078e+06 1613.021494 438.201240 6.122772 -1.792336 1.378694 0.104389 0.351562 0.046616 0.356027 0.110219
34330 34330.0 3.512188e+07 3307.888927 1069.775894 8.160698 1.746472 1.442643 0.142226 0.351562 0.049424 0.371163 0.178742
34331 34331.0 1.181350e+08 6648.535487 2743.478312 5.654373 0.018587 1.558314 0.136803 0.351562 0.048211 0.366386 0.153351
34332 34332.0 1.252829e+09 20165.525309 6452.244225 6.805051 -1.222184 1.310088 0.666123 0.351562 0.037812 0.359105 0.154910
34333 34333.0 1.008217e+10 42271.328020 12552.656437 23.756877 4.167326 1.212179 0.044647 0.585938 0.113283 0.632749 0.093221

34334 rows × 12 columns

## 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). 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` (documentation can be found [here](https://github.com/biomarkersParkinson/paradigma/blob/main/src/paradigma/pipelines/pulse_rate_pipeline.py#:~:text=signal_quality_classification)). ```python from importlib.resources import files from paradigma.pipelines.pulse_rate_pipeline import signal_quality_classification ppg_quality_classifier_package_filename = 'ppg_quality_clf_package.pkl' full_path_to_classifier_package = files('paradigma') / 'assets' / ppg_quality_classifier_package_filename config = PulseRateConfig() df_sqa = signal_quality_classification( df=df_features, config=config, full_path_to_classifier_package=full_path_to_classifier_package ) df_sqa ```
time pred_sqa_proba pred_sqa_acc_label
0 0.0 1.121315e-02 1
1 1.0 7.126135e-03 1
2 2.0 7.017555e-03 1
3 3.0 4.134224e-03 1
4 4.0 9.195340e-04 1
... ... ... ...
34329 34329.0 1.782669e-08 0
34330 34330.0 2.078262e-06 0
34331 34331.0 1.190223e-07 0
34332 34332.0 1.383614e-08 0
34333 34333.0 7.587516e-07 1

34334 rows × 3 columns

#### Store as TSDF The predicted probabilities (and optionally other features) can be stored and loaded in TSDF as demonstrated below. ```python import tsdf from paradigma.util import write_df_data # 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) ``` ```python df_sqa, _, _ = load_tsdf_dataframe(path_to_prepared_data, prefix=f'segment{segment_nr}') df_sqa.head() ```
time pred_sqa_proba pred_sqa_acc_label
0 0.0 0.011213 1.0
1 1.0 0.007126 1.0
2 2.0 0.007018 1.0
3 3.0 0.004134 1.0
4 4.0 0.000920 1.0
## Step 4: Pulse rate estimation For pulse rate estimation, we extract segments of `config.tfd_length` using [estimate_pulse_rate](https://github.com/biomarkersParkinson/paradigma/blob/main/src/paradigma/pipelines/pulse_rate_pipeline.py#:~:text=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. ```python from paradigma.pipelines.pulse_rate_pipeline import estimate_pulse_rate print("The standard default minimal window length for the pulse rate extraction is set to", config.tfd_length, "seconds.") df_pr = estimate_pulse_rate( df_sqa=df_sqa, df_ppg_preprocessed=df_ppg_proc, config=config ) df_pr ``` The standard default minimal window length for the pulse rate extraction is set to 30 seconds.
time pulse_rate
0 47.0 80.372915
1 49.0 79.769382
2 51.0 79.136408
3 53.0 78.606477
4 55.0 77.870461
... ... ...
801 32876.0 78.220133
802 32878.0 78.047301
803 32880.0 78.047301
804 32882.0 78.238326
805 32884.0 78.556701

806 rows × 2 columns

### 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. ```python import pandas as pd from pathlib import Path from importlib.resources import files from paradigma.util import load_tsdf_dataframe from paradigma.config import PPGConfig, IMUConfig, PulseRateConfig from paradigma.preprocessing import preprocess_ppg_data from paradigma.pipelines.pulse_rate_pipeline import extract_signal_quality_features, signal_quality_classification, estimate_pulse_rate # Set the path to where the prepared data is saved path_to_prepared_data = Path('../../example_data') 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('ppg') acc_config = PulseRateConfig('imu') 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 config = PulseRateConfig() df_sqa = signal_quality_classification( df=df_features, config=config, full_path_to_classifier_package=full_path_to_classifier_package ) # 4: Estimate pulse rate df_pr = estimate_pulse_rate( df_sqa=df_sqa, df_ppg_preprocessed=df_ppg_proc, config=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_hr = 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](https://github.com/biomarkersParkinson/paradigma/blob/main/src/paradigma/pipelines/pulse_rate_pipeline.py#:~:text=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. ```python import pprint from paradigma.pipelines.pulse_rate_pipeline import aggregate_pulse_rate pr_values = df_pr['pulse_rate'].values df_pr_agg = aggregate_pulse_rate( pr_values=pr_values, aggregates = ['mode', '99p'] ) pprint.pprint(df_pr_agg) ``` {'metadata': {'nr_pr_est': 806}, 'pr_aggregates': {'99p_pulse_rate': 87.65865011636926, 'mode_pulse_rate': 81.25613346418058}}