berndporr / iirj

An efficient IIR filter library written in JAVA
Apache License 2.0
136 stars 36 forks source link

Unusual results when using Bessel filters of higher orders #25

Closed anleu closed 1 year ago

anleu commented 1 year ago

We've observed that the results from IIRJ align closely with Scipy for most filters. However, for Bessel Filters of order greater than 1, there are significant differences. Visually, the Bessel results from IIRJ appear inconsistent. Given that other filters produce similar results to Scipy, we suspect this isn't a configuration issue, though we can't be entirely certain. Any insights or guidance on this matter would be highly appreciated.

Compare.zip

berndporr commented 1 year ago

Bessel filters have no analytical solution and the code relies on finding the roots of the polynomials. org.apache.commons.math3.analysis.solvers.LaguerreSolver is used for this. As you know it's based on C++ code and converted to java. The C++ code had a buggy root finder and Bessel has been removed a while ago. Looks like the root finder of apache commons has an issue as well (or the Bessel function has been buggy but worked till recently as the roots were not affected). I guess the best plan of action is to remove Bessel completely and rather import the coefficients calculated by Scipy. That's also the way to do it for C++. I certainly have no time to get to the bottom of it as it's non-trivial and takes a lot of time.

berndporr commented 1 year ago

Of course I'd be more than happy to accept pull requests which make sure the impulse responses match. Ideally also the unit tests should receive coefficients from python for comparison.

berndporr commented 1 year ago

Thanks for cross checking though!

anleu commented 1 year ago

It seems that Scipy doesn't use a solver but instead has fixed results for the poles. https://github.com/scipy/scipy/blob/f2ec91c4908f9d67b5445fbfacce7f47518b35d1/scipy/signal/filter_design.py#L2351

Perhaps this is also an approach for IIRJ? I tried it with the given example of order 5. It looks slightly better, but there are still differences. I don't have time at the moment to investigate further. I'll try to take another look in two weeks.

Bessel Order 5 Compare

anleu commented 1 year ago

I will prepare the PR after https://github.com/berndporr/iirj/issues/27 is merged.