swresample: add exact_rational option

give high quality resampling
as good as with linear_interp=on
as fast as without linear_interp=on
tested visually with ffplay
ffplay -f lavfi "aevalsrc='sin(10000*t*t)', aresample=osr=48000, showcqt=gamma=5"
ffplay -f lavfi "aevalsrc='sin(10000*t*t)', aresample=osr=48000:linear_interp=on, showcqt=gamma=5"
ffplay -f lavfi "aevalsrc='sin(10000*t*t)', aresample=osr=48000:exact_rational=on, showcqt=gamma=5"

slightly speed improvement
for fair comparison with -cpuflags 0
audio.wav is ~ 1 hour 44100 stereo 16bit wav file
ffmpeg -i audio.wav -af aresample=osr=48000 -f null -
        old         new
real    13.498s     13.121s
user    13.364s     12.987s
sys      0.131s      0.129s

linear_interp=on
        old         new
real    23.035s     23.050s
user    22.907s     22.917s
sys      0.119s     0.125s

exact_rational=on
real    12.418s
user    12.298s
sys      0.114s

possibility to decrease memory usage if soft compensation is ignored

Signed-off-by: Muhammad Faiz <mfcc64@gmail.com>

10 files changed, 78 insertions, 30 deletions

diff --git a/libswresample/arm/resample_init.c b/libswresample/arm/resample_init.c
index 797c53094a..003fafd29b 100644
--- a/libswresample/arm/resample_init.c
+++ b/libswresample/arm/resample_init.c
@@ -44,11 +44,15 @@ static int ff_resample_common_##TYPE##_neon(ResampleContext *c, void *dest, cons
     int dst_index;                                                                                \
     int index= c->index;                                                                          \
     int frac= c->frac;                                                                            \
-    int sample_index = index >> c->phase_shift;                                                   \
+    int sample_index = 0;                                                                         \
     int x4_aligned_filter_length = c->filter_length & ~3;                                         \
     int x8_aligned_filter_length = c->filter_length & ~7;                                         \
                                                                                                   \
-    index &= c->phase_mask;                                                                       \
+    while (index >= c->phase_count) {                                                             \
+        sample_index++;                                                                           \
+        index -= c->phase_count;                                                                  \
+    }                                                                                             \
+                                                                                                  \
     for (dst_index = 0; dst_index < n; dst_index++) {                                             \
         FELEM *filter = ((FELEM *) c->filter_bank) + c->filter_alloc * index;                     \
                                                                                                   \
@@ -75,8 +79,11 @@ static int ff_resample_common_##TYPE##_neon(ResampleContext *c, void *dest, cons
             frac -= c->src_incr;                                                                  \
             index++;                                                                              \
         }                                                                                         \
-        sample_index += index >> c->phase_shift;                                                  \
-        index &= c->phase_mask;                                                                   \
+                                                                                                  \
+        while (index >= c->phase_count) {                                                         \
+            sample_index++;                                                                       \
+            index -= c->phase_count;                                                              \
+        }                                                                                         \
     }                                                                                             \
                                                                                                   \
     if(update_ctx){                                                                               \
diff --git a/libswresample/options.c b/libswresample/options.c
index eac1a7c91f..816ce47750 100644
--- a/libswresample/options.c
+++ b/libswresample/options.c
@@ -85,6 +85,7 @@ static const AVOption options[]={
 {"filter_size"          , "set swr resampling filter size", OFFSET(filter_size)  , AV_OPT_TYPE_INT  , {.i64=32                    }, 0      , INT_MAX   , PARAM },
 {"phase_shift"          , "set swr resampling phase shift", OFFSET(phase_shift)  , AV_OPT_TYPE_INT  , {.i64=10                    }, 0      , 24        , PARAM },
 {"linear_interp"        , "enable linear interpolation" , OFFSET(linear_interp)  , AV_OPT_TYPE_BOOL , {.i64=0                     }, 0      , 1         , PARAM },
+{"exact_rational"       , "enable exact rational"       , OFFSET(exact_rational) , AV_OPT_TYPE_BOOL , {.i64=0                     }, 0      , 1         , PARAM },
 {"cutoff"               , "set cutoff frequency ratio"  , OFFSET(cutoff)         , AV_OPT_TYPE_DOUBLE,{.dbl=0.                    }, 0      , 1         , PARAM },
 
 /* duplicate option in order to work with avconv */
diff --git a/libswresample/resample.c b/libswresample/resample.c
index d410432658..1b1d83e596 100644
--- a/libswresample/resample.c
+++ b/libswresample/resample.c
@@ -297,13 +297,28 @@ fail:
 
 static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
                                     double cutoff0, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta,
-                                    double precision, int cheby)
+                                    double precision, int cheby, int exact_rational)
 {
     double cutoff = cutoff0? cutoff0 : 0.97;
     double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
     int phase_count= 1<<phase_shift;
 
-    if (!c || c->phase_shift != phase_shift || c->linear!=linear || c->factor != factor
+    if (exact_rational) {
+        int phase_count_exact, phase_count_exact_den;
+
+        av_reduce(&phase_count_exact, &phase_count_exact_den, out_rate, in_rate, INT_MAX);
+        /* FIXME this is not required, but build_filter needs even phase_count */
+        if (phase_count_exact & 1 && phase_count_exact > 1 && phase_count_exact < INT_MAX/2)
+            phase_count_exact *= 2;
+
+        if (phase_count_exact <= phase_count) {
+            /* FIXME this is not required when soft compensation is disabled */
+            phase_count_exact *= phase_count / phase_count_exact;
+            phase_count = phase_count_exact;
+        }
+    }
+
+    if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
            || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
            || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
         c = av_mallocz(sizeof(*c));
@@ -337,6 +352,7 @@ static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_r
 
         c->phase_shift   = phase_shift;
         c->phase_mask    = phase_count - 1;
+        c->phase_count   = phase_count;
         c->linear        = linear;
         c->factor        = factor;
         c->filter_length = FFMAX((int)ceil(filter_size/factor), 1);
@@ -399,7 +415,7 @@ static int swri_resample(ResampleContext *c,
                          uint8_t *dst, const uint8_t *src, int *consumed,
                          int src_size, int dst_size, int update_ctx)
 {
-    if (c->filter_length == 1 && c->phase_shift == 0) {
+    if (c->filter_length == 1 && c->phase_count == 1) {
         int index= c->index;
         int frac= c->frac;
         int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*index;
@@ -418,7 +434,7 @@ static int swri_resample(ResampleContext *c,
             c->index = 0;
         }
     } else {
-        int64_t end_index = (1LL + src_size - c->filter_length) << c->phase_shift;
+        int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
         int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
         int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
 
@@ -438,7 +454,7 @@ static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, A
     int av_unused mm_flags = av_get_cpu_flags();
     int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
                     (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
-    int64_t max_src_size = (INT64_MAX >> (c->phase_shift+1)) / c->src_incr;
+    int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;
 
     if (c->compensation_distance)
         dst_size = FFMIN(dst_size, c->compensation_distance);
@@ -466,11 +482,11 @@ static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, A
 static int64_t get_delay(struct SwrContext *s, int64_t base){
     ResampleContext *c = s->resample;
     int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
-    num *= 1 << c->phase_shift;
+    num *= c->phase_count;
     num -= c->index;
     num *= c->src_incr;
     num -= c->frac;
-    return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr << c->phase_shift);
+    return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr * c->phase_count);
 }
 
 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
@@ -479,9 +495,9 @@ static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
     // They also make it easier to proof that changes and optimizations do not
     // break the upper bound.
     int64_t num = s->in_buffer_count + 2LL + in_samples;
-    num *= 1 << c->phase_shift;
+    num *= c->phase_count;
     num -= c->index;
-    num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) << c->phase_shift, AV_ROUND_UP) + 2;
+    num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) * c->phase_count, AV_ROUND_UP) + 2;
 
     if (c->compensation_distance) {
         if (num > INT_MAX)
@@ -545,10 +561,13 @@ static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const Audio
     }
 
     res = num - *out_sz;
-    *out_idx = c->filter_length + (c->index >> c->phase_shift);
+    *out_idx = c->filter_length;
+    while (c->index < 0) {
+        --*out_idx;
+        c->index += c->phase_count;
+    }
     *out_sz = FFMAX(*out_sz + c->filter_length,
                     1 + c->filter_length * 2) - *out_idx;
-    c->index &= c->phase_mask;
 
     return FFMAX(res, 0);
 }
diff --git a/libswresample/resample.h b/libswresample/resample.h
index ca0cf2b114..53788c484f 100644
--- a/libswresample/resample.h
+++ b/libswresample/resample.h
@@ -40,8 +40,10 @@ typedef struct ResampleContext {
     int frac;
     int src_incr;
     int compensation_distance;
-    int phase_shift;
-    int phase_mask;
+    /* TODO remove phase_shift and phase_mask */
+    attribute_deprecated int phase_shift;
+    attribute_deprecated int phase_mask;
+    int phase_count;
     int linear;
     enum SwrFilterType filter_type;
     double kaiser_beta;
diff --git a/libswresample/resample_template.c b/libswresample/resample_template.c
index d71efd63c0..1636f4e95d 100644
--- a/libswresample/resample_template.c
+++ b/libswresample/resample_template.c
@@ -92,9 +92,13 @@ static int RENAME(resample_common)(ResampleContext *c,
     int dst_index;
     int index= c->index;
     int frac= c->frac;
-    int sample_index = index >> c->phase_shift;
+    int sample_index = 0;
+
+    while (index >= c->phase_count) {
+        sample_index++;
+        index -= c->phase_count;
+    }
 
-    index &= c->phase_mask;
     for (dst_index = 0; dst_index < n; dst_index++) {
         FELEM *filter = ((FELEM *) c->filter_bank) + c->filter_alloc * index;
 
@@ -111,8 +115,11 @@ static int RENAME(resample_common)(ResampleContext *c,
             frac -= c->src_incr;
             index++;
         }
-        sample_index += index >> c->phase_shift;
-        index &= c->phase_mask;
+
+        while (index >= c->phase_count) {
+            sample_index++;
+            index -= c->phase_count;
+        }
     }
 
     if(update_ctx){
@@ -132,12 +139,16 @@ static int RENAME(resample_linear)(ResampleContext *c,
     int dst_index;
     int index= c->index;
     int frac= c->frac;
-    int sample_index = index >> c->phase_shift;
+    int sample_index = 0;
 #if FILTER_SHIFT == 0
     double inv_src_incr = 1.0 / c->src_incr;
 #endif
 
-    index &= c->phase_mask;
+    while (index >= c->phase_count) {
+        sample_index++;
+        index -= c->phase_count;
+    }
+
     for (dst_index = 0; dst_index < n; dst_index++) {
         FELEM *filter = ((FELEM *) c->filter_bank) + c->filter_alloc * index;
         FELEM2 val=0, v2 = 0;
@@ -164,8 +175,11 @@ static int RENAME(resample_linear)(ResampleContext *c,
             frac -= c->src_incr;
             index++;
         }
-        sample_index += index >> c->phase_shift;
-        index &= c->phase_mask;
+
+        while (index >= c->phase_count) {
+            sample_index++;
+            index -= c->phase_count;
+        }
     }
 
     if(update_ctx){
diff --git a/libswresample/soxr_resample.c b/libswresample/soxr_resample.c
index 807f508405..b9c6735028 100644
--- a/libswresample/soxr_resample.c
+++ b/libswresample/soxr_resample.c
@@ -30,7 +30,7 @@
 #include <soxr.h>
 
 static struct ResampleContext *create(struct ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
-        double cutoff, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta, double precision, int cheby){
+        double cutoff, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta, double precision, int cheby, int exact_rational){
     soxr_error_t error;
 
     soxr_datatype_t type =
diff --git a/libswresample/swresample.c b/libswresample/swresample.c
index e02ec81db6..351623b9d6 100644
--- a/libswresample/swresample.c
+++ b/libswresample/swresample.c
@@ -262,7 +262,7 @@ av_cold int swr_init(struct SwrContext *s){
     }
 
     if (s->out_sample_rate!=s->in_sample_rate || (s->flags & SWR_FLAG_RESAMPLE)){
-        s->resample = s->resampler->init(s->resample, s->out_sample_rate, s->in_sample_rate, s->filter_size, s->phase_shift, s->linear_interp, s->cutoff, s->int_sample_fmt, s->filter_type, s->kaiser_beta, s->precision, s->cheby);
+        s->resample = s->resampler->init(s->resample, s->out_sample_rate, s->in_sample_rate, s->filter_size, s->phase_shift, s->linear_interp, s->cutoff, s->int_sample_fmt, s->filter_type, s->kaiser_beta, s->precision, s->cheby, s->exact_rational);
         if (!s->resample) {
             av_log(s, AV_LOG_ERROR, "Failed to initialize resampler\n");
             return AVERROR(ENOMEM);
diff --git a/libswresample/swresample_internal.h b/libswresample/swresample_internal.h
index 5d61de47af..3828b722cc 100644
--- a/libswresample/swresample_internal.h
+++ b/libswresample/swresample_internal.h
@@ -69,7 +69,7 @@ struct DitherContext {
 };
 
 typedef struct ResampleContext * (* resample_init_func)(struct ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
-                                    double cutoff, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta, double precision, int cheby);
+                                    double cutoff, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta, double precision, int cheby, int exact_rational);
 typedef void    (* resample_free_func)(struct ResampleContext **c);
 typedef int     (* multiple_resample_func)(struct ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed);
 typedef int     (* resample_flush_func)(struct SwrContext *c);
@@ -126,6 +126,7 @@ struct SwrContext {
     int filter_size;                                /**< length of each FIR filter in the resampling filterbank relative to the cutoff frequency */
     int phase_shift;                                /**< log2 of the number of entries in the resampling polyphase filterbank */
     int linear_interp;                              /**< if 1 then the resampling FIR filter will be linearly interpolated */
+    int exact_rational;                             /**< if 1 then enable non power of 2 phase_count */
     double cutoff;                                  /**< resampling cutoff frequency (swr: 6dB point; soxr: 0dB point). 1.0 corresponds to half the output sample rate */
     int filter_type;                                /**< swr resampling filter type */
     double kaiser_beta;                                /**< swr beta value for Kaiser window (only applicable if filter_type == AV_FILTER_TYPE_KAISER) */
diff --git a/libswresample/version.h b/libswresample/version.h
index 830f00e430..b8e32c0e41 100644
--- a/libswresample/version.h
+++ b/libswresample/version.h
@@ -29,8 +29,8 @@
 #include "libavutil/avutil.h"
 
 #define LIBSWRESAMPLE_VERSION_MAJOR   2
-#define LIBSWRESAMPLE_VERSION_MINOR   0
-#define LIBSWRESAMPLE_VERSION_MICRO 101
+#define LIBSWRESAMPLE_VERSION_MINOR   1
+#define LIBSWRESAMPLE_VERSION_MICRO 100
 
 #define LIBSWRESAMPLE_VERSION_INT  AV_VERSION_INT(LIBSWRESAMPLE_VERSION_MAJOR, \
                                                   LIBSWRESAMPLE_VERSION_MINOR, \
diff --git a/libswresample/x86/resample_init.c b/libswresample/x86/resample_init.c
index 9d7d5cf89e..f05397fecb 100644
--- a/libswresample/x86/resample_init.c
+++ b/libswresample/x86/resample_init.c
@@ -47,6 +47,10 @@ av_cold void swri_resample_dsp_x86_init(ResampleContext *c)
 {
     int av_unused mm_flags = av_get_cpu_flags();
 
+    /* FIXME use phase_count on asm */
+    if (c->phase_count != 1 << c->phase_shift)
+        return;
+
     switch(c->format){
     case AV_SAMPLE_FMT_S16P:
         if (ARCH_X86_32 && EXTERNAL_MMXEXT(mm_flags)) {