Elliptic filter

[meiyasan]

Hello,

Is such elliptic filter implemented in KFR DSP ?
I couldn't find any hint from source code, if so is there some plan for that ?

[dancazarin]

Hello,

It is not implemented as of now and yes, there are plans to implement it in future releases (KFR 5.x).
Elliptic filter's special cases, the butterworth filter as well as Chebyshev I and II are all implemented since 4.0.

[meiyasan]

Hello !

I have been doing some work around this elliptic filter as I needed a generalized version of it.
If it is fine for you I would like to share the outcome in a side branch to be improved.

Clang++ is not provided the special mathematical functions like elliptic integral (std::comp_ellip1, etc..).
I tried to implement two versions, but there are some double precision limitations in both cases. I hope this might be improved in terms of kfr style implementation and in term of dependencies.

[meiyasan]

I just pushed some thinking. The mainline is based on SciPy package (elllipk, ellipkm1) methods.
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipk.html
Clang++ do not implement any elliptic integral at this time, AFAIK.

I pick up the calculation of angles from : https://people.math.sc.edu/Burkardt/cpp_src/elliptic_integral/elliptic_integral.html
This library looked good, but it seems it have an issue with m~1.
So I included the real calculation of the integral from CEPHES library (elliptic_k, elliptic_km1).
This is using a polynomial approximation with precomputed coefficients. So I think this should be quite fast.

Nevertheless, my implementation still lag of calculation for the inverse jacobian integral due to double precision limitation
(https://github.com/xKZL/kfr/blob/d315de7daa3889de0ef763daacc8409191aa2dcc/include/kfr/math/elliptic.hpp#L253)
I just point where the issue start to appear, but I think we might try to replace inv_jacobi_sn to fix the issue
I have no more time to dedicate, but please feel free to try something.. this would be extremely useful in the future for my usage.

To put this in a nutshell, here is the outcome of my test below:

• This is comparing the k coefficient from zpk representation.
• I have compared different orders ( 0 and 1 are special cases so they are safe)
• As the order increase, the precision limitation starts to increase and the result is very completely limited in the last orders

Test #0 : 0.89125093813374556273032567332848 (python SciPy/stdc++) <> 0.89125093813374556273032567332848 = 1 (1e-0)
Test #1 : 1.96522672836027156861860021308530 (python SciPy/stdc++) <> 1.96522672836027156861860021308530 = 1 (1e-0)
Test #2 : 0.50118723362727279901918109317194 (python SciPy/stdc++) <> 0.50118723362727179981845893053105 = 0 (1e-16)
Test #3 : 0.98794734469230960360874860270997 (python SciPy/stdc++) <> 0.98794734469230882645263136510038 = 0 (1e-16)
Test #4 : 0.50118723362727934933502638159553 (python SciPy/stdc++) <> 0.50118723362725281500473784035421 = 0 (1e-15)
Test #5 : 0.97538440598394993141795339397504 (python SciPy/stdc++) <> 0.97538440598390829805452995060477 = 0 (1e-15)
Test #6 : 0.50118723362745221105996051846887 (python SciPy/stdc++) <> 0.50118723362739070470439628479653 = 0 (1e-14)
Test #7 : 0.97511086984454220516482791936141 (python SciPy/stdc++) <> 0.97511086985267536597632442862959 = 0 (1e-12)
Test #8 : 0.50118723362753969663430098080426 (python SciPy/stdc++) <> 0.50118723359427075347838353991392 = 0 (1e-11)
Test #9 : 0.97510485666884338940008092322387 (python SciPy/stdc++) <> 0.97510484997981838883873706436133 = 0 (1e-9)
Test #10: 0.50118723375847462619958605500869 (python SciPy/stdc++) <> 0.50118672696715338421569185811677 = 0 (1e-7)
Test #11: 0.97510467362126196366745034538326 (python SciPy/stdc++) <> 0.97502524988276029205280792666599 = 0 (1e-5)
Test #12: 0.50118498126436306083775207298459 (python SciPy/stdc++) <> 0.49903983562735582113489840594411 = 0 (1e-2)
Test #13: 0.97482310868854160634811023555812 (python SciPy/stdc++) <> 0.88735287300246623587440808478277 = 0 (1e-2)
Test #14: 0.49486370015702924041178789593687 (python SciPy/stdc++) <> 0.32953759545438848777010321100533 = 0 (1e-2)
Test #15: 0.83465109606081722137815859241528 (python SciPy/stdc++) <> 0.49342000704387722898047741182381 = 0 (1e-2)
Test #16: 0.30887740092250570711485124775209 (python SciPy/stdc++) <> 0.14213035309011415319169202575722 = 0 (1e-2)
Test #17: 0.48361023745374215332404332912120 (python SciPy/stdc++) <> 0.24555920485643306649947703590442 = 0 (1e-2)
Test #18: 0.14515935341003305403262402251130 (python SciPy/stdc++) <> 0.05455482844115393248340950549391 = 0 (1e-2)
Test #19: 0.26007573487724366945172960186028 (python SciPy/stdc++) <> 0.11009230039696740743870861933828 = 0 (1e-2)
Test #20: 0.06214721652055405637371521265777 (python SciPy/stdc++) <> 0.01915810306404193072427055710704 = 0 (1e-3)
Test #21: 0.12856256785699293754277050538803 (python SciPy/stdc++) <> 0.04456440191168813819144745025369 = 0 (1e-2)
Test #22: 0.02471231445727731929062898075244 (python SciPy/stdc++) <> 0.00621417411672309152187443359594 = 0 (1e-3)
Test #23: 0.05842982129840767341333318540819 (python SciPy/stdc++) <> 0.01641719101767708313688309829103 = 0 (1e-3)
Test #24: 0.00919341192048158013794267873209 (python SciPy/stdc++) <> 0.00187271832797217328663019753065 = 0 (1e-4)
Test #25: 0.02453838324842558954452798047896 (python SciPy/stdc++) <> 0.00554806695055968385199562931120 = 0 (1e-3)
Test #26: 0.00321322302422841340682757582670 (python SciPy/stdc++) <> 0.00052696377910764501985296792696 = 0 (1e-4)
Test #27: 0.00957488847253192346120620470629 (python SciPy/stdc++) <> 0.00173220837320768204434240367106 = 0 (1e-4)
Test #28: 0.00105876470543555894175680176517 (python SciPy/stdc++) <> 0.00013908931851696425111876431746 = 0 (1e-4)
Test #29: 0.00348895553891140223004563303277 (python SciPy/stdc++) <> 0.00050281167427838478929669197015 = 0 (1e-4)
Test #30: 0.00032993230118648272106499086398 (python SciPy/stdc++) <> 0.00003458250633305212775941636649 = 0 (1e-5)
Test #31: 0.00119263844932737412413148447854 (python SciPy/stdc++) <> 0.00013645609690712717548995158711 = 0 (1e-4)

[meiyasan]

Hello, I haven't managed to get rid of this precision issue.
Any help from an KFR expert would be great ! @dancazarin

That would be really great if this could be implemented in a near future.

[dancazarin]

Hello,
Original CEPHES library has no license specified, I don't think it's possible to use it with KFR without licensing issues. Author has explicitly granted rights to some projects using CEPHES but still hasn't added license details to the source code available on his website.
Boost has elliptic functions too, as well as libstdc++ 7 and msvc 19.14.
Maybe it's a good option to implement elliptic filters in KFR but require implementation of the special math functions to be available under std or boost namespace.

[meiyasan]

It seems to be MIT license (original?)
https://en.smath.com/view/CephesMathLibrary/license

In the scope of my work, I would like to avoid duplicating libraries. Especially boost if I "just" need to compute elliptic. A native solution would be better for sure

[meiyasan]

@dancazarin ?

[meiyasan]

Sorry for the spam, but do you have any opinion about my previous message ? @dancazarin

[dancazarin]

smath.com is a C# wrapper so MIT license is not the license of Cephes library. The original is at https://www.netlib.org/cephes/ and still doesn't mention any license. In strict words, lack of license means no rights to use the code without explicit permission for every project. Scipy team has asked explicit permission from the author to use Cephes in scipy, so did other developers. As an option you may add Cephes to your own fork of KFR if you feel comfortable with using unlicensed (private) code. Also, as I already wrote, it's possible to implement that way: use elliptic functions from boost or std (actual STL implementations already have it) or disable elliptic filters if they cannot be found. Are your toolchain support elliptic functions in std namespace?

[meiyasan]

ChatGPT says: "As of my last knowledge update in September 2021, the Cephes library, which provides mathematical functions including elliptic integrals, was typically distributed under a permissive open-source license called the "Cephes Mathematical Library License." This license is quite permissive and allows for both commercial and non-commercial use with minimal restrictions."

I am looking for the references now.. ahah

[meiyasan]

Perhaps we might ask its author :) (steve@moshier.net)
https://github.com/google-deepmind/torch-cephes/blob/master/LICENSE.txt

[dancazarin]

The author definitely has better knowledge of the license details of his library compared to ChatGPT.