libavfilter/asrc_sinc.c - manifest_repos/ffmpeg - Git at Google

 /*
  * Copyright (c) 2008-2009 Rob Sykes <robs@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 <= 0)
         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 = ff_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.f / (.5f * m + rho);

     av_assert0(Fc >= 0 && Fc <= 1);

     for (i = 0; i <= m / 2; i++) {
         float z = i - .5f * m, x = z * M_PI, y = z * mult1;
         h[i] = x ? sinf(Fc * x) / x : Fc;
         sum += h[i] *= bessel_I_0(beta * sqrtf(1.f - 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.f) {
         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 = logf(tr_bw / .0005f) / logf(2.f);
         float const *c0 = coefs[av_clip((int)realm, 0, FF_ARRAY_ELEMS(coefs) - 1)];
         float const *c1 = coefs[av_clip(1 + (int)realm, 0, 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.f)
         return .1102f * (att - 8.7f);
     if (att > 20.96f)
         return .58417f * powf(att - 20.96f, .4f) + .07886f * (att - 20.96f);
     return 0;
 }

 static void kaiser_params(float att, float Fc, float tr_bw, float *beta, int *num_taps)
 {
     *beta = *beta < 0.f ? kaiser_beta(att, tr_bw * .5f / Fc): *beta;
     att = att < 60.f ? (att - 7.95f) / (2.285f * M_PI * 2.f) :
         ((.0007528358f-1.577737e-05 * *beta) * *beta + 0.6248022f) * *beta + .06186902f;
     *num_taps = !*num_taps ? ceilf(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.f || Fc >= 1.f) {
         *num_taps = 0;
         return NULL;
     }

     att = att ? att : 120.f;

     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)
 {
     for (int 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 logf(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.f ? 100.f - phase : phase) / 50.f;
     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);

     /* The first part is for work (+2 for (UN)PACK), the latter for pi_wraps. */
     work = av_calloc((work_len + 2) + (work_len / 2 + 1), sizeof(float));
     if (!work)
         return AVERROR(ENOMEM);
     pi_wraps = &work[work_len + 2];

     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) {
         av_free(work);
         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 * .7f) - (delta > detect * .7f));

         prev_angle2 = angle;
         cum_2pi += adjust;
         angle += cum_2pi;
         detect = M_PI;
         delta = angle - prev_angle1;
         adjust = detect * ((delta < -detect * .7f) - (delta > detect * .7f));
         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 = (.997f - (2 - phase1) * .22f) * *len + .5f;
         end = (.997f + (0 - phase1) * .22f) * *len + .5f;
         begin = peak - (begin & ~3);
         end = peak + 1 + ((end + 3) & ~3);
         *len = end - begin;
         *h = av_realloc_f(*h, *len, sizeof(**h));
         if (!*h) {
             av_free(work);
             return AVERROR(ENOMEM);
         }
     }

     for (i = 0; i < *len; i++) {
         (*h)[i] = work[(begin + (phase > 50.f ? *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.f - 100.f * *post_len / (*len - 1));

     av_free(work);

     return 0;
 }

 static int config_output(AVFilterLink *outlink)
 {
     AVFilterContext *ctx = outlink->src;
     SincContext *s = ctx->priv;
     float Fn = s->sample_rate * .5f;
     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.f) {
         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,
 };
	/*
	* Copyright (c) 2008-2009 Rob Sykes <robs@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 <= 0)
	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 = ff_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.f / (.5f * m + rho);

	av_assert0(Fc >= 0 && Fc <= 1);

	for (i = 0; i <= m / 2; i++) {
	float z = i - .5f * m, x = z * M_PI, y = z * mult1;
	h[i] = x ? sinf(Fc * x) / x : Fc;
	sum += h[i] = bessel_I_0(beta sqrtf(1.f - 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.f) {
	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 = logf(tr_bw / .0005f) / logf(2.f);
	float const *c0 = coefs[av_clip((int)realm, 0, FF_ARRAY_ELEMS(coefs) - 1)];
	float const *c1 = coefs[av_clip(1 + (int)realm, 0, 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.f)
	return .1102f * (att - 8.7f);
	if (att > 20.96f)
	return .58417f * powf(att - 20.96f, .4f) + .07886f * (att - 20.96f);
	return 0;
	}

	static void kaiser_params(float att, float Fc, float tr_bw, float beta, int num_taps)
	{
	beta = beta < 0.f ? kaiser_beta(att, tr_bw * .5f / Fc): *beta;
	att = att < 60.f ? (att - 7.95f) / (2.285f * M_PI * 2.f) :
	((.0007528358f-1.577737e-05 * beta) beta + 0.6248022f) *beta + .06186902f;
	num_taps = !num_taps ? ceilf(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.f \|\| Fc >= 1.f) {
	*num_taps = 0;
	return NULL;
	}

	att = att ? att : 120.f;

	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)
	{
	for (int 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 logf(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.f ? 100.f - phase : phase) / 50.f;
	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);

	/* The first part is for work (+2 for (UN)PACK), the latter for pi_wraps. */
	work = av_calloc((work_len + 2) + (work_len / 2 + 1), sizeof(float));
	if (!work)
	return AVERROR(ENOMEM);
	pi_wraps = &work[work_len + 2];

	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) {
	av_free(work);
	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 * .7f) - (delta > detect * .7f));

	prev_angle2 = angle;
	cum_2pi += adjust;
	angle += cum_2pi;
	detect = M_PI;
	delta = angle - prev_angle1;
	adjust = detect * ((delta < -detect * .7f) - (delta > detect * .7f));
	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 = (.997f - (2 - phase1) * .22f) * *len + .5f;
	end = (.997f + (0 - phase1) * .22f) * *len + .5f;
	begin = peak - (begin & ~3);
	end = peak + 1 + ((end + 3) & ~3);
	*len = end - begin;
	h = av_realloc_f(h, len, sizeof(*h));
	if (!*h) {
	av_free(work);
	return AVERROR(ENOMEM);
	}
	}

	for (i = 0; i < *len; i++) {
	(h)[i] = work[(begin + (phase > 50.f ? 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.f - 100.f * post_len / (len - 1));

	av_free(work);

	return 0;
	}

	static int config_output(AVFilterLink *outlink)
	{
	AVFilterContext *ctx = outlink->src;
	SincContext *s = ctx->priv;
	float Fn = s->sample_rate * .5f;
	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.f) {
	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,
	};