From 41705eaa320d12e619224fac2e13f0ed014541e0 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 2 Mar 2016 11:45:20 -0500 Subject: [PATCH 1/3] ENH: Faster raw resampling --- mne/cuda.py | 52 ++++++++++++++++++------------ mne/filter.py | 33 +++++++++++++++---- mne/io/base.py | 10 ++++-- mne/io/fiff/tests/test_raw_fiff.py | 30 ++++++++--------- mne/tests/test_filter.py | 17 ++++++---- 5 files changed, 91 insertions(+), 51 deletions(-) diff --git a/mne/cuda.py b/mne/cuda.py index 130f8bc38e5..02ae6263dde 100644 --- a/mne/cuda.py +++ b/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]) diff --git a/mne/filter.py b/mne/filter.py index 93a054797a8..b0000d97702 100644 --- a/mne/filter.py +++ b/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,30 @@ 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 pad that gets us to a power of 2 with at least + # 100 on each end + npad = 2 ** int(np.ceil(np.log2(x_len + 200))) - x_len + npad, extra = divmod(npad, 2) + npads = np.array([npad, npad + extra]) # XXX + 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) + assert np.abs(to_removes[1] - to_removes[0]) <= int(np.ceil(ratio)) # figure out windowing function if window is not None: @@ -1344,13 +1363,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) diff --git a/mne/io/base.py b/mne/io/base.py index d39fac51c66..4a5d6c6c29d 100644 --- a/mne/io/base.py +++ b/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,10 @@ 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.') if not self.preload: raise RuntimeError('Can only resample preloaded data') diff --git a/mne/io/fiff/tests/test_raw_fiff.py b/mne/io/fiff/tests/test_raw_fiff.py index 951bf58bd42..d0496fdd825 100644 --- a/mne/io/fiff/tests/test_raw_fiff.py +++ b/mne/io/fiff/tests/test_raw_fiff.py @@ -860,7 +860,7 @@ def test_resample(): raw_resamp = raw.copy() sfreq = raw.info['sfreq'] # test parallel on upsample - raw_resamp.resample(sfreq * 2, n_jobs=2) + raw_resamp.resample(sfreq * 2, n_jobs=2, npad='auto') assert_equal(raw_resamp.n_times, len(raw_resamp.times)) raw_resamp.save(op.join(tempdir, 'raw_resamp-raw.fif')) raw_resamp = Raw(op.join(tempdir, 'raw_resamp-raw.fif'), preload=True) @@ -869,7 +869,7 @@ def test_resample(): assert_equal(raw_resamp._data.shape[1], raw_resamp.n_times) assert_equal(raw._data.shape[0], raw_resamp._data.shape[0]) # test non-parallel on downsample - raw_resamp.resample(sfreq, n_jobs=1) + raw_resamp.resample(sfreq, n_jobs=1, npad='auto') assert_equal(raw_resamp.info['sfreq'], sfreq) assert_equal(raw._data.shape, raw_resamp._data.shape) assert_equal(raw.first_samp, raw_resamp.first_samp) @@ -892,9 +892,9 @@ def test_resample(): raw3 = raw.copy() raw4 = raw.copy() raw1 = concatenate_raws([raw1, raw2]) - raw1.resample(10.) - raw3.resample(10.) - raw4.resample(10.) + raw1.resample(10., npad='auto') + raw3.resample(10., npad='auto') + raw4.resample(10., npad='auto') raw3 = concatenate_raws([raw3, raw4]) assert_array_equal(raw1._data, raw3._data) assert_array_equal(raw1._first_samps, raw3._first_samps) @@ -909,12 +909,12 @@ def test_resample(): # basic decimation stim = [1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0] raw = RawArray([stim], create_info(1, len(stim), ['stim'])) - assert_allclose(raw.resample(8.)._data, + assert_allclose(raw.resample(8., npad='auto')._data, [[1, 1, 0, 0, 1, 1, 0, 0]]) # decimation of multiple stim channels raw = RawArray(2 * [stim], create_info(2, len(stim), 2 * ['stim'])) - assert_allclose(raw.resample(8.)._data, + assert_allclose(raw.resample(8., npad='auto')._data, [[1, 1, 0, 0, 1, 1, 0, 0], [1, 1, 0, 0, 1, 1, 0, 0]]) @@ -922,7 +922,7 @@ def test_resample(): # done naively stim = [0, 0, 0, 1, 1, 0, 0, 0] raw = RawArray([stim], create_info(1, len(stim), ['stim'])) - assert_allclose(raw.resample(4.)._data, + assert_allclose(raw.resample(4., npad='auto')._data, [[0, 1, 1, 0]]) # two events are merged in this case (warning) @@ -930,7 +930,7 @@ def test_resample(): raw = RawArray([stim], create_info(1, len(stim), ['stim'])) with warnings.catch_warnings(record=True) as w: warnings.simplefilter('always') - raw.resample(8.) + raw.resample(8., npad='auto') assert_true(len(w) == 1) # events are dropped in this case (warning) @@ -938,28 +938,28 @@ def test_resample(): raw = RawArray([stim], create_info(1, len(stim), ['stim'])) with warnings.catch_warnings(record=True) as w: warnings.simplefilter('always') - raw.resample(4.) + raw.resample(4., npad='auto') assert_true(len(w) == 1) # test resampling events: this should no longer give a warning stim = [0, 1, 1, 0, 0, 1, 1, 0] raw = RawArray([stim], create_info(1, len(stim), ['stim'])) events = find_events(raw) - raw, events = raw.resample(4., events=events) + raw, events = raw.resample(4., events=events, npad='auto') assert_equal(events, np.array([[0, 0, 1], [2, 0, 1]])) # test copy flag stim = [1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0] raw = RawArray([stim], create_info(1, len(stim), ['stim'])) - raw_resampled = raw.resample(4., copy=True) + raw_resampled = raw.resample(4., npad='auto', copy=True) assert_true(raw_resampled is not raw) - raw_resampled = raw.resample(4., copy=False) + raw_resampled = raw.resample(4., npad='auto', copy=False) assert_true(raw_resampled is raw) # resample should still work even when no stim channel is present raw = RawArray(np.random.randn(1, 100), create_info(1, 100, ['eeg'])) - raw.resample(10) - assert_true(len(raw) == 10) + raw.resample(10, npad='auto') + assert_equal(len(raw), 10) @testing.requires_testing_data diff --git a/mne/tests/test_filter.py b/mne/tests/test_filter.py index 216f01760b3..fc897526e5e 100644 --- a/mne/tests/test_filter.py +++ b/mne/tests/test_filter.py @@ -32,7 +32,7 @@ def test_1d_filter(): h = np.concatenate([[1.], np.zeros(n_filter - 1)]) # ensure we pad the signal the same way for both filters n_pad = max(min(n_filter, n_signal - 1), 0) - x_pad = _smart_pad(x, n_pad) + x_pad = _smart_pad(x, np.array([n_pad, n_pad])) for zero_phase in (True, False): # compute our expected result the slow way if zero_phase: @@ -342,12 +342,15 @@ def test_cuda(): for o in out]) == tot) # check resampling - a = rng.randn(3, sig_len_secs * sfreq) - a1 = resample(a, 1, 2, n_jobs=2, npad=0) - a2 = resample(a, 1, 2, n_jobs='cuda', npad=0) - a3 = resample(a, 2, 1, n_jobs=2, npad=0) - a4 = resample(a, 2, 1, n_jobs='cuda', npad=0) - assert_array_almost_equal(a3, a4, 14) + for window in ('boxcar', 'triang'): + for N in (997, 1000): # one prime, one even + a = rng.randn(2, N) + for fro, to in ((1, 2), (2, 1), (1, 3), (3, 1)): + a1 = resample(a, fro, to, n_jobs=1, npad='auto', + window=window) + a2 = resample(a, fro, to, n_jobs='cuda', npad='auto', + window=window) + assert_allclose(a1, a2, rtol=1e-7, atol=1e-14) assert_array_almost_equal(a1, a2, 14) assert_array_equal(resample([0, 0], 2, 1, n_jobs='cuda'), [0., 0., 0., 0.]) assert_array_equal(resample(np.zeros(2, np.float32), 2, 1, n_jobs='cuda'), From 9b6555bad2779d0425a80de4bd265ccc86fd1d92 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 2 Mar 2016 11:51:10 -0500 Subject: [PATCH 2/3] FIX: Fix warning [ci skip] --- mne/io/base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mne/io/base.py b/mne/io/base.py index 4a5d6c6c29d..6623b4021e4 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -1011,7 +1011,8 @@ def resample(self, sfreq, npad=None, window='boxcar', stim_picks=None, 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.') + '"auto" in 0.12. Please set the value explicitly.', + DeprecationWarning) if not self.preload: raise RuntimeError('Can only resample preloaded data') From 954fe8a235f3fb5cc8f35166ff51719494b35e55 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 2 Mar 2016 12:30:48 -0500 Subject: [PATCH 3/3] FIX: minor fixes --- examples/preprocessing/plot_resample.py | 9 +++++---- mne/epochs.py | 11 +++++++++-- mne/evoked.py | 11 +++++++++-- mne/filter.py | 11 ++++++----- mne/inverse_sparse/tests/test_gamma_map.py | 2 +- mne/source_estimate.py | 11 +++++++++-- mne/tests/test_cov.py | 2 +- mne/tests/test_epochs.py | 4 ++-- mne/tests/test_evoked.py | 4 ++-- mne/tests/test_proj.py | 2 +- tutorials/plot_cluster_stats_spatio_temporal.py | 4 ++-- tutorials/plot_cluster_stats_spatio_temporal_2samp.py | 2 +- ...r_stats_spatio_temporal_repeated_measures_anova.py | 2 +- 13 files changed, 49 insertions(+), 26 deletions(-) diff --git a/examples/preprocessing/plot_resample.py b/examples/preprocessing/plot_resample.py index d58f9b54d61..026a67849b5 100644 --- a/examples/preprocessing/plot_resample.py +++ b/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)) diff --git a/mne/epochs.py b/mne/epochs.py index 8b6e0ba0f60..edebb526e77 100644 --- a/mne/epochs.py +++ b/mne/epochs.py @@ -1438,7 +1438,7 @@ 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 @@ -1446,8 +1446,10 @@ def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1, ---------- 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 diff --git a/mne/evoked.py b/mne/evoked.py index bff55a96b4f..fac24fb5736 100644 --- a/mne/evoked.py +++ b/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) diff --git a/mne/filter.py b/mne/filter.py index b0000d97702..4b569d4b0c7 100644 --- a/mne/filter.py +++ b/mne/filter.py @@ -1322,11 +1322,11 @@ def resample(x, up, down, npad=100, axis=-1, window='boxcar', n_jobs=1, if isinstance(npad, string_types): if npad != 'auto': raise ValueError(bad_msg) - # Figure out pad that gets us to a power of 2 with at least - # 100 on each end - npad = 2 ** int(np.ceil(np.log2(x_len + 200))) - x_len + # 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]) # XXX + npads = np.array([npad, npad + extra], int) else: if npad != int(npad): raise ValueError(bad_msg) @@ -1341,7 +1341,8 @@ def resample(x, up, down, npad=100, axis=-1, window='boxcar', n_jobs=1, to_removes = [int(round(ratio * npads[0]))] to_removes.append(new_len - final_len - to_removes[0]) to_removes = np.array(to_removes) - assert np.abs(to_removes[1] - to_removes[0]) <= int(np.ceil(ratio)) + # 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: diff --git a/mne/inverse_sparse/tests/test_gamma_map.py b/mne/inverse_sparse/tests/test_gamma_map.py index f28b9170df3..98f20d9e73c 100644 --- a/mne/inverse_sparse/tests/test_gamma_map.py +++ b/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) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 7f18a59553a..91d9cf3fb9e 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -505,7 +505,7 @@ def crop(self, tmin=None, tmax=None): return self # return self for chaining methods @verbose - def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1, + def resample(self, sfreq, npad=None, window='boxcar', n_jobs=1, verbose=None): """Resample data @@ -513,8 +513,10 @@ def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1, ---------- 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 @@ -530,6 +532,11 @@ def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1, Note that the sample rate of the original data is inferred from tstep. """ + 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) # resampling in sensor instead of source space gives a somewhat # different result, so we don't allow it self._remove_kernel_sens_data_() diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py index 106f68d73a9..0edb7b75fc6 100644 --- a/mne/tests/test_cov.py +++ b/mne/tests/test_cov.py @@ -448,7 +448,7 @@ def test_compute_covariance_auto_reg(): """Test automated regularization""" raw = Raw(raw_fname, preload=True) - raw.resample(100) # much faster estimation + raw.resample(100, npad='auto') # much faster estimation events = find_events(raw, stim_channel='STI 014') event_ids = [1, 2, 3, 4] reject = dict(mag=4e-12) diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py index 8686287f8e2..aab5cb63882 100644 --- a/mne/tests/test_epochs.py +++ b/mne/tests/test_epochs.py @@ -1123,8 +1123,8 @@ def test_resample(): epochs1 = EpochsArray(data, deepcopy(info), events) epochs2 = EpochsArray(data, deepcopy(info), events) epochs = concatenate_epochs([epochs1, epochs2]) - epochs1.resample(epochs1.info['sfreq'] // 2) - epochs2.resample(epochs2.info['sfreq'] // 2) + epochs1.resample(epochs1.info['sfreq'] // 2, npad='auto') + epochs2.resample(epochs2.info['sfreq'] // 2, npad='auto') epochs = concatenate_epochs([epochs1, epochs2]) for e in epochs1, epochs2, epochs: assert_equal(e.times[0], epochs.tmin) diff --git a/mne/tests/test_evoked.py b/mne/tests/test_evoked.py index 44ff9fa7b6a..c882d1a24d0 100644 --- a/mne/tests/test_evoked.py +++ b/mne/tests/test_evoked.py @@ -175,7 +175,7 @@ def test_evoked_resample(): # upsample, write it out, read it in ave = read_evokeds(fname, 0) sfreq_normal = ave.info['sfreq'] - ave.resample(2 * sfreq_normal) + ave.resample(2 * sfreq_normal, npad=100) write_evokeds(op.join(tempdir, 'evoked-ave.fif'), ave) ave_up = read_evokeds(op.join(tempdir, 'evoked-ave.fif'), 0) @@ -184,7 +184,7 @@ def test_evoked_resample(): # and compare the original to the downsampled upsampled version ave_new = read_evokeds(op.join(tempdir, 'evoked-ave.fif'), 0) - ave_new.resample(sfreq_normal) + ave_new.resample(sfreq_normal, npad=100) assert_array_almost_equal(ave_normal.data, ave_new.data, 2) assert_array_almost_equal(ave_normal.times, ave_new.times) diff --git a/mne/tests/test_proj.py b/mne/tests/test_proj.py index 23700afd0d5..05a0bb3a256 100644 --- a/mne/tests/test_proj.py +++ b/mne/tests/test_proj.py @@ -210,7 +210,7 @@ def test_compute_proj_raw(): # test resampled-data projector, upsampling instead of downsampling # here to save an extra filtering (raw would have to be LP'ed to be equiv) raw_resamp = cp.deepcopy(raw) - raw_resamp.resample(raw.info['sfreq'] * 2, n_jobs=2) + raw_resamp.resample(raw.info['sfreq'] * 2, n_jobs=2, npad='auto') with warnings.catch_warnings(record=True) as w: warnings.simplefilter('always') projs = compute_proj_raw(raw_resamp, duration=None, stop=raw_time, diff --git a/tutorials/plot_cluster_stats_spatio_temporal.py b/tutorials/plot_cluster_stats_spatio_temporal.py index ef2a7df5ee7..9e5c051a4a2 100644 --- a/tutorials/plot_cluster_stats_spatio_temporal.py +++ b/tutorials/plot_cluster_stats_spatio_temporal.py @@ -75,10 +75,10 @@ # Let's average and compute inverse, resampling to speed things up evoked1 = epochs1.average() -evoked1.resample(50) +evoked1.resample(50, npad='auto') condition1 = apply_inverse(evoked1, inverse_operator, lambda2, method) evoked2 = epochs2.average() -evoked2.resample(50) +evoked2.resample(50, npad='auto') condition2 = apply_inverse(evoked2, inverse_operator, lambda2, method) # Let's only deal with t > 0, cropping to reduce multiple comparisons diff --git a/tutorials/plot_cluster_stats_spatio_temporal_2samp.py b/tutorials/plot_cluster_stats_spatio_temporal_2samp.py index 19d774c2138..69e512fbb81 100644 --- a/tutorials/plot_cluster_stats_spatio_temporal_2samp.py +++ b/tutorials/plot_cluster_stats_spatio_temporal_2samp.py @@ -34,7 +34,7 @@ # Load stc to in common cortical space (fsaverage) stc = mne.read_source_estimate(stc_fname) -stc.resample(50) +stc.resample(50, npad='auto') stc = mne.morph_data('sample', 'fsaverage', stc, grade=5, smooth=20, subjects_dir=subjects_dir) diff --git a/tutorials/plot_cluster_stats_spatio_temporal_repeated_measures_anova.py b/tutorials/plot_cluster_stats_spatio_temporal_repeated_measures_anova.py index 1a5ecf1215c..c0f851fe19a 100644 --- a/tutorials/plot_cluster_stats_spatio_temporal_repeated_measures_anova.py +++ b/tutorials/plot_cluster_stats_spatio_temporal_repeated_measures_anova.py @@ -83,7 +83,7 @@ conditions = [] for cond in ['l_aud', 'r_aud', 'l_vis', 'r_vis']: # order is important evoked = epochs[cond].average() - evoked.resample(50) + evoked.resample(50, npad='auto') condition = apply_inverse(evoked, inverse_operator, lambda2, method) # Let's only deal with t > 0, cropping to reduce multiple comparisons condition.crop(0, None)