Zero phase digital filtering

[rob-luke]

filtfilt function and associated tests are included.

This code is merely a wrapper for the work done by @mbauman

[rob-luke]

Code coverage decreases on coveralls because I only do a cursory test of the filt! and filt_stepstate

[coveralls]

Coverage decreased (-3.28%) when pulling c82c630 on codles:master into 782776f on JuliaDSP:master.

[coveralls]

Coverage decreased (-2.82%) when pulling 1e24269 on codles:master into 782776f on JuliaDSP:master.

[mbauman]

Thanks! This is great. It's awesome to see that we match Matlab here… and good on you for digging through my dirty trash to pull out the filt_stepstate function — that won't become a part of base Julia, and will stay in DSP.jl for good.

I just updated my PR over at Julia to improve the performance of filt. Even when that gets merged, though, I suppose we'll need an implementation here to allow folks with older versions of Julia to use DSP.jl. I can also help you make sure you hit all the code paths in the testing — the biggest one that's missing could be solved with a filter that has a scalar a.

[simonster]

reverse! should be better here since it doesn't need to allocate.

[rob-luke]

I only managed to get this to work for the 1d case.

[simonster]

This may be microoptimizing, but is it possible that we could do this without having to make two copies of x (one padded and one with the padding removed) if we could keep the filter state?

[rob-luke]

Thank you both for the feedback.

@mbauman I will change the tilt_stepstate to reflect the fact its staying in DSP.jl. Unless you want to create a separate file for it? Also I will put filt in a separate backwards compatibility file with your updates (again unless you want to do it?).

@simonster Thanks for the comments. I will make the changes you suggested.

Do we want this to work on non vector data? Perhaps columns of matrices?

[mbauman]

Do we want this to work on non vector data? Perhaps columns of matrices?

Yes, I think that'd be nice to allow. It's a pretty simple modification to the base filt! function to allow this and it incurs almost no overhead in the 1-d case… I played with adding it while I was working on the PR but I figured I'd submit that change separately to keep things simple. I have a stashed version that I'll try to submit later tonight.

It makes the initial states a bit trickier, but Julia has a nice advantage over SciPy in its array shapes - filt_stepstate(…) * x[1, :] will broadcast the initial states into a matrix (whereas SciPy needs some shape-wrangling to make it work with row-major arrays).

[mbauman]

Alright, I'm done futzing with the base implementation of filt now! Sorry to keep changing the API on you, but the latest version should be stable now. I hope. It allows for arrays and heterogeneous argument types:

filt(b, a, x, [si])
filt!(out, b, a, x, [si])  # so if you want to filter in-place, you now need: filt!(x, b, a, x)

[coveralls]

Coverage increased (+0.45%) when pulling d7e7674 on codles:master into 782776f on JuliaDSP:master.

[coveralls]

Coverage increased (+0.43%) when pulling 894cc14 on codles:master into e23dc1a on JuliaDSP:master.

[coveralls]

Coverage increased (+0.43%) when pulling 92ceb61 on codles:master into e23dc1a on JuliaDSP:master.

[simonster]

A few more thoughts on this:

• Julia 0.3 should be released soon (there is an rc now!). Perhaps rather than copy filt! into the package, we could just drop support for 0.2. @jfsantos, what do you think? If we leave the filt! implementation in, there is a typo in the name: it should be BackwardsCompatibility (or maybe Julia02?). Would also be great to only load it if VERSION < v"0.3.0-".
• filtfilt should be able to accept a Filter object as designed by FilterDesign. These two lines with filtfilt in place of filt should be sufficient.
• It would be nice to support filters designed as second-order sections without first converting them back to transfer functions. Maybe we should be designing our IIR filters as second-order sections by default, since they are more numerically stable. (This is what MATLAB does if you use its object-oriented filter API.) I can implement filt! for SOSFilter. We can determine the initial states using the same strategy, but individually for each biquad.

[jfsantos]

I think it's safe to drop support for 0.2.

I also agree about designing SOS filters directly instead of converting back and forth.

[coveralls]

Coverage increased (+0.15%) when pulling f3a9b2e on codles:master into e23dc1a on JuliaDSP:master.

[rob-luke]

Sorry for the commit noise. Should I develop this elsewhere then commit once its ready?

I will drop support for 0.2
Add filter support
Add 2d filter support

[coveralls]

Coverage increased (+0.16%) when pulling d376372 on codles:master into e23dc1a on JuliaDSP:master.

[rob-luke]

This only fails on the release not on the nightlies. Because of the lack of filt!

I think I have taken everyone's suggestion in to account. Have I missed anything?

cc: @mbauman @simonster @jfsantos

[coveralls]

Coverage increased (+0.16%) when pulling e7288f0 on codles:master into e23dc1a on JuliaDSP:master.

[simonster]

I merged this code with some tweaks to reduce the number of temporaries it needs to allocate down to the bare minimum (at least until we have a way to preserve the state vector when filtering). I'll work on SOS filtfilt next, and write up some docs.

[mbauman]

This is wonderful to have. Thanks @codles and @simonster!