[FFmpeg-devel] [PATCH] avfilter: add sinc source filter
Paul B Mahol
onemda at gmail.com
Sun Oct 14 00:14:59 EEST 2018
Signed-off-by: Paul B Mahol <onemda at gmail.com>
---
libavfilter/Makefile | 1 +
libavfilter/allfilters.c | 1 +
libavfilter/asrc_sinc.c | 457 +++++++++++++++++++++++++++++++++++++++
3 files changed, 459 insertions(+)
create mode 100644 libavfilter/asrc_sinc.c
diff --git a/libavfilter/Makefile b/libavfilter/Makefile
index 62cc2f561f..b03b2457eb 100644
--- a/libavfilter/Makefile
+++ b/libavfilter/Makefile
@@ -141,6 +141,7 @@ OBJS-$(CONFIG_ANOISESRC_FILTER) += asrc_anoisesrc.o
OBJS-$(CONFIG_ANULLSRC_FILTER) += asrc_anullsrc.o
OBJS-$(CONFIG_FLITE_FILTER) += asrc_flite.o
OBJS-$(CONFIG_HILBERT_FILTER) += asrc_hilbert.o
+OBJS-$(CONFIG_SINC_FILTER) += asrc_sinc.o
OBJS-$(CONFIG_SINE_FILTER) += asrc_sine.o
OBJS-$(CONFIG_ANULLSINK_FILTER) += asink_anullsink.o
diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
index 5e72803b13..725bac94a0 100644
--- a/libavfilter/allfilters.c
+++ b/libavfilter/allfilters.c
@@ -134,6 +134,7 @@ extern AVFilter ff_asrc_anoisesrc;
extern AVFilter ff_asrc_anullsrc;
extern AVFilter ff_asrc_flite;
extern AVFilter ff_asrc_hilbert;
+extern AVFilter ff_asrc_sinc;
extern AVFilter ff_asrc_sine;
extern AVFilter ff_asink_anullsink;
diff --git a/libavfilter/asrc_sinc.c b/libavfilter/asrc_sinc.c
new file mode 100644
index 0000000000..0b2714f04c
--- /dev/null
+++ b/libavfilter/asrc_sinc.c
@@ -0,0 +1,457 @@
+/*
+ * Copyright (c) 2008-2009 Rob Sykes <robs at users.sourceforge.net>
+ * Copyright (c) 2017 Paul B Mahol
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "libavutil/avassert.h"
+#include "libavutil/opt.h"
+
+#include "libavcodec/avfft.h"
+
+#include "audio.h"
+#include "avfilter.h"
+#include "internal.h"
+
+typedef struct SincContext {
+ const AVClass *class;
+
+ int sample_rate, nb_samples;
+ float att, beta, phase, Fc0, Fc1, tbw0, tbw1;
+ int num_taps[2];
+ int round;
+
+ int n, rdft_len;
+ float *coeffs;
+ int64_t pts;
+
+ RDFTContext *rdft, *irdft;
+} SincContext;
+
+static int request_frame(AVFilterLink *outlink)
+{
+ AVFilterContext *ctx = outlink->src;
+ SincContext *s = ctx->priv;
+ const float *coeffs = s->coeffs;
+ AVFrame *frame = NULL;
+ int nb_samples;
+
+ nb_samples = FFMIN(s->nb_samples, s->n - s->pts);
+ if (!nb_samples)
+ return AVERROR_EOF;
+
+ if (!(frame = ff_get_audio_buffer(outlink, nb_samples)))
+ return AVERROR(ENOMEM);
+
+ memcpy(frame->data[0], coeffs + s->pts, nb_samples * sizeof(float));
+
+ frame->pts = s->pts;
+ s->pts += nb_samples;
+
+ return ff_filter_frame(outlink, frame);
+}
+
+static int query_formats(AVFilterContext *ctx)
+{
+ SincContext *s = ctx->priv;
+ static const int64_t chlayouts[] = { AV_CH_LAYOUT_MONO, -1 };
+ int sample_rates[] = { s->sample_rate, -1 };
+ static const enum AVSampleFormat sample_fmts[] = { AV_SAMPLE_FMT_FLT,
+ AV_SAMPLE_FMT_NONE };
+ AVFilterFormats *formats;
+ AVFilterChannelLayouts *layouts;
+ int ret;
+
+ formats = ff_make_format_list(sample_fmts);
+ if (!formats)
+ return AVERROR(ENOMEM);
+ ret = ff_set_common_formats (ctx, formats);
+ if (ret < 0)
+ return ret;
+
+ layouts = avfilter_make_format64_list(chlayouts);
+ if (!layouts)
+ return AVERROR(ENOMEM);
+ ret = ff_set_common_channel_layouts(ctx, layouts);
+ if (ret < 0)
+ return ret;
+
+ formats = ff_make_format_list(sample_rates);
+ if (!formats)
+ return AVERROR(ENOMEM);
+ return ff_set_common_samplerates(ctx, formats);
+}
+
+static float bessel_I_0(float x)
+{
+ float term = 1, sum = 1, last_sum, x2 = x / 2;
+ int i = 1;
+
+ do {
+ float y = x2 / i++;
+
+ last_sum = sum;
+ sum += term *= y * y;
+ } while (sum != last_sum);
+
+ return sum;
+}
+
+static float *make_lpf(int num_taps, float Fc, float beta, float rho,
+ float scale, int dc_norm)
+{
+ int i, m = num_taps - 1;
+ float *h = av_calloc(num_taps, sizeof(*h)), sum = 0;
+ float mult = scale / bessel_I_0(beta), mult1 = 1 / (.5 * m + rho);
+
+ av_assert0(Fc >= 0 && Fc <= 1);
+
+ for (i = 0; i <= m / 2; i++) {
+ float z = i - .5 * m, x = z * M_PI, y = z * mult1;
+ h[i] = x? sin(Fc * x) / x : Fc;
+ sum += h[i] *= bessel_I_0(beta * sqrt(1 - y * y)) * mult;
+ if (m - i != i) {
+ h[m - i] = h[i];
+ sum += h[i];
+ }
+ }
+
+ for (i = 0; dc_norm && i < num_taps; i++)
+ h[i] *= scale / sum;
+
+ return h;
+}
+
+static float kaiser_beta(float att, float tr_bw)
+{
+ if (att >= 60) {
+ static const float coefs[][4] = {
+ {-6.784957e-10, 1.02856e-05, 0.1087556, -0.8988365 + .001},
+ {-6.897885e-10, 1.027433e-05, 0.10876, -0.8994658 + .002},
+ {-1.000683e-09, 1.030092e-05, 0.1087677, -0.9007898 + .003},
+ {-3.654474e-10, 1.040631e-05, 0.1087085, -0.8977766 + .006},
+ {8.106988e-09, 6.983091e-06, 0.1091387, -0.9172048 + .015},
+ {9.519571e-09, 7.272678e-06, 0.1090068, -0.9140768 + .025},
+ {-5.626821e-09, 1.342186e-05, 0.1083999, -0.9065452 + .05},
+ {-9.965946e-08, 5.073548e-05, 0.1040967, -0.7672778 + .085},
+ {1.604808e-07, -5.856462e-05, 0.1185998, -1.34824 + .1},
+ {-1.511964e-07, 6.363034e-05, 0.1064627, -0.9876665 + .18},
+ };
+ float realm = log(tr_bw / .0005) / log(2.);
+ float const *c0 = coefs[av_clip((int) realm, 0, (int) FF_ARRAY_ELEMS(coefs) - 1)];
+ float const *c1 = coefs[av_clip(1 + (int) realm, 0, (int) FF_ARRAY_ELEMS(coefs) - 1)];
+ float b0 = ((c0[0] * att + c0[1]) * att + c0[2]) * att + c0[3];
+ float b1 = ((c1[0] * att + c1[1]) * att + c1[2]) * att + c1[3];
+
+ return b0 + (b1 - b0) * (realm - (int) realm);
+ }
+ if (att > 50)
+ return .1102 * (att - 8.7);
+ if (att > 20.96)
+ return .58417 * powf(att - 20.96, .4) + .07886 * (att - 20.96);
+ return 0;
+}
+
+static void kaiser_params(float att, float Fc, float tr_bw, float *beta, int *num_taps)
+{
+ *beta = *beta < 0 ? kaiser_beta(att, tr_bw * .5 / Fc): *beta;
+ att = att < 60 ? (att - 7.95) / (2.285 * M_PI * 2) :
+ ((.0007528358-1.577737e-05**beta)**beta+.6248022)**beta+.06186902;
+ *num_taps = !*num_taps ? ceil(att/tr_bw + 1) : *num_taps;
+}
+
+static float *lpf(float Fn, float Fc, float tbw, int *num_taps, float att, float *beta, int round)
+{
+ int n = *num_taps;
+
+ if ((Fc /= Fn) <= 0 || Fc >= 1) {
+ *num_taps = 0;
+ return NULL;
+ }
+
+ att = att ? att : 120;
+
+ kaiser_params(att, Fc, (tbw ? tbw / Fn : .05f) * .5f, beta, num_taps);
+
+ if (!n) {
+ n = *num_taps;
+ *num_taps = av_clip(n, 11, 32767);
+ if (round)
+ *num_taps = 1 + 2 * (int)((int)((*num_taps / 2) * Fc + .5f) / Fc + .5f);
+ }
+
+ return make_lpf(*num_taps |= 1, Fc, *beta, 0.f, 1.f, 0);
+}
+
+static void invert(float *h, int n)
+{
+ int i;
+
+ for (i = 0; i < n; i++)
+ h[i] = -h[i];
+
+ h[(n - 1) / 2] += 1;
+}
+
+#define PACK(h, n) h[1] = h[n]
+#define UNPACK(h, n) h[n] = h[1], h[n + 1] = h[1] = 0;
+#define SQR(a) ((a) * (a))
+
+static float safe_log(float x)
+{
+ av_assert0(x >= 0);
+ if (x)
+ return log(x);
+ return -26;
+}
+
+static int fir_to_phase(SincContext *s, float **h, int *len, int *post_len, float phase)
+{
+ float *pi_wraps, *work, phase1 = (phase > 50 ? 100 - phase : phase) / 50;
+ int i, work_len, begin, end, imp_peak = 0, peak = 0;
+ float imp_sum = 0, peak_imp_sum = 0;
+ float prev_angle2 = 0, cum_2pi = 0, prev_angle1 = 0, cum_1pi = 0;
+
+ for (i = *len, work_len = 2 * 2 * 8; i > 1; work_len <<= 1, i >>= 1);
+
+ work = av_calloc(work_len + 2, sizeof(*work)); /* +2: (UN)PACK */
+ pi_wraps = av_calloc(((work_len + 2) / 2), sizeof(*pi_wraps));
+ if (!work || !pi_wraps)
+ return AVERROR(ENOMEM);
+
+ memcpy(work, *h, *len * sizeof(*work));
+
+ av_rdft_end(s->rdft);
+ av_rdft_end(s->irdft);
+ s->rdft = s->irdft = NULL;
+ s->rdft = av_rdft_init(av_log2(work_len), DFT_R2C);
+ s->irdft = av_rdft_init(av_log2(work_len), IDFT_C2R);
+ if (!s->rdft || !s->irdft)
+ return AVERROR(ENOMEM);
+
+ av_rdft_calc(s->rdft, work); /* Cepstral: */
+ UNPACK(work, work_len);
+
+ for (i = 0; i <= work_len; i += 2) {
+ float angle = atan2f(work[i + 1], work[i]);
+ float detect = 2 * M_PI;
+ float delta = angle - prev_angle2;
+ float adjust = detect * ((delta < -detect * .7) - (delta > detect * .7));
+
+ prev_angle2 = angle;
+ cum_2pi += adjust;
+ angle += cum_2pi;
+ detect = M_PI;
+ delta = angle - prev_angle1;
+ adjust = detect * ((delta < -detect * .7) - (delta > detect * .7));
+ prev_angle1 = angle;
+ cum_1pi += fabsf(adjust); /* fabs for when 2pi and 1pi have combined */
+ pi_wraps[i >> 1] = cum_1pi;
+
+ work[i] = safe_log(sqrtf(SQR(work[i]) + SQR(work[i + 1])));
+ work[i + 1] = 0;
+ }
+
+ PACK(work, work_len);
+ av_rdft_calc(s->irdft, work);
+
+ for (i = 0; i < work_len; i++)
+ work[i] *= 2.f / work_len;
+
+ for (i = 1; i < work_len / 2; i++) { /* Window to reject acausal components */
+ work[i] *= 2;
+ work[i + work_len / 2] = 0;
+ }
+ av_rdft_calc(s->rdft, work);
+
+ for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
+ work[i + 1] = phase1 * i / work_len * pi_wraps[work_len >> 1] + (1 - phase1) * (work[i + 1] + pi_wraps[i >> 1]) - pi_wraps[i >> 1];
+
+ work[0] = exp(work[0]);
+ work[1] = exp(work[1]);
+ for (i = 2; i < work_len; i += 2) {
+ float x = expf(work[i]);
+
+ work[i ] = x * cosf(work[i + 1]);
+ work[i + 1] = x * sinf(work[i + 1]);
+ }
+
+ av_rdft_calc(s->irdft, work);
+ for (i = 0; i < work_len; i++)
+ work[i] *= 2.f / work_len;
+
+ /* Find peak pos. */
+ for (i = 0; i <= (int) (pi_wraps[work_len >> 1] / M_PI + .5f); i++) {
+ imp_sum += work[i];
+ if (fabs(imp_sum) > fabs(peak_imp_sum)) {
+ peak_imp_sum = imp_sum;
+ peak = i;
+ }
+ if (work[i] > work[imp_peak]) /* For debug check only */
+ imp_peak = i;
+ }
+
+ while (peak && fabsf(work[peak - 1]) > fabsf(work[peak]) && (work[peak - 1] * work[peak] > 0)) {
+ peak--;
+ }
+
+ if (!phase1) {
+ begin = 0;
+ } else if (phase1 == 1) {
+ begin = peak - *len / 2;
+ } else {
+ begin = (.997 - (2 - phase1) * .22) * *len + .5;
+ end = (.997 + (0 - phase1) * .22) * *len + .5;
+ begin = peak - (begin & ~3);
+ end = peak + 1 + ((end + 3) & ~3);
+ *len = end - begin;
+ *h = av_realloc(*h, *len * sizeof(**h));
+ if (!*h) {
+ av_free(pi_wraps);
+ av_free(work);
+ return AVERROR(ENOMEM);
+ }
+ }
+ for (i = 0; i < *len; i++) {
+ (*h)[i] = work[(begin + (phase > 50 ? *len - 1 - i : i) + work_len) & (work_len - 1)];
+ }
+ *post_len = phase > 50 ? peak - begin : begin + *len - (peak + 1);
+
+ av_log(s, AV_LOG_DEBUG, "%d nPI=%g peak-sum@%i=%g (val@%i=%g); len=%i post=%i (%g%%)\n",
+ work_len, pi_wraps[work_len >> 1] / M_PI, peak, peak_imp_sum, imp_peak,
+ work[imp_peak], *len, *post_len, 100 - 100. * *post_len / (*len - 1));
+
+ av_free(pi_wraps);
+ av_free(work);
+
+ return 0;
+}
+
+static int config_output(AVFilterLink *outlink)
+{
+ AVFilterContext *ctx = outlink->src;
+ SincContext *s = ctx->priv;
+ float Fn = s->sample_rate * .5;
+ float *h[2];
+ int i, n, post_peak, longer;
+
+ outlink->sample_rate = s->sample_rate;
+ s->pts = 0;
+
+ if (s->Fc0 >= Fn || s->Fc1 >= Fn) {
+ av_log(ctx, AV_LOG_ERROR,
+ "filter frequency must be less than %d/2.\n", s->sample_rate);
+ return AVERROR(EINVAL);
+ }
+
+ h[0] = lpf(Fn, s->Fc0, s->tbw0, &s->num_taps[0], s->att, &s->beta, s->round);
+ h[1] = lpf(Fn, s->Fc1, s->tbw1, &s->num_taps[1], s->att, &s->beta, s->round);
+
+ if (h[0])
+ invert(h[0], s->num_taps[0]);
+
+ longer = s->num_taps[1] > s->num_taps[0];
+ n = s->num_taps[longer];
+
+ if (h[0] && h[1]) {
+ for (i = 0; i < s->num_taps[!longer]; i++)
+ h[longer][i + (n - s->num_taps[!longer]) / 2] += h[!longer][i];
+
+ if (s->Fc0 < s->Fc1)
+ invert(h[longer], n);
+
+ av_free(h[!longer]);
+ }
+
+ if (s->phase != 50) {
+ int ret = fir_to_phase(s, &h[longer], &n, &post_peak, s->phase);
+
+ if (ret < 0)
+ return ret;
+ } else {
+ post_peak = n >> 1;
+ }
+
+ s->n = 1 << (av_log2(n) + 1);
+ s->rdft_len = 1 << av_log2(n);
+ s->coeffs = av_calloc(s->n, sizeof(*s->coeffs));
+ if (!s->coeffs)
+ return AVERROR(ENOMEM);
+
+ for (i = 0; i < n; i++)
+ s->coeffs[i] = h[longer][i];
+ av_free(h[longer]);
+
+ av_rdft_end(s->rdft);
+ av_rdft_end(s->irdft);
+ s->rdft = s->irdft = NULL;
+
+ return 0;
+}
+
+static av_cold void uninit(AVFilterContext *ctx)
+{
+ SincContext *s = ctx->priv;
+
+ av_freep(&s->coeffs);
+ av_rdft_end(s->rdft);
+ av_rdft_end(s->irdft);
+ s->rdft = s->irdft = NULL;
+}
+
+static const AVFilterPad sinc_outputs[] = {
+ {
+ .name = "default",
+ .type = AVMEDIA_TYPE_AUDIO,
+ .config_props = config_output,
+ .request_frame = request_frame,
+ },
+ { NULL }
+};
+
+#define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
+#define OFFSET(x) offsetof(SincContext, x)
+
+static const AVOption sinc_options[] = {
+ { "sample_rate", "set sample rate", OFFSET(sample_rate), AV_OPT_TYPE_INT, {.i64=44100}, 1, INT_MAX, AF },
+ { "r", "set sample rate", OFFSET(sample_rate), AV_OPT_TYPE_INT, {.i64=44100}, 1, INT_MAX, AF },
+ { "nb_samples", "set the number of samples per requested frame", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=1024}, 1, INT_MAX, AF },
+ { "n", "set the number of samples per requested frame", OFFSET(nb_samples), AV_OPT_TYPE_INT, {.i64=1024}, 1, INT_MAX, AF },
+ { "hp", "set high-pass filter frequency", OFFSET(Fc0), AV_OPT_TYPE_FLOAT, {.dbl=0}, 0, INT_MAX, AF },
+ { "lp", "set low-pass filter frequency", OFFSET(Fc1), AV_OPT_TYPE_FLOAT, {.dbl=0}, 0, INT_MAX, AF },
+ { "phase", "set filter phase response", OFFSET(phase), AV_OPT_TYPE_FLOAT, {.dbl=50}, 0, 100, AF },
+ { "beta", "set kaiser window beta", OFFSET(beta), AV_OPT_TYPE_FLOAT, {.dbl=-1}, -1, 256, AF },
+ { "att", "set stop-band attenuation", OFFSET(att), AV_OPT_TYPE_FLOAT, {.dbl=120}, 40, 180, AF },
+ { "round", "enable rounding", OFFSET(round), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, AF },
+ { "hptaps", "set number of taps for high-pass filter", OFFSET(num_taps[0]), AV_OPT_TYPE_INT, {.i64=0}, 0, 32768, AF },
+ { "lptaps", "set number of taps for low-pass filter", OFFSET(num_taps[1]), AV_OPT_TYPE_INT, {.i64=0}, 0, 32768, AF },
+ { NULL }
+};
+
+AVFILTER_DEFINE_CLASS(sinc);
+
+AVFilter ff_asrc_sinc = {
+ .name = "sinc",
+ .description = NULL_IF_CONFIG_SMALL("Generate a sinc kaiser-windowed low-pass, high-pass, band-pass, or band-reject FIR coefficients."),
+ .priv_size = sizeof(SincContext),
+ .priv_class = &sinc_class,
+ .query_formats = query_formats,
+ .uninit = uninit,
+ .inputs = NULL,
+ .outputs = sinc_outputs,
+};
--
2.17.1
More information about the ffmpeg-devel
mailing list