Merge pull request #155 from ErlendHaa/rolling-median

Add rolling median filter processor

camille/process/__init__.py

@@ -14,6 +14,7 @@
 * :func:`~camille.process.high_pass`
 * :func:`~camille.process.band_pass`
 * :func:`~camille.process.mooring_fatigue`
+* :func:`~camille.process.rolling_median`
 """
 
 from .atm_stb import process as atm_stb
@@ -24,6 +25,7 @@
 from .pass_filter import high_pass
 from .pass_filter import band_pass
 from .mooring_fatigue import process as mooring_fatigue
+from .rolling_median import process as rolling_median
 
 __all__ = [
     'atm_stb',
@@ -34,4 +36,5 @@
     'high_pass',
     'band_pass',
     'mooring_fatigue',
+    'rolling_median'
 ]

camille/process/rolling_median.py

@@ -0,0 +1,41 @@
+import pandas as pd
+import numpy as np
+
+def process(signal, wsize, tolerance=None):
+    """ Process rolling median
+
+    Rolling window calculations
+
+    Parameters
+    ----------
+    signal : pd.Series
+    wsize : int or pd.Timedelta
+        size of the rolling window
+    tolerance : float, optional
+        threshold for outliers, default to the standard deviation of signal
+
+    Returns
+    -------
+    pd.Series
+        filtered signal
+
+    Examples
+    --------
+    >>> s = pd.Series(signal, index=t)
+    >>> processed = process.rolling_median(s, wsize=20, tolerance=2.0)
+    """
+
+    dt = signal.index[1] - signal.index[0]
+    if isinstance(wsize, pd.Timedelta):
+        if wsize < dt:
+            problem = "wsize is smaller than the samplerate of signal: {} < {}".format(wsize, dt)
+            raise ValueError(problem)
+
+        wsize = int(np.ceil(wsize / dt))
+
+    if tolerance == None: tolerance = 1.0 * signal.std()
+
+    rolling = signal.rolling(window=wsize, min_periods=1, center=True).median()
+    outliers = signal[(signal - rolling).abs() >= tolerance]
+    cleaned = signal.drop(outliers.index).iloc[wsize:-wsize]
+    return cleaned

tests/process/test_rolling_median.py

@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+import numpy as np
+import pandas as pd
+
+from camille import process
+
+def test_process():
+    index = pd.date_range('1/1/2018', periods=100, freq='S')
+
+    t = np.linspace(-np.pi, np.pi, num=100)
+    signal = np.sin(t)
+    signal[40] *= 5   #create an outlier
+
+    s = pd.Series(signal, index=index)
+    dt = pd.Timedelta(seconds=5)
+
+    processed = process.rolling_median(s, wsize=dt)
+
+    expected = s.drop(s.index[40]).iloc[5:-5]
+    pd.testing.assert_series_equal(expected, processed)