[FFmpeg-devel] avfilter: add arbitrary audio FIR filter

Signed-off-by: Paul B Mahol <onemda@gmail.com>
---
 doc/filters.texi            |   7 +
 libavfilter/Makefile        |   1 +
 libavfilter/af_afirfilter.c | 411 ++++++++++++++++++++++++++++++++++++++++++++
 libavfilter/allfilters.c    |   1 +
 4 files changed, 420 insertions(+)
 create mode 100644 libavfilter/af_afirfilter.c

On Mon, May 1, 2017 at 5:02 AM, Paul B Mahol <onemda@gmail.com> wrote:
> Signed-off-by: Paul B Mahol <onemda@gmail.com>
> ---
>  doc/filters.texi            |   7 +
>  libavfilter/Makefile        |   1 +
>  libavfilter/af_afirfilter.c | 411 ++++++++++++++++++++++++++++++++++++++++++++
>  libavfilter/allfilters.c    |   1 +
>  4 files changed, 420 insertions(+)
>  create mode 100644 libavfilter/af_afirfilter.c
>
> diff --git a/doc/filters.texi b/doc/filters.texi
> index 119e747..d0f6cc4 100644
> --- a/doc/filters.texi
> +++ b/doc/filters.texi
> @@ -878,6 +878,13 @@ afftfilt="1-clip((b/nb)*b,0,1)"
>  @end example
>  @end itemize
>
> +@section afirfilter
> +
> +Apply an Arbitary Frequency Impulse Response filter.
> +
> +This filter uses second stream as FIR coefficients.
> +Even channels hold real and odds channels hold imaginary coefficients.
> +
>  @anchor{aformat}
>  @section aformat
>
> diff --git a/libavfilter/Makefile b/libavfilter/Makefile
> index 66c36e4..1a0f24b 100644
> --- a/libavfilter/Makefile
> +++ b/libavfilter/Makefile
> @@ -38,6 +38,7 @@ OBJS-$(CONFIG_AEMPHASIS_FILTER)              += af_aemphasis.o
>  OBJS-$(CONFIG_AEVAL_FILTER)                  += aeval.o
>  OBJS-$(CONFIG_AFADE_FILTER)                  += af_afade.o
>  OBJS-$(CONFIG_AFFTFILT_FILTER)               += af_afftfilt.o window_func.o
> +OBJS-$(CONFIG_AFIRFILTER_FILTER)             += af_afirfilter.o
>  OBJS-$(CONFIG_AFORMAT_FILTER)                += af_aformat.o
>  OBJS-$(CONFIG_AGATE_FILTER)                  += af_agate.o
>  OBJS-$(CONFIG_AINTERLEAVE_FILTER)            += f_interleave.o
> diff --git a/libavfilter/af_afirfilter.c b/libavfilter/af_afirfilter.c
> new file mode 100644
> index 0000000..e127579
> --- /dev/null
> +++ b/libavfilter/af_afirfilter.c
> @@ -0,0 +1,411 @@
> +/*
> + * 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
> + */
> +
> +/**
> + * @file
> + * An arbitrary audio FIR filter
> + */
> +
> +#include "libavutil/audio_fifo.h"
> +#include "libavutil/avassert.h"
> +#include "libavutil/channel_layout.h"
> +#include "libavutil/common.h"
> +#include "libavutil/opt.h"
> +#include "libavcodec/avfft.h"
> +
> +#include "audio.h"
> +#include "avfilter.h"
> +#include "formats.h"
> +#include "hermite.h"
> +#include "internal.h"
> +
> +typedef struct FIRContext {
> +    const AVClass *class;
> +
> +    int n;
> +    int eof_coeffs;
> +    int have_coeffs;
> +    int nb_taps;
> +    int fft_length;
> +    int nb_channels;
> +    int one2many;
> +
> +    FFTContext *fft, *ifft;
> +    FFTComplex **fft_data;
> +    FFTComplex **fft_coef;
> +
> +    AVAudioFifo *fifo[2];
> +    AVFrame *in[2];
> +    AVFrame *buffer;
> +    int64_t pts;
> +    int hop_size;
> +    int start, end;
> +} FIRContext;
> +
> +static int fir_filter(FIRContext *s, AVFilterLink *outlink)
> +{
> +    AVFilterContext *ctx = outlink->src;
> +    int start = s->start, end = s->end;
> +    int ret = 0, n, ch, j, k;
> +    int nb_samples;
> +    AVFrame *out;
> +
> +    nb_samples = FFMIN(s->fft_length, av_audio_fifo_size(s->fifo[0]));
> +
> +    s->in[0] = ff_get_audio_buffer(ctx->inputs[0], nb_samples);
> +    if (!s->in[0])
> +        return AVERROR(ENOMEM);
> +    av_audio_fifo_peek(s->fifo[0], (void **)s->in[0]->data, nb_samples);
> +
> +    for (ch = 0; ch < outlink->channels; ch++) {
> +        const float *src = (float *)s->in[0]->extended_data[ch];
> +        float *buf = (float *)s->buffer->extended_data[ch];
> +        FFTComplex *fft_data = s->fft_data[ch];
> +        FFTComplex *fft_coef = s->fft_coef[ch];
> +
> +        memset(fft_data, 0, sizeof(*fft_data) * s->fft_length);
> +        for (n = 0; n < s->fft_length; n++) {
> +            fft_data[n].re = src[n];
> +            fft_data[n].im = 0;
> +        }
> +
> +        av_fft_permute(s->fft, fft_data);
> +        av_fft_calc(s->fft, fft_data);
> +
> +        fft_data[0].re *= fft_coef[0].re;
> +        fft_data[0].im *= fft_coef[0].im;
> +        for (n = 1; n < s->fft_length; n++) {
> +            const float re = fft_data[n].re;
> +            const float im = fft_data[n].im;
> +
> +            fft_data[n].re = re * fft_coef[n].re - im * fft_coef[n].im;
> +            fft_data[n].im = re * fft_coef[n].im + im * fft_coef[n].re;
> +        }
> +
> +        av_fft_permute(s->ifft, fft_data);
> +        av_fft_calc(s->ifft, fft_data);
> +
> +        start = s->start;
> +        end = s->end;
> +        k = end;
> +
> +        for (n = 0, j = start; j < k && n < s->fft_length; n++, j++) {
> +            buf[j] = fft_data[n].re;
> +        }
> +
> +        for (; n < s->fft_length; n++, j++) {
> +            buf[j] = fft_data[n].re;
> +        }
> +
> +        start += s->hop_size;
> +        end = j;
> +    }
> +
> +    s->start = start;
> +    s->end   = end;
> +
> +    if (start >= s->fft_length) {
> +        float *dst, *buf;
> +
> +        start -= s->fft_length;
> +        end   -= s->fft_length;
> +
> +        s->start = start;
> +        s->end = end;
> +
> +        out = ff_get_audio_buffer(outlink, s->fft_length);
> +        if (!out)
> +            return AVERROR(ENOMEM);
> +
> +        out->pts = s->pts;
> +        s->pts += s->fft_length;
> +
> +        for (ch = 0; ch < s->nb_channels; ch++) {
> +            dst = (float *)out->extended_data[ch];
> +            buf = (float *)s->buffer->extended_data[ch];
> +
> +            for (n = 0; n < s->fft_length; n++)
> +                dst[n] = buf[n];
> +            memmove(buf, buf + s->fft_length, s->fft_length * 4);
> +        }

Is this overlap-save?



> +
> +        ret = ff_filter_frame(outlink, out);
> +    }
> +
> +    av_audio_fifo_drain(s->fifo[0], s->hop_size);
> +    av_frame_free(&s->in[0]);
> +
> +    return ret;
> +}
> +
> +static int convert_coeffs(AVFilterContext *ctx)
> +{
> +    FIRContext *s = ctx->priv;
> +    int ch, n;
> +
> +    s->nb_taps = av_audio_fifo_size(s->fifo[1]);
> +    if (s->nb_taps > 131072) {
> +        av_log(ctx, AV_LOG_ERROR, "Too big number of taps: %d > 131072.\n", s->nb_taps);
> +        return AVERROR(EINVAL);
> +    }
> +
> +    for (n = 1; (1 << n) < s->nb_taps; n++);
> +    s->n = n;
> +    s->fft_length = 1 << s->n;
> +
> +    for (ch = 0; ch < ctx->inputs[0]->channels; ch++) {
> +        s->fft_data[ch] = av_calloc(s->fft_length, sizeof(**s->fft_data));
> +        if (!s->fft_data[ch])
> +            return AVERROR(ENOMEM);
> +    }
> +
> +    for (ch = 0; ch < ctx->inputs[1]->channels; ch++) {
> +        s->fft_coef[ch] = av_calloc(s->fft_length, sizeof(**s->fft_coef));
> +        if (!s->fft_coef[ch])
> +            return AVERROR(ENOMEM);
> +    }
> +
> +    s->hop_size = s->nb_taps / 4;

In theory, hop_size should be <= fft_length - nb_taps + 1



> +    if (s->hop_size <= 0) {
> +        av_log(ctx, AV_LOG_ERROR, "Too small number of taps: %d < 4.\n", s->nb_taps);
> +        return AVERROR(EINVAL);
> +    }
> +
> +    s->buffer = ff_get_audio_buffer(ctx->inputs[0], s->fft_length * 2);
> +    if (!s->buffer)
> +        return AVERROR(ENOMEM);
> +
> +    s->fft  = av_fft_init(s->n, 0);
> +    s->ifft = av_fft_init(s->n, 1);
> +    if (!s->fft || !s->ifft)
> +        return AVERROR(ENOMEM);
> +
> +    s->in[1] = ff_get_audio_buffer(ctx->inputs[1], s->nb_taps);
> +    if (!s->in[1])
> +        return AVERROR(ENOMEM);
> +
> +    av_audio_fifo_read(s->fifo[1], (void **)s->in[1]->data, s->nb_taps);
> +    for (ch = 0; ch < ctx->inputs[1]->channels; ch++) {
> +        FFTComplex *fft_coef = s->fft_coef[ch];
> +        const float *re = (const float *)s->in[1]->extended_data[0 + !s->one2many * ch];
> +        const float *im = (const float *)s->in[1]->extended_data[1 + !s->one2many * ch];

What is the meaning of imaginary coeffs in time domain?



> +        const float scale = 1.f / s->fft_length;
> +        const int offset = (s->fft_length - s->nb_taps) / 2;
> +
> +        memset(fft_coef, 0, sizeof(*fft_coef) * s->fft_length);
> +        for (n = 0; n < s->nb_taps; n++) {
> +            fft_coef[n + offset].re = re[n] * scale;
> +            fft_coef[n + offset].im = im[n] * scale;
> +        }
> +        av_fft_permute(s->fft, fft_coef);
> +        av_fft_calc(s->fft, fft_coef);
> +    }
> +
> +    av_frame_free(&s->in[1]);
> +    s->have_coeffs = 1;
> +
> +    return 0;
> +}
> +
> +static int read_coeffs(AVFilterLink *link, AVFrame *frame)
> +{
> +    AVFilterContext *ctx = link->dst;
> +    FIRContext *s = ctx->priv;
> +
> +    av_audio_fifo_write(s->fifo[1], (void **)frame->extended_data,
> +                        frame->nb_samples);
> +    av_frame_free(&frame);
> +
> +    return 0;
> +}
> +
> +static int filter_frame(AVFilterLink *link, AVFrame *frame)
> +{
> +    AVFilterContext *ctx = link->dst;
> +    FIRContext *s = ctx->priv;
> +    AVFilterLink *outlink = ctx->outputs[0];
> +    int ret = 0;
> +
> +    av_audio_fifo_write(s->fifo[0], (void **)frame->extended_data,
> +                        frame->nb_samples);
> +    av_frame_free(&frame);
> +
> +    if (!s->have_coeffs && s->eof_coeffs) {
> +        ret = convert_coeffs(ctx);
> +        if (ret < 0)
> +            return ret;
> +    }
> +
> +    if (s->have_coeffs) {
> +        while (av_audio_fifo_size(s->fifo[0]) >= s->fft_length) {
> +            ret = fir_filter(s, outlink);
> +            if (ret < 0)
> +                break;
> +        }
> +    }
> +    return ret;
> +}
> +
> +static int request_frame(AVFilterLink *outlink)
> +{
> +    AVFilterContext *ctx = outlink->src;
> +    FIRContext *s = ctx->priv;
> +    int ret;
> +
> +    if (!s->eof_coeffs) {
> +        ret = ff_request_frame(ctx->inputs[1]);
> +        if (ret == AVERROR_EOF) {
> +            s->eof_coeffs = 1;
> +            ret = 0;
> +        }
> +        return ret;
> +    }
> +    ret = ff_request_frame(ctx->inputs[0]);
> +    if (ret == AVERROR_EOF) {
> +        while (av_audio_fifo_size(s->fifo[0]) > 0) {
> +            ret = fir_filter(s, outlink);
> +            if (ret < 0)
> +                return ret;
> +        }
> +        ret = AVERROR_EOF;
> +    }
> +    return ret;
> +}
> +
> +static int query_formats(AVFilterContext *ctx)
> +{
> +    AVFilterFormats *formats;
> +    AVFilterChannelLayouts *layouts = NULL;
> +    static const enum AVSampleFormat sample_fmts[] = {
> +        AV_SAMPLE_FMT_FLTP,
> +        AV_SAMPLE_FMT_NONE
> +    };
> +    int ret, i;
> +
> +    layouts = ff_all_channel_counts();
> +    if ((ret = ff_channel_layouts_ref(layouts, &ctx->outputs[0]->in_channel_layouts)) < 0)
> +        return ret;
> +
> +    for (i = 0; i < 2; i++) {
> +        layouts = ff_all_channel_counts();
> +        if ((ret = ff_channel_layouts_ref(layouts, &ctx->inputs[i]->out_channel_layouts)) < 0)
> +            return ret;
> +    }
> +
> +    formats = ff_make_format_list(sample_fmts);
> +    if ((ret = ff_set_common_formats(ctx, formats)) < 0)
> +        return ret;
> +
> +    formats = ff_all_samplerates();
> +    return ff_set_common_samplerates(ctx, formats);
> +}
> +
> +static int config_output(AVFilterLink *outlink)
> +{
> +    AVFilterContext *ctx = outlink->src;
> +    FIRContext *s = ctx->priv;
> +
> +    if (ctx->inputs[0]->channels * 2 != ctx->inputs[1]->channels &&
> +        ctx->inputs[1]->channels != 2) {
> +        av_log(ctx, AV_LOG_ERROR,
> +               "Second input must have double number of channels as first input or "
> +               "exactly 2 channels.\n");
> +        return AVERROR(EINVAL);
> +    }
> +
> +    s->one2many = ctx->inputs[1]->channels == 2;
> +    outlink->sample_rate = ctx->inputs[0]->sample_rate;
> +    outlink->time_base   = ctx->inputs[0]->time_base;
> +    outlink->channel_layout = ctx->inputs[0]->channel_layout;
> +    outlink->channels = ctx->inputs[0]->channels;
> +
> +    s->fifo[0] = av_audio_fifo_alloc(ctx->inputs[0]->format, ctx->inputs[0]->channels, 1024);
> +    s->fifo[1] = av_audio_fifo_alloc(ctx->inputs[1]->format, ctx->inputs[1]->channels, 1024);
> +    if (!s->fifo[0] || !s->fifo[1])
> +        return AVERROR(ENOMEM);
> +
> +    s->fft_data = av_calloc(outlink->channels, sizeof(*s->fft_data));
> +    s->fft_coef = av_calloc(ctx->inputs[1]->channels, sizeof(*s->fft_coef));
> +    if (!s->fft_data || !s->fft_coef)
> +        return AVERROR(ENOMEM);
> +    s->nb_channels = outlink->channels;
> +
> +    return 0;
> +}
> +
> +static av_cold void uninit(AVFilterContext *ctx)
> +{
> +    FIRContext *s = ctx->priv;
> +    int ch;
> +
> +    for (ch = 0; ch < s->nb_channels; ch++) {
> +        if (s->fft_data)
> +            av_freep(&s->fft_data[ch]);
> +    }
> +    av_freep(&s->fft_data);
> +
> +    for (ch = 0; ch < s->nb_channels; ch++) {
> +        if (s->fft_coef)
> +            av_freep(&s->fft_coef[ch]);
> +    }
> +    av_freep(&s->fft_coef);
> +
> +    av_fft_end(s->fft);
> +    av_fft_end(s->ifft);
> +
> +    av_frame_free(&s->in[0]);
> +    av_frame_free(&s->in[1]);
> +
> +    av_audio_fifo_free(s->fifo[0]);
> +    av_audio_fifo_free(s->fifo[1]);
> +}
> +
> +static const AVFilterPad afirfilter_inputs[] = {
> +    {
> +        .name           = "main",
> +        .type           = AVMEDIA_TYPE_AUDIO,
> +        .filter_frame   = filter_frame,
> +    },{
> +        .name           = "coefficients",
> +        .type           = AVMEDIA_TYPE_AUDIO,
> +        .filter_frame   = read_coeffs,
> +    },
> +    { NULL }
> +};
> +
> +static const AVFilterPad afirfilter_outputs[] = {
> +    {
> +        .name          = "default",
> +        .type          = AVMEDIA_TYPE_AUDIO,
> +        .config_props  = config_output,
> +        .request_frame = request_frame,
> +    },
> +    { NULL }
> +};
> +
> +AVFilter ff_af_afirfilter = {
> +    .name           = "afirfilter",
> +    .description    = NULL_IF_CONFIG_SMALL("Apply Finite Impulse Response filter with supplied coefficients in 2nd stream."),
> +    .priv_size      = sizeof(FIRContext),
> +    .query_formats  = query_formats,
> +    .uninit         = uninit,
> +    .inputs         = afirfilter_inputs,
> +    .outputs        = afirfilter_outputs,
> +};
> diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
> index 8fb87eb..8bfe1ae 100644
> --- a/libavfilter/allfilters.c
> +++ b/libavfilter/allfilters.c
> @@ -50,6 +50,7 @@ static void register_all(void)
>      REGISTER_FILTER(AEVAL,          aeval,          af);
>      REGISTER_FILTER(AFADE,          afade,          af);
>      REGISTER_FILTER(AFFTFILT,       afftfilt,       af);
> +    REGISTER_FILTER(AFIRFILTER,     afirfilter,     af);
>      REGISTER_FILTER(AFORMAT,        aformat,        af);
>      REGISTER_FILTER(AGATE,          agate,          af);
>      REGISTER_FILTER(AINTERLEAVE,    ainterleave,    af);
> --
> 2.9.3
>
> _______________________________________________
> ffmpeg-devel mailing list
> ffmpeg-devel@ffmpeg.org
> http://ffmpeg.org/mailman/listinfo/ffmpeg-devel

On 5/1/17, Muhammad Faiz <mfcc64@gmail.com> wrote:
> On Mon, May 1, 2017 at 5:02 AM, Paul B Mahol <onemda@gmail.com> wrote:

[...]

>> +
>> +        for (ch = 0; ch < s->nb_channels; ch++) {
>> +            dst = (float *)out->extended_data[ch];
>> +            buf = (float *)s->buffer->extended_data[ch];
>> +
>> +            for (n = 0; n < s->fft_length; n++)
>> +                dst[n] = buf[n];
>> +            memmove(buf, buf + s->fft_length, s->fft_length * 4);
>> +        }
>
> Is this overlap-save?
>

Yes.

>
>
>> +
>> +        ret = ff_filter_frame(outlink, out);
>> +    }
>> +
>> +    av_audio_fifo_drain(s->fifo[0], s->hop_size);
>> +    av_frame_free(&s->in[0]);
>> +
>> +    return ret;
>> +}
>> +
>> +static int convert_coeffs(AVFilterContext *ctx)
>> +{
>> +    FIRContext *s = ctx->priv;
>> +    int ch, n;
>> +
>> +    s->nb_taps = av_audio_fifo_size(s->fifo[1]);
>> +    if (s->nb_taps > 131072) {
>> +        av_log(ctx, AV_LOG_ERROR, "Too big number of taps: %d >
>> 131072.\n", s->nb_taps);
>> +        return AVERROR(EINVAL);
>> +    }
>> +
>> +    for (n = 1; (1 << n) < s->nb_taps; n++);
>> +    s->n = n;
>> +    s->fft_length = 1 << s->n;
>> +
>> +    for (ch = 0; ch < ctx->inputs[0]->channels; ch++) {
>> +        s->fft_data[ch] = av_calloc(s->fft_length,
>> sizeof(**s->fft_data));
>> +        if (!s->fft_data[ch])
>> +            return AVERROR(ENOMEM);
>> +    }
>> +
>> +    for (ch = 0; ch < ctx->inputs[1]->channels; ch++) {
>> +        s->fft_coef[ch] = av_calloc(s->fft_length,
>> sizeof(**s->fft_coef));
>> +        if (!s->fft_coef[ch])
>> +            return AVERROR(ENOMEM);
>> +    }
>> +
>> +    s->hop_size = s->nb_taps / 4;
>
> In theory, hop_size should be <= fft_length - nb_taps + 1

Fixed.

>
>
>
>> +    if (s->hop_size <= 0) {
>> +        av_log(ctx, AV_LOG_ERROR, "Too small number of taps: %d < 4.\n",
>> s->nb_taps);
>> +        return AVERROR(EINVAL);
>> +    }
>> +
>> +    s->buffer = ff_get_audio_buffer(ctx->inputs[0], s->fft_length * 2);
>> +    if (!s->buffer)
>> +        return AVERROR(ENOMEM);
>> +
>> +    s->fft  = av_fft_init(s->n, 0);
>> +    s->ifft = av_fft_init(s->n, 1);
>> +    if (!s->fft || !s->ifft)
>> +        return AVERROR(ENOMEM);
>> +
>> +    s->in[1] = ff_get_audio_buffer(ctx->inputs[1], s->nb_taps);
>> +    if (!s->in[1])
>> +        return AVERROR(ENOMEM);
>> +
>> +    av_audio_fifo_read(s->fifo[1], (void **)s->in[1]->data, s->nb_taps);
>> +    for (ch = 0; ch < ctx->inputs[1]->channels; ch++) {
>> +        FFTComplex *fft_coef = s->fft_coef[ch];
>> +        const float *re = (const float *)s->in[1]->extended_data[0 +
>> !s->one2many * ch];
>> +        const float *im = (const float *)s->in[1]->extended_data[1 +
>> !s->one2many * ch];
>
> What is the meaning of imaginary coeffs in time domain?

Removed.


See new patch.