add reduce_ringing() function

Author: nbara

examples/example_dering.py

@@ -0,0 +1,42 @@
+"""
+Ringing artifact reduction example
+==================================
+
+This example shows how to subtract the impulse response from a filter to
+reduce ringing artifacts.
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+from scipy.signal import butter, lfilter
+
+from meegkit.detrend import reduce_ringing
+
+# import config  # plotting utils
+
+np.random.seed(9)
+
+###############################################################################
+# Detrending
+# =============================================================================
+
+###############################################################################
+# Basic example with a linear trend
+# -----------------------------------------------------------------------------
+# Simulate the effect of filtering a signal containing a discontinuity, and try
+# to remove the resulting ringing artifact by subtracing the opposite of the
+# impulse response.
+
+x = np.arange(1000) < 1
+[b, a] = butter(6, 0.2)     # Butterworth filter design
+x = lfilter(b, a, x) * 50    # Filter data using above filter
+x = np.roll(x, 500)
+x = x[:, None] + np.random.randn(1000, 2)
+y = reduce_ringing(x, samples=np.array([500]))
+
+plt.figure()
+plt.plot(x + np.array([-10, 10]), 'C0', label='before')
+plt.plot(y + np.array([-10, 10]), 'C1:', label='after')
+plt.legend()
+plt.show()

meegkit/__init__.py

@@ -1,3 +1,5 @@
 """M/EEG denoising utilities in python."""
 from . import utils
-from . import dss, sns, tspca, star, cca
+from . import dss, sns, tspca, star, cca, detrend
+
+__all__ = ['cca', 'detrend', 'dss', 'sns', 'star']

meegkit/cca.py

@@ -263,10 +263,6 @@ def nt_cca(X=None, Y=None, lags=None, C=None, m=None, thresh=1e-12, sfreq=1):
     .. warning:: A, B are scaled so that (X * A)^2 and (Y * B)^2 are identity
                  matrices (differs from sklearn).
 
-    See Also
-    --------
-    nt_cov_lags, nt_relshift, nt_cov, nt_pca in NoiseTools.
-
     """
     if (X is None and Y is not None) or (Y is None and X is not None):
         raise AttributeError('Either *both* X and Y should be defined, or C!')

meegkit/detrend.py

@@ -1,8 +1,11 @@
 """Robust detrending."""
 import numpy as np
 
+from scipy.signal import lfilter
+
 from .utils import demean, mrdivide, pca, unfold
 from .utils.matrix import _check_weights
+from .utils.sig import stmcb
 
 
 def detrend(x, order, w=None, basis='polynomials', threshold=3, n_iter=4,
@@ -209,6 +212,66 @@ def regress(x, r, w=None, threshold=1e-7, return_mean=False):
     return b, z
 
 
+def reduce_ringing(X, samples, order=10, n_samples=100, extra=50, threshold=3,
+                   show=False):
+    """Subtract filter impulse response from signal at given samples.
+
+    Parameters
+    ----------
+    X: ndarray, shape=(n_times, n_chans[, n_trials])
+        Data containing ringing artifacts.
+    samples : list of ints
+        Sample indices where to find ringing artifacts.
+    order : int
+        Order of polynomial trend (default=10).
+    n_samples = 100
+        Number of samples over which to estimate impulse response
+        (default=100).
+    extra : int
+        Samples before stimulus to anchor trend (default=50).
+    threshold: float
+        Threshold for robust detrending (default=3).
+
+    Returns
+    -------
+    y : ndarray, shape=(n_times, n_chans[, n_trials])
+        Clean data.
+
+    """
+    NNUM = 8
+    NDEN = 8  # number of filter coeffs
+
+    # remove samples too close to beginning or end
+    samples = samples[samples > extra]
+    samples = samples[samples < X.shape[0] - n_samples]
+
+    y = X.copy()
+    for i, s in enumerate(samples):
+        for c in range(X.shape[1]):
+            # select portion to fit filter response, remove polynomial trend
+            response = X[s - extra:s + n_samples, c]
+            # response = detrend(response, order, threshold)
+            response = response[extra:]
+
+            # estimate filter parameters - helps ensure stable filter
+            response = np.r_[(response, np.zeros(response.shape))]
+            [B, A] = stmcb(response, q=NNUM, p=NDEN, niter=20)
+
+            # estimate filter response to event
+            pulse = np.arange(n_samples) < 1
+            model = lfilter(B, A, pulse)
+            idx = s + np.arange(model.shape[0])
+            y[idx, c] = X[idx, c] - model
+
+    if show:
+        w = np.zeros((X.shape[0], X.shape[1]))
+        for s in samples:
+            w[s:s+n_samples, :] = 1
+        _plot_detrend(X, y, w)
+
+    return y
+
+
 def _plot_detrend(x, y, w):
     """Plot detrending results."""
     import matplotlib.pyplot as plt
@@ -228,8 +291,9 @@ def _plot_detrend(x, y, w):
     ax1.legend()
 
     ax2 = f.add_subplot(gs[3, 0])
-    ax2.imshow(w.T, aspect='auto', cmap='Greys')
-    ax2.set_yticks(np.arange(1, n_chans + 1, 1))
+    ax2.pcolormesh(w.T, cmap='Greys')
+    ax2.set_yticks(np.arange(0, n_chans) + 0.5)
+    ax2.set_yticklabels(['ch{}'.format(i) for i in np.arange(n_chans)])
     ax2.set_xlim(0, n_times)
     ax2.set_ylabel('ch. weights')
     ax2.set_xlabel('samples')

meegkit/dss.py

@@ -1,3 +1,4 @@
+"""Denoising source separation."""
 import numpy as np
 from scipy import linalg
 

meegkit/sns.py

@@ -1,3 +1,4 @@
+"""Sensor noise suppression."""
 import numpy as np
 
 from .tspca import tsr

meegkit/star.py

@@ -1,3 +1,4 @@
+"""Sparse time-artefact removal."""
 import numpy as np
 from scipy.signal import filtfilt
 

meegkit/tspca.py

@@ -1,3 +1,4 @@
+"""Time-shift PCA."""
 import numpy as np
 
 from .utils import (demean, fold, multishift, normcol, pca, regcov, tscov,

meegkit/utils/denoise.py

@@ -1,9 +1,7 @@
-from __future__ import division, print_function
-
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib import gridspec
-from scipy import linalg
+from scipy import linalg, signal
 
 from .matrix import demean, fold, theshapeof, unfold
 

meegkit/utils/sig.py

@@ -603,19 +603,14 @@ def prony(h, nb, na):
     If the largest order specified is greater than the length of H, H is padded
     with zeros.
 
-    Examples
-    --------
-    Fit an IIR model to an impulse response of a lowpass filter.
-
-    >>> [b,a] = butter(4,0.2);
-    >>> impulseResp = impz(b,a);                % obtain impulse response
-    >>> denOrder=4; numOrder=4;                 % system function of order 4
-    >>> [Num,Den]=prony(impulseResp,numOrder,denOrder);
-    >>> subplot(211);                           % impulse response and input
-    >>> stem(impz(Num,Den,length(impulseResp)));
-    >>> title('Impulse Response with Prony Design');
-    >>> subplot(212);
-    >>> stem(impulseResp); title('Input Impulse Response');
+    Parameters
+    ----------
+    h : array
+        Impulse response.
+    nb : int
+        Numerator order.
+    na : int
+        Denominator order.
 
     References
     ----------