MRG: Faster raw resampling

examples/preprocessing/plot_resample.py

@@ -37,7 +37,7 @@
 
 # Downsample to 100 Hz
 print('Original sampling rate:', epochs.info['sfreq'], 'Hz')
-epochs_resampled = epochs.resample(100, copy=True)
+epochs_resampled = epochs.resample(100, npad='auto', copy=True)
 print('New sampling rate:', epochs_resampled.info['sfreq'], 'Hz')
 
 # Plot a piece of data to see the effects of downsampling
@@ -65,7 +65,7 @@
 # can also be done on non-preloaded data.
 
 # Resample to 300 Hz
-raw_resampled = raw.resample(300, copy=True)
+raw_resampled = raw.resample(300, npad='auto', copy=True)
 
 ###############################################################################
 # Because resampling also affects the stim channels, some trigger onsets might
@@ -75,11 +75,12 @@
 print('Number of events before resampling:', len(mne.find_events(raw)))
 
 # Resample to 100 Hz (generates warning)
-raw_resampled = raw.resample(100, copy=True)
+raw_resampled = raw.resample(100, npad='auto', copy=True)
 print('Number of events after resampling:',
       len(mne.find_events(raw_resampled)))
 
 # To avoid losing events, jointly resample the data and event matrix
 events = mne.find_events(raw)
-raw_resampled, events_resampled = raw.resample(100, events=events, copy=True)
+raw_resampled, events_resampled = raw.resample(100, npad='auto',
+                                               events=events, copy=True)
 print('Number of events after resampling:', len(events_resampled))

mne/cuda.py

@@ -3,7 +3,7 @@
 # License: BSD (3-clause)
 
 import numpy as np
-from scipy.fftpack import fft, ifft
+from scipy.fftpack import fft, ifft, rfft, irfft
 
 from .utils import sizeof_fmt, logger, get_config, warn
 
@@ -295,7 +295,7 @@ def setup_cuda_fft_resample(n_jobs, W, new_len):
     return n_jobs, cuda_dict, W
 
 
-def fft_resample(x, W, new_len, npad, to_remove,
+def fft_resample(x, W, new_len, npads, to_removes,
                  cuda_dict=dict(use_cuda=False)):
     """Do FFT resampling with a filter function (possibly using CUDA)
 
@@ -307,9 +307,10 @@ def fft_resample(x, W, new_len, npad, to_remove,
         The filtering function to apply.
     new_len : int
         The size of the output array (before removing padding).
-    npad : int
-        Amount of padding to apply before resampling.
-    to_remove : int
+    npads : tuple of int
+        Amount of padding to apply to the start and end of the
+        signal before resampling.
+    to_removes : tuple of int
         Number of samples to remove after resampling.
     cuda_dict : dict
         Dictionary constructed using setup_cuda_multiply_repeated().
@@ -322,18 +323,28 @@ def fft_resample(x, W, new_len, npad, to_remove,
     # add some padding at beginning and end to make this work a little cleaner
     if x.dtype != np.float64:
         x = x.astype(np.float64)
-    x = _smart_pad(x, npad)
+    x = _smart_pad(x, npads)
     old_len = len(x)
     shorter = new_len < old_len
     if not cuda_dict['use_cuda']:
         N = int(min(new_len, old_len))
-        sl_1 = slice((N + 1) // 2)
-        y_fft = np.zeros(new_len, np.complex128)
-        x_fft = fft(x).ravel() * W
+        # The below is equivalent to this, but faster
+        # sl_1 = slice((N + 1) // 2)
+        # y_fft = np.zeros(new_len, np.complex128)
+        # x_fft = fft(x).ravel() * W
+        # y_fft[sl_1] = x_fft[sl_1]
+        # sl_2 = slice(-(N - 1) // 2, None)
+        # y_fft[sl_2] = x_fft[sl_2]
+        # y = np.real(ifft(y_fft, overwrite_x=True)).ravel()
+        x_fft = rfft(x).ravel()
+        x_fft *= W[np.arange(1, len(x) + 1) // 2].real
+        y_fft = np.zeros(new_len, np.float64)
+        sl_1 = slice(N)
         y_fft[sl_1] = x_fft[sl_1]
-        sl_2 = slice(-(N - 1) // 2, None)
-        y_fft[sl_2] = x_fft[sl_2]
-        y = np.real(ifft(y_fft, overwrite_x=True)).ravel()
+        if min(new_len, old_len) % 2 == 0:
+            if new_len > old_len:
+                y_fft[N - 1] /= 2.
+        y = irfft(y_fft, overwrite_x=True).ravel()
     else:
         cudafft = _get_cudafft()
         cuda_dict['x'].set(np.concatenate((x, np.zeros(max(new_len - old_len,
@@ -356,10 +367,10 @@ def fft_resample(x, W, new_len, npad, to_remove,
         y = cuda_dict['x'].get()[:new_len if shorter else None]
 
     # now let's trim it back to the correct size (if there was padding)
-    if to_remove > 0:
+    if (to_removes > 0).any():
         keep = np.ones((new_len), dtype='bool')
-        keep[:to_remove] = False
-        keep[-to_remove:] = False
+        keep[:to_removes[0]] = False
+        keep[-to_removes[1]:] = False
         y = np.compress(keep, y)
 
     return y
@@ -372,11 +383,12 @@ def fft_resample(x, W, new_len, npad, to_remove,
 def _smart_pad(x, n_pad):
     """Pad vector x
     """
-    if n_pad == 0:
+    if (n_pad == 0).all():
         return x
-    elif n_pad < 0:
+    elif (n_pad < 0).any():
         raise RuntimeError('n_pad must be non-negative')
     # need to pad with zeros if len(x) <= npad
-    z_pad = np.zeros(max(n_pad - len(x) + 1, 0), dtype=x.dtype)
-    return np.concatenate([z_pad, 2 * x[0] - x[n_pad:0:-1], x,
-                           2 * x[-1] - x[-2:-n_pad - 2:-1], z_pad])
+    l_z_pad = np.zeros(max(n_pad[0] - len(x) + 1, 0), dtype=x.dtype)
+    r_z_pad = np.zeros(max(n_pad[0] - len(x) + 1, 0), dtype=x.dtype)
+    return np.concatenate([l_z_pad, 2 * x[0] - x[n_pad[0]:0:-1], x,
+                           2 * x[-1] - x[-2:-n_pad[1] - 2:-1], r_z_pad])

mne/epochs.py

@@ -1438,16 +1438,18 @@ def crop(self, tmin=None, tmax=None, copy=False):
         return this_epochs
 
     @verbose
-    def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1,
+    def resample(self, sfreq, npad=None, window='boxcar', n_jobs=1,
                  copy=False, verbose=None):
         """Resample preloaded data
 
         Parameters
         ----------
         sfreq : float
             New sample rate to use
-        npad : int
+        npad : int | str
             Amount to pad the start and end of the data.
+            Can also be "auto" to use a padding that will result in
+            a power-of-two size (can be much faster).
         window : string or tuple
             Window to use in resampling. See scipy.signal.resample.
         n_jobs : int
@@ -1472,6 +1474,11 @@ def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1,
         # XXX this could operate on non-preloaded data, too
         if not self.preload:
             raise RuntimeError('Can only resample preloaded data')
+        if npad is None:
+            npad = 100
+            warn('npad is currently taken to be 100, but will be changed to '
+                 '"auto" in 0.12. Please set the value explicitly.',
+                 DeprecationWarning)
 
         inst = self.copy() if copy else self
 

mne/evoked.py

@@ -705,7 +705,7 @@ def as_type(self, ch_type='grad', mode='fast'):
         from .forward import _as_meg_type_evoked
         return _as_meg_type_evoked(self, ch_type=ch_type, mode=mode)
 
-    def resample(self, sfreq, npad=100, window='boxcar'):
+    def resample(self, sfreq, npad=None, window='boxcar'):
         """Resample data
 
         This function operates in-place.
@@ -714,11 +714,18 @@ def resample(self, sfreq, npad=100, window='boxcar'):
         ----------
         sfreq : float
             New sample rate to use
-        npad : int
+        npad : int | str
             Amount to pad the start and end of the data.
+            Can also be "auto" to use a padding that will result in
+            a power-of-two size (can be much faster).
         window : string or tuple
             Window to use in resampling. See scipy.signal.resample.
         """
+        if npad is None:
+            npad = 100
+            warn('npad is currently taken to be 100, but will be changed to '
+                 '"auto" in 0.12. Please set the value explicitly.',
+                 DeprecationWarning)
         sfreq = float(sfreq)
         o_sfreq = self.info['sfreq']
         self.data = resample(self.data, sfreq, o_sfreq, npad, -1, window)

mne/filter.py

@@ -159,7 +159,7 @@ def _1d_overlap_filter(x, h_fft, n_h, n_edge, zero_phase, cuda_dict):
         n_fft = cuda_dict['x'].size  # account for CUDA's modification of h_fft
     else:
         n_fft = len(h_fft)
-    x_ext = _smart_pad(x, n_edge)
+    x_ext = _smart_pad(x, np.array([n_edge, n_edge]))
     n_x = len(x_ext)
     x_filtered = np.zeros_like(x_ext)
 
@@ -1265,8 +1265,9 @@ def resample(x, up, down, npad=100, axis=-1, window='boxcar', n_jobs=1,
         Factor to upsample by.
     down : float
         Factor to downsample by.
-    npad : integer
+    npad : int | str
         Number of samples to use at the beginning and end for padding.
+        Can be "auto" to pad to the next highest power of 2.
     axis : int
         Axis along which to resample (default is the last axis).
     window : string or tuple
@@ -1317,12 +1318,31 @@ def resample(x, up, down, npad=100, axis=-1, window='boxcar', n_jobs=1,
     if x_len == 0:
         warn('x has zero length along last axis, returning a copy of x')
         return x.copy()
+    bad_msg = 'npad must be "auto" or an integer'
+    if isinstance(npad, string_types):
+        if npad != 'auto':
+            raise ValueError(bad_msg)
+        # Figure out reasonable pad that gets us to a power of 2
+        min_add = min(x_len // 8, 100) * 2
+        npad = 2 ** int(np.ceil(np.log2(x_len + min_add))) - x_len
+        npad, extra = divmod(npad, 2)
+        npads = np.array([npad, npad + extra], int)
+    else:
+        if npad != int(npad):
+            raise ValueError(bad_msg)
+        npads = np.array([npad, npad], int)
+    del npad
 
     # prep for resampling now
     x_flat = x.reshape((-1, x_len))
-    orig_len = x_len + 2 * npad  # length after padding
+    orig_len = x_len + npads.sum()  # length after padding
     new_len = int(round(ratio * orig_len))  # length after resampling
-    to_remove = np.round(ratio * npad).astype(int)
+    final_len = int(round(ratio * x_len))
+    to_removes = [int(round(ratio * npads[0]))]
+    to_removes.append(new_len - final_len - to_removes[0])
+    to_removes = np.array(to_removes)
+    # This should hold:
+    # assert np.abs(to_removes[1] - to_removes[0]) <= int(np.ceil(ratio))
 
     # figure out windowing function
     if window is not None:
@@ -1344,13 +1364,13 @@ def resample(x, up, down, npad=100, axis=-1, window='boxcar', n_jobs=1,
     # do the resampling using an adaptation of scipy's FFT-based resample()
     # use of the 'flat' window is recommended for minimal ringing
     if n_jobs == 1:
-        y = np.zeros((len(x_flat), new_len - 2 * to_remove), dtype=x.dtype)
+        y = np.zeros((len(x_flat), new_len - to_removes.sum()), dtype=x.dtype)
         for xi, x_ in enumerate(x_flat):
-            y[xi] = fft_resample(x_, W, new_len, npad, to_remove,
+            y[xi] = fft_resample(x_, W, new_len, npads, to_removes,
                                  cuda_dict)
     else:
         parallel, p_fun, _ = parallel_func(fft_resample, n_jobs)
-        y = parallel(p_fun(x_, W, new_len, npad, to_remove, cuda_dict)
+        y = parallel(p_fun(x_, W, new_len, npads, to_removes, cuda_dict)
                      for x_ in x_flat)
         y = np.array(y)
 

mne/inverse_sparse/tests/test_gamma_map.py

@@ -44,7 +44,7 @@ def test_gamma_map():
     forward = pick_types_forward(forward, meg=False, eeg=True)
     evoked = read_evokeds(fname_evoked, condition=0, baseline=(None, 0),
                           proj=False)
-    evoked.resample(50)
+    evoked.resample(50, npad=100)
     evoked.crop(tmin=0.1, tmax=0.16)  # crop to nice window near samp border
 
     cov = read_cov(fname_cov)

mne/io/base.py

@@ -945,7 +945,7 @@ def notch_filter(self, freqs, picks=None, filter_length='10s',
                                   picks=picks, n_jobs=n_jobs, copy=False)
 
     @verbose
-    def resample(self, sfreq, npad=100, window='boxcar', stim_picks=None,
+    def resample(self, sfreq, npad=None, window='boxcar', stim_picks=None,
                  n_jobs=1, events=None, copy=False, verbose=None):
         """Resample data channels.
 
@@ -973,8 +973,10 @@ def resample(self, sfreq, npad=100, window='boxcar', stim_picks=None,
         ----------
         sfreq : float
             New sample rate to use.
-        npad : int
+        npad : int | str
             Amount to pad the start and end of the data.
+            Can also be "auto" to use a padding that will result in
+            a power-of-two size (can be much faster).
         window : string or tuple
             Window to use in resampling. See scipy.signal.resample.
         stim_picks : array of int | None
@@ -1006,6 +1008,11 @@ def resample(self, sfreq, npad=100, window='boxcar', stim_picks=None,
         For some data, it may be more accurate to use npad=0 to reduce
         artifacts. This is dataset dependent -- check your data!
         """
+        if npad is None:
+            npad = 100
+            warn('npad is currently taken to be 100, but will be changed to '
+                 '"auto" in 0.12. Please set the value explicitly.',
+                 DeprecationWarning)
         if not self.preload:
             raise RuntimeError('Can only resample preloaded data')