Resample function optimizations

[choldgraf]

Currently it takes me quite a long time to carry out a resample function. I've got about 120 channels with about 10 minutes of data, sampled at 5Khz, which equals about 120 x 3,000,000 data points. Obviously my first step in this process is to resample the data such that it's not so densely-sampled. However, this can take a really long time and use a lot of memory.

I have two questions regarding this:

1. Is the resampling optimized at all in terms of making sure that the signal is a ^2 in length? In the past I've tried to pad the signal with 0's to the next pow2, but this can be a memory issue if there are enough datapoints. Another option is to find a close number with lots of factors, no?
2. What kind of speedup do you all see when you use CUDA? Is it worth the time to set up in terms of the performance boost?

Anyway just trying to see if anyone has thoughts on improving the speed of these functions. If more people use MNE that do ECoG research (and thus have little control over the recording parameters), it may prove useful.

[agramfort]

the culprit is the resample method from scipy.signal that honestly sucks... it does a dummy fft over the entire signal length... the problem has to be fixed in scipy. you might want to try: http://scikits.appspot.com/samplerate

[larsoner]

In my experince CUDA is about an order of magnitude faster. I'm pretty sure that was both for FIR filtering and resampling. What would really help the CPU case is using upfirdn as discussed in #1814.

[choldgraf]

Yeah, the scipy signal processing is the one thing that makes me miss matlab (just a little bit). I've tried creating functions like make_len_pow_2 that just dumbly the next power of 2 and appends 0's to the signal. That definitely speeds things up, but can explode the memory. What if something like that was added to the resample function as an option to npad? So if you put in a string like 'pow2' then it would find how many datapoints to add to make it a power of 2, then split it in half and add it to the front/back of the signal?

[larsoner]

+1 for the string idea, we already do something similar for filter lengths

[choldgraf]

Did you already begin an upfirdn contribution to scipy? That may be outside of my signal processing chops :/

[larsoner]

No, and it's going to take some time to get right. So for now the padding to power of 2 should help.

[larsoner]

Someone actually has a SWIG'ed version available they gave me permission to relicense as BSD for scipy, but I'm not sure how complex that problem is going to be.

[choldgraf]

Ah ok - for the resample string thing, do you know if there is a function
in python for "give me a number, and I will tell you the next highest
number with greater than N factors"? I've tried padding to a power of two,
but once you're in the higher regions it becomes memory-limited. That said,
this sounds like it may be too complicated to include natively :P.

Actually another thought might be to allow someone to specify the total
length that they want the signal to be, rather than the amount of padding.
E.g.: have some parameter length, take the diff from current length,
divide that in half, and add zeros to the front/back accordingly.

On Wed, Apr 29, 2015 at 7:00 AM, Eric Larson notifications@github.com
wrote:

Someone actually has a SWIG'ed version available they gave me permission
to relicense as BSD for scipy, but I'm not sure how complex that problem is
going to be.

—
Reply to this email directly or view it on GitHub
#2035 (comment)
.

[larsoner]

No, I don't know of one. If memory is an issue, one thing that might fix that is iterating over channels instead of operating on them as a contiguous block (e.g., as axis=1). The overhead for a few hundred channels should be pretty negligible, but you'd want to test it to make sure.

[choldgraf]

Ya that is true - if you do it for one channel at a time it's not a big
deal, so in that case you would basically do:

  all_res_chans = []
  for chan in channels:
    chan = pad_to_power_2(chan)
    chan_res = resample(chan)
    all_res_chans.append(remove_padding(chan_res))
  make_back_into_Raw(all_res_chans)

or something like this

On Wed, Apr 29, 2015 at 9:55 AM, Eric Larson notifications@github.com
wrote:

No, I don't know of one. If memory is an issue, one thing that might fix
that is iterating over channels instead of operating on them as a
contiguous block (e.g., as axis=1). The overhead for a few hundred
channels should be pretty negligible, but you'd want to test it to make
sure.

—
Reply to this email directly or view it on GitHub
#2035 (comment)
.

[larsoner]

Yeah, that's the pseudocode anyway. Getting it to work with the current resampling might take some work, but should be doable, at least for the n_jobs=1 case. Sometimes we branch based on whether n_jobs == 1 or not, and this would be a case where it could be useful to do so. You can't really avoid the memcpy-type operation with multiple jobs, unfortunately, so it will only work for a single-job case.

You're welcome to take a stab at it if you have time. I think if n_jobs=1 you'd want to preallocate the new array, and fill it in row-by-row by resampling.

[choldgraf]

mm that's a good idea. I will try to include this if I finish my stuff for
this week early. Have a plane flight tomorrow so maybe then.

On Wed, Apr 29, 2015 at 10:05 AM, Eric Larson notifications@github.com
wrote:

Yeah, that's the pseudocode anyway. Getting it to work with the current
resampling might take some work, but should be doable, at least for the
n_jobs=1 case. Sometimes we branch based on whether n_jobs == 1 or not,
and this would be a case where it could be useful to do so. You can't
really avoid the memcpy-type operation with multiple jobs, unfortunately,
so it will only work for a single-job case.

You're welcome to take a stab at it if you have time. I think if n_jobs=1
you'd want to preallocate the new array, and fill it in row-by-row by
resampling.

—
Reply to this email directly or view it on GitHub
#2035 (comment)
.

[larsoner]

Hopefully I'll have time to actually make this PR work and get it into scipy:

scipy/scipy#5186

[larsoner]

An update -- this should be made ~2x faster in the near term (months) by scipy/scipy#5592 and the related scipy/scipy#5610. By my estimate, using upfirdn through scipy/scipy#5186 will require more time (~6-12 months).

[choldgraf]

That's great - looking forward to not dreading my resampling steps :)

[larsoner]

FYI scipy/scipy#5610 (upfirdn) has been merged and I've started a upfirdn-based resample discussion in scipy/scipy#5746.

[choldgraf]

wohoo! wheels are a-turning :)

On Wed, Jan 20, 2016 at 10:43 AM, Eric Larson notifications@github.com
wrote:

FYI scipy/scipy#5610 scipy/scipy#5610 (upfirdn)
has been merged and I've started a upfirdn-based resample discussion in
scipy/scipy#5746 scipy/scipy#5746.

—
Reply to this email directly or view it on GitHub
#2035 (comment)
.

[jona-sassenhagen]

Indeed, v. cool.