diff --git a/neurokit2/misc/report.py b/neurokit2/misc/report.py index 2da99ea3cf..14f3421972 100644 --- a/neurokit2/misc/report.py +++ b/neurokit2/misc/report.py @@ -100,7 +100,7 @@ def summarize_table(signals): def text_combine(info): """Reformat dictionary describing processing methods as strings to be inserted into HTML file.""" preprocessing = '

Preprocessing

' - for key in ["text_cleaning", "text_peaks"]: + for key in ["text_cleaning", "text_peaks", "text_quality"]: if key in info.keys(): preprocessing += info[key] + "
" ref = '

References

' diff --git a/neurokit2/ppg/__init__.py b/neurokit2/ppg/__init__.py index 6da006006c..c56ef7911f 100644 --- a/neurokit2/ppg/__init__.py +++ b/neurokit2/ppg/__init__.py @@ -11,6 +11,7 @@ from .ppg_peaks import ppg_peaks from .ppg_plot import ppg_plot from .ppg_process import ppg_process +from .ppg_quality import ppg_quality from .ppg_segment import ppg_segment from .ppg_simulate import ppg_simulate @@ -23,6 +24,7 @@ "ppg_rate", "ppg_process", "ppg_plot", + "ppg_quality", "ppg_methods", "ppg_intervalrelated", "ppg_eventrelated", diff --git a/neurokit2/ppg/ppg_methods.py b/neurokit2/ppg/ppg_methods.py index 57554c7cd3..88eeea3597 100644 --- a/neurokit2/ppg/ppg_methods.py +++ b/neurokit2/ppg/ppg_methods.py @@ -4,6 +4,7 @@ from ..misc.report import get_kwargs from .ppg_clean import ppg_clean from .ppg_findpeaks import ppg_findpeaks +from .ppg_quality import ppg_quality def ppg_methods( @@ -11,6 +12,7 @@ def ppg_methods( method="elgendi", method_cleaning="default", method_peaks="default", + method_quality="default", **kwargs, ): """**PPG Preprocessing Methods** @@ -38,6 +40,11 @@ def ppg_methods( will be set to the value of ``"method"``. Defaults to ``"default"``. For more information, see the ``"method"`` argument of :func:`.ppg_findpeaks`. + method_quality: str + The method used to assess PPG signal quality. If ``"default"``, + will be set to the value of ``"templatematch"``. Defaults to ``"templatematch"``. + For more information, see the ``"method"`` argument + of :func:`.ppg_quality`. **kwargs Other arguments to be passed to :func:`.ppg_clean` and :func:`.ppg_findpeaks`. @@ -51,7 +58,7 @@ def ppg_methods( See Also -------- - ppg_process, ppg_clean, ppg_findpeaks + ppg_process, ppg_clean, ppg_findpeaks, ppg_quality Examples -------- @@ -59,7 +66,9 @@ def ppg_methods( import neurokit2 as nk - methods = nk.ppg_methods(sampling_rate=100, method="elgendi", method_cleaning="nabian2018") + methods = nk.ppg_methods( + sampling_rate=100, method="elgendi", + method_cleaning="nabian2018", method_quality="templatematch") print(methods["text_cleaning"]) print(methods["references"][0]) @@ -71,7 +80,12 @@ def ppg_methods( else str(method_cleaning).lower() ) method_peaks = ( - str(method).lower() if method_peaks == "default" else str(method_peaks).lower() + str(method).lower() + if method_peaks == "default" + else str(method_peaks).lower() + ) + method_quality = ( + str(method_quality).lower() ) # Create dictionary with all inputs @@ -80,16 +94,19 @@ def ppg_methods( "method": method, "method_cleaning": method_cleaning, "method_peaks": method_peaks, + "method_quality": method_quality, **kwargs, } - # Get arguments to be passed to cleaning and peak finding functions + # Get arguments to be passed to cleaning, peak finding, and quality assessment functions kwargs_cleaning, report_info = get_kwargs(report_info, ppg_clean) kwargs_peaks, report_info = get_kwargs(report_info, ppg_findpeaks) + kwargs_quality, report_info = get_kwargs(report_info, ppg_quality) # Save keyword arguments in dictionary report_info["kwargs_cleaning"] = kwargs_cleaning report_info["kwargs_peaks"] = kwargs_peaks + report_info["kwargs_quality"] = kwargs_quality # Initialize refs list with NeuroKit2 reference refs = ["""Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., @@ -158,5 +175,40 @@ def ppg_methods( "text_peaks" ] = f"The peak detection was carried out using the method {method_peaks}." + # 2. Quality + # ---------- + if method_quality in ["templatematch"]: + report_info[ + "text_quality" + ] = ( + "The quality assessment was carried out using template-matching, approximately as described " + + "in Orphanidou et al. (2015)." + ) + refs.append( + """Orphanidou C, Bonnici T, Charlton P, Clifton D, Vallance D, Tarassenko L (2015) + Signal-quality indices for the electrocardiogram and photoplethysmogram: Derivation + and applications to wireless monitoring + IEEE Journal of Biomedical and Health Informatics 19(3): 832–838. doi:10.1109/JBHI.2014.2338351.""" + ) + elif method_quality in ["disimilarity"]: + report_info[ + "text_quality" + ] = ( + "The quality assessment was carried out using a disimilarity measure of positive-peaked beats, " + + "approximately as described in Sabeti et al. (2019)." + ) + refs.append( + """Sabeti E, Reamaroon N, Mathis M, Gryak J, Sjoding M, Najarian K (2019) + Signal quality measure for pulsatile physiological signals using + morphological features: Applications in reliability measure for pulse oximetry + Informatics in Medicine Unlocked 16: 100222. doi:10.1016/j.imu.2019.100222.""" + ) + elif method_quality in ["none"]: + report_info["text_quality"] = "There was no quality assessment carried out." + else: + report_info[ + "text_quality" + ] = f"The quality assessment was carried out using the method {method_quality}." + report_info["references"] = list(np.unique(refs)) return report_info diff --git a/neurokit2/ppg/ppg_peaks.py b/neurokit2/ppg/ppg_peaks.py index 6b0eaf7240..e7ec31bddf 100644 --- a/neurokit2/ppg/ppg_peaks.py +++ b/neurokit2/ppg/ppg_peaks.py @@ -3,6 +3,7 @@ from ..ecg.ecg_peaks import _ecg_peaks_plot_artefacts from ..signal import signal_fixpeaks, signal_formatpeaks +from ..stats import rescale from .ppg_findpeaks import ppg_findpeaks @@ -143,6 +144,7 @@ def _ppg_peaks_plot( info=None, sampling_rate=1000, raw=None, + quality=None, ax=None, ): x_axis = np.linspace(0, len(ppg_cleaned) / sampling_rate, len(ppg_cleaned)) @@ -154,6 +156,29 @@ def _ppg_peaks_plot( ax.set_xlabel("Time (seconds)") ax.set_title("PPG signal and peaks") + # Quality Area ------------------------------------------------------------- + if quality is not None: + quality = rescale( + quality, + to=[ + np.min([np.min(raw), np.min(ppg_cleaned)]), + np.max([np.max(raw), np.max(ppg_cleaned)]), + ], + ) + minimum_line = np.full(len(x_axis), quality.min()) + + # Plot quality area first + ax.fill_between( + x_axis, + minimum_line, + quality, + alpha=0.12, + zorder=0, + interpolate=True, + facecolor="#4CAF50", + label="Signal quality", + ) + # Raw Signal --------------------------------------------------------------- if raw is not None: ax.plot(x_axis, raw, color="#B0BEC5", label="Raw signal", zorder=1) diff --git a/neurokit2/ppg/ppg_plot.py b/neurokit2/ppg/ppg_plot.py index b77905e4d9..cffbec8509 100644 --- a/neurokit2/ppg/ppg_plot.py +++ b/neurokit2/ppg/ppg_plot.py @@ -92,6 +92,7 @@ def ppg_plot(ppg_signals, info=None, static=True): info=info, sampling_rate=info["sampling_rate"], raw=ppg_signals["PPG_Raw"].values, + quality=ppg_signals["PPG_Quality"].values, ax=ax0, ) diff --git a/neurokit2/ppg/ppg_process.py b/neurokit2/ppg/ppg_process.py index e90a5d5880..21a94993ba 100644 --- a/neurokit2/ppg/ppg_process.py +++ b/neurokit2/ppg/ppg_process.py @@ -8,10 +8,11 @@ from .ppg_methods import ppg_methods from .ppg_peaks import ppg_peaks from .ppg_plot import ppg_plot +from .ppg_quality import ppg_quality def ppg_process( - ppg_signal, sampling_rate=1000, method="elgendi", report=None, **kwargs + ppg_signal, sampling_rate=1000, method="elgendi", method_quality="templatematch", report=None, **kwargs ): """**Process a photoplethysmogram (PPG) signal** @@ -26,6 +27,9 @@ def ppg_process( method : str The processing pipeline to apply. Can be one of ``"elgendi"``. Defaults to ``"elgendi"``. + method_quality : str + The quality assessment approach to use. Can be one of ``"templatematch"``, ``"disimilarity"``. + Defaults to ``"templatematch"``. report : str The filename of a report containing description and figures of processing (e.g. ``"myreport.html"``). Needs to be supplied if a report file @@ -69,7 +73,7 @@ def ppg_process( """ # Sanitize input ppg_signal = as_vector(ppg_signal) - methods = ppg_methods(sampling_rate=sampling_rate, method=method, **kwargs) + methods = ppg_methods(sampling_rate=sampling_rate, method=method, method_quality=method_quality, **kwargs) # Clean signal ppg_cleaned = ppg_clean( @@ -94,12 +98,22 @@ def ppg_process( info["PPG_Peaks"], sampling_rate=sampling_rate, desired_length=len(ppg_cleaned) ) + # Assess signal quality + quality = ppg_quality( + ppg_cleaned, + ppg_pw_peaks=info["PPG_Peaks"], + sampling_rate=sampling_rate, + method=methods["method_quality"], + **methods["kwargs_quality"] + ) + # Prepare output signals = pd.DataFrame( { "PPG_Raw": ppg_signal, "PPG_Clean": ppg_cleaned, "PPG_Rate": rate, + "PPG_Quality": quality, "PPG_Peaks": peaks_signal["PPG_Peaks"].values, } ) diff --git a/neurokit2/ppg/ppg_quality.py b/neurokit2/ppg/ppg_quality.py new file mode 100644 index 0000000000..c53a609f94 --- /dev/null +++ b/neurokit2/ppg/ppg_quality.py @@ -0,0 +1,196 @@ +# - * - coding: utf-8 - * - + +import numpy as np + +from ..epochs import epochs_to_df +from ..signal import signal_interpolate +from .ppg_peaks import ppg_peaks +from .ppg_segment import ppg_segment + + +def ppg_quality( + ppg_cleaned, ppg_pw_peaks=None, sampling_rate=1000, method="templatematch", approach=None +): + """**PPG Signal Quality Assessment** + + Assess the quality of the PPG Signal using various methods: + + * The ``"templatematch"`` method (loosely based on Orphanidou et al., 2015) computes a continuous + index of quality of the PPG signal, by calculating the correlation coefficient between each + individual pulse wave and an average (template) pulse wave shape. This index is therefore + relative: 1 corresponds to pulse waves that are closest to the average pulse wave shape (i.e. + correlate exactly with it) and 0 corresponds to there being no correlation with the average + pulse wave shape. Note that 1 does not necessarily mean "good": use this index with care and + plot it alongside your PPG signal to see if it makes sense. + + * The ``"disimilarity"`` method (loosely based on Sabeti et al., 2019) computes a continuous index + of quality of the PPG signal, by calculating the level of disimilarity between each individual + pulse wave and an average (template) pulse wave shape (after they are normalised). A value of + zero indicates no disimilarity (i.e. equivalent pulse wave shapes), whereas values above or below + indicate increasing disimilarity. The original method used dynamic time-warping to align the pulse + waves prior to calculating the level of dsimilarity, whereas this implementation does not currently + include this step. + + Parameters + ---------- + ppg_cleaned : Union[list, np.array, pd.Series] + The cleaned PPG signal in the form of a vector of values. + ppg_pw_peaks : tuple or list + The list of PPG pulse wave peak samples returned by ``ppg_peaks()``. If None, peaks is computed from + the signal input. + sampling_rate : int + The sampling frequency of the signal (in Hz, i.e., samples/second). + method : str + The method for computing PPG signal quality, can be ``"templatematch"`` (default). + + Returns + ------- + array + Vector containing the quality index ranging from 0 to 1 for ``"templatematch"`` method, + or an unbounded value (where 0 indicates high quality) for ``"disimilarity"`` method. + + See Also + -------- + ppg_segment + + References + ---------- + * Orphanidou, C. et al. (2015). "Signal-quality indices for the electrocardiogram and photoplethysmogram: + derivation and applications to wireless monitoring". IEEE Journal of Biomedical and Health Informatics, 19(3), 832-8. + + Examples + -------- + * **Example 1:** 'templatematch' method + + .. ipython:: python + + import neurokit2 as nk + + ppg = nk.ppg_simulate(duration=30, sampling_rate=300, heart_rate=80) + ppg_cleaned = nk.ppg_clean(ppg, sampling_rate=300) + quality = nk.ppg_quality(ppg_cleaned, sampling_rate=300, method="templatematch") + + @savefig p_ppg_quality.png scale=100% + nk.signal_plot([ppg_cleaned, quality], standardize=True) + @suppress + plt.close() + + """ + + method = method.lower() # remove capitalised letters + + # Run selected quality assessment method + if method in ["templatematch"]: + quality = _ppg_quality_templatematch( + ppg_cleaned, ppg_pw_peaks=ppg_pw_peaks, sampling_rate=sampling_rate + ) + elif method in ["disimilarity"]: + quality = _ppg_quality_disimilarity( + ppg_cleaned, ppg_pw_peaks=ppg_pw_peaks, sampling_rate=sampling_rate + ) + + return quality + +# ============================================================================= +# Calculate a template pulse wave +# ============================================================================= +def _calc_template_pw(ppg_cleaned, ppg_pw_peaks=None, sampling_rate=1000): + + # Sanitize inputs + if ppg_pw_peaks is None: + _, ppg_pw_peaks = ppg_peaks(ppg_cleaned, sampling_rate=sampling_rate) + ppg_pw_peaks = ppg_pw_peaks["PPG_Peaks"] + + # Get heartbeats + heartbeats = ppg_segment(ppg_cleaned, ppg_pw_peaks, sampling_rate) + pw_data = epochs_to_df(heartbeats).pivot( + index="Label", columns="Time", values="Signal" + ) + pw_data.index = pw_data.index.astype(int) + pw_data = pw_data.sort_index() + + # Filter Nans + missing = pw_data.T.isnull().sum().values + nonmissing = np.where(missing == 0)[0] + pw_data = pw_data.iloc[nonmissing, :] + + # Find template pulse wave + templ_pw = pw_data.mean() + + return templ_pw, pw_data, ppg_pw_peaks + + +# ============================================================================= +# Template-matching method +# ============================================================================= +def _ppg_quality_templatematch(ppg_cleaned, ppg_pw_peaks=None, sampling_rate=1000): + + # Obtain individual pulse waves and template pulse wave + templ_pw, pw_data, ppg_pw_peaks = _calc_template_pw( + ppg_cleaned, ppg_pw_peaks=ppg_pw_peaks, sampling_rate=sampling_rate + ) + + # Find individual correlation coefficients (CCs) + cc = np.zeros(len(ppg_pw_peaks)-1) + for beat_no in range(0,len(ppg_pw_peaks)-1): + temp = np.corrcoef(pw_data.iloc[beat_no], templ_pw) + cc[beat_no] = temp[0,1] + + # Interpolate beat-by-beat CCs + quality = signal_interpolate( + ppg_pw_peaks[0:-1], cc, x_new=np.arange(len(ppg_cleaned)), method="quadratic" + ) + + return quality + +# ============================================================================= +# Disimilarity measure method +# ============================================================================= +def _norm_sum_one(pw): + + # ensure all values are positive + pw = pw - pw.min() + 1 + + # normalise pulse wave to sum to one + pw = [x / sum(pw) for x in pw] + + return pw + +def _calc_dis(pw1, pw2): + # following the methodology in https://doi.org/10.1016/j.imu.2019.100222 (Sec. 3.1.2.5) + + # convert to numpy arrays + pw1 = np.array(pw1) + pw2 = np.array(pw2) + + # normalise to sum to one + pw1 = _norm_sum_one(pw1) + pw2 = _norm_sum_one(pw2) + + # ignore any elements which are zero because log(0) is -inf + rel_els = (pw1 != 0) & (pw2 != 0) + + # calculate disimilarity measure (using pw2 as the template) + dis = np.sum(pw2[rel_els] * np.log(pw2[rel_els] / pw1[rel_els])) + + return dis + + +def _ppg_quality_disimilarity(ppg_cleaned, ppg_pw_peaks=None, sampling_rate=1000): + + # Obtain individual pulse waves and template pulse wave + templ_pw, pw_data, ppg_pw_peaks = _calc_template_pw( + ppg_cleaned, ppg_pw_peaks=ppg_pw_peaks, sampling_rate=sampling_rate + ) + + # Find individual disimilarity measures + dis = np.zeros(len(ppg_pw_peaks)-1) + for beat_no in range(0,len(ppg_pw_peaks)-1): + dis[beat_no] = _calc_dis(pw_data.iloc[beat_no], templ_pw) + + # Interpolate beat-by-beat dis's + quality = signal_interpolate( + ppg_pw_peaks[0:-1], dis, x_new=np.arange(len(ppg_cleaned)), method="previous" + ) + + return quality \ No newline at end of file