MRG: Add polyphase resampling

[larsoner]

scipy.signal.resample_poly has been in SciPy since 0.18. Time to expose it.

Originally I was going to make it the new default -- it does not assume signal periodicity (an implicit assumption of FFT-based resampling) and could be faster for integer sample rate changes (especially if the signal of interest has a non-easily-factorable number of samples so the FFT/IFFT will be slow). But given how long we've been using fft and the fact that sample and somato use wacky non-integral sample rates, I guess we can leave it as fft.

[codecov]

Codecov Report

Merging #5136 into master will increase coverage by <.01%.
The diff coverage is 84.31%.

@@            Coverage Diff             @@
##           master    #5136      +/-   ##
==========================================
+ Coverage   88.13%   88.13%   +<.01%     
==========================================
  Files         358      358              
  Lines       65784    65843      +59     
  Branches    11191    11204      +13     
==========================================
+ Hits        57976    58031      +55     
- Misses       4983     4987       +4     
  Partials     2825     2825

[agramfort]

I am not against this feature but i'd like to be sure it brings something useful on some settings. Faster? more accurate?

[larsoner]

See PR description about assumptions, integer values, and primes. Do you want code examples?

[agramfort]

yes. What is in the impact on a final ERP/ERF for example? is it significantly faster? I imagine it does not change the results of an ERF/ERP study no?

[larsoner]

It depends on options. By default npad chooses something that is designed to make the FFT option not too slow (though it's not always perfect). If you want to eliminate padding or set it to some fixed value rather than have it vary based on what makes a good FFT length, then you can hit cases like:

import time
import numpy as np
import mne
raw = mne.io.RawArray(np.ones((1, 26594)),
                      mne.create_info(1, 1000., 'eeg'))
for method in ('poly', 'fft'):
    t0 = time.time()
    raw.copy().resample(500, npad=0, method=method)
    print(time.time() - t0)

Getting 0.0035 sec for poly and 0.2329 sec for FFT (because the resulting length of 13297 is prime).

[jona-sassenhagen]

Wouldn't it make sense to make functions out of all three of these?

[larsoner]

Sorry I meant to put "what do you mean" here. The triaging needs to be done somewhere and this is the first level where mode actually matters (e.g. upstream callers pad before and remove padding after the same way regardless of method, etc.). So I don't see a split elsewhere being cleaner

[larsoner]

Ah you must mean move the other two bits into functions. I guess we could, you think it should be more readable?

[jona-sassenhagen]

Yes.

Also, looking at this .... the function is called fft_resample. That would be a bit misleading, right?

[larsoner]

Yes but probably not worth the effort to rename

[jona-sassenhagen]

Can't we add the benefits/drawbacks to the docs of the resample method?

[larsoner]

What do you mean?

…

[jona-sassenhagen]

method: str
    Can be ...
    `poly` is faster and more accurate if the resulting length is prime ... 

Something like that? Or at least a reference to a section of the docs explaining it.

[larsoner]

@jona-sassenhagen can you see if the documentation is better now?

[jona-sassenhagen]

Yes, and I think factoring out the fft resampling into a function looks better, too.

So there are no expected differences in quality? Otherwise, that should also probably be documented.

If not, lgtm

[larsoner]

There can be slight differences in behavior in some cases, e.g.

https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.resample_poly.html

But for our data I doubt it would make much difference, especially with a non-zero pad (which is the default anyway).

[larsoner]

@jona-sassenhagen okay for you?

[jona-sassenhagen]

Sure.

[jona-sassenhagen]

Shall I merge? Or do you indeed want to wait for .17?

[agramfort]

please don't merge. give me time to test and look I am doing the 0.16 release now.

[jona-sassenhagen]

great!

[cbrnr]

@agramfort 0.16.1 has been on PyPI since at least yesterday - is this intentional?

[agramfort]

yes I screwed up 0.16 so we are already at 0.16.1 I just sent the announcement email

[jona-sassenhagen]

Rebase and merge right?

[jona-sassenhagen]

Anything standing in the way of merging this?

[jona-sassenhagen]

(I mean, besides for the merge conflict etc of course :D)

[agramfort]

sorry this is not high on my todo list now. However if I can make a comment I think we should have method="auto" that defaults to fft or poly depending on what will be more efficient. Otherwise people will not discover this option.

[larsoner]

Since the assumptions of the two methods (and resulting properties like ringing) are not equivalent for all signals between the two methods, I would not be comfortable switching between them automatically :(

We could always kill this PR and let experts use polyphase manually by getting data and using *Array constructors.

[agramfort]

I am just uneasy to add more complex code for stuff that already work and will only benefit to a tiny fraction of people.

[larsoner]

Okay let's leave it out for now at least then