ECG tutorial

import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from pprint import pprint

import physio

Detect ECG R peaks: quick way  —————————–

The fastest way is to use compute_ecg() using a predefined parameter preset. Here is a simple example.

raw_ecg = np.load('ecg1.npy')
srate = 1000.
times = np.arange(raw_ecg.size) / srate

ecg, ecg_peaks = physio.compute_ecg(raw_ecg, srate, parameter_preset='human_ecg')

r_peak_ind = ecg_peaks['peak_index'].values

fig, axs = plt.subplots(nrows=2, sharex=True)
ax = axs[0]
ax.plot(times, raw_ecg)
ax.scatter(times[r_peak_ind], raw_ecg[r_peak_ind], marker='o', color='magenta')
ax.set_ylabel('raw ecg')

ax = axs[1]
ax.plot(times, ecg)
ax.scatter(times[r_peak_ind], ecg[r_peak_ind], marker='o', color='magenta')
ax.set_ylabel('ecg')

ax.set_xlim(95, 125)

(95.0, 125.0)

Detect ECG R peaks: Parameters tuning  ————————————-

 Here is a simple recipe to change the default parameters.  We change here some filtering parameters (frequency band, filter type, order of the filter).

# get paramseters predefined set for 'human_ecg'
# this is a nested dict of parameter of every step
parameters = physio.get_ecg_parameters('human_ecg')
pprint(parameters)
# lets change on parameter in the structure
parameters['preprocess']['band'] = [2., 40.]
parameters['preprocess']['ftype'] = 'bessel'
parameters['preprocess']['order'] = 4
pprint(parameters)

ecg, ecg_peaks = physio.compute_ecg(raw_ecg, srate, parameters=parameters)

r_peak_ind = ecg_peaks['peak_index'].values

fig, ax = plt.subplots()
ax.plot(times, ecg)
ax.scatter(times[r_peak_ind], ecg[r_peak_ind], marker='o', color='magenta')
ax.set_ylabel('ecg')
ax.set_xlim(95, 125)

{'peak_clean': {'min_interval_ms': 400.0},
 'peak_detection': {'exclude_sweep_ms': 4.0, 'thresh': 'auto'},
 'preprocess': {'band': [5.0, 45.0],
                'ftype': 'bessel',
                'normalize': True,
                'order': 5}}
{'peak_clean': {'min_interval_ms': 400.0},
 'peak_detection': {'exclude_sweep_ms': 4.0, 'thresh': 'auto'},
 'preprocess': {'band': [2.0, 40.0],
                'ftype': 'bessel',
                'normalize': True,
                'order': 4}}

(95.0, 125.0)

ECG: compute time-domain heart rate variability (HRV) metrics  ——————–

metrics = physio.compute_ecg_metrics(ecg_peaks, min_interval_ms=500., max_interval_ms=2000.)
print(metrics)

HRV_Mean         734.442260
HRV_SD           145.928935
HRV_Median       678.000000
HRV_Mad          115.642973
HRV_CV             0.198694
HRV_MCV            0.170565
HRV_Asymmetry    -56.442260
HRV_RMSSD        103.591943
dtype: float64

ECG : compute instantaneous rate  ——————————–

The RR-interval (aka rri) time series is a common tool to analyse the heart rate variability (HRV).

 This is equivalent to compute the instantaneous heart rate.

Heart rate [bpm] = 1 / rri * 60

Most people use rri in ms, we feel that the use of heart rate in bpm is more intuitive. With bpm unit, an increase in the curve means heart rate acceleration. With ms unit, an increase in the curve means heart rate deceleration.

Feel free to use the units you prefer (bpm or ms)

new_times = times[::10]
instantaneous_rate = physio.compute_instantaneous_rate(
    ecg_peaks,
    new_times,
    limits=None,
    units='bpm',
    interpolation_kind='linear',
)
rri = physio.compute_instantaneous_rate(
    ecg_peaks,
    new_times,
    limits=None,
    units='ms',
    interpolation_kind='linear',
)

fig, axs = plt.subplots(nrows=2, sharex=True)
ax = axs[0]
ax.plot(new_times, instantaneous_rate)
ax.set_ylabel('heart rate [bpm]')
ax = axs[1]
ax.plot(new_times, rri)
ax.set_ylabel('rri [ms]')
ax.set_xlabel('time [s]')
ax.set_xlim(100, 150)

(100.0, 150.0)

ECG: compute frequency-domain heart rate variability (HRV) metrics  ————————-

freqency_bands = {'lf': (0.04, .15), 'hf' : (0.15, .4)}
psd_freqs, psd, psd_metrics = physio.compute_hrv_psd(
    ecg_peaks,
    sample_rate=100.,
    limits=None,
    units='bpm',
    freqency_bands=freqency_bands,
    window_s=250.,
    interpolation_kind='cubic',
)

print(psd_metrics)
fig, ax = plt.subplots()
# ax.semilogy(psd_freqs, psd)
ax.plot(psd_freqs, psd)
colors = {'lf': '#B8860B', 'hf' : '#D2691E'}
for name, freq_band in freqency_bands.items():
    ax.axvspan(*freq_band, alpha=0.1, color=colors[name], label=f'{name} : {psd_metrics[name]}')
ax.set_xlim(0, 0.6)
ax.set_xlabel('freq [Hz]')
ax.set_ylabel('HRV PSD')
ax.legend()

plt.show()

/home/docs/checkouts/readthedocs.org/user_builds/physio/envs/latest/lib/python3.11/site-packages/physio/ecg.py:282: UserWarning: The window is not optimal 250.0s compared to the lowest frequency 0.04Hz
  warnings.warn(f'The window is not optimal {window_s}s compared to the lowest frequency {min_freq}Hz')
/home/docs/checkouts/readthedocs.org/user_builds/physio/envs/latest/lib/python3.11/site-packages/physio/ecg.py:284: UserWarning: The duration is not optimal 299.236s compared to the lowest frequency 0.04Hz
  warnings.warn(f'The duration is not optimal {ecg_duration_s}s compared to the lowest frequency {min_freq}Hz')
lf    189.865574
hf     33.799386
dtype: float64