ENH: Add sosfiltfilt

[larsoner]

This one is something I meant to add from back when I worked on #3717. Since #3717 (EDIT: #6137) is some ways off presumably it's probably better to get this in than wait. It would be great to have for 0.18 if possible, but not strictly necessary.

[larsoner]

@WarrenWeckesser @endolith you might be interested in this since you were involved in sosfilt and zpk2sos way back when :)

[endolith]

"cascaded"

[e-q]

This function could also check that sos[:,3] == np.ones(sos.shape[0])

[larsoner]

Thanks for the reviews, comments addressed. I think (hope?) this is actually a pretty straightforward change so it seems like merging for 0.18 is doable. @e-q do you agree?

[e-q]

Looks like the current failures are a missing entry in __init__.py, and SkipTest not importing in python 2.7 and 3.4.

[larsoner]

@rgommers do you think it's reasonable to try for 0.18 for this?

[rgommers]

@rgommers do you think it's reasonable to try for 0.18 for this?

Why not? It's not a blocker so I don't think we should delay the branching for it (up to @ev-br of course), but it seems doable to get this merged by Sunday.

[larsoner]

Great. Tests are passing now, @e-q are you happy with it now?

[WarrenWeckesser]

lfilter, filtfilt and sosfilt can filter and return complex arrays. I haven't tried it, but I don't see why this function shouldn't do the same.

[larsoner]

I copied this from filtfilt, maybe there they are restricted by the Gustafsson implementation or something. I'll add some simple tests. Maybe both this and filtfilt need to be updated...

[WarrenWeckesser]

This will be great to get into scipy.signal. I've been trying it out, and ran into an issue with odd order filters. I'm not sure this is actually a problem, and if it is a problem, I'm not sure if it is a problem with sosfiltfilt or something else.

Here's a script. When order is odd, the outputs of filtfilt and sosfiltfilt don't match.

from __future__ import print_function

import numpy as np
from scipy.signal import butter, filtfilt, sosfiltfilt, zpk2tf, zpk2sos
import matplotlib.pyplot as plt


order = 5

z, p, k = butter(order, 0.15, output='zpk')
b, a = zpk2tf(z, p, k)
sos = zpk2sos(z, p, k)

x = np.random.randn(80)
y1 = filtfilt(b, a, x)
y2 = sosfiltfilt(sos, x)

if not np.allclose(y1, y2):
    print("Filtered signals are not close.")
    print("Max difference:", np.max(np.abs(y1 - y2)))
    plt.plot(y1 - y2)
    plt.show()
else:
    print("Filtered signals are approximately the same.")

Here's the plot of y1 - y2 when order is 5:

Because the order is odd, the SOS representation has an extra zero and pole at the origin. I suspect that is the source of the problem. This might be unavoidable. I'll continue experimenting and try to get a better handle on the issue this evening.

[larsoner]

Do they match for sosfilt and lfilter? If so, what if there are initial conditions?

[WarrenWeckesser]

Do they match for sosfilt and lfilter?

Good question. When I replaced filtfilt and sosfiltfilt with lfilter and sosfilt (resp.), the filtered signals match for even and odd orders.

[larsoner]

I bet it's the pad length

[WarrenWeckesser]

Yes, it looks like it is something to do with padding. The outputs of y1 = filtfilt(b, a, x, padtype=None) and y2 = sosfiltfilt(sos, x, padtype=None) agree for even and odd orders.

[WarrenWeckesser]

filtfilt and sosfiltfilt also agree if I explicitly give them the same padlen. So the difference comes from using different default pad lengths.

I also noticed that the description of the default value of padlen in the docstring needs to be fixed. Currently it says The default value is``3 * max(len(a), len(b))``.

[WarrenWeckesser]

The filtered outputs y1 and y2 in the script agree for even and odd orders if I use

y1 = filtfilt(b, a, x)
y2 = sosfiltfilt(sos, x, padlen=6*sos.shape[0] + (0 if (order & 1) else 3))

In general, sosfiltfilt won't know if one of the sections is actually a first order section (with a pole and zero at the origin that cancel), so that formula won't work as the default. However, if we do b, a = sos2tf(sos), then len(a) is 2*sos.shape[0] + 1, so the best default padlen (for consistency with filtfilt) is probably 3*(2*sos.shape[0] + 1).

(By the way, I have no idea where the original default of 3*max(len(a), len(b)) for filtfilt comes from. It seems arbitrary. Why should the default pad length depend on the order like that?)

[WarrenWeckesser]

...so the best default padlen (for consistency with filtfilt) is probably 3_(2_sos.shape[0] + 1).

And now that I've taken another look at the code, I see that that is exactly what the default padlen in sosfiltfilt is. Sorry for the noise!

[WarrenWeckesser]

So now my only suggestions are for some tweaks to the documentation. In the description of the return value, just drop the reference to the data type of the result. The function does, in fact, return a complex array when the input is complex. And in the description of padlen, fix the description of the default value. Other than that, this looks ready to go.

[larsoner]

Sorry @WarrenWeckesser I should have saved you some time by telling you what the current behavior was, but was on mobile :(

Without looking at the math, I suspect that we get the actual order by looking at whether the third coefficient is zero. I'll have to check. Then we can get exact equivalence for odd-order systems.

[larsoner]

(in terms of padlen I mean)

[larsoner]

@WarrenWeckesser added an equivalence test for odd order (it passes, hooray!) and updated the computation and docstrings. If you have time, see if you agree.

[larsoner]

Why should the default pad length depend on the order like that?

I think this comes from MATLAB, since I know they do something similar. Octave has a 3 * (max(length(a), length(b)) - 1) in there. Unfortunately it's not a great choice, since what (I think) really matters is the effectively length of the impulse response, which depends on the specific locations of the poles and zeros rather than their count. I know people who have been burned by this heuristic, so maybe we should add a .. note:: about it, but I'd have to think about the phrasing.

[larsoner]

Although thinking about it more, having that many values ensures that your system at least has all of its delay taps depend some aspect of the input signal, beyond the single edge value (which we otherwise use to set the initial conditions), which might also be part of the motivation.

[WarrenWeckesser]

Your updated default padlen looks reasonable, so other than the comment I made about the description of the default padlen in the docstring (i.e. it is 3 times the stated amount), this looks ready to go.

[larsoner]

Thanks for another short-notice look @WarrenWeckesser, docstring fixed.

@e-q feel free to have a look in the next few hours if you have time, otherwise since people seem happy with it I'll go ahead and merge so it makes it to 0.18.

[WarrenWeckesser]

How about squashing the commits? Then I'll merge. (In the "squash commits to keep the history relative simple and clean" vs "keep every little committed tweak" debate, I'm in the "squash" camp.)

[WarrenWeckesser]

Wait, there is one more thing. Add a note about sosfiltfilt to the 0.18 release notes.

[larsoner]

Add a note about sosfiltfilt to the 0.18 release notes.

Hah! I was just doing that when you posted it. Not sure if great minds think alike or fools seldom differ here :)

I'll add it and squash and ping you when the CIs are done. Thanks for the quick reviews.

[WarrenWeckesser]

Merged. Thanks, Eric!