RFC: filt improvements

[mbauman]

I've made a handful of improvements and extensions to the filt function on the way to adding filtfilt to DSP.jl.

• This redefines filt as filt! with a few modifications to perform its filtering in-place. filt is defined as a simple wrapper, copying the vectors that filt! modifies.
• Both filt and filt! now accept an optional fourth argument argument for the initial filter state vector. This vector is modified in-place in filt!
• filt! takes an optional fifth argument for its output array. This can alias the input, and it does so by default.
• The default filter state is now automatically computed such that its response to a step function is steady state.

Here's an example of why allowing that third part is important. This is a 5 pole low-pass butterworth filter from DSP.jl. In red is the old implementation an initial state of 0s, and in blue is the new default a matched initial state:

A few questions:

• This changes the output of filt. I think that this is the "right" default. But it will change results. Is that ok? I am not a digital filter guru; I just use them a lot. Does this make sense to do in other applications? Edit: Hm. Thinking about this more, I don't think this should happen in the base filt function. I can roll back the second commit, and it can become part of a more domain-specific package, like DSP.jl.
• Licensing: I was looking at SciPy's implementation as I worked on this. But the only part that actually comes from SciPy is the 3-line algorithm for solving the initial state. SciPy is 3-clause BSD, but surely that algorithm is published somewhere and/or not copyrightable? (only applies to 19ba34c.)
• See my one inline question on copy.
• If approved, I'll add NEWS and update the doc.

[timholy]

I haven't looked this over in detail, but just a couple of comments:

• I like the ability to filter in-place. Another alternative, however, is to allow a specific output array to be passed in. When the user wants to filter in-place, one could use the same array for the input and the output. A note in the documentation for filt! could point out that this is permissible.
• These changes in the initialization are cool, but I agree that changing the default results of filt would be problematic. Better as a separate function or placed in a package. (I kind of wonder if filt should just live in DSP.)

[mbauman]

Ah, yes, I got a little overzealous in setting the initial conditions. Setting them like this is a common thing to do in filtfilt, but both SciPy's lfilterand Matlab's filter seem to behave as though they start from an initial zero state. Which makes sense to me — the filter need not be stable if you're only going one direction. I think the first commit is still very much worth doing as it'd enable changing the filter state for filtfilt.

[StefanKarpinski]

@mbauman, this is off-topic, but what are your in-terminal Gadfly settings? These plots look great. I have the iTerm plotting setup but I haven't tweaked the settings at all so it still has a light background.

[timholy]

@mbauman, absolutely, setting filter state would be a big improvement. I would just use a default that preserves current behavior.

[mbauman]

PR updated, should be good to review at this point.

@StefanKarpinski Here's my .juliarc.jl. The Gadfly theme and its swapping are pilfered from @Keno. It's a little messy, but it works for now.

[simonster]

This looks good to me, but I think there may be some value to @timholy's suggestion of making filt! accept separate input and output arrays (which could potentially be the same array). My benchmarks show that the implementation in this PR is actually already slightly faster than the current implementation because val can be stored in a register (otherwise LLVM needs to consider the possibility that y and si are aliased), but if I make the same change to the current implementation, that is faster still (by ~10% with a long vector), because there is no need to copy x to the output array before filtering. You can probably also gain some performance by putting x[i] in a variable since that would allow it to be stored in a register.

[timholy]

Exactly. The point being that calling copy forces you to make one traversal through the input, with all the consequent cache-misses, before you make the "real" traversal through the input to do the filtering (incurring a second round of cache-misses). In contrast, if filt allocates a blank output array and passes it to filt!, then you only make one traversal through the input.

When the output array is different from the input array, you might think there should be two cache misses because you're traversing two arrays. At least for the experiments I've run, that doesn't seem to cost you any extra time---I think hardware prefetching is getting pretty good 😄, and it fetches the two regions in parallel, so you only pay the price of a cache miss once. Might be worth running some timing experiments to check this explicitly, though, before you go to a lot of effort to rework this.

[mbauman]

Ah, very clever. I didn't quite catch on to the power of using two arguments for this from the first comment. Thank you both for running the performance tests and your advice.

I have updated the PR with your suggestions, and it works wonderfully. I no longer modify the state vector in-place to minimize surprises (and it makes for a simpler API/doc).

[timholy]

If we get it along with this PR, I don't see a need to make it a separate PR.

[simonster]

Can we move out to be the first argument? I understand that then it can't be optional, but I think it's worth it for consistency, since most mutating functions in Base take the output as their first argument.

[mbauman]

Sure thing. This strikes me as a good use-case for keyword arguments as their overhead will likely be insignificant compared to the actual filtering. I'd like to preserve the ability to call filt!(b,a,x) (which I realize bucks that convention… we can ax it if you like). See my update — does that look good to you?

[mbauman]

Thanks for the guidance, @simonster and @timholy.

[simonster]

Thanks for implementing this!

[timholy]

👍

[daddesio]

@mbauman

I no longer modify the state vector in-place to minimize surprises (and it makes for a simpler API/doc).

Could we bring that feature back? It would be useful to be able to chain function calls like such:

filter_state = zeros(filter_order);

# First audio frame
decode_frame!(b, a, excitation)
filt!(sub(output, :, 1), b, a, excitation, filter_state)

# Second audio frame
decode_frame!(b, a, excitation)
filt!(sub(output, :, 2), b, a, excitation, filter_state)

# Third audio frame
decode_frame!(b, a, excitation)
filt!(sub(output, :, 3), b, a, excitation, filter_state)

# ...

In linear predictive coding, one typically varies the vocal tract parameters (and the resulting IIR filter) every 10-20ms.

If you don't tell me what to put in filter_state, I have to reproduce it from scratch by refiltering x[n] for the last filter_order samples. Alternatively, instead of having filt accept the state of the circular buffer, it could accept the last inputs and outputs; or (if we don't like that) we could introduce a new function which creates a state for you to pass to filt, a la Matlab's filtic.

By the way, excellent work. :)

[mbauman]

Interesting, and sounds reasonable. Do you mind opening a new issue with this request? I'm not sure if anyone is depending upon the initial state vector not getting modified, and I've not done much with DSP lately. DSP.jl has filt_stepstate, and I bet they'd be happy to add a filtic-like function there.