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 double x, y, w, t, s;
148 double *tab = av_malloc_array(tap_count+1, sizeof(*tab));
149 double *sin_lut = av_malloc_array(phase_count / 2 + 1, sizeof(*sin_lut));
150 const int center= (tap_count-1)/2;
152 if (!tab || !sin_lut)
155 /* if upsampling, only need to interpolate, no filter */
159 av_assert0(phase_count == 1 || phase_count % 2 == 0);
162 for (ph = 0; ph <= phase_count / 2; ph++)
163 sin_lut[ph] = sin(M_PI * ph / phase_count);
165 for(ph = 0; ph <= phase_count / 2; 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 (tap_count % 2 == 0) {
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 (tap_count % 2 == 0) {
220 for (i = 0; i < tap_count; i++)
221 ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int32_t*)filter)[ph * alloc + i];
224 for (i = 1; i <= tap_count; i++)
225 ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
226 av_clipl_int32(llrint(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
229 case AV_SAMPLE_FMT_FLTP:
230 for(i=0;i<tap_count;i++)
231 ((float*)filter)[ph * alloc + i] = tab[i] * scale / norm;
232 if (tap_count % 2 == 0) {
233 for (i = 0; i < tap_count; i++)
234 ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((float*)filter)[ph * alloc + i];
237 for (i = 1; i <= tap_count; i++)
238 ((float*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
241 case AV_SAMPLE_FMT_DBLP:
242 for(i=0;i<tap_count;i++)
243 ((double*)filter)[ph * alloc + i] = tab[i] * scale / norm;
244 if (tap_count % 2 == 0) {
245 for (i = 0; i < tap_count; i++)
246 ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((double*)filter)[ph * alloc + i];
249 for (i = 1; i <= tap_count; i++)
250 ((double*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
259 double sine[LEN + tap_count];
260 double filtered[LEN];
261 double maxff=-2, minff=2, maxsf=-2, minsf=2;
262 for(i=0; i<LEN; i++){
263 double ss=0, sf=0, ff=0;
264 for(j=0; j<LEN+tap_count; j++)
265 sine[j]= cos(i*j*M_PI/LEN);
266 for(j=0; j<LEN; j++){
269 for(k=0; k<tap_count; k++)
270 sum += filter[ph * tap_count + k] * sine[k+j];
271 filtered[j]= sum / (1<<FILTER_SHIFT);
272 ss+= sine[j + center] * sine[j + center];
273 ff+= filtered[j] * filtered[j];
274 sf+= sine[j + center] * filtered[j];
279 maxff= FFMAX(maxff, ff);
280 minff= FFMIN(minff, ff);
281 maxsf= FFMAX(maxsf, sf);
282 minsf= FFMIN(minsf, sf);
284 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);
298 static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
299 double cutoff0, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta,
300 double precision, int cheby)
302 double cutoff = cutoff0? cutoff0 : 0.97;
303 double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
304 int phase_count= 1<<phase_shift;
306 if (!c || c->phase_shift != phase_shift || c->linear!=linear || c->factor != factor
307 || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
308 || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
309 c = av_mallocz(sizeof(*c));
315 c->felem_size= av_get_bytes_per_sample(c->format);
318 case AV_SAMPLE_FMT_S16P:
319 c->filter_shift = 15;
321 case AV_SAMPLE_FMT_S32P:
322 c->filter_shift = 30;
324 case AV_SAMPLE_FMT_FLTP:
325 case AV_SAMPLE_FMT_DBLP:
329 av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
333 if (filter_size/factor > INT32_MAX/256) {
334 av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
338 c->phase_shift = phase_shift;
339 c->phase_mask = phase_count - 1;
342 c->filter_length = FFMAX((int)ceil(filter_size/factor), 1);
343 c->filter_alloc = FFALIGN(c->filter_length, 8);
344 c->filter_bank = av_calloc(c->filter_alloc, (phase_count+1)*c->felem_size);
345 c->filter_type = filter_type;
346 c->kaiser_beta = kaiser_beta;
349 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))
351 memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, c->filter_bank, (c->filter_alloc-1)*c->felem_size);
352 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);
355 c->compensation_distance= 0;
356 if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
358 while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
362 c->ideal_dst_incr = c->dst_incr;
363 c->dst_incr_div = c->dst_incr / c->src_incr;
364 c->dst_incr_mod = c->dst_incr % c->src_incr;
366 c->index= -phase_count*((c->filter_length-1)/2);
369 swri_resample_dsp_init(c);
373 av_freep(&c->filter_bank);
378 static void resample_free(ResampleContext **c){
381 av_freep(&(*c)->filter_bank);
385 static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
386 c->compensation_distance= compensation_distance;
387 if (compensation_distance)
388 c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
390 c->dst_incr = c->ideal_dst_incr;
392 c->dst_incr_div = c->dst_incr / c->src_incr;
393 c->dst_incr_mod = c->dst_incr % c->src_incr;
398 static int swri_resample(ResampleContext *c,
399 uint8_t *dst, const uint8_t *src, int *consumed,
400 int src_size, int dst_size, int update_ctx)
402 if (c->filter_length == 1 && c->phase_shift == 0) {
405 int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*index;
406 int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
407 int new_size = (src_size * (int64_t)c->src_incr - frac + c->dst_incr - 1) / c->dst_incr;
409 dst_size= FFMIN(dst_size, new_size);
410 c->dsp.resample_one(dst, src, dst_size, index2, incr);
412 index += dst_size * c->dst_incr_div;
413 index += (frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
414 av_assert2(index >= 0);
417 c->frac = (frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
421 int64_t end_index = (1LL + src_size - c->filter_length) << c->phase_shift;
422 int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
423 int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
425 dst_size = FFMIN(dst_size, delta_n);
427 *consumed = c->dsp.resample(c, dst, src, dst_size, update_ctx);
436 static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
438 int av_unused mm_flags = av_get_cpu_flags();
439 int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
440 (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
441 int64_t max_src_size = (INT64_MAX >> (c->phase_shift+1)) / c->src_incr;
443 if (c->compensation_distance)
444 dst_size = FFMIN(dst_size, c->compensation_distance);
445 src_size = FFMIN(src_size, max_src_size);
447 for(i=0; i<dst->ch_count; i++){
448 ret= swri_resample(c, dst->ch[i], src->ch[i],
449 consumed, src_size, dst_size, i+1==dst->ch_count);
454 if (c->compensation_distance) {
455 c->compensation_distance -= ret;
456 if (!c->compensation_distance) {
457 c->dst_incr = c->ideal_dst_incr;
458 c->dst_incr_div = c->dst_incr / c->src_incr;
459 c->dst_incr_mod = c->dst_incr % c->src_incr;
466 static int64_t get_delay(struct SwrContext *s, int64_t base){
467 ResampleContext *c = s->resample;
468 int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
469 num *= 1 << c->phase_shift;
473 return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr << c->phase_shift);
476 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
477 ResampleContext *c = s->resample;
478 // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
479 // They also make it easier to proof that changes and optimizations do not
480 // break the upper bound.
481 int64_t num = s->in_buffer_count + 2LL + in_samples;
482 num *= 1 << c->phase_shift;
484 num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) << c->phase_shift, AV_ROUND_UP) + 2;
486 if (c->compensation_distance) {
488 return AVERROR(EINVAL);
490 num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
495 static int resample_flush(struct SwrContext *s) {
496 AudioData *a= &s->in_buffer;
498 if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
500 av_assert0(a->planar);
501 for(i=0; i<a->ch_count; i++){
502 for(j=0; j<s->in_buffer_count; j++){
503 memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j )*a->bps,
504 a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
507 s->in_buffer_count += (s->in_buffer_count+1)/2;
511 // in fact the whole handle multiple ridiculously small buffers might need more thinking...
512 static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
513 int in_count, int *out_idx, int *out_sz)
515 int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;
520 if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
524 for (n = *out_sz; n < num; n++) {
525 for (ch = 0; ch < src->ch_count; ch++) {
526 memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
527 src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
531 // if not enough data is in, return and wait for more
532 if (num < c->filter_length + 1) {
534 *out_idx = c->filter_length;
539 for (n = 1; n <= c->filter_length; n++) {
540 for (ch = 0; ch < src->ch_count; ch++) {
541 memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
542 dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
548 *out_idx = c->filter_length + (c->index >> c->phase_shift);
549 *out_sz = FFMAX(*out_sz + c->filter_length,
550 1 + c->filter_length * 2) - *out_idx;
551 c->index &= c->phase_mask;
553 return FFMAX(res, 0);
556 struct Resampler const swri_resampler={
563 invert_initial_buffer,