From 60e8cfeaab6f90fc27fcf1cd5a17910c2cccf7a7 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 6 Jan 2021 20:26:39 -0500 Subject: [PATCH 1/4] Initial CPM translation from MATLAB. --- phys2denoise/metrics/cardiac.py | 116 ++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 phys2denoise/metrics/cardiac.py diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py new file mode 100644 index 0000000..4addd89 --- /dev/null +++ b/phys2denoise/metrics/cardiac.py @@ -0,0 +1,116 @@ +"""Denoising metrics for cardio recordings.""" +import numpy as np +from scipy import signal + + +def cpm(Fs, PPGlocs, HR, cardiac): + """Calculate CPM regressors.""" + Ts = 1 / Fs + HPF_f = 0.008 + shift_CPM = 0.5 + shift_CPM_VA = 0.5 + shift_RETR = 0 + M_order = 8 + filt_b, filt_a = signal.butter(2, HPF_f * 2 * Ts, "high") + + T_del = 5 # delete first 5 seconds from output + nDel = int(T_del * Fs) + + HRmean = np.mean(HR) + + voxel = cardiac[nDel:-nDel] + NV = len(voxel) + voxel = signal.filtfilt(filt_b, filt_a, voxel) + N = len(cardiac) + time = np.arange(0, (N - 1) * Ts, Ts) + + memory = 60 / HRmean + CPM_IR = func_CPM_cos(Ts, memory, M_order) + + u = np.zeros(time.shape) + uA = u.copy() + nPeaks = len(PPGlocs) + for i in range(nPeaks): + t = PPGlocs[i] + val, loc = np.min(np.abs(time - t)) + u[loc] = 1 + uA[loc] = cardiac[loc] + + CPM_regr_all = np.zeros((N, M_order)) + CPM_Amp_regr_all = np.zeros((N, M_order)) + for m in range(M_order * 2): + u_conv = np.convolve(u, CPM_IR[:, m]) + u_conv = u_conv[:N] + CPM_regr_all[:, m] = u_conv[:] + x = np.convolve(uA, CPM_IR[:, m]) + x = x[:N] + CPM_Amp_regr_all[:, m] = x + + CPM_regr_all = signal.filtfilt(filt_b, filt_a, CPM_regr_all) + CPM_Amp_regr_all = signal.filtfilt(filt_b, filt_a, CPM_Amp_regr_all) + + RETR_regr_all = RETR_Card_regressors_v2(time, PPGlocs, M_order) + RETR_regr_all = signal.filtfilt(filt_b, filt_a, RETR_regr_all) + + ind = np.arange(NV) + ind = np.round(ind + nDel + (shift_CPM * Fs)) + CPM_regr = CPM_regr_all[ind, :] + + ind = np.arange(NV) + ind = np.round(ind + nDel + (shift_CPM_VA * Fs)) + CPM_Amp_regr = CPM_Amp_regr_all[ind, :] + + ind = np.arange(NV) + ind = np.round(ind + nDel + (shift_RETR * Fs)) + RETR_regr = RETR_regr_all[ind, :] + + regr_CPM = np.hstack(CPM_regr, np.ones((NV, 1))) + regr_CPM_Amp = np.hstack(CPM_Amp_regr, np.ones((NV, 1))) + RETR_regr = np.hstack(RETR_regr, np.ones((NV, 1))) + return regr_CPM, regr_CPM_Amp, RETR_regr + + +def RETR_Card_regressors_v2(time, locsECG, M): + """Calculate RETROICOR cardiac regressors.""" + NV = len(time) + Phi = np.zeros((NV, 1)) + for i in range(NV): + t = time[i] + _, minI = np.min(np.abs(locsECG - t)) + + minOnLeft = t - locsECG[minI] > 0 + if (minI == 1) and not minOnLeft: + t2 = locsECG[minI] + t1 = t2 - 1 + elif (minI == len(locsECG)) and minOnLeft: + t1 = locsECG[minI] + t2 = t1 + 1 + elif minOnLeft: + t1 = locsECG[minI] + t2 = locsECG[minI + 1] + else: + t1 = locsECG[minI - 1] + t2 = locsECG[minI] + + Phi[i] = 2 * np.pi * (t - t1) / (t2 - t1) + + Regr = np.zeros((NV, M * 2)) + for i in range(M): + Regr[:, ((i - 1) * 2) + 1] = np.cos(i * Phi) + Regr[:, i * 2] = np.sin(i * Phi) + + # Regr=[zeros(1,M*2);diff(Regr)]; + return Regr + + +def func_CPM_cos(Ts, memory, M): + """Calculate CPM.""" + t_win = np.arange(0, memory, Ts) + nIR = len(t_win) + + IR_all = np.zeros((nIR, M * 2)) + for m in range(M): + IR_all[:, ((m - 1) * 2) + 1] = np.cos(m * 2 * np.pi * t_win / memory) - 1 + IR_all[:, ((m - 1) * 2) + 2] = np.sin(m * 2 * np.pi * t_win / memory) + + return IR_all From 466be3ba69500936a3393cf610c5fb0153cede57 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 7 Jan 2021 14:38:02 -0500 Subject: [PATCH 2/4] More work. --- phys2denoise/metrics/__init__.py | 4 + phys2denoise/metrics/cardiac.py | 233 ++++++++++++++++++++----------- 2 files changed, 156 insertions(+), 81 deletions(-) diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py index e69de29..c63a621 100644 --- a/phys2denoise/metrics/__init__.py +++ b/phys2denoise/metrics/__init__.py @@ -0,0 +1,4 @@ +"""Denoising metrics from cardiac recordings.""" +from .cardiac import cpm + +__all__ = ["cpm"] diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index 4addd89..fff4b37 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -3,86 +3,143 @@ from scipy import signal -def cpm(Fs, PPGlocs, HR, cardiac): - """Calculate CPM regressors.""" - Ts = 1 / Fs - HPF_f = 0.008 - shift_CPM = 0.5 - shift_CPM_VA = 0.5 - shift_RETR = 0 - M_order = 8 - filt_b, filt_a = signal.butter(2, HPF_f * 2 * Ts, "high") - - T_del = 5 # delete first 5 seconds from output - nDel = int(T_del * Fs) - - HRmean = np.mean(HR) - - voxel = cardiac[nDel:-nDel] - NV = len(voxel) - voxel = signal.filtfilt(filt_b, filt_a, voxel) - N = len(cardiac) - time = np.arange(0, (N - 1) * Ts, Ts) - - memory = 60 / HRmean - CPM_IR = func_CPM_cos(Ts, memory, M_order) - - u = np.zeros(time.shape) - uA = u.copy() - nPeaks = len(PPGlocs) - for i in range(nPeaks): - t = PPGlocs[i] - val, loc = np.min(np.abs(time - t)) - u[loc] = 1 - uA[loc] = cardiac[loc] - - CPM_regr_all = np.zeros((N, M_order)) - CPM_Amp_regr_all = np.zeros((N, M_order)) - for m in range(M_order * 2): - u_conv = np.convolve(u, CPM_IR[:, m]) - u_conv = u_conv[:N] - CPM_regr_all[:, m] = u_conv[:] - x = np.convolve(uA, CPM_IR[:, m]) - x = x[:N] - CPM_Amp_regr_all[:, m] = x - - CPM_regr_all = signal.filtfilt(filt_b, filt_a, CPM_regr_all) - CPM_Amp_regr_all = signal.filtfilt(filt_b, filt_a, CPM_Amp_regr_all) - - RETR_regr_all = RETR_Card_regressors_v2(time, PPGlocs, M_order) - RETR_regr_all = signal.filtfilt(filt_b, filt_a, RETR_regr_all) - - ind = np.arange(NV) - ind = np.round(ind + nDel + (shift_CPM * Fs)) - CPM_regr = CPM_regr_all[ind, :] - - ind = np.arange(NV) - ind = np.round(ind + nDel + (shift_CPM_VA * Fs)) - CPM_Amp_regr = CPM_Amp_regr_all[ind, :] - - ind = np.arange(NV) - ind = np.round(ind + nDel + (shift_RETR * Fs)) - RETR_regr = RETR_regr_all[ind, :] - - regr_CPM = np.hstack(CPM_regr, np.ones((NV, 1))) - regr_CPM_Amp = np.hstack(CPM_Amp_regr, np.ones((NV, 1))) - RETR_regr = np.hstack(RETR_regr, np.ones((NV, 1))) - return regr_CPM, regr_CPM_Amp, RETR_regr +def cpm(samplerate, PPGlocs, HR, cardiac): + """Calculate cardiac pulsatility model (CPM) regressors. + + Parameters + ---------- + samplerate : float + Sample rate in Hertz. + PPGlocs : 1D numpy.ndarray + Index of PPG peaks. + HR : numpy.ndarray + Instantaneous heart rate? + cardiac : 1D numpy.ndarray + Raw PPG signal. + + Returns + ------- + cpm_regressors + cpm_amplitude_regressors + retroicor_regressors + + Notes + ----- + CPM was first developed in [1]_ and the original code is implemented in MATLAB. + The original implementation, from which this code has been adapted, is released with an Apache + 2.0 license. + More information about the CPM license can be found in the ``phys2denoise`` LICENSE. + + Here we log all meaningful changes made to this code from the original CPM code: + 1. Translation from MATLAB to Python. + + References + ---------- + .. [1]: Kassinopoulos, M., & Mitsis, G. D. (2020). + Physiological noise modeling in fMRI based on the pulsatile + component of photoplethysmograph. bioRxiv. + https://doi.org/10.1101/2020.06.01.128306 + """ + timestep = 1 / samplerate + + # number of seconds by which to shift each regressor + SHIFT_CPM = 0.5 + SHIFT_CPM_VA = 0.5 + SHIFT_RETR = 0 + + MODEL_ORDER = 8 + + # Get filter parameters + HIGH_PASS_FILTER_FREQ = 0.008 # in Hertz + filt_b, filt_a = signal.butter(2, HIGH_PASS_FILTER_FREQ * 2 * timestep, "high") + + N_SECS_TO_REMOVE = 5 # delete first and last 5 seconds from output + n_vals_to_remove = int(N_SECS_TO_REMOVE * samplerate) + + HRmean = np.mean(HR) # mean heart-rate in beats-per-second? + memory = 60 / HRmean # average minutes-per-beat? + + cardiac_reduced = cardiac[n_vals_to_remove:-n_vals_to_remove] + n_timepoints_reduced = len(cardiac_reduced) + # high-pass filter the reduced cardiac signal + cardiac_reduced = signal.filtfilt(filt_b, filt_a, cardiac_reduced) + + n_timepoints = len(cardiac) + time = np.arange(0, n_timepoints * timestep, timestep) + + cpm_ir = func_CPM_cos(timestep, memory, MODEL_ORDER) + + # peaks in timeseries form (peaks are 1, all other timepoints are 0) + peak_timeseries = np.zeros(time.shape) + # peak amplitudes in timeseries form (peaks are amplitude value, all other timepoints are 0) + peak_amplitudes = peak_timeseries.copy() + for peak_time in PPGlocs: + time_distance = np.abs(time - peak_time) + closest_time_idx = np.argmin(time_distance) + peak_timeseries[closest_time_idx] = 1 + peak_amplitudes[closest_time_idx] = cardiac[closest_time_idx] + + cpm_regressors = np.zeros((n_timepoints, MODEL_ORDER * 2)) + cpm_amplitude_regressors = np.zeros((n_timepoints, MODEL_ORDER * 2)) + for i_col in range(MODEL_ORDER * 2): + cpm_nonamplitudes = np.convolve(peak_timeseries, cpm_ir[:, i_col]) + cpm_nonamplitudes = cpm_nonamplitudes[:n_timepoints] + cpm_regressors[:, i_col] = cpm_nonamplitudes + cpm_amplitudes = np.convolve(peak_amplitudes, cpm_ir[:, i_col]) + cpm_amplitudes = cpm_amplitudes[:n_timepoints] + cpm_amplitude_regressors[:, i_col] = cpm_amplitudes + + cpm_regressors = signal.filtfilt(filt_b, filt_a, cpm_regressors) + cpm_amplitude_regressors = signal.filtfilt(filt_b, filt_a, cpm_amplitude_regressors) + + retroicor_regressors = RETR_Card_regressors_v2(time, PPGlocs, MODEL_ORDER) + retroicor_regressors = signal.filtfilt(filt_b, filt_a, retroicor_regressors) + + # Select relevant timepoints from regressors. + idx = np.arange(n_timepoints_reduced) + temp_idx = np.round(idx + n_vals_to_remove + (SHIFT_CPM * samplerate)).astype(int) + cpm_regressors = cpm_regressors[temp_idx, :] + + temp_idx = np.round(idx + n_vals_to_remove + (SHIFT_CPM_VA * samplerate)).astype(int) + cpm_amplitude_regressors = cpm_amplitude_regressors[temp_idx, :] + + temp_idx = np.round(idx + n_vals_to_remove + (SHIFT_RETR * samplerate)).astype(int) + retroicor_regressors = retroicor_regressors[temp_idx, :] + + # Append ones? + cpm_regressors = np.hstack((cpm_regressors, np.ones((n_timepoints_reduced, 1)))) + cpm_amplitude_regressors = np.hstack( + (cpm_amplitude_regressors, np.ones((n_timepoints_reduced, 1))) + ) + retroicor_regressors = np.hstack( + (retroicor_regressors, np.ones((n_timepoints_reduced, 1))) + ) + return cpm_regressors, cpm_amplitude_regressors, retroicor_regressors def RETR_Card_regressors_v2(time, locsECG, M): - """Calculate RETROICOR cardiac regressors.""" - NV = len(time) - Phi = np.zeros((NV, 1)) - for i in range(NV): - t = time[i] - _, minI = np.min(np.abs(locsECG - t)) - - minOnLeft = t - locsECG[minI] > 0 - if (minI == 1) and not minOnLeft: + """Calculate RETROICOR cardiac regressors. + + Parameters + ---------- + time + locsECG + M + + Returns + ------- + Regr + """ + n_timepoints = len(time) + Phi = np.zeros((n_timepoints)) + for i_time, timepoint in enumerate(time): + minI = np.argmin(np.abs(locsECG - timepoint)) + + minOnLeft = timepoint - locsECG[minI] > 0 + if (minI == 0) and not minOnLeft: t2 = locsECG[minI] t1 = t2 - 1 - elif (minI == len(locsECG)) and minOnLeft: + elif (minI == (len(locsECG) - 1)) and minOnLeft: t1 = locsECG[minI] t2 = t1 + 1 elif minOnLeft: @@ -92,20 +149,34 @@ def RETR_Card_regressors_v2(time, locsECG, M): t1 = locsECG[minI - 1] t2 = locsECG[minI] - Phi[i] = 2 * np.pi * (t - t1) / (t2 - t1) + Phi[i_time] = 2 * np.pi * (timepoint - t1) / (t2 - t1) - Regr = np.zeros((NV, M * 2)) + Regr = np.zeros((n_timepoints, M * 2)) for i in range(M): Regr[:, ((i - 1) * 2) + 1] = np.cos(i * Phi) Regr[:, i * 2] = np.sin(i * Phi) - # Regr=[zeros(1,M*2);diff(Regr)]; return Regr -def func_CPM_cos(Ts, memory, M): - """Calculate CPM.""" - t_win = np.arange(0, memory, Ts) +def func_CPM_cos(timestep, memory, M): + """Calculate CPM? + + Parameters + ---------- + timestep : float + Sample rate of cardiac recording, in seconds. + memory : float + ??? + M : int + Model order? + + Returns + ------- + IR_all : numpy.ndarray + CPM metric at different model orders? + """ + t_win = np.arange(0, memory, timestep) nIR = len(t_win) IR_all = np.zeros((nIR, M * 2)) From 16771ef2615edf60a11627f7c858f5763a5724bd Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 7 Jan 2021 17:46:56 -0500 Subject: [PATCH 3/4] Reorganize arguments. --- phys2denoise/metrics/cardiac.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index fff4b37..ffa56ed 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -3,19 +3,19 @@ from scipy import signal -def cpm(samplerate, PPGlocs, HR, cardiac): +def cpm(cardiac, i_hr, peaks, samplerate): """Calculate cardiac pulsatility model (CPM) regressors. Parameters ---------- - samplerate : float - Sample rate in Hertz. - PPGlocs : 1D numpy.ndarray - Index of PPG peaks. - HR : numpy.ndarray - Instantaneous heart rate? cardiac : 1D numpy.ndarray Raw PPG signal. + i_hr : numpy.ndarray + Instantaneous heart rate? + peaks : 1D numpy.ndarray + Index of PPG peaks. + samplerate : float + Sample rate in Hertz. Returns ------- @@ -56,7 +56,7 @@ def cpm(samplerate, PPGlocs, HR, cardiac): N_SECS_TO_REMOVE = 5 # delete first and last 5 seconds from output n_vals_to_remove = int(N_SECS_TO_REMOVE * samplerate) - HRmean = np.mean(HR) # mean heart-rate in beats-per-second? + HRmean = np.mean(i_hr) # mean heart-rate in beats-per-second? memory = 60 / HRmean # average minutes-per-beat? cardiac_reduced = cardiac[n_vals_to_remove:-n_vals_to_remove] @@ -73,7 +73,7 @@ def cpm(samplerate, PPGlocs, HR, cardiac): peak_timeseries = np.zeros(time.shape) # peak amplitudes in timeseries form (peaks are amplitude value, all other timepoints are 0) peak_amplitudes = peak_timeseries.copy() - for peak_time in PPGlocs: + for peak_time in peaks: time_distance = np.abs(time - peak_time) closest_time_idx = np.argmin(time_distance) peak_timeseries[closest_time_idx] = 1 @@ -92,7 +92,7 @@ def cpm(samplerate, PPGlocs, HR, cardiac): cpm_regressors = signal.filtfilt(filt_b, filt_a, cpm_regressors) cpm_amplitude_regressors = signal.filtfilt(filt_b, filt_a, cpm_amplitude_regressors) - retroicor_regressors = RETR_Card_regressors_v2(time, PPGlocs, MODEL_ORDER) + retroicor_regressors = RETR_Card_regressors_v2(time, peaks, MODEL_ORDER) retroicor_regressors = signal.filtfilt(filt_b, filt_a, retroicor_regressors) # Select relevant timepoints from regressors. From fce3ccf1b2494815135717fda66d3fbae89ad19e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 11 Feb 2021 10:58:31 -0500 Subject: [PATCH 4/4] Fix style issue. --- phys2denoise/metrics/cardiac.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index ffa56ed..3eb6006 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -160,7 +160,7 @@ def RETR_Card_regressors_v2(time, locsECG, M): def func_CPM_cos(timestep, memory, M): - """Calculate CPM? + """Calculate CPM cosine values. Parameters ----------