From patchwork Mon May 22 23:35:56 2023 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Michael Niedermayer X-Patchwork-Id: 41783 Delivered-To: ffmpegpatchwork2@gmail.com Received: by 2002:a05:6a20:c51c:b0:10c:5e6f:955f with SMTP id gm28csp175116pzb; Mon, 22 May 2023 16:36:14 -0700 (PDT) X-Google-Smtp-Source: ACHHUZ5Qt3DuZhjMC+9vM+ra6DOkP9HyFoptJAvqjDPEWsrdBq8vI+xXq3EqsMO3AeB9CcEhItku X-Received: by 2002:a17:907:6e29:b0:969:813c:9868 with SMTP id sd41-20020a1709076e2900b00969813c9868mr13861695ejc.18.1684798574616; Mon, 22 May 2023 16:36:14 -0700 (PDT) ARC-Seal: i=1; a=rsa-sha256; t=1684798574; cv=none; d=google.com; s=arc-20160816; b=TDPgJERP8/6UeJK8lGumt2VVeiEKeMe6CozERToSG8UywGCzmccSUJwDSpWjixAR8H iu077oTrUREMGUp6WKRekexXgBDxrHy+aWyZpW4aKOvtYcYRjt8yH3+5sTUGPkDCKtr5 lx1XYjc0nc2g5Fjil7PhfYBJST+IemDvF+0W0TIb8XToeH5C/uXT7ij3ePNmuBpqp69T UYJwFWaNs3klWq3e+I3MhoolZ+6Ba34vJTSmAB3bnRIv6o2b9IIkQoD1ky5CiJ/paHfH gTIlIKAxGe9cXCNlSu6QhDwPfvxTaRerpL9zFT3RPpO+7UuZMmCdEu5j2Tg/bzh+ur/R OuqQ== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=arc-20160816; h=sender:errors-to:content-transfer-encoding:mime-version:reply-to :list-subscribe:list-help:list-post:list-archive:list-unsubscribe :list-id:precedence:subject:message-id:date:to:from:delivered-to; bh=iVnkDBe8eft+DkXE5geWU2REdazaAxSI0GA6XYFso98=; b=m3ZpBUk6GHV7dBvknLjr9Swrya4hIA68ilxuOKuRS3+YiuCFRu+koxR66cYarQ9Klx uK5pgYEckEb1FvcLpd38yzjPjOX3yGrO/9jMxzGqQHtmVFw6nAxxiM6yPq913oI0gbtP 9+E4OkxHSoFnfhkHTmn+119ic6A0dmGAibF/QdgKvJjHlslbSVgU9DI5TNSrDVlVbyJd el0kpgsnzEMcuAFxqvvxDyagUyIvpWDvnsaPnyxfmwrqGkfGHbJV7+wwfeVJ4tPidp1V ee9jNiQCLsr4jJ/ieCOs2NRKHgDoxjMfU/ZGoi02zDd4VNTZ0W2lI12+dtvuu71eT56u S57w== ARC-Authentication-Results: i=1; mx.google.com; spf=pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) smtp.mailfrom=ffmpeg-devel-bounces@ffmpeg.org Return-Path: Received: from ffbox0-bg.mplayerhq.hu (ffbox0-bg.ffmpeg.org. [79.124.17.100]) by mx.google.com with ESMTP id w22-20020a170907271600b0096f231529f1si1209599ejk.870.2023.05.22.16.36.13; Mon, 22 May 2023 16:36:14 -0700 (PDT) Received-SPF: pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) client-ip=79.124.17.100; Authentication-Results: mx.google.com; spf=pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) smtp.mailfrom=ffmpeg-devel-bounces@ffmpeg.org Received: from [127.0.1.1] (localhost [127.0.0.1]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTP id 9CB4D68C08B; Tue, 23 May 2023 02:36:09 +0300 (EEST) X-Original-To: ffmpeg-devel@ffmpeg.org Delivered-To: ffmpeg-devel@ffmpeg.org Received: from relay4-d.mail.gandi.net (relay4-d.mail.gandi.net [217.70.183.196]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTPS id E066A68BBA2 for ; Tue, 23 May 2023 02:36:02 +0300 (EEST) Received: (Authenticated sender: michael@niedermayer.cc) by mail.gandi.net (Postfix) with ESMTPSA id 1691DE0003 for ; Mon, 22 May 2023 23:36:01 +0000 (UTC) From: Michael Niedermayer To: FFmpeg development discussions and patches Date: Tue, 23 May 2023 01:35:56 +0200 Message-Id: <20230522233601.20552-1-michael@niedermayer.cc> X-Mailer: git-send-email 2.17.1 Subject: [FFmpeg-devel] [PATCH 1/6] Move bessel_i0() from swresample/resample to avutil/mathematics X-BeenThere: ffmpeg-devel@ffmpeg.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: FFmpeg development discussions and patches List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Reply-To: FFmpeg development discussions and patches MIME-Version: 1.0 Errors-To: ffmpeg-devel-bounces@ffmpeg.org Sender: "ffmpeg-devel" X-TUID: vDgXRIkBvJG/ 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 ++ libswresample/resample.c | 106 +-------------------------------------- 4 files changed, 112 insertions(+), 105 deletions(-) diff --git a/doc/APIchanges b/doc/APIchanges index 130955ef02..b256df9c54 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-xx-xx - xxxxxxxxxx - lavu 58.x.100 - mathematics.h + Add av_bessel_i0() + 2023-05-xx - xxxxxxxxxx - lavu 58.8.100 - frame.h Add av_frame_replace(). 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/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);