From 75918016ab16cd69973e818759f027e78022dee4 Mon Sep 17 00:00:00 2001 From: Michael Niedermayer Date: Mon, 22 May 2023 23:34:37 +0200 Subject: [PATCH] Move bessel_i0() from swresample/resample to avutil/mathematics 0th order modified bessel function of the first kind are used in multiple places, lets avoid having 3+ different implementations I picked this one as its accurate and quite fast, it can be replaced if a better one is found Signed-off-by: Michael Niedermayer --- doc/APIchanges | 3 ++ libavutil/mathematics.c | 104 ++++++++++++++++++++++++++++++++++++++ libavutil/mathematics.h | 4 ++ libavutil/version.h | 2 +- libswresample/resample.c | 106 +-------------------------------------- 5 files changed, 113 insertions(+), 106 deletions(-) diff --git a/doc/APIchanges b/doc/APIchanges index 2706afcf9d..0232d05b10 100644 --- a/doc/APIchanges +++ b/doc/APIchanges @@ -2,6 +2,9 @@ The last version increases of all libraries were on 2023-02-09 API changes, most recent first: +2023-05-29 - xxxxxxxxxx - lavu 58.12.100 - mathematics.h + Add av_bessel_i0() + 2023-05-xx - xxxxxxxxxx - lavc 60.15.100 - avcodec.h Add AVHWAccel.update_thread_context, AVHWAccel.free_frame_priv, AVHWAccel.flush. diff --git a/libavutil/mathematics.c b/libavutil/mathematics.c index b878317d63..61aeb7c029 100644 --- a/libavutil/mathematics.c +++ b/libavutil/mathematics.c @@ -213,3 +213,107 @@ int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t i return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts); } } + +static inline double eval_poly(const double *coeff, int size, double x) { + double sum = coeff[size-1]; + int i; + for (i = size-2; i >= 0; --i) { + sum *= x; + sum += coeff[i]; + } + return sum; +} + +/** + * 0th order modified bessel function of the first kind. + * Algorithm taken from the Boost project, source: + * https://searchcode.com/codesearch/view/14918379/ + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0 (see notice below). + * Boost Software License - Version 1.0 - August 17th, 2003 +Permission is hereby granted, free of charge, to any person or organization +obtaining a copy of the software and accompanying documentation covered by +this license (the "Software") to use, reproduce, display, distribute, +execute, and transmit the Software, and to prepare derivative works of the +Software, and to permit third-parties to whom the Software is furnished to +do so, all subject to the following: + +The copyright notices in the Software and this entire statement, including +the above license grant, this restriction and the following disclaimer, +must be included in all copies of the Software, in whole or in part, and +all derivative works of the Software, unless such copies or derivative +works are solely in the form of machine-executable object code generated by +a source language processor. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. + */ + +double av_bessel_i0(double x) { +// Modified Bessel function of the first kind of order zero +// minimax rational approximations on intervals, see +// Blair and Edwards, Chalk River Report AECL-4928, 1974 + static const double p1[] = { + -2.2335582639474375249e+15, + -5.5050369673018427753e+14, + -3.2940087627407749166e+13, + -8.4925101247114157499e+11, + -1.1912746104985237192e+10, + -1.0313066708737980747e+08, + -5.9545626019847898221e+05, + -2.4125195876041896775e+03, + -7.0935347449210549190e+00, + -1.5453977791786851041e-02, + -2.5172644670688975051e-05, + -3.0517226450451067446e-08, + -2.6843448573468483278e-11, + -1.5982226675653184646e-14, + -5.2487866627945699800e-18, + }; + static const double q1[] = { + -2.2335582639474375245e+15, + 7.8858692566751002988e+12, + -1.2207067397808979846e+10, + 1.0377081058062166144e+07, + -4.8527560179962773045e+03, + 1.0, + }; + static const double p2[] = { + -2.2210262233306573296e-04, + 1.3067392038106924055e-02, + -4.4700805721174453923e-01, + 5.5674518371240761397e+00, + -2.3517945679239481621e+01, + 3.1611322818701131207e+01, + -9.6090021968656180000e+00, + }; + static const double q2[] = { + -5.5194330231005480228e-04, + 3.2547697594819615062e-02, + -1.1151759188741312645e+00, + 1.3982595353892851542e+01, + -6.0228002066743340583e+01, + 8.5539563258012929600e+01, + -3.1446690275135491500e+01, + 1.0, + }; + double y, r, factor; + if (x == 0) + return 1.0; + x = fabs(x); + if (x <= 15) { + y = x * x; + return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y); + } + else { + y = 1 / x - 1.0 / 15; + r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y); + factor = exp(x) / sqrt(x); + return factor * r; + } +} diff --git a/libavutil/mathematics.h b/libavutil/mathematics.h index 3e7b52eed6..e213bab68c 100644 --- a/libavutil/mathematics.h +++ b/libavutil/mathematics.h @@ -288,6 +288,10 @@ int64_t av_rescale_delta(AVRational in_tb, int64_t in_ts, AVRational fs_tb, int */ int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t inc); +/** + * 0th order modified bessel function of the first kind. + */ +double av_bessel_i0(double x); /** * @} diff --git a/libavutil/version.h b/libavutil/version.h index 032d0c294d..0d8434493f 100644 --- a/libavutil/version.h +++ b/libavutil/version.h @@ -79,7 +79,7 @@ */ #define LIBAVUTIL_VERSION_MAJOR 58 -#define LIBAVUTIL_VERSION_MINOR 11 +#define LIBAVUTIL_VERSION_MINOR 12 #define LIBAVUTIL_VERSION_MICRO 100 #define LIBAVUTIL_VERSION_INT AV_VERSION_INT(LIBAVUTIL_VERSION_MAJOR, \ diff --git a/libswresample/resample.c b/libswresample/resample.c index 8f9efc3f21..bd54a7002f 100644 --- a/libswresample/resample.c +++ b/libswresample/resample.c @@ -30,110 +30,6 @@ #include "libavutil/cpu.h" #include "resample.h" -static inline double eval_poly(const double *coeff, int size, double x) { - double sum = coeff[size-1]; - int i; - for (i = size-2; i >= 0; --i) { - sum *= x; - sum += coeff[i]; - } - return sum; -} - -/** - * 0th order modified bessel function of the first kind. - * Algorithm taken from the Boost project, source: - * https://searchcode.com/codesearch/view/14918379/ - * Use, modification and distribution are subject to the - * Boost Software License, Version 1.0 (see notice below). - * Boost Software License - Version 1.0 - August 17th, 2003 -Permission is hereby granted, free of charge, to any person or organization -obtaining a copy of the software and accompanying documentation covered by -this license (the "Software") to use, reproduce, display, distribute, -execute, and transmit the Software, and to prepare derivative works of the -Software, and to permit third-parties to whom the Software is furnished to -do so, all subject to the following: - -The copyright notices in the Software and this entire statement, including -the above license grant, this restriction and the following disclaimer, -must be included in all copies of the Software, in whole or in part, and -all derivative works of the Software, unless such copies or derivative -works are solely in the form of machine-executable object code generated by -a source language processor. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT -SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE -FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, -ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -DEALINGS IN THE SOFTWARE. - */ - -static double bessel(double x) { -// Modified Bessel function of the first kind of order zero -// minimax rational approximations on intervals, see -// Blair and Edwards, Chalk River Report AECL-4928, 1974 - static const double p1[] = { - -2.2335582639474375249e+15, - -5.5050369673018427753e+14, - -3.2940087627407749166e+13, - -8.4925101247114157499e+11, - -1.1912746104985237192e+10, - -1.0313066708737980747e+08, - -5.9545626019847898221e+05, - -2.4125195876041896775e+03, - -7.0935347449210549190e+00, - -1.5453977791786851041e-02, - -2.5172644670688975051e-05, - -3.0517226450451067446e-08, - -2.6843448573468483278e-11, - -1.5982226675653184646e-14, - -5.2487866627945699800e-18, - }; - static const double q1[] = { - -2.2335582639474375245e+15, - 7.8858692566751002988e+12, - -1.2207067397808979846e+10, - 1.0377081058062166144e+07, - -4.8527560179962773045e+03, - 1.0, - }; - static const double p2[] = { - -2.2210262233306573296e-04, - 1.3067392038106924055e-02, - -4.4700805721174453923e-01, - 5.5674518371240761397e+00, - -2.3517945679239481621e+01, - 3.1611322818701131207e+01, - -9.6090021968656180000e+00, - }; - static const double q2[] = { - -5.5194330231005480228e-04, - 3.2547697594819615062e-02, - -1.1151759188741312645e+00, - 1.3982595353892851542e+01, - -6.0228002066743340583e+01, - 8.5539563258012929600e+01, - -3.1446690275135491500e+01, - 1.0, - }; - double y, r, factor; - if (x == 0) - return 1.0; - x = fabs(x); - if (x <= 15) { - y = x * x; - return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y); - } - else { - y = 1 / x - 1.0 / 15; - r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y); - factor = exp(x) / sqrt(x); - return factor * r; - } -} - /** * builds a polyphase filterbank. * @param factor resampling factor @@ -189,7 +85,7 @@ static int build_filter(ResampleContext *c, void *filter, double factor, int tap break; case SWR_FILTER_TYPE_KAISER: w = 2.0*x / (factor*tap_count*M_PI); - y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0))); + y *= av_bessel_i0(kaiser_beta*sqrt(FFMAX(1-w*w, 0))); break; default: av_assert0(0);