ENH: signal: Add Gustafsson's method as an option for the filtfilt function

[WarrenWeckesser]

The paper on which the enhancement is based is available here: http://users.isy.liu.se/en/rt/fredrik/pub.html

[coveralls]

Coverage remained the same when pulling 772dbc4 on WarrenWeckesser:gust into 473dca9 on scipy:master.

[rgommers]

This looks promising at first glance. Could you rebase it? I'll attempt to finally review it soon.

[WarrenWeckesser]

Thanks, Ralf. I'm on the road this week, so I probably won't get back to this until next week at the earliest.

[rgommers]

Here's a rebased branch: https://github.com/rgommers/scipy/tree/gust. It contains one extra commit with fixes for reST issues and for the plot in the example not rendering. Could you take that over?

[rgommers]

Main issue with reST was the difference between single and double backticks. You were consistently using single backticks for code fragments (as well as for referring to functions/variables of course).

[rgommers]

Why a new file for the tests? There are already filtfilt tests in test_signaltools.py. The gust tests should be added there.

[rgommers]

If x0 and x1 are not returned from filtfilt, why are they returned from the private function and tested? If there's a reason, may want to add that to the Notes section. Now it only says that they're not needed.

[rgommers]

In the example, maybe it would help to show the differences at the edges graphically with something like:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


# generate some signal data
t = np.linspace(0, 20 * np.pi, num=1000)
sig = np.cos(t) + 0.3 * np.sin(5*t) + 2 * np.sin(0.2 * t)
sig_noisy = sig + 0.2 * np.random.randn(t.size) + 0.1 * np.cos(25*t)

# Perform low-pass filtering
b, a = signal.butter(8, 0.125)
low1 = signal.filtfilt(b, a, sig_noisy, padlen=50)
low2 = signal.filtfilt(b, a, sig_noisy, method="gust")

# Plot the results
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, sig, 'b-', lw=3, alpha=0.25, label='true signal')
ax.plot(t, sig_noisy, 'b-', lw=0.5, label='noisy signal')
ax.plot(t, low1, 'r-', label='pad')
ax.plot(t, low2, 'g-', label='gust')
ax.legend()

plt.show()

[WarrenWeckesser]

@rgommers: Thanks. I'll grab your branch in a couple days.

I originally developed this as utility outside of scipy, which is why there is a separate test file. I can merge the tests into test_signaltools.py. The private function returns x0 and x1 for testability; I'll add a note about that to the docstring of the private function.

I think I've been using single backticks for years as a styling markup to separate Python text from English text, and used it with names, code snippets, etc.. @rkern recently posted this in a comment on the differential evolution pull request:

Specifically, single-backticks are for cross-referencing names of functions/classes/etc. that are defined elsewhere. Here, x and args should also get double-backticks since they aren't cross-references even though they are just single identifiers.

I think that is just a summary of how Sphinx handles backticks. But in fact, the examples in the "Parameters" section of the guide to documentation (https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt) show single ticks around the variable names. So if single ticks are OK for names that we don't expect to be cross-referenced by Sphinx, why not use single ticks for small code snippets such as x + 1 or a[0]? Well, whatever the "correct" use of backticks is, it would be helpful to get it written down in the numpy/scipy documentation guide .

On the use of from scipy import signal, and referring to the function as signal.filtfilt: I've been working under the assumption that we don't need an import for the function being documented. Requiring an explicit import is fine with me, but this too should be part of the documentation guide.

I guess I'll start working on a pull request for the documentation guide. :)

[rgommers]

_filtfilt_gust is kind of non-trivial to understand in detail, so I made sure that I understood filtfilt_opt in the tests and that the test coverage is complete. That looks good to me.

[rgommers]

IIRC the parameters that get referenced with single backticks were supposed to be linked by numpydoc, therefore we don't treat them as code. Support in numpydoc wasn't completed I think. It's clear though that code snippets should be treated as code (so double backticks) - grey background looks better than italics anyway.

Re signal.filtfilt: the module.funcname style has been preferred for a long time. Looking forward to your doc PR:)

[WarrenWeckesser]

I removed the 0.15.0 milestone. There are other enhancements and bugs with a higher priority than this to work on in the short term. It can wait until 0.16.

[rgommers]

@WarrenWeckesser I've updated https://github.com/rgommers/scipy/tree/gust with a rebase on current master and moved the release notes additions to 0.16.0. I think this PR needs only a few tweaks to complete (move tests and extend example like I suggested above). So maybe good to do that sometime soon?

[WarrenWeckesser]

Rebased with @rgommers's suggested changes.

[rgommers]

thanks to the high commit rate in the holiday season, needs another rebase....

[rgommers]

These three little dots are invalid syntax in this section; interpreted as a doctest line continuation. Removing them makes the plot render.

[rgommers]

nice plot by the way

[rgommers]

This is good to go after the ... is fixed I think.

[ewmoore]

Should there be a versionadded somewhere? Say in the description of the choices to method?

[rgommers]

@ewmoore good point

[WarrenWeckesser]

Should there be a versionadded somewhere? Say in the description of the choices to method?

Yes, thanks for the reminder. I like keeping "version added" notes in the "Notes" section, but I've seen them added in the parameter descriptions, too. We don't have a standard for this.

[WarrenWeckesser]

I added a Notes section to the filtfilt docstring with a comment about adding Gustafsson's method in 0.16.0.

[rgommers]

Looks good, in it goes. Thanks Warren!