3 * Copyright (c) 2004-2012 Michael Niedermayer <michaelni@gmx.at>
4 * bessel function: Copyright (c) 2006 Xiaogang Zhang
6 * This file is part of FFmpeg.
8 * FFmpeg is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
13 * FFmpeg is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with FFmpeg; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 * @author Michael Niedermayer <michaelni@gmx.at>
29 #include "libavutil/avassert.h"
32 static inline double eval_poly(const double *coeff, int size, double x) {
33 double sum = coeff[size-1];
35 for (i = size-2; i >= 0; --i) {
43 * 0th order modified bessel function of the first kind.
44 * Algorithm taken from the Boost project, source:
45 * https://searchcode.com/codesearch/view/14918379/
46 * Use, modification and distribution are subject to the
47 * Boost Software License, Version 1.0 (see notice below).
48 * Boost Software License - Version 1.0 - August 17th, 2003
49 Permission is hereby granted, free of charge, to any person or organization
50 obtaining a copy of the software and accompanying documentation covered by
51 this license (the "Software") to use, reproduce, display, distribute,
52 execute, and transmit the Software, and to prepare derivative works of the
53 Software, and to permit third-parties to whom the Software is furnished to
54 do so, all subject to the following:
56 The copyright notices in the Software and this entire statement, including
57 the above license grant, this restriction and the following disclaimer,
58 must be included in all copies of the Software, in whole or in part, and
59 all derivative works of the Software, unless such copies or derivative
60 works are solely in the form of machine-executable object code generated by
61 a source language processor.
63 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
64 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
65 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
66 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
67 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
68 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
69 DEALINGS IN THE SOFTWARE.
72 static double bessel(double x) {
73 // Modified Bessel function of the first kind of order zero
74 // minimax rational approximations on intervals, see
75 // Blair and Edwards, Chalk River Report AECL-4928, 1974
76 static const double p1[] = {
77 -2.2335582639474375249e+15,
78 -5.5050369673018427753e+14,
79 -3.2940087627407749166e+13,
80 -8.4925101247114157499e+11,
81 -1.1912746104985237192e+10,
82 -1.0313066708737980747e+08,
83 -5.9545626019847898221e+05,
84 -2.4125195876041896775e+03,
85 -7.0935347449210549190e+00,
86 -1.5453977791786851041e-02,
87 -2.5172644670688975051e-05,
88 -3.0517226450451067446e-08,
89 -2.6843448573468483278e-11,
90 -1.5982226675653184646e-14,
91 -5.2487866627945699800e-18,
93 static const double q1[] = {
94 -2.2335582639474375245e+15,
95 7.8858692566751002988e+12,
96 -1.2207067397808979846e+10,
97 1.0377081058062166144e+07,
98 -4.8527560179962773045e+03,
101 static const double p2[] = {
102 -2.2210262233306573296e-04,
103 1.3067392038106924055e-02,
104 -4.4700805721174453923e-01,
105 5.5674518371240761397e+00,
106 -2.3517945679239481621e+01,
107 3.1611322818701131207e+01,
108 -9.6090021968656180000e+00,
110 static const double q2[] = {
111 -5.5194330231005480228e-04,
112 3.2547697594819615062e-02,
113 -1.1151759188741312645e+00,
114 1.3982595353892851542e+01,
115 -6.0228002066743340583e+01,
116 8.5539563258012929600e+01,
117 -3.1446690275135491500e+01,
126 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
129 y = 1 / x - 1.0 / 15;
130 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
131 factor = exp(x) / sqrt(x);
137 * builds a polyphase filterbank.
138 * @param factor resampling factor
139 * @param scale wanted sum of coefficients for each filter
140 * @param filter_type filter type
141 * @param kaiser_beta kaiser window beta
142 * @return 0 on success, negative on error
144 static int build_filter(ResampleContext *c, void *filter, double factor, int tap_count, int alloc, int phase_count, int scale,
145 int filter_type, double kaiser_beta){
147 int ph_nb = phase_count % 2 ? phase_count : phase_count / 2 + 1;
148 double x, y, w, t, s;
149 double *tab = av_malloc_array(tap_count+1, sizeof(*tab));
150 double *sin_lut = av_malloc_array(ph_nb, sizeof(*sin_lut));
151 const int center= (tap_count-1)/2;
152 int ret = AVERROR(ENOMEM);
154 if (!tab || !sin_lut)
157 /* if upsampling, only need to interpolate, no filter */
162 for (ph = 0; ph < ph_nb; ph++)
163 sin_lut[ph] = sin(M_PI * ph / phase_count);
165 for(ph = 0; ph < ph_nb; ph++) {
168 for(i=0;i<=tap_count;i++) {
169 x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
171 else if (factor == 1.0)
176 case SWR_FILTER_TYPE_CUBIC:{
177 const float d= -0.5; //first order derivative = -0.5
178 x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
179 if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*( -x*x + x*x*x);
180 else y= d*(-4 + 8*x - 5*x*x + x*x*x);
182 case SWR_FILTER_TYPE_BLACKMAN_NUTTALL:
183 w = 2.0*x / (factor*tap_count);
185 y *= 0.3635819 - 0.4891775 * t + 0.1365995 * (2*t*t-1) - 0.0106411 * (4*t*t*t - 3*t);
187 case SWR_FILTER_TYPE_KAISER:
188 w = 2.0*x / (factor*tap_count*M_PI);
189 y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
201 /* normalize so that an uniform color remains the same */
203 case AV_SAMPLE_FMT_S16P:
204 for(i=0;i<tap_count;i++)
205 ((int16_t*)filter)[ph * alloc + i] = av_clip_int16(lrintf(tab[i] * scale / norm));
206 if (phase_count % 2) break;
207 if (tap_count % 2 == 0 || tap_count == 1) {
208 for (i = 0; i < tap_count; i++)
209 ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int16_t*)filter)[ph * alloc + i];
212 for (i = 1; i <= tap_count; i++)
213 ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
214 av_clip_int16(lrintf(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
217 case AV_SAMPLE_FMT_S32P:
218 for(i=0;i<tap_count;i++)
219 ((int32_t*)filter)[ph * alloc + i] = av_clipl_int32(llrint(tab[i] * scale / norm));
220 if (phase_count % 2) break;
221 if (tap_count % 2 == 0 || tap_count == 1) {
222 for (i = 0; i < tap_count; i++)
223 ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int32_t*)filter)[ph * alloc + i];
226 for (i = 1; i <= tap_count; i++)
227 ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
228 av_clipl_int32(llrint(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
231 case AV_SAMPLE_FMT_FLTP:
232 for(i=0;i<tap_count;i++)
233 ((float*)filter)[ph * alloc + i] = tab[i] * scale / norm;
234 if (phase_count % 2) break;
235 if (tap_count % 2 == 0 || tap_count == 1) {
236 for (i = 0; i < tap_count; i++)
237 ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((float*)filter)[ph * alloc + i];
240 for (i = 1; i <= tap_count; i++)
241 ((float*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
244 case AV_SAMPLE_FMT_DBLP:
245 for(i=0;i<tap_count;i++)
246 ((double*)filter)[ph * alloc + i] = tab[i] * scale / norm;
247 if (phase_count % 2) break;
248 if (tap_count % 2 == 0 || tap_count == 1) {
249 for (i = 0; i < tap_count; i++)
250 ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((double*)filter)[ph * alloc + i];
253 for (i = 1; i <= tap_count; i++)
254 ((double*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
263 double sine[LEN + tap_count];
264 double filtered[LEN];
265 double maxff=-2, minff=2, maxsf=-2, minsf=2;
266 for(i=0; i<LEN; i++){
267 double ss=0, sf=0, ff=0;
268 for(j=0; j<LEN+tap_count; j++)
269 sine[j]= cos(i*j*M_PI/LEN);
270 for(j=0; j<LEN; j++){
273 for(k=0; k<tap_count; k++)
274 sum += filter[ph * tap_count + k] * sine[k+j];
275 filtered[j]= sum / (1<<FILTER_SHIFT);
276 ss+= sine[j + center] * sine[j + center];
277 ff+= filtered[j] * filtered[j];
278 sf+= sine[j + center] * filtered[j];
283 maxff= FFMAX(maxff, ff);
284 minff= FFMIN(minff, ff);
285 maxsf= FFMAX(maxsf, sf);
286 minsf= FFMIN(minsf, sf);
288 av_log(NULL, AV_LOG_ERROR, "i:%4d ss:%f ff:%13.6e-%13.6e sf:%13.6e-%13.6e\n", i, ss, maxff, minff, maxsf, minsf);
303 static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
304 double cutoff0, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta,
305 double precision, int cheby, int exact_rational)
307 double cutoff = cutoff0? cutoff0 : 0.97;
308 double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
309 int phase_count= 1<<phase_shift;
310 int phase_count_compensation = phase_count;
312 if (exact_rational) {
313 int phase_count_exact, phase_count_exact_den;
315 av_reduce(&phase_count_exact, &phase_count_exact_den, out_rate, in_rate, INT_MAX);
316 if (phase_count_exact <= phase_count) {
317 phase_count_compensation = phase_count_exact * (phase_count / phase_count_exact);
318 phase_count = phase_count_exact;
322 if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
323 || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
324 || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
325 c = av_mallocz(sizeof(*c));
331 c->felem_size= av_get_bytes_per_sample(c->format);
334 case AV_SAMPLE_FMT_S16P:
335 c->filter_shift = 15;
337 case AV_SAMPLE_FMT_S32P:
338 c->filter_shift = 30;
340 case AV_SAMPLE_FMT_FLTP:
341 case AV_SAMPLE_FMT_DBLP:
345 av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
349 if (filter_size/factor > INT32_MAX/256) {
350 av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
354 c->phase_count = phase_count;
357 c->filter_length = FFMAX((int)ceil(filter_size/factor), 1);
358 c->filter_alloc = FFALIGN(c->filter_length, 8);
359 c->filter_bank = av_calloc(c->filter_alloc, (phase_count+1)*c->felem_size);
360 c->filter_type = filter_type;
361 c->kaiser_beta = kaiser_beta;
362 c->phase_count_compensation = phase_count_compensation;
365 if (build_filter(c, (void*)c->filter_bank, factor, c->filter_length, c->filter_alloc, phase_count, 1<<c->filter_shift, filter_type, kaiser_beta))
367 memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, c->filter_bank, (c->filter_alloc-1)*c->felem_size);
368 memcpy(c->filter_bank + (c->filter_alloc*phase_count )*c->felem_size, c->filter_bank + (c->filter_alloc - 1)*c->felem_size, c->felem_size);
371 c->compensation_distance= 0;
372 if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
374 while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
378 c->ideal_dst_incr = c->dst_incr;
379 c->dst_incr_div = c->dst_incr / c->src_incr;
380 c->dst_incr_mod = c->dst_incr % c->src_incr;
382 c->index= -phase_count*((c->filter_length-1)/2);
385 swri_resample_dsp_init(c);
389 av_freep(&c->filter_bank);
394 static void resample_free(ResampleContext **c){
397 av_freep(&(*c)->filter_bank);
401 static int rebuild_filter_bank_with_compensation(ResampleContext *c)
403 uint8_t *new_filter_bank;
404 int new_src_incr, new_dst_incr;
405 int phase_count = c->phase_count_compensation;
408 if (phase_count == c->phase_count)
411 av_assert0(!c->frac && !c->dst_incr_mod && !c->compensation_distance);
413 new_filter_bank = av_calloc(c->filter_alloc, (phase_count + 1) * c->felem_size);
414 if (!new_filter_bank)
415 return AVERROR(ENOMEM);
417 ret = build_filter(c, new_filter_bank, c->factor, c->filter_length, c->filter_alloc,
418 phase_count, 1 << c->filter_shift, c->filter_type, c->kaiser_beta);
420 av_freep(&new_filter_bank);
423 memcpy(new_filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, new_filter_bank, (c->filter_alloc-1)*c->felem_size);
424 memcpy(new_filter_bank + (c->filter_alloc*phase_count )*c->felem_size, new_filter_bank + (c->filter_alloc - 1)*c->felem_size, c->felem_size);
426 if (!av_reduce(&new_src_incr, &new_dst_incr, c->src_incr,
427 c->dst_incr * (int64_t)(phase_count/c->phase_count), INT32_MAX/2))
429 av_freep(&new_filter_bank);
430 return AVERROR(EINVAL);
433 c->src_incr = new_src_incr;
434 c->dst_incr = new_dst_incr;
435 while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
439 c->ideal_dst_incr = c->dst_incr;
440 c->dst_incr_div = c->dst_incr / c->src_incr;
441 c->dst_incr_mod = c->dst_incr % c->src_incr;
442 c->index *= phase_count / c->phase_count;
443 c->phase_count = phase_count;
444 av_freep(&c->filter_bank);
445 c->filter_bank = new_filter_bank;
449 static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
452 if (compensation_distance && sample_delta) {
453 ret = rebuild_filter_bank_with_compensation(c);
458 c->compensation_distance= compensation_distance;
459 if (compensation_distance)
460 c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
462 c->dst_incr = c->ideal_dst_incr;
464 c->dst_incr_div = c->dst_incr / c->src_incr;
465 c->dst_incr_mod = c->dst_incr % c->src_incr;
470 static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
472 int av_unused mm_flags = av_get_cpu_flags();
473 int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
474 (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
475 int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;
477 if (c->compensation_distance)
478 dst_size = FFMIN(dst_size, c->compensation_distance);
479 src_size = FFMIN(src_size, max_src_size);
483 if (c->filter_length == 1 && c->phase_count == 1) {
484 int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*c->index;
485 int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
486 int new_size = (src_size * (int64_t)c->src_incr - c->frac + c->dst_incr - 1) / c->dst_incr;
488 dst_size = FFMAX(FFMIN(dst_size, new_size), 0);
490 for (i = 0; i < dst->ch_count; i++) {
491 c->dsp.resample_one(dst->ch[i], src->ch[i], dst_size, index2, incr);
492 if (i+1 == dst->ch_count) {
493 c->index += dst_size * c->dst_incr_div;
494 c->index += (c->frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
495 av_assert2(c->index >= 0);
496 *consumed = c->index;
497 c->frac = (c->frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
503 int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
504 int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
505 int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
506 int (*resample_func)(struct ResampleContext *c, void *dst,
507 const void *src, int n, int update_ctx);
509 dst_size = FFMAX(FFMIN(dst_size, delta_n), 0);
511 /* resample_linear and resample_common should have same behavior
512 * when frac and dst_incr_mod are zero */
513 resample_func = (c->linear && (c->frac || c->dst_incr_mod)) ?
514 c->dsp.resample_linear : c->dsp.resample_common;
515 for (i = 0; i < dst->ch_count; i++)
516 *consumed = resample_func(c, dst->ch[i], src->ch[i], dst_size, i+1 == dst->ch_count);
523 if (c->compensation_distance) {
524 c->compensation_distance -= dst_size;
525 if (!c->compensation_distance) {
526 c->dst_incr = c->ideal_dst_incr;
527 c->dst_incr_div = c->dst_incr / c->src_incr;
528 c->dst_incr_mod = c->dst_incr % c->src_incr;
535 static int64_t get_delay(struct SwrContext *s, int64_t base){
536 ResampleContext *c = s->resample;
537 int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
538 num *= c->phase_count;
542 return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr * c->phase_count);
545 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
546 ResampleContext *c = s->resample;
547 // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
548 // They also make it easier to proof that changes and optimizations do not
549 // break the upper bound.
550 int64_t num = s->in_buffer_count + 2LL + in_samples;
551 num *= c->phase_count;
553 num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) * c->phase_count, AV_ROUND_UP) + 2;
555 if (c->compensation_distance) {
557 return AVERROR(EINVAL);
559 num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
564 static int resample_flush(struct SwrContext *s) {
565 AudioData *a= &s->in_buffer;
567 if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
569 av_assert0(a->planar);
570 for(i=0; i<a->ch_count; i++){
571 for(j=0; j<s->in_buffer_count; j++){
572 memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j )*a->bps,
573 a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
576 s->in_buffer_count += (s->in_buffer_count+1)/2;
580 // in fact the whole handle multiple ridiculously small buffers might need more thinking...
581 static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
582 int in_count, int *out_idx, int *out_sz)
584 int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;
589 if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
593 for (n = *out_sz; n < num; n++) {
594 for (ch = 0; ch < src->ch_count; ch++) {
595 memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
596 src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
600 // if not enough data is in, return and wait for more
601 if (num < c->filter_length + 1) {
603 *out_idx = c->filter_length;
608 for (n = 1; n <= c->filter_length; n++) {
609 for (ch = 0; ch < src->ch_count; ch++) {
610 memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
611 dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
617 *out_idx = c->filter_length;
618 while (c->index < 0) {
620 c->index += c->phase_count;
622 *out_sz = FFMAX(*out_sz + c->filter_length,
623 1 + c->filter_length * 2) - *out_idx;
625 return FFMAX(res, 0);
628 struct Resampler const swri_resampler={
635 invert_initial_buffer,