Add resample and decimate utility methods for Analog and IrregularlysampledSignals (#776)

* Add decimate and resampling method for analog and irregular signal

* Add decimate and resampling method for analog and irregular signal

* Refactor method names and imports and add requirements to unittests

* Refactor method names and imports and add requirements to unittests

* Remove non-implemented resampling method and refactor

* Refactor for PEP8

* Refactor for PEP8

* Improve docstrings and add crossreferences.

Author: JuliaSprenger

neo/core/analogsignal.py

@@ -23,6 +23,13 @@
 
 import logging
 
+try:
+    import scipy.signal
+except ImportError as err:
+    HAVE_SCIPY = False
+else:
+    HAVE_SCIPY = True
+
 import numpy as np
 import quantities as pq
 
@@ -539,3 +546,92 @@ def splice(self, signal, copy=False):
         else:
             self[i:j, :] = signal
             return self
+
+    def downsample(self, downsampling_factor, **kwargs):
+        """
+        Downsample the data of a signal.
+        This method reduces the number of samples of the AnalogSignal to a fraction of the
+        original number of samples, defined by `downsampling_factor`.
+        This method is a wrapper of scipy.signal.decimate and accepts the same set of keyword
+        arguments, except for specifying the axis of resampling, which is fixed to the first axis
+        here.
+
+        Parameters:
+        -----------
+        downsampling_factor: integer
+            Factor used for decimation of samples. Scipy recommends to call decimate multiple times
+            for downsampling factors higher than 13 when using IIR downsampling (default).
+
+        Returns:
+        --------
+        downsampled_signal: :class:`AnalogSignal`
+            New instance of a :class:`AnalogSignal` object containing the resampled data points.
+            The original :class:`AnalogSignal` is not modified.
+
+        Note:
+        -----
+        For resampling the signal with a fixed number of samples, see `resample` method.
+        """
+
+        if not HAVE_SCIPY:
+            raise ImportError('Decimating requires availability of scipy.signal')
+
+        # Resampling is only permitted along the time axis (axis=0)
+        if 'axis' in kwargs:
+            kwargs.pop('axis')
+
+        downsampled_data = scipy.signal.decimate(self.magnitude, downsampling_factor, axis=0,
+                                                 **kwargs)
+        downsampled_signal = self.duplicate_with_new_data(downsampled_data)
+
+        # since the number of channels stays the same, we can also copy array annotations here
+        downsampled_signal.array_annotations = self.array_annotations.copy()
+        downsampled_signal.sampling_rate = self.sampling_rate / downsampling_factor
+
+        return downsampled_signal
+
+    def resample(self, sample_count, **kwargs):
+        """
+        Resample the data points of the signal.
+        This method interpolates the signal and returns a new signal with a fixed number of
+        samples defined by `sample_count`.
+        This method is a wrapper of scipy.signal.resample and accepts the same set of keyword
+        arguments, except for specifying the axis of resampling which is fixed to the first axis
+        here, and the sample positions. .
+
+        Parameters:
+        -----------
+        sample_count: integer
+            Number of desired samples. The resulting signal starts at the same sample as the
+            original and is sampled regularly.
+
+        Returns:
+        --------
+        resampled_signal: :class:`AnalogSignal`
+            New instance of a :class:`AnalogSignal` object containing the resampled data points.
+            The original :class:`AnalogSignal` is not modified.
+
+        Note:
+        -----
+        For reducing the number of samples to a fraction of the original, see `downsample` method
+        """
+
+        if not HAVE_SCIPY:
+            raise ImportError('Resampling requires availability of scipy.signal')
+
+        # Resampling is only permitted along the time axis (axis=0)
+        if 'axis' in kwargs:
+            kwargs.pop('axis')
+        if 't' in kwargs:
+            kwargs.pop('t')
+
+        resampled_data, resampled_times = scipy.signal.resample(self.magnitude, sample_count,
+                                                                t=self.times, axis=0, **kwargs)
+
+        resampled_signal = self.duplicate_with_new_data(resampled_data)
+        resampled_signal.sampling_rate = (sample_count / self.shape[0]) * self.sampling_rate
+
+        # since the number of channels stays the same, we can also copy array annotations here
+        resampled_signal.array_annotations = self.array_annotations.copy()
+
+        return resampled_signal

neo/core/irregularlysampledsignal.py

@@ -24,11 +24,20 @@
 from __future__ import absolute_import, division, print_function
 
 from copy import deepcopy, copy
+
+try:
+    import scipy.signal
+except ImportError as err:
+    HAVE_SCIPY = False
+else:
+    HAVE_SCIPY = True
+
 import numpy as np
 import quantities as pq
 
 from neo.core.baseneo import BaseNeo, MergeError, merge_annotations
 from neo.core.basesignal import BaseSignal
+from neo.core.analogsignal import AnalogSignal
 from neo.core.channelindex import ChannelIndex
 from neo.core.dataobject import DataObject
 
@@ -346,21 +355,51 @@ def mean(self, interpolation=None):
         else:
             raise NotImplementedError
 
-    def resample(self, at=None, interpolation=None):
-        '''
-        Resample the signal, returning either an :class:`AnalogSignal` object
-        or another :class:`IrregularlySampledSignal` object.
+    def resample(self, sample_count, **kwargs):
+        """
+        Resample the data points of the signal.
+        This method interpolates the signal and returns a new signal with a fixed number of
+        samples defined by `sample_count`.
+        This function is a wrapper of scipy.signal.resample and accepts the same set of keyword
+        arguments, except for specifying the axis of resampling which is fixed to the first axis
+        here, and the sample positions. .
 
-        Arguments:
-            :at: either a :class:`Quantity` array containing the times at
-                 which samples should be created (times must be within the
-                 signal duration, there is no extrapolation), a sampling rate
-                 with dimensions (1/Time) or a sampling interval
-                 with dimensions (Time).
-            :interpolation: one of: None, 'linear'
-        '''
-        # further interpolation methods could be added
-        raise NotImplementedError
+        Parameters:
+        -----------
+        sample_count: integer
+            Number of desired samples. The resulting signal starts at the same sample as the
+            original and is sampled regularly.
+
+        Returns:
+        --------
+        resampled_signal: :class:`AnalogSignal`
+            New instance of a :class:`AnalogSignal` object containing the resampled data points.
+            The original :class:`AnalogSignal` is not modified.
+        """
+
+        if not HAVE_SCIPY:
+            raise ImportError('Resampling requires availability of scipy.signal')
+
+        # Resampling is only permitted along the time axis (axis=0)
+        if 'axis' in kwargs:
+            kwargs.pop('axis')
+        if 't' in kwargs:
+            kwargs.pop('t')
+
+        resampled_data, resampled_times = scipy.signal.resample(self.magnitude, sample_count,
+                                                                t=self.times.magnitude,
+                                                                axis=0, **kwargs)
+
+        new_sampling_rate = (sample_count - 1) / self.duration
+        resampled_signal = AnalogSignal(resampled_data, units=self.units, dtype=self.dtype,
+                                        t_start=self.t_start,
+                                        sampling_rate=new_sampling_rate,
+                                        array_annotations=self.array_annotations.copy(),
+                                        **self.annotations.copy())
+
+        # since the number of channels stays the same, we can also copy array annotations here
+        resampled_signal.array_annotations = self.array_annotations.copy()
+        return resampled_signal
 
     def time_slice(self, t_start, t_stop):
         '''

neo/test/coretest/test_analogsignal.py

@@ -24,6 +24,13 @@
 else:
     HAVE_IPYTHON = True
 
+try:
+    import scipy
+except ImportError:
+    HAVE_SCIPY = False
+else:
+    HAVE_SCIPY = True
+
 from numpy.testing import assert_array_equal
 from neo.core.analogsignal import AnalogSignal, _get_sampling_rate
 from neo.core.channelindex import ChannelIndex
@@ -1158,6 +1165,87 @@ def test_array_annotations_getitem(self):
         assert_arrays_equal(result3.array_annotations['label'], np.array(['abc', 'def']))
         self.assertIsInstance(result3.array_annotations, ArrayDict)
 
+    @unittest.skipUnless(HAVE_SCIPY, "requires Scipy")
+    def test_downsample(self):
+        # generate signal long enough for decimating
+        data = np.sin(np.arange(1500) / 30).reshape(500, 3) * pq.mV
+        signal = AnalogSignal(data, sampling_rate=30000 * pq.Hz)
+
+        # test decimation using different decimation factors
+        factors = [1, 9, 10, 11]
+        for factor in factors:
+            desired = signal[::factor].magnitude
+            result = signal.downsample(factor)
+
+            self.assertEqual(np.ceil(signal.shape[0] / factor), result.shape[0])
+            self.assertEqual(signal.shape[-1],
+                             result.shape[-1])  # preserve number of recording traces
+            self.assertAlmostEqual(signal.sampling_rate, factor * result.sampling_rate)
+            # only comparing center values due to border effects
+            np.testing.assert_allclose(desired[3:-3], result.magnitude[3:-3], rtol=0.05, atol=0.1)
+
+    @unittest.skipUnless(HAVE_SCIPY, "requires Scipy")
+    def test_resample_less_samples(self):
+        # generate signal long enough for resampling
+        data = np.sin(np.arange(1500) / 30).reshape(3, 500).T * pq.mV
+        signal = AnalogSignal(data, sampling_rate=30000 * pq.Hz)
+
+        # test resampling using different numbers of desired samples
+        sample_counts = [10, 100, 400]
+        for sample_count in sample_counts:
+            sample_ids = np.linspace(0, signal.shape[0], sample_count, dtype=int, endpoint=False)
+            desired = signal.magnitude[sample_ids]
+            result = signal.resample(sample_count)
+
+            self.assertEqual(sample_count, result.shape[0])
+            self.assertEqual(signal.shape[-1],
+                             result.shape[-1])  # preserve number of recording traces
+            self.assertAlmostEqual(sample_count / signal.shape[0] * signal.sampling_rate,
+                                   result.sampling_rate)
+            # only comparing center values due to border effects
+            np.testing.assert_allclose(desired[3:-3], result.magnitude[3:-3], rtol=0.05, atol=0.1)
+
+    @unittest.skipUnless(HAVE_SCIPY, "requires Scipy")
+    def test_resample_more_samples(self):
+        # generate signal long enough for resampling
+        data = np.sin(np.arange(1500) / 100).T * pq.mV
+        signal = AnalogSignal(data, sampling_rate=30000 * pq.Hz)
+
+        # test resampling using different numbers of desired samples
+        factor = 2
+        sample_count = factor * signal.shape[0]
+        desired = np.interp(np.arange(sample_count) / factor, np.arange(signal.shape[0]),
+                            signal.magnitude.flatten()).reshape(-1, 1)
+        result = signal.resample(sample_count)
+
+        self.assertEqual(sample_count, result.shape[0])
+        self.assertEqual(signal.shape[-1], result.shape[-1])  # preserve number of recording traces
+        self.assertAlmostEqual(sample_count / signal.shape[0] * signal.sampling_rate,
+                               result.sampling_rate)
+        # only comparing center values due to border effects
+        np.testing.assert_allclose(desired[10:-10], result.magnitude[10:-10], rtol=0.0, atol=0.1)
+
+    @unittest.skipUnless(HAVE_SCIPY, "requires Scipy")
+    def test_compare_resample_and_downsample(self):
+        # generate signal long enough for resampling
+        data = np.sin(np.arange(1500) / 30).reshape(3, 500).T * pq.mV
+        signal = AnalogSignal(data, sampling_rate=30000 * pq.Hz)
+
+        # test resampling using different numbers of desired samples
+        sample_counts = [10, 100, 250]
+        for sample_count in sample_counts:
+            downsampling_factor = int(signal.shape[0] / sample_count)
+            desired = signal.downsample(downsampling_factor=downsampling_factor)
+            result = signal.resample(sample_count)
+
+            self.assertEqual(desired.shape[0], result.shape[0])
+            self.assertEqual(desired.shape[-1],
+                             result.shape[-1])  # preserve number of recording traces
+            self.assertAlmostEqual(desired.sampling_rate, result.sampling_rate)
+            # only comparing center values due to border effects
+            np.testing.assert_allclose(desired.magnitude[3:-3], result.magnitude[3:-3], rtol=0.05,
+                                       atol=0.1)
+
 
 class TestAnalogSignalEquality(unittest.TestCase):
     def test__signals_with_different_data_complement_should_be_not_equal(self):

neo/test/coretest/test_irregularysampledsignal.py

@@ -23,6 +23,13 @@
 else:
     HAVE_IPYTHON = True
 
+try:
+    import scipy
+except ImportError:
+    HAVE_SCIPY = False
+else:
+    HAVE_SCIPY = True
+
 from neo.core.irregularlysampledsignal import IrregularlySampledSignal
 from neo.core import Segment, ChannelIndex
 from neo.core.baseneo import MergeError
@@ -352,9 +359,6 @@ def test_simple_statistics(self):
     def test_mean_interpolation_NotImplementedError(self):
         self.assertRaises(NotImplementedError, self.signal1.mean, True)
 
-    def test_resample_NotImplementedError(self):
-        self.assertRaises(NotImplementedError, self.signal1.resample, True)
-
     def test__rescale_same(self):
         result = self.signal1.copy()
         result = result.rescale(pq.mV)
@@ -710,6 +714,14 @@ def test__copy_should_preserve_parent_objects(self):
         self.assertIs(result.segment, self.signal1.segment)
         self.assertIs(result.channel_index, self.signal1.channel_index)
 
+    @unittest.skipUnless(HAVE_SCIPY, "requires Scipy")
+    def test_resample(self):
+        factors = [1, 2, 10]
+        for factor in factors:
+            result = self.signal1.resample(self.signal1.shape[0] * factor)
+            np.testing.assert_allclose(self.signal1.magnitude, result.magnitude[::factor],
+                                       rtol=1e-7, atol=0)
+
 
 class TestIrregularlySampledSignalCombination(unittest.TestCase):
     def setUp(self):