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;
153 if (!tab || !sin_lut)
156 /* if upsampling, only need to interpolate, no filter */
161 for (ph = 0; ph < ph_nb; ph++)
162 sin_lut[ph] = sin(M_PI * ph / phase_count);
164 for(ph = 0; ph < ph_nb; ph++) {
167 for(i=0;i<=tap_count;i++) {
168 x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
170 else if (factor == 1.0)
175 case SWR_FILTER_TYPE_CUBIC:{
176 const float d= -0.5; //first order derivative = -0.5
177 x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
178 if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*( -x*x + x*x*x);
179 else y= d*(-4 + 8*x - 5*x*x + x*x*x);
181 case SWR_FILTER_TYPE_BLACKMAN_NUTTALL:
182 w = 2.0*x / (factor*tap_count);
184 y *= 0.3635819 - 0.4891775 * t + 0.1365995 * (2*t*t-1) - 0.0106411 * (4*t*t*t - 3*t);
186 case SWR_FILTER_TYPE_KAISER:
187 w = 2.0*x / (factor*tap_count*M_PI);
188 y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
200 /* normalize so that an uniform color remains the same */
202 case AV_SAMPLE_FMT_S16P:
203 for(i=0;i<tap_count;i++)
204 ((int16_t*)filter)[ph * alloc + i] = av_clip_int16(lrintf(tab[i] * scale / norm));
205 if (phase_count % 2) break;
206 if (tap_count % 2 == 0 || tap_count == 1) {
207 for (i = 0; i < tap_count; i++)
208 ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int16_t*)filter)[ph * alloc + i];
211 for (i = 1; i <= tap_count; i++)
212 ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
213 av_clip_int16(lrintf(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
216 case AV_SAMPLE_FMT_S32P:
217 for(i=0;i<tap_count;i++)
218 ((int32_t*)filter)[ph * alloc + i] = av_clipl_int32(llrint(tab[i] * scale / norm));
219 if (phase_count % 2) break;
220 if (tap_count % 2 == 0 || tap_count == 1) {
221 for (i = 0; i < tap_count; i++)
222 ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int32_t*)filter)[ph * alloc + i];
225 for (i = 1; i <= tap_count; i++)
226 ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
227 av_clipl_int32(llrint(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
230 case AV_SAMPLE_FMT_FLTP:
231 for(i=0;i<tap_count;i++)
232 ((float*)filter)[ph * alloc + i] = tab[i] * scale / norm;
233 if (phase_count % 2) break;
234 if (tap_count % 2 == 0 || tap_count == 1) {
235 for (i = 0; i < tap_count; i++)
236 ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((float*)filter)[ph * alloc + i];
239 for (i = 1; i <= tap_count; i++)
240 ((float*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
243 case AV_SAMPLE_FMT_DBLP:
244 for(i=0;i<tap_count;i++)
245 ((double*)filter)[ph * alloc + i] = tab[i] * scale / norm;
246 if (phase_count % 2) break;
247 if (tap_count % 2 == 0 || tap_count == 1) {
248 for (i = 0; i < tap_count; i++)
249 ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((double*)filter)[ph * alloc + i];
252 for (i = 1; i <= tap_count; i++)
253 ((double*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
262 double sine[LEN + tap_count];
263 double filtered[LEN];
264 double maxff=-2, minff=2, maxsf=-2, minsf=2;
265 for(i=0; i<LEN; i++){
266 double ss=0, sf=0, ff=0;
267 for(j=0; j<LEN+tap_count; j++)
268 sine[j]= cos(i*j*M_PI/LEN);
269 for(j=0; j<LEN; j++){
272 for(k=0; k<tap_count; k++)
273 sum += filter[ph * tap_count + k] * sine[k+j];
274 filtered[j]= sum / (1<<FILTER_SHIFT);
275 ss+= sine[j + center] * sine[j + center];
276 ff+= filtered[j] * filtered[j];
277 sf+= sine[j + center] * filtered[j];
282 maxff= FFMAX(maxff, ff);
283 minff= FFMIN(minff, ff);
284 maxsf= FFMAX(maxsf, sf);
285 minsf= FFMIN(minsf, sf);
287 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);
301 static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
302 double cutoff0, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta,
303 double precision, int cheby, int exact_rational)
305 double cutoff = cutoff0? cutoff0 : 0.97;
306 double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
307 int phase_count= 1<<phase_shift;
308 int phase_count_compensation = phase_count;
310 if (exact_rational) {
311 int phase_count_exact, phase_count_exact_den;
313 av_reduce(&phase_count_exact, &phase_count_exact_den, out_rate, in_rate, INT_MAX);
314 if (phase_count_exact <= phase_count) {
315 phase_count_compensation = phase_count_exact * (phase_count / phase_count_exact);
316 phase_count = phase_count_exact;
320 if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
321 || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
322 || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
323 c = av_mallocz(sizeof(*c));
329 c->felem_size= av_get_bytes_per_sample(c->format);
332 case AV_SAMPLE_FMT_S16P:
333 c->filter_shift = 15;
335 case AV_SAMPLE_FMT_S32P:
336 c->filter_shift = 30;
338 case AV_SAMPLE_FMT_FLTP:
339 case AV_SAMPLE_FMT_DBLP:
343 av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
347 if (filter_size/factor > INT32_MAX/256) {
348 av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
352 c->phase_count = phase_count;
355 c->filter_length = FFMAX((int)ceil(filter_size/factor), 1);
356 c->filter_alloc = FFALIGN(c->filter_length, 8);
357 c->filter_bank = av_calloc(c->filter_alloc, (phase_count+1)*c->felem_size);
358 c->filter_type = filter_type;
359 c->kaiser_beta = kaiser_beta;
360 c->phase_count_compensation = phase_count_compensation;
363 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))
365 memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, c->filter_bank, (c->filter_alloc-1)*c->felem_size);
366 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);
369 c->compensation_distance= 0;
370 if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
372 while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
376 c->ideal_dst_incr = c->dst_incr;
377 c->dst_incr_div = c->dst_incr / c->src_incr;
378 c->dst_incr_mod = c->dst_incr % c->src_incr;
380 c->index= -phase_count*((c->filter_length-1)/2);
383 swri_resample_dsp_init(c);
387 av_freep(&c->filter_bank);
392 static void resample_free(ResampleContext **c){
395 av_freep(&(*c)->filter_bank);
399 static int rebuild_filter_bank_with_compensation(ResampleContext *c)
401 uint8_t *new_filter_bank;
402 int new_src_incr, new_dst_incr;
403 int phase_count = c->phase_count_compensation;
406 if (phase_count == c->phase_count)
409 av_assert0(!c->frac && !c->dst_incr_mod && !c->compensation_distance);
411 new_filter_bank = av_calloc(c->filter_alloc, (phase_count + 1) * c->felem_size);
412 if (!new_filter_bank)
413 return AVERROR(ENOMEM);
415 ret = build_filter(c, new_filter_bank, c->factor, c->filter_length, c->filter_alloc,
416 phase_count, 1 << c->filter_shift, c->filter_type, c->kaiser_beta);
418 av_freep(&new_filter_bank);
421 memcpy(new_filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, new_filter_bank, (c->filter_alloc-1)*c->felem_size);
422 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);
424 if (!av_reduce(&new_src_incr, &new_dst_incr, c->src_incr,
425 c->dst_incr * (int64_t)(phase_count/c->phase_count), INT32_MAX/2))
427 av_freep(&new_filter_bank);
428 return AVERROR(EINVAL);
431 c->src_incr = new_src_incr;
432 c->dst_incr = new_dst_incr;
433 while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
437 c->ideal_dst_incr = c->dst_incr;
438 c->dst_incr_div = c->dst_incr / c->src_incr;
439 c->dst_incr_mod = c->dst_incr % c->src_incr;
440 c->index *= phase_count / c->phase_count;
441 c->phase_count = phase_count;
442 av_freep(&c->filter_bank);
443 c->filter_bank = new_filter_bank;
447 static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
450 if (compensation_distance) {
451 ret = rebuild_filter_bank_with_compensation(c);
456 c->compensation_distance= compensation_distance;
457 if (compensation_distance)
458 c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
460 c->dst_incr = c->ideal_dst_incr;
462 c->dst_incr_div = c->dst_incr / c->src_incr;
463 c->dst_incr_mod = c->dst_incr % c->src_incr;
468 static int swri_resample(ResampleContext *c,
469 uint8_t *dst, const uint8_t *src, int *consumed,
470 int src_size, int dst_size, int update_ctx)
472 if (c->filter_length == 1 && c->phase_count == 1) {
475 int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*index;
476 int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
477 int new_size = (src_size * (int64_t)c->src_incr - frac + c->dst_incr - 1) / c->dst_incr;
479 dst_size= FFMIN(dst_size, new_size);
480 c->dsp.resample_one(dst, src, dst_size, index2, incr);
482 index += dst_size * c->dst_incr_div;
483 index += (frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
484 av_assert2(index >= 0);
487 c->frac = (frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
491 int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
492 int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
493 int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
495 dst_size = FFMIN(dst_size, delta_n);
497 *consumed = c->dsp.resample(c, dst, src, dst_size, update_ctx);
506 static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
508 int av_unused mm_flags = av_get_cpu_flags();
509 int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
510 (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
511 int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;
513 if (c->compensation_distance)
514 dst_size = FFMIN(dst_size, c->compensation_distance);
515 src_size = FFMIN(src_size, max_src_size);
517 for(i=0; i<dst->ch_count; i++){
518 ret= swri_resample(c, dst->ch[i], src->ch[i],
519 consumed, src_size, dst_size, i+1==dst->ch_count);
524 if (c->compensation_distance) {
525 c->compensation_distance -= ret;
526 if (!c->compensation_distance) {
527 c->dst_incr = c->ideal_dst_incr;
528 c->dst_incr_div = c->dst_incr / c->src_incr;
529 c->dst_incr_mod = c->dst_incr % c->src_incr;
536 static int64_t get_delay(struct SwrContext *s, int64_t base){
537 ResampleContext *c = s->resample;
538 int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
539 num *= c->phase_count;
543 return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr * c->phase_count);
546 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
547 ResampleContext *c = s->resample;
548 // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
549 // They also make it easier to proof that changes and optimizations do not
550 // break the upper bound.
551 int64_t num = s->in_buffer_count + 2LL + in_samples;
552 num *= c->phase_count;
554 num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) * c->phase_count, AV_ROUND_UP) + 2;
556 if (c->compensation_distance) {
558 return AVERROR(EINVAL);
560 num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
565 static int resample_flush(struct SwrContext *s) {
566 AudioData *a= &s->in_buffer;
568 if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
570 av_assert0(a->planar);
571 for(i=0; i<a->ch_count; i++){
572 for(j=0; j<s->in_buffer_count; j++){
573 memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j )*a->bps,
574 a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
577 s->in_buffer_count += (s->in_buffer_count+1)/2;
581 // in fact the whole handle multiple ridiculously small buffers might need more thinking...
582 static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
583 int in_count, int *out_idx, int *out_sz)
585 int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;
590 if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
594 for (n = *out_sz; n < num; n++) {
595 for (ch = 0; ch < src->ch_count; ch++) {
596 memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
597 src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
601 // if not enough data is in, return and wait for more
602 if (num < c->filter_length + 1) {
604 *out_idx = c->filter_length;
609 for (n = 1; n <= c->filter_length; n++) {
610 for (ch = 0; ch < src->ch_count; ch++) {
611 memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
612 dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
618 *out_idx = c->filter_length;
619 while (c->index < 0) {
621 c->index += c->phase_count;
623 *out_sz = FFMAX(*out_sz + c->filter_length,
624 1 + c->filter_length * 2) - *out_idx;
626 return FFMAX(res, 0);
629 struct Resampler const swri_resampler={
636 invert_initial_buffer,