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(lrintf(tab[i] * scale / norm), INT16_MIN, INT16_MAX);
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(lrintf(tab[i] * scale / (norm - tab[0] + tab[tap_count])), INT16_MIN, INT16_MAX);
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 c->ideal_dst_incr = c->dst_incr;
359 c->dst_incr_div = c->dst_incr / c->src_incr;
360 c->dst_incr_mod = c->dst_incr % c->src_incr;
362 c->index= -phase_count*((c->filter_length-1)/2);
365 swri_resample_dsp_init(c);
369 av_freep(&c->filter_bank);
374 static void resample_free(ResampleContext **c){
377 av_freep(&(*c)->filter_bank);
381 static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
382 c->compensation_distance= compensation_distance;
383 if (compensation_distance)
384 c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
386 c->dst_incr = c->ideal_dst_incr;
388 c->dst_incr_div = c->dst_incr / c->src_incr;
389 c->dst_incr_mod = c->dst_incr % c->src_incr;
394 static int swri_resample(ResampleContext *c,
395 uint8_t *dst, const uint8_t *src, int *consumed,
396 int src_size, int dst_size, int update_ctx)
398 if (c->filter_length == 1 && c->phase_shift == 0) {
401 int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*index;
402 int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
403 int new_size = (src_size * (int64_t)c->src_incr - frac + c->dst_incr - 1) / c->dst_incr;
405 dst_size= FFMIN(dst_size, new_size);
406 c->dsp.resample_one(dst, src, dst_size, index2, incr);
408 index += dst_size * c->dst_incr_div;
409 index += (frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
410 av_assert2(index >= 0);
413 c->frac = (frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
417 int64_t end_index = (1LL + src_size - c->filter_length) << c->phase_shift;
418 int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
419 int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
421 dst_size = FFMIN(dst_size, delta_n);
423 *consumed = c->dsp.resample(c, dst, src, dst_size, update_ctx);
432 static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
434 int av_unused mm_flags = av_get_cpu_flags();
435 int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
436 (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
437 int64_t max_src_size = (INT64_MAX >> (c->phase_shift+1)) / c->src_incr;
439 if (c->compensation_distance)
440 dst_size = FFMIN(dst_size, c->compensation_distance);
441 src_size = FFMIN(src_size, max_src_size);
443 for(i=0; i<dst->ch_count; i++){
444 ret= swri_resample(c, dst->ch[i], src->ch[i],
445 consumed, src_size, dst_size, i+1==dst->ch_count);
450 if (c->compensation_distance) {
451 c->compensation_distance -= ret;
452 if (!c->compensation_distance) {
453 c->dst_incr = c->ideal_dst_incr;
454 c->dst_incr_div = c->dst_incr / c->src_incr;
455 c->dst_incr_mod = c->dst_incr % c->src_incr;
462 static int64_t get_delay(struct SwrContext *s, int64_t base){
463 ResampleContext *c = s->resample;
464 int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
465 num *= 1 << c->phase_shift;
469 return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr << c->phase_shift);
472 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
473 ResampleContext *c = s->resample;
474 // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
475 // They also make it easier to proof that changes and optimizations do not
476 // break the upper bound.
477 int64_t num = s->in_buffer_count + 2LL + in_samples;
478 num *= 1 << c->phase_shift;
480 num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) << c->phase_shift, AV_ROUND_UP) + 2;
482 if (c->compensation_distance) {
484 return AVERROR(EINVAL);
486 num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
491 static int resample_flush(struct SwrContext *s) {
492 AudioData *a= &s->in_buffer;
494 if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
496 av_assert0(a->planar);
497 for(i=0; i<a->ch_count; i++){
498 for(j=0; j<s->in_buffer_count; j++){
499 memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j )*a->bps,
500 a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
503 s->in_buffer_count += (s->in_buffer_count+1)/2;
507 // in fact the whole handle multiple ridiculously small buffers might need more thinking...
508 static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
509 int in_count, int *out_idx, int *out_sz)
511 int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;
516 if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
520 for (n = *out_sz; n < num; n++) {
521 for (ch = 0; ch < src->ch_count; ch++) {
522 memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
523 src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
527 // if not enough data is in, return and wait for more
528 if (num < c->filter_length + 1) {
530 *out_idx = c->filter_length;
535 for (n = 1; n <= c->filter_length; n++) {
536 for (ch = 0; ch < src->ch_count; ch++) {
537 memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
538 dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
544 *out_idx = c->filter_length + (c->index >> c->phase_shift);
545 *out_sz = FFMAX(*out_sz + c->filter_length,
546 1 + c->filter_length * 2) - *out_idx;
547 c->index &= c->phase_mask;
549 return FFMAX(res, 0);
552 struct Resampler const swri_resampler={
559 invert_initial_buffer,