diff mbox

[FFmpeg-devel] lavfi: add mbfequalizer filter.

Message ID 20191223111720.87345-1-george@nsup.org
State New
Headers show

Commit Message

Nicolas George Dec. 23, 2019, 11:17 a.m. UTC
TODO changelog and version bump

Signed-off-by: Nicolas George <george@nsup.org>
---
 doc/filters.texi              |  48 +++
 libavfilter/Makefile          |   1 +
 libavfilter/af_mbfequalizer.c | 567 ++++++++++++++++++++++++++++++++++
 libavfilter/allfilters.c      |   1 +
 4 files changed, 617 insertions(+)
 create mode 100644 libavfilter/af_mbfequalizer.c


I think I'll rewrite the bands parser to allow per-channel options and make
the structure of the code clearer. But I can already send the whole patch
for your consideration.

Comments

James Almer Dec. 23, 2019, 1:15 p.m. UTC | #1
On 12/23/2019 8:17 AM, Nicolas George wrote:
> TODO changelog and version bump
> 
> Signed-off-by: Nicolas George <george@nsup.org>
> ---
>  doc/filters.texi              |  48 +++
>  libavfilter/Makefile          |   1 +
>  libavfilter/af_mbfequalizer.c | 567 ++++++++++++++++++++++++++++++++++
>  libavfilter/allfilters.c      |   1 +
>  4 files changed, 617 insertions(+)
>  create mode 100644 libavfilter/af_mbfequalizer.c
> 
> 
> I think I'll rewrite the bands parser to allow per-channel options and make
> the structure of the code clearer. But I can already send the whole patch
> for your consideration.

[...]

> +static int filter_frame(AVFilterLink *link, AVFrame *frame)
> +{
> +    AVFilterContext *ctx = link->dst;
> +    EqContext *s = ctx->priv;
> +    AVFrame *frame_out;
> +
> +    if (av_frame_is_writable(frame)) {
> +        frame_out = frame;
> +    } else {
> +        frame_out = ff_get_audio_buffer(ctx->outputs[0], frame->nb_samples);
> +        if (!frame_out) {
> +            av_frame_free(&frame);
> +            return AVERROR(ENOMEM);
> +        }
> +    }

Can't you use av_frame_make_writable() instead to simplify this?
Nicolas George Dec. 23, 2019, 1:25 p.m. UTC | #2
Thanks for the comment.

James Almer (12019-12-23):
> > +        frame_out = ff_get_audio_buffer(ctx->outputs[0], frame->nb_samples);

> Can't you use av_frame_make_writable() instead to simplify this?

No, because it does not use ff_get_audio_buffer(), and therefore does
not allow direct rendering. Also, it unnecessarily copies the data. I
have already noted that this is a frequent pattern in filters that
should be factored soon.

Regards,
James Almer Dec. 23, 2019, 1:38 p.m. UTC | #3
On 12/23/2019 10:25 AM, Nicolas George wrote:
> Thanks for the comment.
> 
> James Almer (12019-12-23):
>>> +        frame_out = ff_get_audio_buffer(ctx->outputs[0], frame->nb_samples);
> 
>> Can't you use av_frame_make_writable() instead to simplify this?
> 
> No, because it does not use ff_get_audio_buffer(), and therefore does
> not allow direct rendering. Also, it unnecessarily copies the data. I
> have already noted that this is a frequent pattern in filters that
> should be factored soon.
> 
> Regards,

I see, thanks for the clarification.

> 
> 
> _______________________________________________
> ffmpeg-devel mailing list
> ffmpeg-devel@ffmpeg.org
> https://ffmpeg.org/mailman/listinfo/ffmpeg-devel
> 
> To unsubscribe, visit link above, or email
> ffmpeg-devel-request@ffmpeg.org with subject "unsubscribe".
>
Paul B Mahol Jan. 31, 2020, 3:39 p.m. UTC | #4
On 12/23/19, Nicolas George <george@nsup.org> wrote:
> Thanks for the comment.
>
> James Almer (12019-12-23):
>> > +        frame_out = ff_get_audio_buffer(ctx->outputs[0],
>> > frame->nb_samples);
>
>> Can't you use av_frame_make_writable() instead to simplify this?
>
> No, because it does not use ff_get_audio_buffer(), and therefore does
> not allow direct rendering. Also, it unnecessarily copies the data. I
> have already noted that this is a frequent pattern in filters that
> should be factored soon.
>

I do not get it. Which filters are affected by this?

> Regards,
>
> --
>   Nicolas George
>
Paul B Mahol Jan. 31, 2020, 3:42 p.m. UTC | #5
On 12/23/19, Nicolas George <george@nsup.org> wrote:
> TODO changelog and version bump
>
> Signed-off-by: Nicolas George <george@nsup.org>
> ---
>  doc/filters.texi              |  48 +++
>  libavfilter/Makefile          |   1 +
>  libavfilter/af_mbfequalizer.c | 567 ++++++++++++++++++++++++++++++++++
>  libavfilter/allfilters.c      |   1 +
>  4 files changed, 617 insertions(+)
>  create mode 100644 libavfilter/af_mbfequalizer.c
>
>
> I think I'll rewrite the bands parser to allow per-channel options and make
> the structure of the code clearer. But I can already send the whole patch
> for your consideration.
>
>
> diff --git a/doc/filters.texi b/doc/filters.texi
> index 8c5d3a5760..c6917fad38 100644
> --- a/doc/filters.texi
> +++ b/doc/filters.texi
> @@ -4383,6 +4383,54 @@
> lv2=p=http\\\\://www.openavproductions.com/artyfx#bitta:c=crush=0.3
>  @end example
>  @end itemize
>
> +@section mbfequalizer
> +Multiband fast equalizer.
> +
> +As an equalizer, it allows to selectively alter the volume of the different
> +frequencies in the audio signal. This particular equalizer does not allow a
> +very sharp separation of the frequencies, but that gives it other
> +advantages. It is rather fast. It does not add latency. It supports all
> +sample formats, and for integer formats it uses bit-exact integer
> +arithmetic.
> +
> +It accepts the following parameters:
> +
> +@table @option
> +@item defs
> +Definition of bands and gain. The syntax is:
> +"gain/freq/gain[/freq/gain...][|gain/freq/gain...]"
> +where @var{freq} is the cut frequency between bands and @var{gain} the
> +associated gain.
> +Gain and bands for several channels can be defined by separating them with
> +the pipe "|" character; if only one definition is given, it is used for all
> +channels; otherwise, the number of definitions must match the number of
> +input channels.
> +Frequencies should be in ascending order. Gain greater than 1 or negative
> +may result in overflows and ugly noise or other strange results.
> +
> +@item global_gain
> +Gain applied to all channels and all bands. Can be greater than 1 without
> +risk; in case of overflow, soft clipping is applied.
> +
> +@item keep_all_bands
> +By default, bands beyond @var{sample_rate}/4 are discarded. If this option
> +is enabled, they are kept; it can result in strange effects.
> +@end table
> +
> +A Gnuplot formula mapping the frequency to the gain is printed at verbose
> +log level. It can be used to check the result, proably with
> +"set logscale x; db(x) = 20*log(x)/log(10); plot [20:2000] db(...)".
> +
> +@subsection Examples
> +
> +Assuming 3.1 input, reduce high frequencies on left and right, keep a
> narrow
> +band on center and dampen bass on LFE. All gains were lowered by 5dB, then
> +raised again using global gain, and the -2.8dB for the middle band was
> found
> +by looking at the graph.
> +@example
> +mbfequalizer=-5dB/7040/-15dB|-5dB/7040/-15dB|0/220/-2.8dB/880/0|-10dB/110/-5dB:global_gain=+5dB
> +@end example
> +
>  @section mcompand
>  Multiband Compress or expand the audio's dynamic range.
>
> diff --git a/libavfilter/Makefile b/libavfilter/Makefile
> index 37d4eee858..f5690a50db 100644
> --- a/libavfilter/Makefile
> +++ b/libavfilter/Makefile
> @@ -124,6 +124,7 @@ OBJS-$(CONFIG_LOWPASS_FILTER)                +=
> af_biquads.o
>  OBJS-$(CONFIG_LOWSHELF_FILTER)               += af_biquads.o
>  OBJS-$(CONFIG_LV2_FILTER)                    += af_lv2.o
>  OBJS-$(CONFIG_MCOMPAND_FILTER)               += af_mcompand.o
> +OBJS-$(CONFIG_MBFEQUALIZER_FILTER)           += af_mbfequalizer.o
>  OBJS-$(CONFIG_PAN_FILTER)                    += af_pan.o
>  OBJS-$(CONFIG_REPLAYGAIN_FILTER)             += af_replaygain.o
>  OBJS-$(CONFIG_RESAMPLE_FILTER)               += af_resample.o
> diff --git a/libavfilter/af_mbfequalizer.c b/libavfilter/af_mbfequalizer.c
> new file mode 100644
> index 0000000000..aadd3cb3da
> --- /dev/null
> +++ b/libavfilter/af_mbfequalizer.c
> @@ -0,0 +1,567 @@
> +/*
> + * Copyright (c) 2019 Nicolas George <george@nsup.org>
> + *
> + * 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
> + */
> +
> +/*
> + * Low-pass-high-pass equalizer using a simple sliding weighed average.
> + *
> + * Input signal: u(t)
> + * Low freqs:    a(t)
> + * High freqs:   b(t)
> + *
> + * a(t) = (1-k) a(t-1) + k u(t)
> + * b(t) = u(t) - a(t)
> + *
> + * F = sample frequency; f = input frequency; ω = 2πf/F
> + * For u(t) = exp(iωt),
> + * a(t) =                 k/(1-exp(-iω)(1-k)) exp(iwt)
> + * b(t) = (1-k)(1-exp(-iω))/(1-exp(-iω)(1-k)) exp(iwt)
> + * hi/low = (1-k)(1-exp(-iω)) / k
> + * R(k,w) = (1-k)*(1-exp(-I*w)) / k
> + *
> + * Cut when |hi/low| = 1, i.e. |1-exp(iω)| = k/(1-k)
> + * Triangle 0, 1, exp(iω) has sides 1, 1, k/(1-k)
> + * w = 2 asin(1/2 k/(1-k))
> + * k/(1-k) = 2 sin(ω/2)
> + * k = 2 sin(ω/2) / (1 + 2 sin(ω/2))
> + *
> + * Gnuplot script: response for cut frequency 110 Hz and 440 Hz.
> + * I = {0,1}
> + * F = 48000
> + * w(f) = 2*pi*f/F
> + * db(x) = 20*log(x)/log(10)
> + * set samples 1000
> + * set grid
> + * set logscale x 2
> + * k(f) = 2*sin(w(f)/2) / (1+2*sin(w(f)/2))
> + * d(f, x) = abs(exp(I*w(x))-(1-k(f)))
> + * al(f, x) = k(f) / d(f, x)
> + * ah(f, x) = abs((1-k(f))*(1-exp(I*w(x)))) / d(f, x)
> + * plot [20:16000] [-20:0] db(al(110, x)), db(ah(110, x)), db(al(880, x)),
> db(ah(880, x))
> + *
> + * Approximation without sin():
> + * ω - ω³/6  ≤  sin(ω)   ≤ ω - ω³/6  + ω⁵/120
> + * ω - ω³/24 ≤ 2sin(ω/2) ≤ ω - ω³/24 + ω⁵/1920
> + * Accurate enough with 0 ≤ ω ≤ π/2.
> + *
> + * Clip to [-A;A] with soft decay on [A(1-k); +∞[:
> + * s(x) = (Ak)² / (x - A(1-2k))
> + */
> +
> +/*
> +
> +aevalsrc=sin(13.75*4.35*exp(0.7*t)):d=10
> +arealtime,asplit=3[a][b][c];
> +[a]showwaves=s=1764x128[ap];
> +[b]mbfequalizer=1:0,showwaves=s=1764x128[bp];
> +[c]mbfequalizer=0:1,showwaves=s=1764x128[cp];
> +[ap][bp][cp]vstack=3
> +
> +*/
> +
> +#include "libavutil/avassert.h"
> +#include "libavutil/bprint.h"
> +#include "libavutil/eval.h"
> +#include "libavutil/opt.h"
> +#include "avfilter.h"
> +#include "audio.h"
> +#include "internal.h"
> +
> +typedef struct EqBand {
> +    int k; /* k(2πf) */
> +    int g; /* gain of the high band */
> +} EqBand;
> +
> +typedef struct EqContext EqContext;
> +struct EqContext {
> +    const AVClass *class;
> +    char *defs;
> +    EqBand *bands;
> +    void *vals;
> +    void (*process_samples)(EqContext *s, unsigned channels, unsigned
> samples,
> +                            AVFrame *frame, AVFrame *frame_out);
> +    double global_gain_opt;
> +    unsigned channels;
> +    int global_gain;
> +    int keep_all_bands;
> +};
> +
> +#define FPOINT 15
> +#define SOFTCLIP_SHIFT 2 /* k = 1/(1<<SHIFT) */
> +
> +#if 1&&HAVE_FAST_64BIT
> +
> +static inline int fpmul_int(int a, int b)
> +{
> +    return ((int64_t)a * b) >> FPOINT;
> +}
> +
> +#else
> +
> +static inline int fpmul_int(int a, int b)
> +{
> +#define MASK ((1 << FPOINT) - 1)
> +    return (((a >> FPOINT) * (b >> FPOINT)) << FPOINT)
> +         + (a >> FPOINT) * (b & MASK)
> +         + (b >> FPOINT) * (a & MASK)
> +         + (((a & MASK) * (b & MASK)) >> FPOINT);
> +#undef MASK
> +}
> +
> +#endif
> +
> +static int64_t fpmul_int64_t(int a, int64_t b)
> +{
> +    uint64_t a0 = a & 0xFFFFFFFF;
> +    uint64_t b0 = b & 0xFFFFFFFF;
> +    /*int64_t  a1 = a >> 32;*/ /* a1 = 0 */
> +    int64_t  b1 = b >> 32;
> +    uint64_t r0 = a0 * b0;
> +    int64_t r1a = (int64_t)a0 * b1;
> +    /*int64_t r1b = (int64_t)b0 * a1;*/
> +    int64_t r2  = 0 /*a1 * b1*/;
> +    r2 += (r1a >> 32) /*+ (r1b >> 32)*/;
> +    r1a &= 0xFFFFFFFF;
> +    /*r1b &= 0xFFFFFFFF;*/
> +    r2 += (r1a /*+ r1b*/ + (r0 >> 32)) >> 32;
> +    r0 += (r1a /*+ r1b*/) << 32;
> +    return (int64_t)(r0 >> FPOINT) + (r2 << (64 - FPOINT));
> +}
> +
> +static inline float fpmul_float(int a, float b)
> +{
> +    return a * b / (1 << FPOINT);
> +}
> +
> +static inline double fpmul_double(int a, double b)
> +{
> +    return a * b / (1 << FPOINT);
> +}
> +
> +static int approx_factor(int f, int sf)
> +{
> +    /* 31089 / 151 =~ 2*Pi*32768 / 1000 */
> +    int w = av_rescale(f, 31089, 151 * sf);
> +    int s = w - fpmul_int(w, fpmul_int(w, w)) / 24;
> +    int r = (s << FPOINT) / ((1 << FPOINT) + s);
> +    return r;
> +}
> +
> +static void reorder_bands(EqBand *bands, size_t nb_bands, unsigned
> sample_rate)
> +{
> +    size_t i, j;
> +
> +    if (nb_bands == 1)
> +        return;
> +    /* bands arrive in ascending order with kf = gain below,
> +       and k = 0 on the final band;
> +       we want them in descending order with kf = gain above;
> +       therefore, reverse the order of k except the final 0,
> +       and reverse the order of kf including the final 0 */
> +    for (i = 0, j = nb_bands - 2; i < j; i++, j--)
> +        FFSWAP(int, bands[i].k, bands[j].k);
> +    for (i = 0, j = nb_bands - 1; i < j; i++, j--)
> +        FFSWAP(int, bands[i].g, bands[j].g);
> +    for (i = 0; i < nb_bands - 1; i++)
> +        bands[i].k = approx_factor(bands[i].k, sample_rate);
> +}
> +
> +static void gain_to_bprint(AVBPrint *bp, unsigned rate, EqBand *bands)
> +{
> +    double k = bands->k / (double)(1 << FPOINT);
> +    double g = bands->g / (double)(1 << FPOINT);
> +
> +    if (bands->k == 0) {
> +        av_bprintf(bp, "%.5f", g);
> +        return;
> +    }
> +    av_bprintf(bp, "(%.5f*%.5f*(1-exp(-{0,1}*2*pi*x/%d))+%.5f*", g, 1 - k,
> rate, k);
> +    gain_to_bprint(bp, rate, bands + 1);
> +    av_bprintf(bp, ")");
> +    av_bprintf(bp, "/(1-exp(-{0,1}*2*pi*x/%d)*%.5f)", rate, 1 - k);
> +}
> +
> +static void print_gain(AVFilterContext *ctx, unsigned channel, EqBand
> *bands)
> +{
> +    EqContext *s = ctx->priv;
> +    AVBPrint bp;
> +    unsigned rate = ctx->outputs[0]->sample_rate;
> +
> +    av_bprint_init(&bp, 0, AV_BPRINT_SIZE_AUTOMATIC);
> +    if (s->global_gain != 1 << FPOINT)
> +        av_bprintf(&bp, "%.5f*", s->global_gain / (double)(1 << FPOINT));
> +    av_bprintf(&bp, "abs(");
> +    gain_to_bprint(&bp, rate, bands);
> +    av_bprintf(&bp, ")");
> +    if (av_bprint_is_complete(&bp))
> +        av_log(ctx, AV_LOG_VERBOSE, "Channel %d: gain: %s\n", channel,
> bp.str);
> +    else
> +        av_log(ctx, AV_LOG_WARNING, "Not enough memory to build gain
> formula; no consequence on output.\n");
> +    av_bprint_finalize(&bp, NULL);
> +}
> +
> +static inline int unpack_sample_uint8_t(uint8_t s)
> +{
> +    return ((int)s - 0x80) << FPOINT;
> +}
> +
> +static inline int unpack_sample_int16_t(int16_t s)
> +{
> +    return ((int)s) << FPOINT;
> +}
> +
> +static inline int64_t unpack_sample_int32_t(int32_t s)
> +{
> +    return ((int64_t)s) << FPOINT;
> +}
> +
> +static inline float unpack_sample_float(float s)
> +{
> +    return s;
> +}
> +
> +static inline double unpack_sample_double(double s)
> +{
> +    return s;
> +}
> +
> +#define SOFTCLIP_PLUS(A, D, x) (A - D * D / ((x) - (A - 2 * D)))
> +#define SOFTCLIP_INT(A, D, x) ((x) < -A + D ? -SOFTCLIP_PLUS(A, D, -(x)) :
> \
> +                               (x) > +A - D ? +SOFTCLIP_PLUS(A, D, +(x)) :
> (x))
> +#define SOFTCLIP(A, x) SOFTCLIP_INT(A, (A / 8), x)
> +
> +static inline uint8_t pack_sample_uint8_t(int gain, int v)
> +{
> +    v >>= FPOINT;
> +    v = fpmul_int(gain, v);
> +    return 0x80 + SOFTCLIP(((int64_t)0x80), v);
> +}
> +
> +static inline int16_t pack_sample_int16_t(int gain, int v)
> +{
> +    v >>= FPOINT;
> +    v = fpmul_int(gain, v);
> +    return SOFTCLIP(0x8000, v);
> +}
> +
> +static inline int32_t pack_sample_int32_t(int gain, int64_t v)
> +{
> +    v >>= FPOINT;
> +    v = fpmul_int64_t(gain, v);
> +    return SOFTCLIP(((int64_t)0x80000000), v);
> +}
> +
> +static inline float pack_sample_float(int gain, float v)
> +{
> +    v = fpmul_float(gain, v);
> +    return SOFTCLIP(1.0, v);
> +}
> +
> +static inline double pack_sample_double(int gain, double v)
> +{
> +    v = fpmul_double(gain, v);
> +    return SOFTCLIP(1.0, v);
> +}
> +
> +#define PROCESS_SAMPLE(type_smp, type_calc) do { \
> +    type_calc iv = unpack_sample_##type_smp(*(samples_in++)); \
> +    type_calc outval = 0; \
> +    for (; band->k; band++) { \
> +        type_calc ov1 = *val += fpmul_##type_calc(band->k, iv - *val); \
> +        type_calc ov2 = iv - ov1; \
> +        val++; \
> +        outval += fpmul_##type_calc(band->g, ov2); \
> +        iv = ov1; \
> +    } \
> +    outval += fpmul_##type_calc(band->g, iv); \
> +    band++; \
> +    *(samples_out++) = pack_sample_##type_smp(s->global_gain, outval); \
> +} while (0);
> +
> +#define DEFINE_PROCESS_SAMPLES(fmt, type_smp, type_calc) \
> +static void process_samples_##fmt(EqContext *s, unsigned channels, unsigned
> samples, \
> +                                  AVFrame *frame, AVFrame *frame_out) \
> +{ \
> +    type_smp *samples_in = (type_smp *)frame->data[0]; \
> +    type_smp *samples_out = (type_smp *)frame_out->data[0]; \
> +    unsigned sample, ch; \
> +    for (sample = 0; sample < samples; sample++) { \
> +        type_calc *val = s->vals; \
> +        EqBand *band = s->bands; \
> +        for (ch = 0; ch < channels; ch++) { \
> +            if (s->channels == 1) \
> +                band = s->bands; \
> +            PROCESS_SAMPLE(type_smp, type_calc); \
> +        } \
> +    } \
> +} \
> +static void process_samples_##fmt##p(EqContext *s, unsigned channels,
> unsigned samples, \
> +                                     AVFrame *frame, AVFrame *frame_out) \
> +{ \
> +    EqBand *band0, *band; \
> +    unsigned sample, ch; \
> +    type_calc *val0, *val; \
> +    band0 = s->bands; \
> +    val0 = s->vals; \
> +    for (ch = 0; ch < channels; ch++) { \
> +        type_smp *samples_in = (type_smp *)frame->data[ch]; \
> +        type_smp *samples_out = (type_smp *)frame_out->data[ch]; \
> +        for (sample = 0; sample < samples; sample++) { \
> +            band = band0; \
> +            val = val0; \
> +            PROCESS_SAMPLE(type_smp, type_calc); \
> +        } \
> +        if (s->channels > 1) \
> +            band0 = band; \
> +        val0 = val; \
> +    } \
> +}
> +
> +DEFINE_PROCESS_SAMPLES(u8,  uint8_t, int)
> +DEFINE_PROCESS_SAMPLES(s16, int16_t, int)
> +DEFINE_PROCESS_SAMPLES(s32, int32_t, int64_t)
> +DEFINE_PROCESS_SAMPLES(flt, float,   float)
> +DEFINE_PROCESS_SAMPLES(dbl, double,  double)
> +
> +static int filter_frame(AVFilterLink *link, AVFrame *frame)
> +{
> +    AVFilterContext *ctx = link->dst;
> +    EqContext *s = ctx->priv;
> +    AVFrame *frame_out;
> +
> +    if (av_frame_is_writable(frame)) {
> +        frame_out = frame;
> +    } else {
> +        frame_out = ff_get_audio_buffer(ctx->outputs[0],
> frame->nb_samples);
> +        if (!frame_out) {
> +            av_frame_free(&frame);
> +            return AVERROR(ENOMEM);
> +        }
> +    }
> +    s->process_samples(s, link->channels, frame->nb_samples, frame,
> frame_out);
> +    if (frame_out != frame)
> +        av_frame_free(&frame);
> +    return ff_filter_frame(ctx->outputs[0], frame_out);
> +}
> +
> +static int init(AVFilterContext *ctx)
> +{
> +    EqContext *s = ctx->priv;
> +    EqBand *bands = NULL, *band;
> +    unsigned nb_bands = 0, nb_channels = 1;
> +    const char *defs = s->defs ? s->defs : "1";
> +    char *tail;
> +    double gain, freq;
> +    int ret = AVERROR(EINVAL);
> +
> +    s->global_gain = round(s->global_gain_opt * (1 << FPOINT));
> +    while (1) {
> +        gain = av_strtod(defs, &tail);
> +        if (tail == defs) {
> +            av_log(ctx, AV_LOG_ERROR, "Error parsing defs: expected gain,
> got \"%s\"\n", defs);
> +            goto error;
> +        }
> +        defs = tail;
> +        band = av_dynarray2_add((void **)&bands, &nb_bands, sizeof(*bands),
> NULL);
> +        if (!band)
> +            return AVERROR(ENOMEM);
> +        band->g = round(gain * (1 << FPOINT));
> +        band->k = 0;
> +        if (*defs == 0)
> +            break;
> +        if (*defs == '/') {
> +            defs++;
> +            freq = av_strtod(defs, &tail);
> +            if (tail == defs || *tail != '/') {
> +                av_log(ctx, AV_LOG_ERROR, "Error parsing defs: expected
> /freq/, got \"%s\"\n", defs - 1);
> +                goto error;
> +            }
> +            defs = tail + 1;
> +            band->k = round(1000 * freq);
> +        } else if (*defs == '|') {
> +            defs++;
> +            nb_channels++;
> +        }
> +    }
> +    av_assert0(nb_bands > 0);
> +    s->bands = av_realloc_f(bands, nb_bands, sizeof(*bands));
> +    if (!s->bands)
> +        return AVERROR(ENOMEM);
> +    s->channels = nb_channels;
> +    return 0;
> +
> +error:
> +    av_free(bands);
> +    return ret;
> +}
> +
> +static void uninit(AVFilterContext *ctx)
> +{
> +    EqContext *s = ctx->priv;
> +    av_freep(&s->vals);
> +    av_freep(&s->bands);
> +}
> +
> +static int query_formats(AVFilterContext *ctx)
> +{
> +    EqContext *s = ctx->priv;
> +    static const enum AVSampleFormat sample_fmts[] = {
> +        AV_SAMPLE_FMT_U8,
> +        AV_SAMPLE_FMT_U8P,
> +        AV_SAMPLE_FMT_S16,
> +        AV_SAMPLE_FMT_S16P,
> +        AV_SAMPLE_FMT_S32,
> +        AV_SAMPLE_FMT_S32P,
> +        AV_SAMPLE_FMT_FLT,
> +        AV_SAMPLE_FMT_FLTP,
> +        AV_SAMPLE_FMT_DBL,
> +        AV_SAMPLE_FMT_DBLP,

So all sample formats expect S64(P), you could use less lines to list them.


> +        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;
> +    if (s->channels == 1) {
> +        layouts = ff_all_channel_counts();
> +    } else {
> +        ret = ff_add_channel_layout(&layouts,
> FF_COUNT2LAYOUT(s->channels));
> +        if (ret < 0)
> +            return ret;
> +    }
> +    if (!layouts)
> +        return AVERROR(ENOMEM);
> +    ret = ff_set_common_channel_layouts(ctx, layouts);
> +    if (ret < 0)
> +        return ret;
> +    return 0;
> +}
> +
> +static int config_props(AVFilterLink *link)
> +{
> +    AVFilterContext *ctx = link->dst;
> +    EqContext *s = ctx->priv;
> +    EqBand *band0, *band;
> +    unsigned nb_vals, ch, val_size;
> +
> +    av_assert0(ctx->outputs[0]->channels == ctx->inputs[0]->channels);
> +    av_assert0(s->channels == 1 || link->channels == s->channels);
> +    nb_vals = 0;
> +    if (!s->keep_all_bands) {
> +        band0 = band = s->bands;
> +        for (ch = 0; ch < s->channels; ch++) {
> +            while (band->k && band->k <= link->sample_rate * 250)
> +                *(band0++) = *(band++);
> +            while (band->k)
> +                av_log(ctx, AV_LOG_VERBOSE, "Channel %d: skipping band %f
> Hz\n",
> +                       ch, (band++)->k / 1000.0);
> +            (band0++)->k = 0;
> +            band++;
> +        }
> +    }
> +    band0 = s->bands;
> +    for (ch = 0; ch < s->channels; ch++) {
> +        band = band0;
> +        while ((band++)->k)
> +            nb_vals++;
> +        reorder_bands(band0, band - band0, link->sample_rate);
> +        print_gain(ctx, ch, band0);
> +        band0 = band;
> +    }
> +    switch (link->format) {
> +    case AV_SAMPLE_FMT_U8:   s->process_samples = process_samples_u8;
> break;
> +    case AV_SAMPLE_FMT_U8P:  s->process_samples = process_samples_u8p;
> break;
> +    case AV_SAMPLE_FMT_S16:  s->process_samples = process_samples_s16;
> break;
> +    case AV_SAMPLE_FMT_S16P: s->process_samples = process_samples_s16p;
> break;
> +    case AV_SAMPLE_FMT_S32:  s->process_samples = process_samples_s32;
> break;
> +    case AV_SAMPLE_FMT_S32P: s->process_samples = process_samples_s32p;
> break;
> +    case AV_SAMPLE_FMT_FLT:  s->process_samples = process_samples_flt;
> break;
> +    case AV_SAMPLE_FMT_FLTP: s->process_samples = process_samples_fltp;
> break;
> +    case AV_SAMPLE_FMT_DBL:  s->process_samples = process_samples_dbl;
> break;
> +    case AV_SAMPLE_FMT_DBLP: s->process_samples = process_samples_dblp;
> break;
> +    default: av_assert0(0);
> +    }
> +    switch (link->format) {
> +    case AV_SAMPLE_FMT_U8:
> +    case AV_SAMPLE_FMT_U8P:
> +    case AV_SAMPLE_FMT_S16:
> +    case AV_SAMPLE_FMT_S16P: val_size = sizeof(int);     break;
> +    case AV_SAMPLE_FMT_S32:
> +    case AV_SAMPLE_FMT_S32P: val_size = sizeof(int64_t); break;
> +    case AV_SAMPLE_FMT_FLT:
> +    case AV_SAMPLE_FMT_FLTP: val_size = sizeof(float);   break;
> +    case AV_SAMPLE_FMT_DBL:
> +    case AV_SAMPLE_FMT_DBLP: val_size = sizeof(double);  break;
> +    default: av_assert0(0);
> +    }
> +    s->vals = av_calloc(nb_vals, val_size);
> +    if (!s->vals)
> +        return AVERROR(ENOMEM);
> +    return 0;
> +}
> +
> +#define OFFSET(x) offsetof(EqContext, x)
> +#define FLAGS AV_OPT_FLAG_AUDIO_PARAM | AV_OPT_FLAG_FILTERING_PARAM
> +static const AVOption mbfequalizer_options[] = {
> +    { "defs", "definitions of the bands", OFFSET(defs), AV_OPT_TYPE_STRING,
> { .str = NULL }, 0, 0, FLAGS },
> +    { "global_gain", "gain applied to all channels",
> OFFSET(global_gain_opt), AV_OPT_TYPE_DOUBLE,
> +                                                     { .dbl = 1 }, 0,
> INT_MAX >> FPOINT, FLAGS },
> +    { "keep_all_bands", "keep bands beyond sample_rate/4",
> OFFSET(keep_all_bands), AV_OPT_TYPE_BOOL,
> +                                                           { .i64 = 0 }, 0,
> 1, FLAGS },
> +    { NULL }
> +};
> +
> +AVFILTER_DEFINE_CLASS(mbfequalizer);
> +
> +static const AVFilterPad mbfequalizer_inputs[] = {
> +    {
> +        .name         = "default",
> +        .type         = AVMEDIA_TYPE_AUDIO,
> +        .config_props = config_props,
> +        .filter_frame = filter_frame,
> +    },
> +    { NULL }
> +};
> +
> +static const AVFilterPad mbfequalizer_outputs[] = {
> +    {
> +        .name = "default",
> +        .type = AVMEDIA_TYPE_AUDIO,
> +    },
> +    { NULL }
> +};
> +
> +AVFilter ff_af_mbfequalizer = {
> +    .name          = "mbfequalizer",
> +    .description   = NULL_IF_CONFIG_SMALL("Multiband fast equalizer"),
> +    .init          = init,
> +    .uninit        = uninit,
> +    .query_formats = query_formats,
> +    .priv_size     = sizeof(EqContext),
> +    .priv_class    = &mbfequalizer_class,
> +    .inputs        = mbfequalizer_inputs,
> +    .outputs       = mbfequalizer_outputs,
> +};
> diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
> index c295f8e403..1f3a1a7827 100644
> --- a/libavfilter/allfilters.c
> +++ b/libavfilter/allfilters.c
> @@ -117,6 +117,7 @@ extern AVFilter ff_af_lowpass;
>  extern AVFilter ff_af_lowshelf;
>  extern AVFilter ff_af_lv2;
>  extern AVFilter ff_af_mcompand;
> +extern AVFilter ff_af_mbfequalizer;
>  extern AVFilter ff_af_pan;
>  extern AVFilter ff_af_replaygain;
>  extern AVFilter ff_af_resample;
> --
> 2.24.0
>
> _______________________________________________
> ffmpeg-devel mailing list
> ffmpeg-devel@ffmpeg.org
> https://ffmpeg.org/mailman/listinfo/ffmpeg-devel
>
> To unsubscribe, visit link above, or email
> ffmpeg-devel-request@ffmpeg.org with subject "unsubscribe".
diff mbox

Patch

diff --git a/doc/filters.texi b/doc/filters.texi
index 8c5d3a5760..c6917fad38 100644
--- a/doc/filters.texi
+++ b/doc/filters.texi
@@ -4383,6 +4383,54 @@  lv2=p=http\\\\://www.openavproductions.com/artyfx#bitta:c=crush=0.3
 @end example
 @end itemize
 
+@section mbfequalizer
+Multiband fast equalizer.
+
+As an equalizer, it allows to selectively alter the volume of the different
+frequencies in the audio signal. This particular equalizer does not allow a
+very sharp separation of the frequencies, but that gives it other
+advantages. It is rather fast. It does not add latency. It supports all
+sample formats, and for integer formats it uses bit-exact integer
+arithmetic.
+
+It accepts the following parameters:
+
+@table @option
+@item defs
+Definition of bands and gain. The syntax is:
+"gain/freq/gain[/freq/gain...][|gain/freq/gain...]"
+where @var{freq} is the cut frequency between bands and @var{gain} the
+associated gain.
+Gain and bands for several channels can be defined by separating them with
+the pipe "|" character; if only one definition is given, it is used for all
+channels; otherwise, the number of definitions must match the number of
+input channels.
+Frequencies should be in ascending order. Gain greater than 1 or negative
+may result in overflows and ugly noise or other strange results.
+
+@item global_gain
+Gain applied to all channels and all bands. Can be greater than 1 without
+risk; in case of overflow, soft clipping is applied.
+
+@item keep_all_bands
+By default, bands beyond @var{sample_rate}/4 are discarded. If this option
+is enabled, they are kept; it can result in strange effects.
+@end table
+
+A Gnuplot formula mapping the frequency to the gain is printed at verbose
+log level. It can be used to check the result, proably with
+"set logscale x; db(x) = 20*log(x)/log(10); plot [20:2000] db(...)".
+
+@subsection Examples
+
+Assuming 3.1 input, reduce high frequencies on left and right, keep a narrow
+band on center and dampen bass on LFE. All gains were lowered by 5dB, then
+raised again using global gain, and the -2.8dB for the middle band was found
+by looking at the graph.
+@example
+mbfequalizer=-5dB/7040/-15dB|-5dB/7040/-15dB|0/220/-2.8dB/880/0|-10dB/110/-5dB:global_gain=+5dB
+@end example
+
 @section mcompand
 Multiband Compress or expand the audio's dynamic range.
 
diff --git a/libavfilter/Makefile b/libavfilter/Makefile
index 37d4eee858..f5690a50db 100644
--- a/libavfilter/Makefile
+++ b/libavfilter/Makefile
@@ -124,6 +124,7 @@  OBJS-$(CONFIG_LOWPASS_FILTER)                += af_biquads.o
 OBJS-$(CONFIG_LOWSHELF_FILTER)               += af_biquads.o
 OBJS-$(CONFIG_LV2_FILTER)                    += af_lv2.o
 OBJS-$(CONFIG_MCOMPAND_FILTER)               += af_mcompand.o
+OBJS-$(CONFIG_MBFEQUALIZER_FILTER)           += af_mbfequalizer.o
 OBJS-$(CONFIG_PAN_FILTER)                    += af_pan.o
 OBJS-$(CONFIG_REPLAYGAIN_FILTER)             += af_replaygain.o
 OBJS-$(CONFIG_RESAMPLE_FILTER)               += af_resample.o
diff --git a/libavfilter/af_mbfequalizer.c b/libavfilter/af_mbfequalizer.c
new file mode 100644
index 0000000000..aadd3cb3da
--- /dev/null
+++ b/libavfilter/af_mbfequalizer.c
@@ -0,0 +1,567 @@ 
+/*
+ * Copyright (c) 2019 Nicolas George <george@nsup.org>
+ *
+ * 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
+ */
+
+/*
+ * Low-pass-high-pass equalizer using a simple sliding weighed average.
+ *
+ * Input signal: u(t)
+ * Low freqs:    a(t)
+ * High freqs:   b(t)
+ *
+ * a(t) = (1-k) a(t-1) + k u(t)
+ * b(t) = u(t) - a(t)
+ *
+ * F = sample frequency; f = input frequency; ω = 2πf/F
+ * For u(t) = exp(iωt),
+ * a(t) =                 k/(1-exp(-iω)(1-k)) exp(iwt)
+ * b(t) = (1-k)(1-exp(-iω))/(1-exp(-iω)(1-k)) exp(iwt)
+ * hi/low = (1-k)(1-exp(-iω)) / k
+ * R(k,w) = (1-k)*(1-exp(-I*w)) / k
+ *
+ * Cut when |hi/low| = 1, i.e. |1-exp(iω)| = k/(1-k)
+ * Triangle 0, 1, exp(iω) has sides 1, 1, k/(1-k)
+ * w = 2 asin(1/2 k/(1-k))
+ * k/(1-k) = 2 sin(ω/2)
+ * k = 2 sin(ω/2) / (1 + 2 sin(ω/2))
+ *
+ * Gnuplot script: response for cut frequency 110 Hz and 440 Hz.
+ * I = {0,1}
+ * F = 48000
+ * w(f) = 2*pi*f/F
+ * db(x) = 20*log(x)/log(10)
+ * set samples 1000
+ * set grid
+ * set logscale x 2
+ * k(f) = 2*sin(w(f)/2) / (1+2*sin(w(f)/2))
+ * d(f, x) = abs(exp(I*w(x))-(1-k(f)))
+ * al(f, x) = k(f) / d(f, x)
+ * ah(f, x) = abs((1-k(f))*(1-exp(I*w(x)))) / d(f, x)
+ * plot [20:16000] [-20:0] db(al(110, x)), db(ah(110, x)), db(al(880, x)), db(ah(880, x))
+ *
+ * Approximation without sin():
+ * ω - ω³/6  ≤  sin(ω)   ≤ ω - ω³/6  + ω⁵/120
+ * ω - ω³/24 ≤ 2sin(ω/2) ≤ ω - ω³/24 + ω⁵/1920
+ * Accurate enough with 0 ≤ ω ≤ π/2.
+ *
+ * Clip to [-A;A] with soft decay on [A(1-k); +∞[:
+ * s(x) = (Ak)² / (x - A(1-2k))
+ */
+
+/*
+
+aevalsrc=sin(13.75*4.35*exp(0.7*t)):d=10
+arealtime,asplit=3[a][b][c];
+[a]showwaves=s=1764x128[ap];
+[b]mbfequalizer=1:0,showwaves=s=1764x128[bp];
+[c]mbfequalizer=0:1,showwaves=s=1764x128[cp];
+[ap][bp][cp]vstack=3
+
+*/
+
+#include "libavutil/avassert.h"
+#include "libavutil/bprint.h"
+#include "libavutil/eval.h"
+#include "libavutil/opt.h"
+#include "avfilter.h"
+#include "audio.h"
+#include "internal.h"
+
+typedef struct EqBand {
+    int k; /* k(2πf) */
+    int g; /* gain of the high band */
+} EqBand;
+
+typedef struct EqContext EqContext;
+struct EqContext {
+    const AVClass *class;
+    char *defs;
+    EqBand *bands;
+    void *vals;
+    void (*process_samples)(EqContext *s, unsigned channels, unsigned samples,
+                            AVFrame *frame, AVFrame *frame_out);
+    double global_gain_opt;
+    unsigned channels;
+    int global_gain;
+    int keep_all_bands;
+};
+
+#define FPOINT 15
+#define SOFTCLIP_SHIFT 2 /* k = 1/(1<<SHIFT) */
+
+#if 1&&HAVE_FAST_64BIT
+
+static inline int fpmul_int(int a, int b)
+{
+    return ((int64_t)a * b) >> FPOINT;
+}
+
+#else
+
+static inline int fpmul_int(int a, int b)
+{
+#define MASK ((1 << FPOINT) - 1)
+    return (((a >> FPOINT) * (b >> FPOINT)) << FPOINT)
+         + (a >> FPOINT) * (b & MASK)
+         + (b >> FPOINT) * (a & MASK)
+         + (((a & MASK) * (b & MASK)) >> FPOINT);
+#undef MASK
+}
+
+#endif
+
+static int64_t fpmul_int64_t(int a, int64_t b)
+{
+    uint64_t a0 = a & 0xFFFFFFFF;
+    uint64_t b0 = b & 0xFFFFFFFF;
+    /*int64_t  a1 = a >> 32;*/ /* a1 = 0 */
+    int64_t  b1 = b >> 32;
+    uint64_t r0 = a0 * b0;
+    int64_t r1a = (int64_t)a0 * b1;
+    /*int64_t r1b = (int64_t)b0 * a1;*/
+    int64_t r2  = 0 /*a1 * b1*/;
+    r2 += (r1a >> 32) /*+ (r1b >> 32)*/;
+    r1a &= 0xFFFFFFFF;
+    /*r1b &= 0xFFFFFFFF;*/
+    r2 += (r1a /*+ r1b*/ + (r0 >> 32)) >> 32;
+    r0 += (r1a /*+ r1b*/) << 32;
+    return (int64_t)(r0 >> FPOINT) + (r2 << (64 - FPOINT));
+}
+
+static inline float fpmul_float(int a, float b)
+{
+    return a * b / (1 << FPOINT);
+}
+
+static inline double fpmul_double(int a, double b)
+{
+    return a * b / (1 << FPOINT);
+}
+
+static int approx_factor(int f, int sf)
+{
+    /* 31089 / 151 =~ 2*Pi*32768 / 1000 */
+    int w = av_rescale(f, 31089, 151 * sf);
+    int s = w - fpmul_int(w, fpmul_int(w, w)) / 24;
+    int r = (s << FPOINT) / ((1 << FPOINT) + s);
+    return r;
+}
+
+static void reorder_bands(EqBand *bands, size_t nb_bands, unsigned sample_rate)
+{
+    size_t i, j;
+
+    if (nb_bands == 1)
+        return;
+    /* bands arrive in ascending order with kf = gain below,
+       and k = 0 on the final band;
+       we want them in descending order with kf = gain above;
+       therefore, reverse the order of k except the final 0,
+       and reverse the order of kf including the final 0 */
+    for (i = 0, j = nb_bands - 2; i < j; i++, j--)
+        FFSWAP(int, bands[i].k, bands[j].k);
+    for (i = 0, j = nb_bands - 1; i < j; i++, j--)
+        FFSWAP(int, bands[i].g, bands[j].g);
+    for (i = 0; i < nb_bands - 1; i++)
+        bands[i].k = approx_factor(bands[i].k, sample_rate);
+}
+
+static void gain_to_bprint(AVBPrint *bp, unsigned rate, EqBand *bands)
+{
+    double k = bands->k / (double)(1 << FPOINT);
+    double g = bands->g / (double)(1 << FPOINT);
+
+    if (bands->k == 0) {
+        av_bprintf(bp, "%.5f", g);
+        return;
+    }
+    av_bprintf(bp, "(%.5f*%.5f*(1-exp(-{0,1}*2*pi*x/%d))+%.5f*", g, 1 - k, rate, k);
+    gain_to_bprint(bp, rate, bands + 1);
+    av_bprintf(bp, ")");
+    av_bprintf(bp, "/(1-exp(-{0,1}*2*pi*x/%d)*%.5f)", rate, 1 - k);
+}
+
+static void print_gain(AVFilterContext *ctx, unsigned channel, EqBand *bands)
+{
+    EqContext *s = ctx->priv;
+    AVBPrint bp;
+    unsigned rate = ctx->outputs[0]->sample_rate;
+
+    av_bprint_init(&bp, 0, AV_BPRINT_SIZE_AUTOMATIC);
+    if (s->global_gain != 1 << FPOINT)
+        av_bprintf(&bp, "%.5f*", s->global_gain / (double)(1 << FPOINT));
+    av_bprintf(&bp, "abs(");
+    gain_to_bprint(&bp, rate, bands);
+    av_bprintf(&bp, ")");
+    if (av_bprint_is_complete(&bp))
+        av_log(ctx, AV_LOG_VERBOSE, "Channel %d: gain: %s\n", channel, bp.str);
+    else
+        av_log(ctx, AV_LOG_WARNING, "Not enough memory to build gain formula; no consequence on output.\n");
+    av_bprint_finalize(&bp, NULL);
+}
+
+static inline int unpack_sample_uint8_t(uint8_t s)
+{
+    return ((int)s - 0x80) << FPOINT;
+}
+
+static inline int unpack_sample_int16_t(int16_t s)
+{
+    return ((int)s) << FPOINT;
+}
+
+static inline int64_t unpack_sample_int32_t(int32_t s)
+{
+    return ((int64_t)s) << FPOINT;
+}
+
+static inline float unpack_sample_float(float s)
+{
+    return s;
+}
+
+static inline double unpack_sample_double(double s)
+{
+    return s;
+}
+
+#define SOFTCLIP_PLUS(A, D, x) (A - D * D / ((x) - (A - 2 * D)))
+#define SOFTCLIP_INT(A, D, x) ((x) < -A + D ? -SOFTCLIP_PLUS(A, D, -(x)) : \
+                               (x) > +A - D ? +SOFTCLIP_PLUS(A, D, +(x)) : (x))
+#define SOFTCLIP(A, x) SOFTCLIP_INT(A, (A / 8), x)
+
+static inline uint8_t pack_sample_uint8_t(int gain, int v)
+{
+    v >>= FPOINT;
+    v = fpmul_int(gain, v);
+    return 0x80 + SOFTCLIP(((int64_t)0x80), v);
+}
+
+static inline int16_t pack_sample_int16_t(int gain, int v)
+{
+    v >>= FPOINT;
+    v = fpmul_int(gain, v);
+    return SOFTCLIP(0x8000, v);
+}
+
+static inline int32_t pack_sample_int32_t(int gain, int64_t v)
+{
+    v >>= FPOINT;
+    v = fpmul_int64_t(gain, v);
+    return SOFTCLIP(((int64_t)0x80000000), v);
+}
+
+static inline float pack_sample_float(int gain, float v)
+{
+    v = fpmul_float(gain, v);
+    return SOFTCLIP(1.0, v);
+}
+
+static inline double pack_sample_double(int gain, double v)
+{
+    v = fpmul_double(gain, v);
+    return SOFTCLIP(1.0, v);
+}
+
+#define PROCESS_SAMPLE(type_smp, type_calc) do { \
+    type_calc iv = unpack_sample_##type_smp(*(samples_in++)); \
+    type_calc outval = 0; \
+    for (; band->k; band++) { \
+        type_calc ov1 = *val += fpmul_##type_calc(band->k, iv - *val); \
+        type_calc ov2 = iv - ov1; \
+        val++; \
+        outval += fpmul_##type_calc(band->g, ov2); \
+        iv = ov1; \
+    } \
+    outval += fpmul_##type_calc(band->g, iv); \
+    band++; \
+    *(samples_out++) = pack_sample_##type_smp(s->global_gain, outval); \
+} while (0);
+
+#define DEFINE_PROCESS_SAMPLES(fmt, type_smp, type_calc) \
+static void process_samples_##fmt(EqContext *s, unsigned channels, unsigned samples, \
+                                  AVFrame *frame, AVFrame *frame_out) \
+{ \
+    type_smp *samples_in = (type_smp *)frame->data[0]; \
+    type_smp *samples_out = (type_smp *)frame_out->data[0]; \
+    unsigned sample, ch; \
+    for (sample = 0; sample < samples; sample++) { \
+        type_calc *val = s->vals; \
+        EqBand *band = s->bands; \
+        for (ch = 0; ch < channels; ch++) { \
+            if (s->channels == 1) \
+                band = s->bands; \
+            PROCESS_SAMPLE(type_smp, type_calc); \
+        } \
+    } \
+} \
+static void process_samples_##fmt##p(EqContext *s, unsigned channels, unsigned samples, \
+                                     AVFrame *frame, AVFrame *frame_out) \
+{ \
+    EqBand *band0, *band; \
+    unsigned sample, ch; \
+    type_calc *val0, *val; \
+    band0 = s->bands; \
+    val0 = s->vals; \
+    for (ch = 0; ch < channels; ch++) { \
+        type_smp *samples_in = (type_smp *)frame->data[ch]; \
+        type_smp *samples_out = (type_smp *)frame_out->data[ch]; \
+        for (sample = 0; sample < samples; sample++) { \
+            band = band0; \
+            val = val0; \
+            PROCESS_SAMPLE(type_smp, type_calc); \
+        } \
+        if (s->channels > 1) \
+            band0 = band; \
+        val0 = val; \
+    } \
+}
+
+DEFINE_PROCESS_SAMPLES(u8,  uint8_t, int)
+DEFINE_PROCESS_SAMPLES(s16, int16_t, int)
+DEFINE_PROCESS_SAMPLES(s32, int32_t, int64_t)
+DEFINE_PROCESS_SAMPLES(flt, float,   float)
+DEFINE_PROCESS_SAMPLES(dbl, double,  double)
+
+static int filter_frame(AVFilterLink *link, AVFrame *frame)
+{
+    AVFilterContext *ctx = link->dst;
+    EqContext *s = ctx->priv;
+    AVFrame *frame_out;
+
+    if (av_frame_is_writable(frame)) {
+        frame_out = frame;
+    } else {
+        frame_out = ff_get_audio_buffer(ctx->outputs[0], frame->nb_samples);
+        if (!frame_out) {
+            av_frame_free(&frame);
+            return AVERROR(ENOMEM);
+        }
+    }
+    s->process_samples(s, link->channels, frame->nb_samples, frame, frame_out);
+    if (frame_out != frame)
+        av_frame_free(&frame);
+    return ff_filter_frame(ctx->outputs[0], frame_out);
+}
+
+static int init(AVFilterContext *ctx)
+{
+    EqContext *s = ctx->priv;
+    EqBand *bands = NULL, *band;
+    unsigned nb_bands = 0, nb_channels = 1;
+    const char *defs = s->defs ? s->defs : "1";
+    char *tail;
+    double gain, freq;
+    int ret = AVERROR(EINVAL);
+
+    s->global_gain = round(s->global_gain_opt * (1 << FPOINT));
+    while (1) {
+        gain = av_strtod(defs, &tail);
+        if (tail == defs) {
+            av_log(ctx, AV_LOG_ERROR, "Error parsing defs: expected gain, got \"%s\"\n", defs);
+            goto error;
+        }
+        defs = tail;
+        band = av_dynarray2_add((void **)&bands, &nb_bands, sizeof(*bands), NULL);
+        if (!band)
+            return AVERROR(ENOMEM);
+        band->g = round(gain * (1 << FPOINT));
+        band->k = 0;
+        if (*defs == 0)
+            break;
+        if (*defs == '/') {
+            defs++;
+            freq = av_strtod(defs, &tail);
+            if (tail == defs || *tail != '/') {
+                av_log(ctx, AV_LOG_ERROR, "Error parsing defs: expected /freq/, got \"%s\"\n", defs - 1);
+                goto error;
+            }
+            defs = tail + 1;
+            band->k = round(1000 * freq);
+        } else if (*defs == '|') {
+            defs++;
+            nb_channels++;
+        }
+    }
+    av_assert0(nb_bands > 0);
+    s->bands = av_realloc_f(bands, nb_bands, sizeof(*bands));
+    if (!s->bands)
+        return AVERROR(ENOMEM);
+    s->channels = nb_channels;
+    return 0;
+
+error:
+    av_free(bands);
+    return ret;
+}
+
+static void uninit(AVFilterContext *ctx)
+{
+    EqContext *s = ctx->priv;
+    av_freep(&s->vals);
+    av_freep(&s->bands);
+}
+
+static int query_formats(AVFilterContext *ctx)
+{
+    EqContext *s = ctx->priv;
+    static const enum AVSampleFormat sample_fmts[] = {
+        AV_SAMPLE_FMT_U8,
+        AV_SAMPLE_FMT_U8P,
+        AV_SAMPLE_FMT_S16,
+        AV_SAMPLE_FMT_S16P,
+        AV_SAMPLE_FMT_S32,
+        AV_SAMPLE_FMT_S32P,
+        AV_SAMPLE_FMT_FLT,
+        AV_SAMPLE_FMT_FLTP,
+        AV_SAMPLE_FMT_DBL,
+        AV_SAMPLE_FMT_DBLP,
+        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;
+    if (s->channels == 1) {
+        layouts = ff_all_channel_counts();
+    } else {
+        ret = ff_add_channel_layout(&layouts, FF_COUNT2LAYOUT(s->channels));
+        if (ret < 0)
+            return ret;
+    }
+    if (!layouts)
+        return AVERROR(ENOMEM);
+    ret = ff_set_common_channel_layouts(ctx, layouts);
+    if (ret < 0)
+        return ret;
+    return 0;
+}
+
+static int config_props(AVFilterLink *link)
+{
+    AVFilterContext *ctx = link->dst;
+    EqContext *s = ctx->priv;
+    EqBand *band0, *band;
+    unsigned nb_vals, ch, val_size;
+
+    av_assert0(ctx->outputs[0]->channels == ctx->inputs[0]->channels);
+    av_assert0(s->channels == 1 || link->channels == s->channels);
+    nb_vals = 0;
+    if (!s->keep_all_bands) {
+        band0 = band = s->bands;
+        for (ch = 0; ch < s->channels; ch++) {
+            while (band->k && band->k <= link->sample_rate * 250)
+                *(band0++) = *(band++);
+            while (band->k)
+                av_log(ctx, AV_LOG_VERBOSE, "Channel %d: skipping band %f Hz\n",
+                       ch, (band++)->k / 1000.0);
+            (band0++)->k = 0;
+            band++;
+        }
+    }
+    band0 = s->bands;
+    for (ch = 0; ch < s->channels; ch++) {
+        band = band0;
+        while ((band++)->k)
+            nb_vals++;
+        reorder_bands(band0, band - band0, link->sample_rate);
+        print_gain(ctx, ch, band0);
+        band0 = band;
+    }
+    switch (link->format) {
+    case AV_SAMPLE_FMT_U8:   s->process_samples = process_samples_u8;   break;
+    case AV_SAMPLE_FMT_U8P:  s->process_samples = process_samples_u8p;  break;
+    case AV_SAMPLE_FMT_S16:  s->process_samples = process_samples_s16;  break;
+    case AV_SAMPLE_FMT_S16P: s->process_samples = process_samples_s16p; break;
+    case AV_SAMPLE_FMT_S32:  s->process_samples = process_samples_s32;  break;
+    case AV_SAMPLE_FMT_S32P: s->process_samples = process_samples_s32p; break;
+    case AV_SAMPLE_FMT_FLT:  s->process_samples = process_samples_flt;  break;
+    case AV_SAMPLE_FMT_FLTP: s->process_samples = process_samples_fltp; break;
+    case AV_SAMPLE_FMT_DBL:  s->process_samples = process_samples_dbl;  break;
+    case AV_SAMPLE_FMT_DBLP: s->process_samples = process_samples_dblp; break;
+    default: av_assert0(0);
+    }
+    switch (link->format) {
+    case AV_SAMPLE_FMT_U8:
+    case AV_SAMPLE_FMT_U8P:
+    case AV_SAMPLE_FMT_S16:
+    case AV_SAMPLE_FMT_S16P: val_size = sizeof(int);     break;
+    case AV_SAMPLE_FMT_S32:
+    case AV_SAMPLE_FMT_S32P: val_size = sizeof(int64_t); break;
+    case AV_SAMPLE_FMT_FLT:
+    case AV_SAMPLE_FMT_FLTP: val_size = sizeof(float);   break;
+    case AV_SAMPLE_FMT_DBL:
+    case AV_SAMPLE_FMT_DBLP: val_size = sizeof(double);  break;
+    default: av_assert0(0);
+    }
+    s->vals = av_calloc(nb_vals, val_size);
+    if (!s->vals)
+        return AVERROR(ENOMEM);
+    return 0;
+}
+
+#define OFFSET(x) offsetof(EqContext, x)
+#define FLAGS AV_OPT_FLAG_AUDIO_PARAM | AV_OPT_FLAG_FILTERING_PARAM
+static const AVOption mbfequalizer_options[] = {
+    { "defs", "definitions of the bands", OFFSET(defs), AV_OPT_TYPE_STRING, { .str = NULL }, 0, 0, FLAGS },
+    { "global_gain", "gain applied to all channels", OFFSET(global_gain_opt), AV_OPT_TYPE_DOUBLE,
+                                                     { .dbl = 1 }, 0, INT_MAX >> FPOINT, FLAGS },
+    { "keep_all_bands", "keep bands beyond sample_rate/4", OFFSET(keep_all_bands), AV_OPT_TYPE_BOOL,
+                                                           { .i64 = 0 }, 0, 1, FLAGS },
+    { NULL }
+};
+
+AVFILTER_DEFINE_CLASS(mbfequalizer);
+
+static const AVFilterPad mbfequalizer_inputs[] = {
+    {
+        .name         = "default",
+        .type         = AVMEDIA_TYPE_AUDIO,
+        .config_props = config_props,
+        .filter_frame = filter_frame,
+    },
+    { NULL }
+};
+
+static const AVFilterPad mbfequalizer_outputs[] = {
+    {
+        .name = "default",
+        .type = AVMEDIA_TYPE_AUDIO,
+    },
+    { NULL }
+};
+
+AVFilter ff_af_mbfequalizer = {
+    .name          = "mbfequalizer",
+    .description   = NULL_IF_CONFIG_SMALL("Multiband fast equalizer"),
+    .init          = init,
+    .uninit        = uninit,
+    .query_formats = query_formats,
+    .priv_size     = sizeof(EqContext),
+    .priv_class    = &mbfequalizer_class,
+    .inputs        = mbfequalizer_inputs,
+    .outputs       = mbfequalizer_outputs,
+};
diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
index c295f8e403..1f3a1a7827 100644
--- a/libavfilter/allfilters.c
+++ b/libavfilter/allfilters.c
@@ -117,6 +117,7 @@  extern AVFilter ff_af_lowpass;
 extern AVFilter ff_af_lowshelf;
 extern AVFilter ff_af_lv2;
 extern AVFilter ff_af_mcompand;
+extern AVFilter ff_af_mbfequalizer;
 extern AVFilter ff_af_pan;
 extern AVFilter ff_af_replaygain;
 extern AVFilter ff_af_resample;