ENH: Improve FIR path of signal.decimate #5975

e-q · 2016-03-16T04:18:08Z

This PR aims to relieve the phase shift due to the group delay of an FIR downsampling filter in signal.decimate through the use of resample_poly.

This requires modifying resample_poly to accept arbitrary FIR filter coefficients, as decimate itself advertises that capability. A test for resample_poly has been added for this capability. Hopefully this isn't too much of a hack. (Thoughts, @Eric89GXL ?)

codecov-io · 2016-03-16T04:52:07Z

@@            master   #5975   diff @@
======================================
  Files          238     238       
  Stmts        43803   43814    +11
  Branches      8211    8214     +3
  Methods          0       0       
======================================
+ Hit          34230   34244    +14
+ Partial       2603    2602     -1
+ Missed        6970    6968     -2

Review entire Coverage Diff as of 31f7fbd

Powered by Codecov. Updated on successful CI builds.

larsoner · 2016-03-16T15:22:33Z

This requires modifying resample_poly to accept arbitrary FIR filter coefficients

It's good to add that capability. I didn't do it originally to keep the PR simple, and to offload the work of getting it right on the next developer. Looks like you fell into my trap :)

larsoner · 2016-03-16T15:23:23Z

scipy/signal/signaltools.py

@@ -1896,6 +1896,9 @@ def resample_poly(x, up, down, axis=0, window=('kaiser', 5.0)):
        Desired window to use to design the low-pass filter. See
        `scipy.signal.get_window` for a list of windows and required
        parameters.
+    num : array_like, optional
+         FIR filter coefficients of the low-pass filter to be used. Overrides
+         the `window` argument if not None.


Just roll this into window. See e.g. how resample allows multiple types for window:

http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.signal.resample.html

larsoner · 2016-03-16T15:28:58Z

So we're swapping in (currently) resample_poly for filtfilt and (soon) upfirdn for lfilter to save time. It would be nice in the tests to have the outputs computed using the naive filtfilt and lfilter algorithms, and compare outputs of decimate (which will use resample_poly and upfirdn under the hood). They should hopefully match well using assert_allclose. It helps ensure we are really getting equivalent operations. WDYT?

larsoner · 2016-03-16T15:33:34Z

BTW can you add a comment about the speed gains? We doing 1 / down operations by using upfirdn compared to lfilter. Also, something neat under the hood is that instead of doing a two-pass operation for filtfilt, resample_poly does a single filter operation by taking advantage of associativity (x(t) * h(t)) * h(-t) = (x(t) * (h(t) * h(-t))), i.e. convolving the filter with it's time-reversed version, then dealing with the zero padding properly.

Speaking of which, this is one way that you might need to change the code -- I think that filtfilt would apply the same filter twice -- once forward and once backward -- whereas resample_poly will only apply it once (and it will assume you've constructed your filter in this symmetrical way). So maybe in this function, the window function/array needs to be convolved with it's time-reversed version before being passed to resample_poly? You should see it in the tests if you implement what I mention above.

e-q · 2016-03-24T03:20:29Z

Thanks for the feedback @Eric89GXL! I've incorporated much of your feedback.

The speed gains of upfirdn are quite good for larger downsampling factors, which is great. Where were you thinking of the speedups being mentioned, besides the commit message?

As for the upfirdn vs. lfilter results, I'm opening a separate PR that adds it to the upfirdn tests (which seems a little more natural to me). I haven't yet tried resample_poly vs. filtfilt yet. I think this should work as long as the filtfilt is explicitly given the zero-padding argument...

larsoner · 2016-03-24T13:58:45Z

Where were you thinking of the speedups being mentioned, besides the commit message?

A comment in the code saying why upfirdn is preferred to lfilter (polyphase resampling avoids unnecessary calculations) would be a good place

I haven't yet tried resample_poly vs. filtfilt yet. I think this should work as long as the filtfilt is explicitly given the zero-padding argument...

I don't expect the zero-padding to really matter. resample_poly must zero pad and remove zeros at the end because it only does a single pass, in one direction, and computes only specific subsamples -- this introduces a phase shift, so the data are not centered after the operation. filtfilt on the other hand does two passes, one forward and one backward, and this is a zero-phase operation. So hopefully you don't have to do anything too special, mostly:

Pass h as the filter to filtfilt
Extend resample_poly to allow for arbitrary filters
Pass convolve(h, h[::-1]) as the filter to resample_poly

e-q · 2016-03-26T04:52:33Z

Ok, testing vs filtfilt is included now. You were exactly right about passing the convolved filter, thanks.

I was unclear before, the zero-padding I was referring to was in regards to the input signal, not the filter. filtfilt does an odd signal extension by default, but can do constant padding. Since, according to the docstring, resample_poly assumes zeros outside of the input samples, I set the first and last samples to zero and used filtfilt's constant padding mode. Anyways, no big deal, the test passes on my machine.

larsoner · 2016-03-26T16:49:45Z

That makes sense. Glad they match. I'll look deeper soon. Any practical numbers on what the resample_poly speedup is like over filtfilt? I'm curious how close it gets to being `down` times faster.

e-q · 2016-03-27T19:54:25Z

Decimating x=np.random.randn(10**7) with resample_poly, using a 31 tap FIR filter, is actually faster than filtfilt by more than a factor of down on my machine.

q    filtfilt    resample_poly   speedup
----------------------------------------
2    591ms       211.0ms          2.8
5    577ms        98.7ms          5.8
13   604ms        44.4ms         13.6

This comes from the additional time cost of slicing the output of filtfilt. Without the slicing, the speedup for q=13 is 12.97, for instance.

larsoner · 2016-03-27T19:55:36Z

Excellent. Is this good to go from your end, then?

e-q · 2016-03-27T19:56:37Z

yep!

larsoner · 2016-03-27T20:00:17Z

scipy/signal/signaltools.py

-    f_c = 1. / max_rate  # cutoff of FIR filter (rel. to Nyquist)
-    half_len = 10 * max_rate  # reasonable cutoff for our sinc-like function
-    h = firwin(2 * half_len + 1, f_c, window=window)
+    n_out = int(np.ceil(x.shape[axis] * up / down))


This used to effectively be floor, why the change to ceil?

I assume it's for the extra sample in case it doesn't divide evenly, In that case it would actually be better to use integer arithmetic if possible so we stay exact, like:

n_out = x.shape['axis'] * up n_out = n_out // down + bool(n_out % down)

larsoner · 2016-03-27T20:08:23Z

Other than my minor comments, LGTM.

e-q · 2016-03-27T20:45:28Z

Good points all around, addressed and squashed. Thanks!

The short answer about the change of n_out for resample_poly was to make it match the number of outputs from decimate, which in turn is given by how many outputs you get from slicing x[::q].

Incidentally, you'll see in the zero_phase=False, ftype='fir' path that the output from upfirdn is truncated to that same amount. It seems that upfirdn tacks a few more output points on than I would expect from the decimation point of view. Why is this?

larsoner · 2016-03-27T20:48:08Z

IIRC upfirdn lets the filter ring out to the end, like np.convolve(..., mode='full')

e-q · 2016-04-05T05:14:17Z

Is this good to go?

larsoner · 2016-04-05T12:23:58Z

scipy/signal/signaltools.py

+    step, so it should be designed to operate on a signal at a sampling
+    frequency higher than the original by a factor of `up//gcd(up, down)`.
+    This function's output will be centered with respect to this array, so it
+    is best to pass a symmetric filter with an odd number of samples to if, as


"number of samples to if" -> "number of samples if"

larsoner · 2016-04-05T12:26:33Z

Other than my one new gripe, +1 for merge. After the rewording, let's wait a day or two to see if anyone else wants to comment then merge

This PR aims to relieve the phase shift due to the group delay of an FIR downsampling filter in signal.decimate through the use of `resample_poly`. This requires modifying `resample_poly` to accept arbitrary FIR filter coefficients, as `decimate` itself advertises that capability. A test for `resample_poly` has been added for this capability. Additionally, the "traditional" FIR path has been sped up via the use of `upfirdn` which only calculates every q'th output.

Since `resample_poly` is used instead of `filtfilt` in `decimate` for FIR decimation for speed reasons, this adds a test to ensure equivalent output of these two methods.

e-q · 2016-04-05T21:22:52Z

Thanks for the look! The wording has been addressed.

larsoner · 2016-04-08T17:02:53Z

Thanks @e-q

pv · 2016-04-09T22:41:04Z

benchmarks/benchmarks/signal_filtering.py

@@ -3,12 +3,29 @@
 import numpy as np

 try:
-    from scipy.signal import lfilter, firwin
+    from scipy.signal import lfilter, firwin, decimate


Usually best to put these in a separate import statement, so that not having decimate does not prevent running eg the firwin benchmarks

larsoner reviewed Mar 16, 2016
View reviewed changes

larsoner added enhancement A new feature or improvement scipy.signal needs-work Items that are pending response from the author labels Mar 22, 2016

e-q force-pushed the firdecimate branch from c3c7222 to e1327bc Compare March 24, 2016 03:20

e-q mentioned this pull request Mar 24, 2016

TST: Break up upfirdn tests & compare to lfilter #5997

Merged

e-q force-pushed the firdecimate branch from e1327bc to 2262453 Compare March 26, 2016 04:48

e-q changed the title ~~WIP, ENH: Improve FIR path of signal.decimate~~ ENH: Improve FIR path of signal.decimate Mar 26, 2016

larsoner reviewed Mar 27, 2016
View reviewed changes

e-q force-pushed the firdecimate branch from 2262453 to 52d8bf4 Compare March 27, 2016 20:40

larsoner reviewed Apr 5, 2016
View reviewed changes

e-q added 2 commits April 5, 2016 13:41

TST: Compare signal.resample_poly to filtfilt

17e7ce6

Since `resample_poly` is used instead of `filtfilt` in `decimate` for FIR decimation for speed reasons, this adds a test to ensure equivalent output of these two methods.

e-q force-pushed the firdecimate branch from 52d8bf4 to 17e7ce6 Compare April 5, 2016 21:04

larsoner merged commit b65314b into scipy:master Apr 8, 2016

ev-br added this to the 0.18.0 milestone Apr 8, 2016

pv reviewed Apr 9, 2016
View reviewed changes

e-q deleted the firdecimate branch April 10, 2016 04:30

larsoner mentioned this pull request Oct 12, 2016

MemoryError when decimating large array #6669

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: Improve FIR path of signal.decimate #5975

ENH: Improve FIR path of signal.decimate #5975

e-q commented Mar 16, 2016

codecov-io commented Mar 16, 2016

larsoner commented Mar 16, 2016

larsoner Mar 16, 2016

larsoner commented Mar 16, 2016

larsoner commented Mar 16, 2016

e-q commented Mar 24, 2016

larsoner commented Mar 24, 2016

e-q commented Mar 26, 2016

larsoner commented Mar 26, 2016 via email

e-q commented Mar 27, 2016

larsoner commented Mar 27, 2016

e-q commented Mar 27, 2016

larsoner Mar 27, 2016

larsoner commented Mar 27, 2016

e-q commented Mar 27, 2016

larsoner commented Mar 27, 2016

e-q commented Apr 5, 2016

larsoner Apr 5, 2016

larsoner commented Apr 5, 2016

e-q commented Apr 5, 2016 via email

larsoner commented Apr 8, 2016

pv Apr 9, 2016

ENH: Improve FIR path of signal.decimate #5975

ENH: Improve FIR path of signal.decimate #5975

Conversation

e-q commented Mar 16, 2016

codecov-io commented Mar 16, 2016

larsoner commented Mar 16, 2016

larsoner Mar 16, 2016

Choose a reason for hiding this comment

larsoner commented Mar 16, 2016

larsoner commented Mar 16, 2016

e-q commented Mar 24, 2016

larsoner commented Mar 24, 2016

e-q commented Mar 26, 2016

larsoner commented Mar 26, 2016 via email

e-q commented Mar 27, 2016

larsoner commented Mar 27, 2016

e-q commented Mar 27, 2016

larsoner Mar 27, 2016

Choose a reason for hiding this comment

larsoner commented Mar 27, 2016

e-q commented Mar 27, 2016

larsoner commented Mar 27, 2016

e-q commented Apr 5, 2016

larsoner Apr 5, 2016

Choose a reason for hiding this comment

larsoner commented Apr 5, 2016

e-q commented Apr 5, 2016 via email

larsoner commented Apr 8, 2016

pv Apr 9, 2016

Choose a reason for hiding this comment