Merge pull request #3 from SRSteinkamp/peak_metrics

WIP: Peak metrics

Author: smoia

physioqc/__init__.py

@@ -1,17 +1,6 @@
 """Hopefully importing everything."""
 
-import pkgutil
-
 from ._version import get_versions
 
-SKIP_MODULES = ["tests"]
-
 __version__ = get_versions()["version"]
 del get_versions
-
-__all__ = []
-for loader, module_name, is_pkg in pkgutil.walk_packages(__path__):
-    if "tests" not in module_name:
-        __all__.append(module_name)
-        _module = loader.find_module(module_name).load_module(module_name)
-        globals()[module_name] = _module

physioqc/metrics/multimodal.py

@@ -1,9 +1,8 @@
-"""These functions compute RETROICOR regressors (Glover et al. 2000)."""
-# import file
+"""These functions compute various non-modality dependent signal processing metrics."""
 import numpy as np
+import peakdet as pk
 
-from .. import references
-from ..due import due
+from .utils import physio_or_numpy
 
 
 def std(signal):
@@ -12,38 +11,222 @@ def std(signal):
 
     Parameters
     ----------
-    std : function
-        Calculate standard deviation across input channels.
-    args : signal (X, n) :obj:`numpy.ndarray`
-        ND array with signal of some human biometric data, hopefully from a living human.
-        signal, of shape (n_channels, )
+    signal : np.array
+        Physiological data
 
     Returns
     -------
     N-sized array :obj:`numpy.ndarray`
         Standard deviation of signal.
     """
+    signal = physio_or_numpy(signal)
     std_val = np.std(signal, axis=0)
     return std_val
 
 
+def mean(signal: np.array):
+    """
+    Calculate mean across input channels of signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Physiological data
+
+    Returns
+    -------
+    N-sized array :obj:`numpy.ndarray`
+        Mean of signal.
+    """
+    signal = physio_or_numpy(signal)
+    mean_val = np.mean(signal, axis=0)
+    return mean_val
+
+
 def tSNR(signal):
     """
-    Calculate temporal signal to noise ration of signal.
+    Calculate temporal signal to noise ratio of signal.
 
     Parameters
     ----------
-    tSNR : function
-        Calculate temporal signal to noise ratio.
-    args : signal
-        ND array with signal of some human biometric data, hopefully from a living human.
+    signal : np.array
+        Physiological data
 
     Returns
     -------
     N-sized array :obj:`numpy.ndarray`
         Temporal signal to noise ratio of signal.
     """
-    me = np.mean(signal, axis=0)
-    std = np.std(signal, axis=0)
-    tSNR_val = me / std
+    signal = physio_or_numpy(signal)
+    tSNR_val = np.mean(signal, axis=0) / np.std(signal, axis=0)
     return tSNR_val
+
+
+def CoV(signal):
+    """
+    Calculate coefficient of variation of signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Physiological data
+
+    Returns
+    -------
+    N-sized array :obj:`numpy.ndarray`
+        Temporal signal to noise ratio of signal.
+    """
+    signal = physio_or_numpy(signal)
+    tSNR_val = np.std(signal, axis=0) / np.mean(signal, axis=0)
+    return tSNR_val
+
+
+def min(signal: np.array):
+    """
+    Calculate min across input channels of signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Physiological data
+
+    Returns
+    -------
+    N-sized array :obj:`numpy.ndarray`
+        min of signal.
+    """
+    signal = physio_or_numpy(signal)
+    min_val = np.min(signal, axis=0)
+    return min_val
+
+
+def max(signal: np.array):
+    """
+    Calculate max across input channels of signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Physiological data
+
+    Returns
+    -------
+    N-sized array :obj:`numpy.ndarray`
+        max of signal.
+    """
+    signal = physio_or_numpy(signal)
+    max_val = np.max(signal, axis=0)
+    return max_val
+
+
+def iqr(signal: np.array, q_high: float = 75, q_low: float = 25):
+    """Calculate the Inter Quantile Range (IQR) over the input signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Physiological data
+    q_high : float, optional
+        higher percentile for IQR, by default 75
+    q_low : float, optional
+        lower percentile for IQR, by default 25
+
+    Returns
+    -------
+    np.array
+        iqr of the signal
+    """
+    signal = physio_or_numpy(signal)
+    p_high, p_low = np.percentile(signal, [q_high, q_low], axis=0)
+    iqr_val = p_high - p_low
+
+    return iqr_val
+
+
+def percentile(signal: np.array, perc: float = 2):
+    """Calculate the percentile perc over the signal.
+
+    Parameters
+    ----------
+    signal : np.array
+        Physiological data
+    perc : float, optional
+        percentile for data, by default 2
+
+    Returns
+    -------
+    np.array
+        percentile of the signal
+    """
+    signal = physio_or_numpy(signal)
+    perc_val = np.percentile(signal, axis=0, q=perc)
+
+    return perc_val
+
+
+def peak_detection(
+    data: pk.Physio,
+    peak_threshold: float = 0.1,
+    peak_dist: int = 60,
+):
+    """
+    Perform peak detection for further metric extraction and plotting.
+
+    Parameters
+    ----------
+    data : pk.Physio
+        A peakdet Physio object
+    peak_threshold : float, optional
+        Threshold for peak detection, by default 0.1
+    peak_dist : int, optional
+        Distance for peak detection, by default 60
+
+    Returns
+    -------
+    pk.Physio
+        Updated pk.Physio class with peaks etc.
+    """
+    ph = pk.operations.peakfind_physio(data, thresh=peak_threshold, dist=peak_dist)
+
+    return ph
+
+
+def peak_distance(ph: pk.Physio):
+    """Calculate the time between peaks (first derivative of onsets).
+
+    Parameters
+    ----------
+    ph : pk.Physio
+        A pk.Physio object, that contains peak information.
+
+    Returns
+    -------
+    np.array
+        np.array of shape [npeaks, ]
+    """
+    # TODO Check if peaks have been estimated.
+    diff_peak = np.diff(ph.peaks, axis=0)
+
+    return diff_peak
+
+
+def peak_amplitude(ph: pk.Physio):
+    """Return the amplitude for each peak in the ph.Physio object (peak - trough).
+
+    Parameters
+    ----------
+    ph : pk.Physio
+        pk.Physio object with peak and trough information.
+
+    Returns
+    -------
+    np.array
+        np.array of shape [npeaks - 1, ]
+    """
+    # TODO Check if peaks have been estimated.
+    # Assuming that peaks and troughs are aligned. Last peak has no trough.
+    peak_amp = ph.data[ph.peaks[:-1]]
+    trough_amp = ph.data[ph.troughs]
+    peak_amplitude = peak_amp - trough_amp
+
+    return peak_amplitude

physioqc/metrics/utils.py

@@ -1,8 +1,6 @@
 """Miscellaneous utility functions for metric calculation."""
 import logging
 
-import numpy as np
-
 LGR = logging.getLogger(__name__)
 LGR.setLevel(logging.INFO)
 
@@ -31,3 +29,23 @@ def print_metric_call(metric, args):
     msg = f"{msg}\n"
 
     LGR.info(msg)
+
+
+def physio_or_numpy(signal):
+    """
+    Return data from a peakdet.physio.Physio object or a np.ndarray-like object.
+
+    Parameters
+    ----------
+    data : peakdet.physio.Physio, np.ndarray, or list
+        object to get data from
+
+    Returns
+    -------
+    np.ndarray-like object
+        Either a np.ndarray or a list
+    """
+    if hasattr(signal, "history"):
+        signal = signal.data
+
+    return signal