PanieriLorenzo/rust-poly

Reverse Bessel roots difference v0.2.0/scipy and main

Closed this issue ยท 5 comments

Hi,
Thank you for this library! I am trying to use it for filter design and stumbled over some differences between the main branch and the latest release on crates.io (v0.2.0). The results for finding the roots of reverse Bessel are different between those versions:

On v0.2.0 the following:

dbg!(Poly64::reverse_bessel(2).unwrap().roots());

evaluates to:

 [
    Complex {
        re: -1.5,
        im: -0.8660254037844386,
    },
    Complex {
        re: -1.5,
        im: 0.8660254037844386,
    },
]

Similar to scipy:

from scipy.signal._filter_design import (
    _bessel_zeros,
)
p = 1/_bessel_zeros(2)
print(p)
[-1.5+0.866025403784439j -1.5-0.866025403784439j]

However, on the current main branch for:

    let num_poles = 2;
    let roots = Poly64::reverse_bessel(num_poles)
        .unwrap()
        .try_roots(1E-14, 1000, 1, None, None, None)
        .unwrap();
    dbg!(&roots);

evaluates to:

 &roots = [
    Complex {
        re: -0.7912878474779199,
        im: 0.0,
    },
    Complex {
        re: 3.7912878474779195,
        im: 0.0,
    },
]

Is there something that I am missing?

You are right, that is definitely weird. This is a bug.

I am currently re-writing the root finding because I was not satisfied with it, so keep in mind that the main branch is the development branch and is probably much less stable than the released versions.

But thanks for checking it out and I'll release a fix soon :) I'd suggest you use 0.2.0 for now.

If using 0.2.0 is not an option for you, I would suggest that you use try_n_roots for now, which seems to work, at least with Newton algorithm

Thank you, I'll try that!

I just tried the version from the main branch because I couldn't reproduce the exact same values that I obtained from scipy with version 0.2.0. I am testing my filter designs against python, so this is how I noticed it.

Are you interested in adding tests against the python equivalent functions? I could help here and do a pull request.

That would be awesome, thanks! Go ahead :) Please follow conventional commits as I use it to generate the CHANGELOG automatically. If you npm install it sets up git hooks that check it for you

found a bug in crate::poly::base::Poly::companion where I forgot a minus sign when calculating the companion matrix, should work now