History

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

Michael Niedermayer committed 3 years ago

Commit 75918011 parent 0c78b0d5 files · +113 −106Browse files

Unified · Split

5 files
  1. doc/APIchanges +3 −0
  2. libavutil/mathematics.c +104 −0
  3. libavutil/mathematics.h +4 −0
  4. libavutil/version.h +1 −1
  5. libswresample/resample.c +1 −105
doc/APIchanges+3−0
@@ -2,6 +2,9 @@ The last version increases of all libraries were on 2023-02-09 · 1 unmodified line
2 2
3API changes, most recent first:3API changes, most recent first:
4 4
52023-05-29 - xxxxxxxxxx - lavu 58.12.100 - mathematics.h
6 Add av_bessel_i0()
7
52023-05-xx - xxxxxxxxxx - lavc 60.15.100 - avcodec.h82023-05-xx - xxxxxxxxxx - lavc 60.15.100 - avcodec.h
6 Add AVHWAccel.update_thread_context, AVHWAccel.free_frame_priv,9 Add AVHWAccel.update_thread_context, AVHWAccel.free_frame_priv,
7 AVHWAccel.flush.10 AVHWAccel.flush.
libavutil/mathematics.c+104−0
@@ -213,3 +213,107 @@ int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t i · 212 unmodified lines
213 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);213 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }214 }
215}215}
216
217static inline double eval_poly(const double *coeff, int size, double x) {
218 double sum = coeff[size-1];
219 int i;
220 for (i = size-2; i >= 0; --i) {
221 sum *= x;
222 sum += coeff[i];
223 }
224 return sum;
225}
226
227/**
228 * 0th order modified bessel function of the first kind.
229 * Algorithm taken from the Boost project, source:
230 * https://searchcode.com/codesearch/view/14918379/
231 * Use, modification and distribution are subject to the
232 * Boost Software License, Version 1.0 (see notice below).
233 * Boost Software License - Version 1.0 - August 17th, 2003
234Permission is hereby granted, free of charge, to any person or organization
235obtaining a copy of the software and accompanying documentation covered by
236this license (the "Software") to use, reproduce, display, distribute,
237execute, and transmit the Software, and to prepare derivative works of the
238Software, and to permit third-parties to whom the Software is furnished to
239do so, all subject to the following:
240
241The copyright notices in the Software and this entire statement, including
242the above license grant, this restriction and the following disclaimer,
243must be included in all copies of the Software, in whole or in part, and
244all derivative works of the Software, unless such copies or derivative
245works are solely in the form of machine-executable object code generated by
246a source language processor.
247
248THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
249IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
250FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
251SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
252FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
253ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
254DEALINGS IN THE SOFTWARE.
255 */
256
257double av_bessel_i0(double x) {
258// Modified Bessel function of the first kind of order zero
259// minimax rational approximations on intervals, see
260// Blair and Edwards, Chalk River Report AECL-4928, 1974
261 static const double p1[] = {
262 -2.2335582639474375249e+15,
263 -5.5050369673018427753e+14,
264 -3.2940087627407749166e+13,
265 -8.4925101247114157499e+11,
266 -1.1912746104985237192e+10,
267 -1.0313066708737980747e+08,
268 -5.9545626019847898221e+05,
269 -2.4125195876041896775e+03,
270 -7.0935347449210549190e+00,
271 -1.5453977791786851041e-02,
272 -2.5172644670688975051e-05,
273 -3.0517226450451067446e-08,
274 -2.6843448573468483278e-11,
275 -1.5982226675653184646e-14,
276 -5.2487866627945699800e-18,
277 };
278 static const double q1[] = {
279 -2.2335582639474375245e+15,
280 7.8858692566751002988e+12,
281 -1.2207067397808979846e+10,
282 1.0377081058062166144e+07,
283 -4.8527560179962773045e+03,
284 1.0,
285 };
286 static const double p2[] = {
287 -2.2210262233306573296e-04,
288 1.3067392038106924055e-02,
289 -4.4700805721174453923e-01,
290 5.5674518371240761397e+00,
291 -2.3517945679239481621e+01,
292 3.1611322818701131207e+01,
293 -9.6090021968656180000e+00,
294 };
295 static const double q2[] = {
296 -5.5194330231005480228e-04,
297 3.2547697594819615062e-02,
298 -1.1151759188741312645e+00,
299 1.3982595353892851542e+01,
300 -6.0228002066743340583e+01,
301 8.5539563258012929600e+01,
302 -3.1446690275135491500e+01,
303 1.0,
304 };
305 double y, r, factor;
306 if (x == 0)
307 return 1.0;
308 x = fabs(x);
309 if (x <= 15) {
310 y = x * x;
311 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 y = 1 / x - 1.0 / 15;
315 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 factor = exp(x) / sqrt(x);
317 return factor * r;
318 }
319}
libavutil/mathematics.h+4−0
@@ -288,6 +288,10 @@ int64_t av_rescale_delta(AVRational in_tb, int64_t in_ts, AVRational fs_tb, int · 287 unmodified lines
288 */288 */
289int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t inc);289int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t inc);
290 290
291/**
292 * 0th order modified bessel function of the first kind.
293 */
294double av_bessel_i0(double x);
291 295
292/**296/**
293 * @}297 * @}
libavutil/version.h+1−1
@@ -79,7 +79,7 @@ · 78 unmodified lines
79 */79 */
80 80
81#define LIBAVUTIL_VERSION_MAJOR 5881#define LIBAVUTIL_VERSION_MAJOR 58
82#define LIBAVUTIL_VERSION_MINOR 1182#define LIBAVUTIL_VERSION_MINOR 12
83#define LIBAVUTIL_VERSION_MICRO 10083#define LIBAVUTIL_VERSION_MICRO 100
84 84
85#define LIBAVUTIL_VERSION_INT AV_VERSION_INT(LIBAVUTIL_VERSION_MAJOR, \85#define LIBAVUTIL_VERSION_INT AV_VERSION_INT(LIBAVUTIL_VERSION_MAJOR, \
libswresample/resample.c+1−105
@@ -30,110 +30,6 @@ · 29 unmodified lines
30#include "libavutil/cpu.h"30#include "libavutil/cpu.h"
31#include "resample.h"31#include "resample.h"
32 32
33static inline double eval_poly(const double *coeff, int size, double x) {
34 double sum = coeff[size-1];
35 int i;
36 for (i = size-2; i >= 0; --i) {
37 sum *= x;
38 sum += coeff[i];
39 }
40 return sum;
41}
42
43/**
44 * 0th order modified bessel function of the first kind.
45 * Algorithm taken from the Boost project, source:
46 * https://searchcode.com/codesearch/view/14918379/
47 * Use, modification and distribution are subject to the
48 * Boost Software License, Version 1.0 (see notice below).
49 * Boost Software License - Version 1.0 - August 17th, 2003
50Permission is hereby granted, free of charge, to any person or organization
51obtaining a copy of the software and accompanying documentation covered by
52this license (the "Software") to use, reproduce, display, distribute,
53execute, and transmit the Software, and to prepare derivative works of the
54Software, and to permit third-parties to whom the Software is furnished to
55do so, all subject to the following:
56
57The copyright notices in the Software and this entire statement, including
58the above license grant, this restriction and the following disclaimer,
59must be included in all copies of the Software, in whole or in part, and
60all derivative works of the Software, unless such copies or derivative
61works are solely in the form of machine-executable object code generated by
62a source language processor.
63
64THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
65IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
66FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
67SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
68FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
69ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
70DEALINGS IN THE SOFTWARE.
71 */
72
73static double bessel(double x) {
74// Modified Bessel function of the first kind of order zero
75// minimax rational approximations on intervals, see
76// Blair and Edwards, Chalk River Report AECL-4928, 1974
77 static const double p1[] = {
78 -2.2335582639474375249e+15,
79 -5.5050369673018427753e+14,
80 -3.2940087627407749166e+13,
81 -8.4925101247114157499e+11,
82 -1.1912746104985237192e+10,
83 -1.0313066708737980747e+08,
84 -5.9545626019847898221e+05,
85 -2.4125195876041896775e+03,
86 -7.0935347449210549190e+00,
87 -1.5453977791786851041e-02,
88 -2.5172644670688975051e-05,
89 -3.0517226450451067446e-08,
90 -2.6843448573468483278e-11,
91 -1.5982226675653184646e-14,
92 -5.2487866627945699800e-18,
93 };
94 static const double q1[] = {
95 -2.2335582639474375245e+15,
96 7.8858692566751002988e+12,
97 -1.2207067397808979846e+10,
98 1.0377081058062166144e+07,
99 -4.8527560179962773045e+03,
100 1.0,
101 };
102 static const double p2[] = {
103 -2.2210262233306573296e-04,
104 1.3067392038106924055e-02,
105 -4.4700805721174453923e-01,
106 5.5674518371240761397e+00,
107 -2.3517945679239481621e+01,
108 3.1611322818701131207e+01,
109 -9.6090021968656180000e+00,
110 };
111 static const double q2[] = {
112 -5.5194330231005480228e-04,
113 3.2547697594819615062e-02,
114 -1.1151759188741312645e+00,
115 1.3982595353892851542e+01,
116 -6.0228002066743340583e+01,
117 8.5539563258012929600e+01,
118 -3.1446690275135491500e+01,
119 1.0,
120 };
121 double y, r, factor;
122 if (x == 0)
123 return 1.0;
124 x = fabs(x);
125 if (x <= 15) {
126 y = x * x;
127 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
128 }
129 else {
130 y = 1 / x - 1.0 / 15;
131 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
132 factor = exp(x) / sqrt(x);
133 return factor * r;
134 }
135}
136
137/**33/**
138 * builds a polyphase filterbank.34 * builds a polyphase filterbank.
139 * @param factor resampling factor35 * @param factor resampling factor
@@ -189,7 +85,7 @@ static int build_filter(ResampleContext *c, void *filter, double factor, int tap · 49 unmodified lines
189 break;85 break;
190 case SWR_FILTER_TYPE_KAISER:86 case SWR_FILTER_TYPE_KAISER:
191 w = 2.0*x / (factor*tap_count*M_PI);87 w = 2.0*x / (factor*tap_count*M_PI);
192 y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));88 y *= av_bessel_i0(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
193 break;89 break;
194 default:90 default:
195 av_assert0(0);91 av_assert0(0);

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t