Add Hilbert RVT (Harrison et al. 2020)

phys2denoise/metrics/__init__.py

@@ -0,0 +1,3 @@
+from .chest_belt import hilbert_rvt
+
+__all__ = ["hilbert_rvt"]

phys2denoise/metrics/chest_belt.py

@@ -0,0 +1,69 @@
+"""Metrics derived from chest belt (respiration) recordings."""
+import numpy as np
+from scipy import signal
+
+from .utils import butter_bandpass_filter, interp_nans, split_complex
+
+
+def hilbert_rvt(resp, sampling_freq):
+    """Calculate a Hilbert transform-based version of the respiratory volume per unit time metric.
+
+    Parameters
+    ----------
+    resp : array, shape (n_timepoints,)
+        The respiratory trace signal in a 1D array.
+    sampling_freq : float
+        The sampling frequency, in Hertz.
+
+    Returns
+    -------
+    rvt : array, shape (n_timepoints,)
+        The respiratory volume per unit time metric.
+
+    References
+    ----------
+    * Harrison, S. J., Bianchi, S., Heinzle, J., Stephan, K. E., Iglesias, S.,
+      & Kasper, L. (2020). A Hilbert-based method for processing respiratory timeseries.
+      bioRxiv. https://doi.org/10.1101/2020.09.30.321562
+    """
+    # Remove low frequency drifts (less than 0.01 Hz) from the breathing signal,
+    # and remove high-frequency noise above 2.0 Hz.
+    # Lowpass filter the data again to more aggressively remove high-frequency noise above 0.75 Hz.
+    # Simplified to just bandpass filtering between 0.01-0.75 in one go.
+    resp_filt = butter_bandpass_filter(
+        resp, fs=sampling_freq, lowcut=0.01, highcut=0.75, order=1
+    )
+
+    # Decompose the signal into magnitude and phase components via the Hilbert transform.
+    analytic_signal = signal.hilbert(resp_filt)
+    magn, phas = split_complex(analytic_signal)
+
+    for i in range(10):
+        # Linearly interpolate any periods where the phase time course decreases,
+        # using the procedure in Figure 2, to remove any artefactual negative frequencies.
+        bad_idx = phas < 0
+        phas[bad_idx] = np.nan
+        phas = interp_nans(phas)
+        # Reconstruct the oscillatory portion of the signal, cos(𝜙(𝑡)),
+        # and lowpass filter at 0.75 Hz to remove any resulting artefacts.
+        oscill = np.cos(phas)
+        oscill = butter_bandpass_filter(oscill, fs=sampling_freq, highcut=0.75)
+        phas = np.arccos(oscill)
+        # This procedure is repeated 10 times, with the new phase timecourse
+        # re-estimated from the filtered oscillatory signal.
+
+    # instantaneous breathing rate is temporal derivative of phase
+    # not sure why pi is here but it's in the preprint's formula
+    ibr = np.append(0, np.diff(phas)) / (2 * np.pi)
+
+    # respiratory volume is twice signal amplitude
+    rv = 2 * magn
+
+    # low-pass filter both(?) at 0.2Hz
+    ibr = butter_bandpass_filter(ibr, fs=sampling_freq, highcut=0.75)
+    rv = butter_bandpass_filter(rv, fs=sampling_freq, highcut=0.75)
+
+    # rvt is product of ibr and rv
+    rvt = ibr * rv
+
+    return rvt

phys2denoise/metrics/utils.py

@@ -0,0 +1,167 @@
+"""Miscellaneous utility functions for metric calculations."""
+import numpy as np
+from scipy import signal
+
+
+def get_butter_filter(fs, lowcut=None, highcut=None, order=5):
+    """Calculate the appropriate parameters for a Butterworth filter.
+
+    Parameters
+    ----------
+    fs : float
+        Sampling frequency of data in Hertz.
+    lowcut : float or None, optional
+        Frequency (in Hertz) under which to remove in the filter.
+        If both lowcut and highcut are not None, then a bandpass filter is applied.
+        If lowcut is None and highcut is not, then a low-pass filter is applied.
+        If highcut is None and lowcut is not, then a high-pass filter is applied.
+        Either lowcut, highcut, or both must not be None.
+    highcut : float or None, optional
+        Frequency (in Hertz) above which to remove in the filter.
+        If both lowcut and highcut are not None, then a bandpass filter is applied.
+        If lowcut is None and highcut is not, then a low-pass filter is applied.
+        If highcut is None and lowcut is not, then a high-pass filter is applied.
+        Either lowcut, highcut, or both must not be None.
+    order : int, optional
+        Nonzero positive integer indicating order of the filter. Default is 1.
+
+    Returns
+    -------
+    b, a : float
+        Parameters from the Butterworth filter.
+
+    Notes
+    -----
+    From https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
+    """
+    nyq = 0.5 * fs
+
+    if (lowcut is not None) and (highcut is not None):
+        low = lowcut / nyq
+        high = highcut / nyq
+        window = (low, high)
+        btype = "bandpass"
+    elif lowcut is not None:
+        window = lowcut / nyq
+        btype = "highpass"
+    elif highcut is not None:
+        window = highcut / nyq
+        btype = "lowpass"
+    elif (lowcut is None) and (highcut is None):
+        raise ValueError("Either lowcut or highcut must be specified.")
+
+    b, a = signal.butter(order, window, btype=btype)
+    return b, a
+
+
+def butter_bandpass_filter(data, fs, lowcut=None, highcut=None, order=1):
+    """Band/low/high-pass filter an array based on cut frequencies.
+
+    Parameters
+    ----------
+    data : array, shape (n_timepoints,)
+        Data to filter.
+    fs : float
+        Sampling frequency of data in Hertz.
+    lowcut : float or None, optional
+        Frequency (in Hertz) under which to remove in the filter.
+        If both lowcut and highcut are not None, then a bandpass filter is applied.
+        If lowcut is None and highcut is not, then a low-pass filter is applied.
+        If highcut is None and lowcut is not, then a high-pass filter is applied.
+        Either lowcut, highcut, or both must not be None.
+    highcut : float or None, optional
+        Frequency (in Hertz) above which to remove in the filter.
+        If both lowcut and highcut are not None, then a bandpass filter is applied.
+        If lowcut is None and highcut is not, then a low-pass filter is applied.
+        If highcut is None and lowcut is not, then a high-pass filter is applied.
+        Either lowcut, highcut, or both must not be None.
+    order : int, optional
+        Nonzero positive integer indicating order of the filter. Default is 1.
+
+    Returns
+    -------
+    y : array, shape (n_timepoints,)
+        Filtered data.
+
+    Notes
+    -----
+    From https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
+    """
+    b, a = get_butter_filter(fs, lowcut, highcut, order=order)
+    y = signal.filtfilt(b, a, data)
+    return y
+
+
+def split_complex(complex_signal):
+    """Split a complex-valued nifti image into magnitude and phase images.
+
+    From complex-flow
+    """
+    real = complex_signal.real
+    imag = complex_signal.imag
+    mag = abs(complex_signal)
+    phase = to_phase(real, imag)
+    return mag, phase
+
+
+def to_phase(real, imag):
+    """Convert real and imaginary data to phase data.
+
+    Equivalent to cmath.phase.
+    https://www.eeweb.com/quizzes/convert-between-real-imaginary-and-magnitude-phase
+
+    From complex-flow
+    """
+    phase = np.arctan2(imag, real)
+    return phase
+
+
+def nan_helper(y):
+    """Handle indices and logical indices of NaNs.
+
+    Parameters
+    ----------
+    y : array, shape (n_timepoints,)
+        A 1D array with possible NaNs.
+
+    Returns
+    -------
+    nans : array, shape (n_timepoints,)
+        Logical indices of NaNs in y.
+    index : function
+        A function, with signature indices = index(logical_indices),
+        to convert logical indices of NaNs to 'equivalent' indices.
+
+    Examples
+    --------
+    >>> # linear interpolation of NaNs
+    >>> nans, x= nan_helper(y)
+    >>> y[nans]= np.interp(x(nans), x(~nans), y[~nans])
+
+    Notes
+    -----
+    From https://stackoverflow.com/a/6520696
+    """
+    return np.isnan(y), lambda z: z.nonzero()[0]
+
+
+def interp_nans(y):
+    """Linearly interpolate across NaNs in a 1D array.
+
+    Parameters
+    ----------
+    y : array, shape (n_timepoints,)
+        A 1D array with possible NaNs.
+
+    Returns
+    -------
+    y : array, shape (n_timepoints,)
+        A 1D array with no NaNs.
+
+    Notes
+    -----
+    From https://stackoverflow.com/a/6520696
+    """
+    nans, x = nan_helper(y)
+    y[nans] = np.interp(x(nans), x(~nans), y[~nans])
+    return y