Note
Go to the end to download the full example code.
Getting started tutorial
Here is a quick overview of physio
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import physio
Reading data
physio aims to process physiological signals in NumPy format.
Therefore, the first step is to load your data into a NumPy array.
physio provides utility functions for reading certain formats, such as:
micromed
brainvision
For example, cardiac activity recorded by a BrainVision system with the channel name “ECG” can be read using
read_one_channel() by specifying:
the path to the file (“/path/to/the/file.vhdr”)
the format (“brainvision”)
the channel name (“ECG”).
For this tutorial, we will use internal files already stored in NumPy format for demonstration purposes. Let’s load them and plot a sample.
raw_resp = np.load('resp1_airflow.npy') # load respi
raw_ecg = np.load('ecg1.npy') # load ecg
srate = 1000. # our example signals have been recorded at 1000 Hz
times = np.arange(raw_resp.size) / srate # build time vector
linewidth = 0.8
fig, axs = plt.subplots(nrows=2, sharex=True)
ax = axs[0]
ax.plot(times, raw_resp, lw = linewidth)
ax.set_title('Raw respiratory airflow signal')
ax.set_ylabel('Amplitude (AU)')
ax.set_ylim(-1750, -1500)
ax = axs[1]
ax.plot(times, raw_ecg, lw = linewidth)
ax.set_title('Raw ECG (electrocardiogram)')
ax.set_ylabel('Amplitude (V)')
ax.set_ylim(0.012, 0.015)
ax.set_xlabel('Time (s)')
ax.set_xlim(185, 225)
(185.0, 225.0)
Analyze respiration
compute_respiration() is a high-level wrapper function that simplifies respiratory signal analysis.
To use this function, you must provide:
raw_resp : the raw respiratory signal as a NumPy array.
srate : the sampling rate of the respiratory signal
parameter_preset : a string specifying the type of respiratory data, which determines the set of parameters used for processing. Can be one of: human_airflow, human_co2, human_belt, rat_plethysmo, or rat_etisens_belt.
When called, compute_respiration() performs the following:
Preprocesses the respiratory signal (returns a NumPy array: resp)
Computes cycle-by-cycle features (returns a pd.DataFrame: resp_cycles)
Warning: The orientation of the raw_resp trace is important (multiply it by -1 for reversing it if necessary).
Inspiration must point downward for human_airflow or human_co2,
because downward deflections are interpreted by compute_respiration() as inspiration.
For rat_plethysmo or rat_etisens_belt, the inspiration–expiration transition must point upward.
resp, resp_cycles = physio.compute_respiration(raw_resp, srate, parameter_preset = 'human_airflow') # set 'human_airflow' as preset because example resp is an airflow from humans
linewidth = 0.9
fig, ax = plt.subplots()
ax.plot(times, raw_resp, lw = linewidth, label = 'raw resp signal')
ax.plot(times, resp, lw = linewidth, label = 'preprocessed resp signal')
ax.set_ylabel('Amplitude (AU)')
ax.set_ylim(-1750, -1500)
ax.set_xlim(185, 225)
ax.set_xlabel('Time (s)')
ax.set_title('From raw to preprocessed respiratory airflow signal')
ax.legend(loc = 'upper right')
<matplotlib.legend.Legend object at 0x70afab407690>
Respiration cycles and features
resp_cycles is a pd.DataFrame where each row corresponds to a respiratory cycle. Columns contain cycle-specific features such as duration, amplitude, volume, etc … See Respiration tutorial for a complete description of the columns.
print(resp_cycles.shape)
print(resp_cycles.columns)
columns = ['cycle_duration', 'inspi_volume', 'expi_volume', 'total_amplitude'] # selection of some columns / metrics
colors = ['tab:blue','darkorange','green','red'] # set colors
units = ['s','AU','AU','AU'] # set units
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(8, 6))
fig.suptitle(f'Distribution of some metrics from N = {resp_cycles.shape[0]} detected cycles')
for ax, col, color, unit in zip(axs.flatten(), columns, colors, units):
ax.hist(resp_cycles[col], bins=30, color=color, edgecolor = 'k')
ax.set_title(col)
ax.set_xlabel(f'{col} ({unit})')
ax.set_ylabel("Count")
fig.tight_layout()
resp_cycles
(30, 21)
Index(['inspi_index', 'expi_index', 'next_inspi_index', 'inspi_time',
'expi_time', 'next_inspi_time', 'cycle_duration', 'inspi_duration',
'expi_duration', 'cycle_freq', 'cycle_ratio', 'inspi_volume',
'expi_volume', 'total_amplitude', 'inspi_amplitude', 'expi_amplitude',
'inspi_peak_index', 'expi_peak_index', 'inspi_peak_time',
'expi_peak_time', 'total_volume'],
dtype='str')
All these metrics highly depend on the quality of detection of the start points of inspiration and expiration, which can be verified by visual inspection:
inspi_ind = resp_cycles['inspi_index'].values # get index of inspiration start points
expi_ind = resp_cycles['expi_index'].values # get index of expiration start points
fig, ax = plt.subplots()
ax.plot(times, resp)
ax.scatter(times[inspi_ind], resp[inspi_ind], color='green', label = 'inspiration start')
ax.scatter(times[expi_ind], resp[expi_ind], color='red', label = 'expiration start')
ax.set_xlim(185, 225)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude (AU)')
ax.set_ylim(-1750, -1500)
ax.set_title('Visual inspection of resp cycle detection')
ax.legend(loc = 'upper right')
<matplotlib.legend.Legend object at 0x70afaaf71790>
Analyze ECG
compute_ecg() is a high-level wrapper function that simplifies ECG signal analysis.
To use this function, you must provide:
raw_ecg : the raw ECG signal as a NumPy array
srate : the sampling rate of the ECG signal
parameter_preset : a string specifying the type of ECG data, which determines the set of parameters used for processing. Can be one of: human_ecg or rat_ecg.
When called, compute_ecg() performs the following:
Preprocesses the ECG signal (returns a NumPy array: ecg). By default, ecg is normalized.
Detects R peaks (returns a pd.DataFrame: ecg_peaks)
Warning: The orientation of the raw_ecg trace is important (multiply it by -1 for reversing it if necessary).
R peaks must point upward for the highest probability of detection by compute_ecg().
Sometimes the R and S peaks of the QRS complex are equivalent, so orientation does not matter,
because either the R or S peak will be detected and used to mark heartbeats.
ecg, ecg_peaks = physio.compute_ecg(raw_ecg, srate, parameter_preset = 'human_ecg') # set 'human_ecg' as preset because example ecg is from human
r_peak_ind = ecg_peaks['peak_index'].values # get index of R peaks
linewidth = 0.8
fig, axs = plt.subplots(nrows=2, sharex=True)
ax = axs[0]
ax.plot(times, raw_ecg, lw = linewidth)
ax.scatter(times[r_peak_ind], raw_ecg[r_peak_ind], color="#C3CB22", label = 'R peak')
ax.set_ylim(0.012, 0.015)
ax.set_ylabel('Amplitude (V)')
ax.set_title('R peak detection, plotted on raw ECG')
ax.legend(loc = 'upper right')
ax = axs[1]
ax.plot(times, ecg, lw = linewidth)
ax.scatter(times[r_peak_ind], ecg[r_peak_ind], color='#A8AE27', label = 'R peak')
ax.set_xlim(185, 225)
ax.set_ylim(-35, 50)
ax.set_ylabel('Amplitude (AU)')
ax.set_title('R peak detection, plotted on preprocessed ECG')
ax.legend(loc = 'upper right')
ax.set_xlabel('Time (s)')
Text(0.5, 23.52222222222222, 'Time (s)')
Heart Rate Variability metrics (time-domain)
compute_ecg_metrics() is a high-level wrapper function that computes time-domain Heart Rate Variability (HRV) metrics from previously detected R peaks (ecg_peaks).
To use this function, you must provide the previously detected R peaks, ecg_peaks:
ecg_peaks : pd.DataFrame, output of the function
compute_ecg()
We can visualize these metrics and the RR interval distribution. See ECG tutorial for a complete description of the metrics.
ecg_metrics = physio.compute_ecg_metrics(ecg_peaks) # ecg_metrics = a pd.Series containing HRV time domain results
peak_times = ecg_peaks['peak_time'].values # get R peak times of detection
rri_s = np.diff(peak_times) # compute RR intervals in seconds
rri_ms = rri_s * 1000 # seconds to milliseconds
fig, ax = plt.subplots()
ax.hist(rri_ms, bins=np.arange(500, 1400, 25), edgecolor = 'k') # plot distribution of RR intervals
ax.axvline(ecg_metrics['HRV_Mean'], color='orange', label='Mean RRi') # vertical line at the mean RRi
ax.axvline(ecg_metrics['HRV_Median'], color='violet', label='Median RRi') # vertical line at the median RRi
ax.axvspan(ecg_metrics['HRV_Median'] - 2*ecg_metrics['HRV_Mad'], ecg_metrics['HRV_Median'] + 2*ecg_metrics['HRV_Mad'], color = 'orange', alpha = 0.1, label = 'Median +/- 2 * Median Absolute Deviation') # vertical span at the med + 2*mad
ax.axvspan(ecg_metrics['HRV_Mean'] - 2*ecg_metrics['HRV_SD'], ecg_metrics['HRV_Mean'] + 2*ecg_metrics['HRV_SD'], color = 'violet', alpha = 0.1, label = 'Mean +/- 2 * Standard-Deviation') # vertical span at the mean + 2*sd
ax.set_xlabel('RR interval (ms)')
ax.set_ylabel('Count')
ax.set_ylim(-5, 90)
ax.legend(loc = 'upper center', ncols = 2)
ax.set_title('Distribution or RR intervals + Heart Rate Variability metrics')
print(ecg_metrics)
HRV_Mean 734.444717
HRV_SD 145.936125
HRV_Median 678.000000
HRV_Mad 115.642973
HRV_CV 0.198703
HRV_MCV 0.170565
HRV_Asymmetry -56.444717
HRV_RMSSD 103.553299
dtype: float64
Heart Rate Variability metrics (frequency-domain)
compute_hrv_psd() is a high-level wrapper function that computes frequency-domain Heart Rate Variability (HRV) metrics from previously detected R peaks (ecg_peaks).
To use this function, you must provide the previously detected R peaks, ecg_peaks:
ecg_peaks : pd.DataFrame, output of the function
compute_ecg()
Many others parameters can be set for this function (frequency bands, size of the Welch’s window, etc…). See ECG tutorial for more information about these.
When called, compute_hrv_psd() performs the following:
Compute an instantaneous heart rate vector from ecg peaks and compute a Fourier transform using Welch method (returns two NumPy arrays: psd_freqs and psd giving frequency vector and power vector, respectively).
Compute HRV frequency metrics by getting power for each frequency band using trapezoïdal rule (returns a pd.Series: psd_metrics)
We can visualize these metrics and the RR interval distribution. See ECG tutorial for a complete description of the metrics.
frequency_bands = {'lf': (0.04, .15), 'hf' : (0.15, .4)} # set classical cutoffs of low and high frequency bands, in a dictionnary
psd_freqs, psd, psd_metrics = physio.compute_hrv_psd(ecg_peaks=ecg_peaks, frequency_bands=frequency_bands)
print(psd_metrics)
fig, ax = plt.subplots()
ax.plot(psd_freqs, psd)
colors = {'lf': '#B8860B', 'hf' : '#D2691E'}
for name, freq_band in frequency_bands.items():
ax.axvspan(*freq_band, alpha=0.1, color=colors[name], label=f'{name} : {round(psd_metrics[name], 2)}') # plot one vertical span for each frequency band
ax.set_xlim(0, 0.6)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power Spectral Density')
ax.legend(loc = 'upper right')
ax.set_title('HRV Power Spectral Density\n-> Frequency-Domain metrics')
lf 189.890278
hf 33.775978
dtype: float64
Text(0.5, 1.0, 'HRV Power Spectral Density\n-> Frequency-Domain metrics')
Cyclic deformation
deform_traces_to_cycle_template() is a tool used to deform traces to a cycle template.
It works by stretching a trace using linear resampling based on a fixed number of points per cycle.
This is helpful to explore if features of a signal are driven by a cyclic phenomenon like respiration.
To use this function, you must provide:
data : nd.array. Axis of the time must always be 0, meaning of shape (n_times,…).
times : np.array. Time vector of the data. Shape = (n_times,)
cycle_times : nd.array with shape (n_cycles, n_segments + 1). Typically, for respiration, cycle_times is an array with 3 columns (inspi_time + expi_time + next_inspi_time) that will make deformation with 2 segments. If cycle_times is 1D, then it is converted to shape (size-1, 2). The end of every cycles must match the start of the next cycle.
points_per_cycle : Number of phase points per cycle
segment_ratios : None or float or list of float. None if 1 segment. Float or list of float if 2 segments. List of floats if > 2 segments. This is a ratio between 0 and 1 where cycle is divided.
output_mode : ‘stacked’ / ‘unstacked’ / ‘unstacked_full’. Format of the outputs. Stacked -> 2D matrix : cycles / points per cycle. Unstacked -> 1D matrix : flattened version of the stacked. Unstacked_full returns extra-outputs.
# Here, we deform the respiratory trace by "itself": the respiratory cycle.
# This leads to an average respiratory template. Importantly, this can be done using one or several segments inside the cycle.
# here we have 3 times per cycle so 2 segments :
# segment 1: inspiration to expiration
# segment 2: expiration to next inspiration
cycle_times = resp_cycles[['inspi_time', 'expi_time', 'next_inspi_time']].values # get 3 timepoints per cycle, dividing each cycle in 2 segments
segment_ratios = 0.4 # ratio between 0 and 1 where cycle is divided into 2 segments
deformed_resp_2seg = physio.deform_traces_to_cycle_template(resp, # resp trace to deform
times, # time vector of resp trace to deform
cycle_times, # times of resp cycles, used to strech
points_per_cycle=40, # number of points per cycle used for linear resampling
segment_ratios=segment_ratios, # ratio between 0 and 1 where cycle is divided into 2 segments
output_mode='stacked' # choose a stacked version of the returned matrix (see docstrings)
)
print(deformed_resp_2seg.shape, cycle_times.shape)
# here we have 2 times per cycle so 1 segment:
# segment 1: inspiration to next inspiration
cycle_times = resp_cycles[['inspi_time', 'next_inspi_time']].values # get 2 timepoints per cycle, dividing each cycle in 1 segment (so not dividing)
deformed_resp_1seg = physio.deform_traces_to_cycle_template(resp, # resp trace to deform
times, # time vector of resp trace to deform
cycle_times, # times of resp cycles, used to strech
points_per_cycle=40, # number of points per cycle used for linear resampling
segment_ratios=None, # ratio between 0 and 1 where cycle is divided into 2 segments. None in this case because 1 segment
output_mode='stacked' # choose a stacked version of the returned matrix (see docstrings)
)
print(deformed_resp_1seg.shape, cycle_times.shape)
fig, axs = plt.subplots(ncols=2, sharey=True, figsize = (8,4))
ax = axs[0]
physio.plot_cyclic_deformation(deformed_resp_2seg, segment_ratios=0.4, two_cycles=False, ax=ax) # using function physio.plot_cyclic_deformation to plot
ax.set_title('Deformation: 2 segments')
ax.set_ylabel('Resp Amplitude (AU)')
ax = axs[1]
physio.plot_cyclic_deformation(deformed_resp_1seg, segment_ratios=None, two_cycles=False, ax=ax) # using function physio.plot_cyclic_deformation to plot
ax.set_title('Deformation: 1 segment')
for ax in axs:
ax.set_xlabel('Points per cycle (= respiratory phase)')
(30, 40) (30, 3)
(30, 40) (30, 2)
Cyclic deformation on ECG
Let’s use the same deform_traces_to_cycle_template() tool for ECG trace.
We can use a simple time vector: the R peak times. In this case, the ECG trace is converted as a 1-segment.
cycle_times = ecg_peaks['peak_time'].values # one cycle = R to R -> one segment by R to R period
deformed_ecg = physio.deform_traces_to_cycle_template(ecg, # ecg trace to deform
times, # time vector of ecg trace to deform
cycle_times, # times of ecg cycles = R to R periods, used to strech
points_per_cycle=300, # number of points per cycle used for linear resampling (more than for resp because ecg trace has faster dynamics)
segment_ratios=None, # 1 segment so None
output_mode='stacked'
)
print(deformed_ecg.shape, cycle_times.shape)
fig, ax = plt.subplots()
physio.plot_cyclic_deformation(deformed_ecg, two_cycles=True, ax=ax) # set two_cycles = True to see a concatenation of two cyclical deformations, so that we can see the whole PQRST complex
ax.set_title('Two ECG cyclical deformations concatenated:\nnormalized view of the whole PQRST complex')
ax.set_ylabel('Amplitude (AU)')
ax.set_xlabel('Points per cycle * 2')
(407, 300) (408,)
Text(0.5, 23.52222222222222, 'Points per cycle * 2')
Total running time of the script: (0 minutes 1.306 seconds)