# -*- coding: utf-8 -*-
from warnings import warn

import numpy as np
import pandas as pd
import scipy.signal

from ..misc import NeuroKitWarning, as_vector
from ..signal import signal_detrend, signal_filter
from ..stats import mad



[docs]
def rsp_clean(rsp_signal, sampling_rate=1000, method="khodadad2018", **kwargs):
    """**Preprocess a respiration (RSP) signal**

    Clean a respiration signal using different sets of parameters, such as:

    * **khodadad2018**: Second order 0.05-3 Hz bandpass Butterworth filter. Note that the
      implementation differs from the referenced paper (see issue #950).
    * **BioSPPy**: Second order 0.1-0.35 Hz bandpass Butterworth filter followed by a constant
      detrending).
    * **hampel**: Applies a median-based Hampel filter by replacing values which are 3 (can be
      changed via ``threshold``) :func:`.mad` away from the rolling median.

    Parameters
    ----------
    rsp_signal : Union[list, np.array, pd.Series]
        The raw respiration channel (as measured, for instance, by a respiration belt).
    sampling_rate : int, optional
        The sampling frequency of :func:`.rsp_signal` (in Hz, i.e., samples/second).
    method : str, optional
        The processing pipeline to apply. Can be one of ``"khodadad2018"`` (default),
        ``"biosppy"`` or ``"hampel"``.
    **kwargs
        Other arguments to pass to the cleaning method.

    Returns
    -------
    array
        Vector containing the cleaned respiratory signal.

    See Also
    --------
    rsp_findpeaks, signal_rate, rsp_amplitude, rsp_process, rsp_plot

    Examples
    --------
    .. ipython:: python

      import pandas as pd
      import neurokit2 as nk

      rsp = nk.rsp_simulate(duration=30, sampling_rate=50, noise=0.1)
      signals = pd.DataFrame({
          "RSP_Raw": rsp,
          "RSP_Khodadad2018": nk.rsp_clean(rsp, sampling_rate=50, method="khodadad2018"),
          "RSP_BioSPPy": nk.rsp_clean(rsp, sampling_rate=50, method="biosppy"),
          "RSP_Hampel": nk.rsp_clean(rsp, sampling_rate=50, method="hampel", threshold=3)
      })
      @savefig p_rsp_clean1.png scale=100%
      signals.plot()
      @suppress
      plt.close()

    References
    ----------
    * Khodadad, D., Nordebo, S., Müller, B., Waldmann, A., Yerworth, R., Becher, T., ... & Bayford,
      R. (2018). Optimized breath detection algorithm in electrical impedance tomography.
      Physiological measurement, 39(9), 094001.
    * Power, J., Lynch, C., Dubin, M., Silver, B., Martin, A., Jones, R.,(2020)
      Characteristics of respiratory measures in young adults scanned at rest,
      including systematic changes and “missed” deep breaths.
      NeuroImage, Volume 204, 116234

    """
    rsp_signal = as_vector(rsp_signal)

    # Missing data
    n_missing = np.sum(np.isnan(rsp_signal))
    if n_missing > 0:
        warn(
            f"There are {n_missing} missing data points in your signal."
            " Filling missing values by using the forward filling method.",
            category=NeuroKitWarning,
        )
        rsp_signal = _rsp_clean_missing(rsp_signal)

    method = method.lower()  # remove capitalised letters
    if method in ["khodadad", "khodadad2018"]:
        clean = _rsp_clean_khodadad2018(rsp_signal, sampling_rate)
    elif method == "biosppy":
        clean = _rsp_clean_biosppy(rsp_signal, sampling_rate)
    elif method in ["power", "power2020", "hampel"]:
        clean = _rsp_clean_hampel(
            rsp_signal,
            **kwargs,
        )
    elif method is None or method == "none":
        clean = rsp_signal
    else:
        raise ValueError(
            "NeuroKit error: rsp_clean(): 'method' should be one of 'khodadad2018', 'biosppy' or 'hampel'."
        )

    return clean



# =============================================================================
# Handle missing data
# =============================================================================
def _rsp_clean_missing(rsp_signal):

    rsp_signal = pd.DataFrame.pad(pd.Series(rsp_signal))

    return rsp_signal


# =============================================================================
# Khodadad et al. (2018)
# =============================================================================
def _rsp_clean_khodadad2018(rsp_signal, sampling_rate=1000):
    """The algorithm is based on (but not an exact implementation of) the "Zero-crossing algorithm with amplitude
    threshold" by `Khodadad et al. (2018)

    <https://iopscience.iop.org/article/10.1088/1361-6579/aad7e6/meta>`_.

    """
    # Slow baseline drifts / fluctuations must be removed from the raw
    # breathing signal (i.e., the signal must be centered around zero) in order
    # to be able to reliable detect zero-crossings.

    # Remove baseline by applying a lowcut at .05Hz (preserves breathing rates
    # higher than 3 breath per minute) and high frequency noise by applying a
    # highcut at 3 Hz (preserves breathing rates slower than 180 breath per
    # minute).
    clean = signal_filter(
        rsp_signal,
        sampling_rate=sampling_rate,
        lowcut=0.05,
        highcut=3,
        order=2,
        method="butterworth",
    )

    return clean


# =============================================================================
# BioSPPy
# =============================================================================
def _rsp_clean_biosppy(rsp_signal, sampling_rate=1000):
    """Uses the same defaults as `BioSPPy.

    <https://github.com/PIA-Group/BioSPPy/blob/master/biosppy/signals/resp.py>`_.

    """
    # Parameters
    order = 2
    frequency = [0.1, 0.35]
    # Normalize frequency to Nyquist Frequency (Fs/2).
    frequency = 2 * np.array(frequency) / sampling_rate

    # Filtering
    b, a = scipy.signal.butter(N=order, Wn=frequency, btype="bandpass", analog=False)
    filtered = scipy.signal.filtfilt(b, a, rsp_signal)

    # Baseline detrending
    clean = signal_detrend(filtered, order=0)

    return clean


# =============================================================================
# Hampel filter
# =============================================================================
def _rsp_clean_hampel(rsp_signal, sampling_rate=1000, window_length=0.1, threshold=3, **kwargs):
    """Explanation MatLabs' https://www.mathworks.com/help/dsp/ref/hampelfilter.html. From
    https://stackoverflow.com/a/51731332.

    Parameters
    ----------
    rsp_signal : Union[list, np.array, pd.Series]
        The raw respiration channel (as measured, for instance, by a respiration belt).
    window_length : int, optional
        Window to be considered when cleaning, by default 0.1. In seconds.
    threshold : float, optional
        Threshold of deviations after which a point is considered an outlier, by default 3.

    """
    # Get window length in samples
    window_length = int(window_length * sampling_rate)

    # Convert to Series to use its rolling methods
    rsp_signal = pd.Series(rsp_signal)

    rolling_median = rsp_signal.rolling(window=window_length, center=True).median()
    rolling_MAD = rsp_signal.rolling(window=window_length, center=True).apply(mad)

    threshold = threshold * rolling_MAD
    difference = np.abs(rsp_signal - rolling_median)
    # Find outliers
    outlier_idx = difference > threshold
    # Substitute outliers with rolling median
    rsp_signal[outlier_idx] = rolling_median[outlier_idx]
    return as_vector(rsp_signal)