git.sesse.net Git - ffmpeg/blob - libavfilter/af_acrossover.c

   1 /*
   2  * This file is part of FFmpeg.
   3  *
   4  * FFmpeg is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Lesser General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2.1 of the License, or (at your option) any later version.
   8  *
   9  * FFmpeg is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Lesser General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Lesser General Public
  15  * License along with FFmpeg; if not, write to the Free Software
  16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  17  */
  18
  19 /**
  20  * @file
  21  * Crossover filter
  22  *
  23  * Split an audio stream into several bands.
  24  */
  25
  26 #include "libavutil/attributes.h"
  27 #include "libavutil/avstring.h"
  28 #include "libavutil/channel_layout.h"
  29 #include "libavutil/eval.h"
  30 #include "libavutil/float_dsp.h"
  31 #include "libavutil/internal.h"
  32 #include "libavutil/opt.h"
  33
  34 #include "audio.h"
  35 #include "avfilter.h"
  36 #include "formats.h"
  37 #include "internal.h"
  38
  39 #define MAX_SPLITS 16
  40 #define MAX_BANDS MAX_SPLITS + 1
  41
  42 #define B0 0
  43 #define B1 1
  44 #define B2 2
  45 #define A1 3
  46 #define A2 4
  47
  48 typedef struct BiquadCoeffs {
  49     double cd[5];
  50     float cf[5];
  51 } BiquadCoeffs;
  52
  53 typedef struct AudioCrossoverContext {
  54     const AVClass *class;
  55
  56     char *splits_str;
  57     int order_opt;
  58     float level_in;
  59
  60     int order;
  61     int filter_count;
  62     int first_order;
  63     int ap_filter_count;
  64     int nb_splits;
  65     float splits[MAX_SPLITS];
  66
  67     BiquadCoeffs lp[MAX_BANDS][20];
  68     BiquadCoeffs hp[MAX_BANDS][20];
  69     BiquadCoeffs ap[MAX_BANDS][20];
  70
  71     AVFrame *xover;
  72
  73     AVFrame *input_frame;
  74     AVFrame *frames[MAX_BANDS];
  75
  76     int (*filter_channels)(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs);
  77
  78     AVFloatDSPContext *fdsp;
  79 } AudioCrossoverContext;
  80
  81 #define OFFSET(x) offsetof(AudioCrossoverContext, x)
  82 #define AF AV_OPT_FLAG_AUDIO_PARAM | AV_OPT_FLAG_FILTERING_PARAM
  83
  84 static const AVOption acrossover_options[] = {
  85     { "split", "set split frequencies", OFFSET(splits_str), AV_OPT_TYPE_STRING, {.str="500"}, 0, 0, AF },
  86     { "order", "set filter order",      OFFSET(order_opt),  AV_OPT_TYPE_INT,    {.i64=1},     0, 9, AF, "m" },
  87     { "2nd",   "2nd order (12 dB/8ve)", 0,                  AV_OPT_TYPE_CONST,  {.i64=0},     0, 0, AF, "m" },
  88     { "4th",   "4th order (24 dB/8ve)", 0,                  AV_OPT_TYPE_CONST,  {.i64=1},     0, 0, AF, "m" },
  89     { "6th",   "6th order (36 dB/8ve)", 0,                  AV_OPT_TYPE_CONST,  {.i64=2},     0, 0, AF, "m" },
  90     { "8th",   "8th order (48 dB/8ve)", 0,                  AV_OPT_TYPE_CONST,  {.i64=3},     0, 0, AF, "m" },
  91     { "10th",  "10th order (60 dB/8ve)",0,                  AV_OPT_TYPE_CONST,  {.i64=4},     0, 0, AF, "m" },
  92     { "12th",  "12th order (72 dB/8ve)",0,                  AV_OPT_TYPE_CONST,  {.i64=5},     0, 0, AF, "m" },
  93     { "14th",  "14th order (84 dB/8ve)",0,                  AV_OPT_TYPE_CONST,  {.i64=6},     0, 0, AF, "m" },
  94     { "16th",  "16th order (96 dB/8ve)",0,                  AV_OPT_TYPE_CONST,  {.i64=7},     0, 0, AF, "m" },
  95     { "18th",  "18th order (108 dB/8ve)",0,                 AV_OPT_TYPE_CONST,  {.i64=8},     0, 0, AF, "m" },
  96     { "20th",  "20th order (120 dB/8ve)",0,                 AV_OPT_TYPE_CONST,  {.i64=9},     0, 0, AF, "m" },
  97     { "level", "set input gain",        OFFSET(level_in),   AV_OPT_TYPE_FLOAT,  {.dbl=1},     0, 1, AF },
  98     { NULL }
  99 };
 100
 101 AVFILTER_DEFINE_CLASS(acrossover);
 102
 103 static av_cold int init(AVFilterContext *ctx)
 104 {
 105     AudioCrossoverContext *s = ctx->priv;
 106     char *p, *arg, *saveptr = NULL;
 107     int i, ret = 0;
 108
 109     s->fdsp = avpriv_float_dsp_alloc(0);
 110     if (!s->fdsp)
 111         return AVERROR(ENOMEM);
 112
 113     p = s->splits_str;
 114     for (i = 0; i < MAX_SPLITS; i++) {
 115         float freq;
 116
 117         if (!(arg = av_strtok(p, " |", &saveptr)))
 118             break;
 119
 120         p = NULL;
 121
 122         if (av_sscanf(arg, "%f", &freq) != 1) {
 123             av_log(ctx, AV_LOG_ERROR, "Invalid syntax for frequency[%d].\n", i);
 124             return AVERROR(EINVAL);
 125         }
 126         if (freq <= 0) {
 127             av_log(ctx, AV_LOG_ERROR, "Frequency %f must be positive number.\n", freq);
 128             return AVERROR(EINVAL);
 129         }
 130
 131         if (i > 0 && freq <= s->splits[i-1]) {
 132             av_log(ctx, AV_LOG_ERROR, "Frequency %f must be in increasing order.\n", freq);
 133             return AVERROR(EINVAL);
 134         }
 135
 136         s->splits[i] = freq;
 137     }
 138
 139     s->nb_splits = i;
 140
 141     for (i = 0; i <= s->nb_splits; i++) {
 142         AVFilterPad pad  = { 0 };
 143         char *name;
 144
 145         pad.type = AVMEDIA_TYPE_AUDIO;
 146         name = av_asprintf("out%d", ctx->nb_outputs);
 147         if (!name)
 148             return AVERROR(ENOMEM);
 149         pad.name = name;
 150
 151         if ((ret = ff_insert_outpad(ctx, i, &pad)) < 0) {
 152             av_freep(&pad.name);
 153             return ret;
 154         }
 155     }
 156
 157     return ret;
 158 }
 159
 160 static void set_lp(BiquadCoeffs *b, double fc, double q, double sr)
 161 {
 162     double omega = 2. * M_PI * fc / sr;
 163     double cosine = cos(omega);
 164     double alpha = sin(omega) / (2. * q);
 165
 166     double b0 = (1. - cosine) / 2.;
 167     double b1 = 1. - cosine;
 168     double b2 = (1. - cosine) / 2.;
 169     double a0 = 1. + alpha;
 170     double a1 = -2. * cosine;
 171     double a2 = 1. - alpha;
 172
 173     b->cd[B0] =  b0 / a0;
 174     b->cd[B1] =  b1 / a0;
 175     b->cd[B2] =  b2 / a0;
 176     b->cd[A1] = -a1 / a0;
 177     b->cd[A2] = -a2 / a0;
 178
 179     b->cf[B0] = b->cd[B0];
 180     b->cf[B1] = b->cd[B1];
 181     b->cf[B2] = b->cd[B2];
 182     b->cf[A1] = b->cd[A1];
 183     b->cf[A2] = b->cd[A2];
 184 }
 185
 186 static void set_hp(BiquadCoeffs *b, double fc, double q, double sr)
 187 {
 188     double omega = 2. * M_PI * fc / sr;
 189     double cosine = cos(omega);
 190     double alpha = sin(omega) / (2. * q);
 191
 192     double b0 = (1. + cosine) / 2.;
 193     double b1 = -1. - cosine;
 194     double b2 = (1. + cosine) / 2.;
 195     double a0 = 1. + alpha;
 196     double a1 = -2. * cosine;
 197     double a2 = 1. - alpha;
 198
 199     b->cd[B0] =  b0 / a0;
 200     b->cd[B1] =  b1 / a0;
 201     b->cd[B2] =  b2 / a0;
 202     b->cd[A1] = -a1 / a0;
 203     b->cd[A2] = -a2 / a0;
 204
 205     b->cf[B0] = b->cd[B0];
 206     b->cf[B1] = b->cd[B1];
 207     b->cf[B2] = b->cd[B2];
 208     b->cf[A1] = b->cd[A1];
 209     b->cf[A2] = b->cd[A2];
 210 }
 211
 212 static void set_ap(BiquadCoeffs *b, double fc, double q, double sr)
 213 {
 214     double omega = 2. * M_PI * fc / sr;
 215     double cosine = cos(omega);
 216     double alpha = sin(omega) / (2. * q);
 217
 218     double a0 = 1. + alpha;
 219     double a1 = -2. * cosine;
 220     double a2 = 1. - alpha;
 221     double b0 = a2;
 222     double b1 = a1;
 223     double b2 = a0;
 224
 225     b->cd[B0] =  b0 / a0;
 226     b->cd[B1] =  b1 / a0;
 227     b->cd[B2] =  b2 / a0;
 228     b->cd[A1] = -a1 / a0;
 229     b->cd[A2] = -a2 / a0;
 230
 231     b->cf[B0] = b->cd[B0];
 232     b->cf[B1] = b->cd[B1];
 233     b->cf[B2] = b->cd[B2];
 234     b->cf[A1] = b->cd[A1];
 235     b->cf[A2] = b->cd[A2];
 236 }
 237
 238 static void set_ap1(BiquadCoeffs *b, double fc, double sr)
 239 {
 240     double omega = 2. * M_PI * fc / sr;
 241
 242     b->cd[A1] = exp(-omega);
 243     b->cd[A2] = 0.;
 244     b->cd[B0] = -b->cd[A1];
 245     b->cd[B1] = 1.;
 246     b->cd[B2] = 0.;
 247
 248     b->cf[B0] = b->cd[B0];
 249     b->cf[B1] = b->cd[B1];
 250     b->cf[B2] = b->cd[B2];
 251     b->cf[A1] = b->cd[A1];
 252     b->cf[A2] = b->cd[A2];
 253 }
 254
 255 static void calc_q_factors(int order, double *q)
 256 {
 257     double n = order / 2.;
 258
 259     for (int i = 0; i < n / 2; i++)
 260         q[i] = 1. / (-2. * cos(M_PI * (2. * (i + 1) + n - 1.) / (2. * n)));
 261 }
 262
 263 static int query_formats(AVFilterContext *ctx)
 264 {
 265     AVFilterFormats *formats;
 266     AVFilterChannelLayouts *layouts;
 267     static const enum AVSampleFormat sample_fmts[] = {
 268         AV_SAMPLE_FMT_FLTP, AV_SAMPLE_FMT_DBLP,
 269         AV_SAMPLE_FMT_NONE
 270     };
 271     int ret;
 272
 273     layouts = ff_all_channel_counts();
 274     if (!layouts)
 275         return AVERROR(ENOMEM);
 276     ret = ff_set_common_channel_layouts(ctx, layouts);
 277     if (ret < 0)
 278         return ret;
 279
 280     formats = ff_make_format_list(sample_fmts);
 281     if (!formats)
 282         return AVERROR(ENOMEM);
 283     ret = ff_set_common_formats(ctx, formats);
 284     if (ret < 0)
 285         return ret;
 286
 287     formats = ff_all_samplerates();
 288     if (!formats)
 289         return AVERROR(ENOMEM);
 290     return ff_set_common_samplerates(ctx, formats);
 291 }
 292
 293 #define BIQUAD_PROCESS(name, type)                             \
 294 static void biquad_process_## name(const type *const c,        \
 295                                    type *b,                    \
 296                                    type *dst, const type *src, \
 297                                    int nb_samples)             \
 298 {                                                              \
 299     const type b0 = c[B0];                                     \
 300     const type b1 = c[B1];                                     \
 301     const type b2 = c[B2];                                     \
 302     const type a1 = c[A1];                                     \
 303     const type a2 = c[A2];                                     \
 304     type z1 = b[0];                                            \
 305     type z2 = b[1];                                            \
 306                                                                \
 307     for (int n = 0; n + 1 < nb_samples; n++) {                 \
 308         type in = src[n];                                      \
 309         type out;                                              \
 310                                                                \
 311         out = in * b0 + z1;                                    \
 312         z1 = b1 * in + z2 + a1 * out;                          \
 313         z2 = b2 * in + a2 * out;                               \
 314         dst[n] = out;                                          \
 315                                                                \
 316         n++;                                                   \
 317         in = src[n];                                           \
 318         out = in * b0 + z1;                                    \
 319         z1 = b1 * in + z2 + a1 * out;                          \
 320         z2 = b2 * in + a2 * out;                               \
 321         dst[n] = out;                                          \
 322     }                                                          \
 323                                                                \
 324     if (nb_samples & 1) {                                      \
 325         const int n = nb_samples - 1;                          \
 326         const type in = src[n];                                \
 327         type out;                                              \
 328                                                                \
 329         out = in * b0 + z1;                                    \
 330         z1 = b1 * in + z2 + a1 * out;                          \
 331         z2 = b2 * in + a2 * out;                               \
 332         dst[n] = out;                                          \
 333     }                                                          \
 334                                                                \
 335     b[0] = z1;                                                 \
 336     b[1] = z2;                                                 \
 337 }
 338
 339 BIQUAD_PROCESS(fltp, float)
 340 BIQUAD_PROCESS(dblp, double)
 341
 342 #define XOVER_PROCESS(name, type, one, ff)                                                  \
 343 static int filter_channels_## name(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs) \
 344 {                                                                                           \
 345     AudioCrossoverContext *s = ctx->priv;                                                   \
 346     AVFrame *in = s->input_frame;                                                           \
 347     AVFrame **frames = s->frames;                                                           \
 348     const int start = (in->channels * jobnr) / nb_jobs;                                     \
 349     const int end = (in->channels * (jobnr+1)) / nb_jobs;                                   \
 350     const int nb_samples = in->nb_samples;                                                  \
 351     const int nb_outs = ctx->nb_outputs;                                                    \
 352                                                                                             \
 353     for (int ch = start; ch < end; ch++) {                                                  \
 354         const type *src = (const type *)in->extended_data[ch];                              \
 355         type *xover = (type *)s->xover->extended_data[ch];                                  \
 356                                                                                             \
 357         s->fdsp->vector_## ff ##mul_scalar((type *)frames[0]->extended_data[ch], src,       \
 358                                     s->level_in, FFALIGN(nb_samples, sizeof(type)));        \
 359                                                                                             \
 360         for (int band = 0; band < nb_outs; band++) {                                        \
 361             for (int f = 0; band + 1 < nb_outs && f < s->filter_count; f++) {               \
 362                 const type *prv = (const type *)frames[band]->extended_data[ch];            \
 363                 type *dst = (type *)frames[band + 1]->extended_data[ch];                    \
 364                 const type *hsrc = f == 0 ? prv : dst;                                      \
 365                 type *hp = xover + nb_outs * 20 + band * 20 + f * 2;                        \
 366                 const type *const hpc = (type *)&s->hp[band][f].c ## ff;                    \
 367                                                                                             \
 368                 biquad_process_## name(hpc, hp, dst, hsrc, nb_samples);                     \
 369             }                                                                               \
 370                                                                                             \
 371             for (int f = 0; band + 1 < nb_outs && f < s->filter_count; f++) {               \
 372                 type *dst = (type *)frames[band]->extended_data[ch];                        \
 373                 const type *lsrc = dst;                                                     \
 374                 type *lp = xover + band * 20 + f * 2;                                       \
 375                 const type *const lpc = (type *)&s->lp[band][f].c ## ff;                    \
 376                                                                                             \
 377                 biquad_process_## name(lpc, lp, dst, lsrc, nb_samples);                     \
 378             }                                                                               \
 379                                                                                             \
 380             for (int aband = band + 1; aband + 1 < nb_outs; aband++) {                      \
 381                 if (s->first_order) {                                                       \
 382                     const type *asrc = (const type *)frames[band]->extended_data[ch];       \
 383                     type *dst = (type *)frames[band]->extended_data[ch];                    \
 384                     type *ap = xover + nb_outs * 40 + (aband * nb_outs + band) * 20;        \
 385                     const type *const apc = (type *)&s->ap[aband][0].c ## ff;               \
 386                                                                                             \
 387                     biquad_process_## name(apc, ap, dst, asrc, nb_samples);                 \
 388                 }                                                                           \
 389                                                                                             \
 390                 for (int f = s->first_order; f < s->ap_filter_count; f++) {                 \
 391                     const type *asrc = (const type *)frames[band]->extended_data[ch];       \
 392                     type *dst = (type *)frames[band]->extended_data[ch];                    \
 393                     type *ap = xover + nb_outs * 40 + (aband * nb_outs + band) * 20 + f * 2;\
 394                     const type *const apc = (type *)&s->ap[aband][f].c ## ff;               \
 395                                                                                             \
 396                     biquad_process_## name(apc, ap, dst, asrc, nb_samples);                 \
 397                 }                                                                           \
 398             }                                                                               \
 399         }                                                                                   \
 400                                                                                             \
 401         for (int band = 0; band < nb_outs && s->first_order; band++) {                      \
 402             if (band & 1) {                                                                 \
 403                 type *dst = (type *)frames[band]->extended_data[ch];                        \
 404                 s->fdsp->vector_## ff ##mul_scalar(dst, dst, -one,                          \
 405                                                    FFALIGN(nb_samples, sizeof(type)));      \
 406             }                                                                               \
 407         }                                                                                   \
 408     }                                                                                       \
 409                                                                                             \
 410     return 0;                                                                               \
 411 }
 412
 413 XOVER_PROCESS(fltp, float, 1.f, f)
 414 XOVER_PROCESS(dblp, double, 1.0, d)
 415
 416 static int config_input(AVFilterLink *inlink)
 417 {
 418     AVFilterContext *ctx = inlink->dst;
 419     AudioCrossoverContext *s = ctx->priv;
 420     int sample_rate = inlink->sample_rate;
 421     double q[16];
 422
 423     s->order = (s->order_opt + 1) * 2;
 424     s->filter_count = s->order / 2;
 425     s->first_order = s->filter_count & 1;
 426     s->ap_filter_count = s->filter_count / 2 + s->first_order;
 427     calc_q_factors(s->order, q);
 428
 429     for (int band = 0; band <= s->nb_splits; band++) {
 430         if (s->first_order) {
 431             set_lp(&s->lp[band][0], s->splits[band], 0.5, sample_rate);
 432             set_hp(&s->hp[band][0], s->splits[band], 0.5, sample_rate);
 433         }
 434
 435         for (int n = s->first_order; n < s->filter_count; n++) {
 436             const int idx = s->filter_count / 2 - ((n + s->first_order) / 2 - s->first_order) - 1;
 437
 438             set_lp(&s->lp[band][n], s->splits[band], q[idx], sample_rate);
 439             set_hp(&s->hp[band][n], s->splits[band], q[idx], sample_rate);
 440         }
 441
 442         if (s->first_order)
 443             set_ap1(&s->ap[band][0], s->splits[band], sample_rate);
 444
 445         for (int n = s->first_order; n < s->ap_filter_count; n++) {
 446             const int idx = (s->filter_count / 2 - ((n * 2 + s->first_order) / 2 - s->first_order) - 1);
 447
 448             set_ap(&s->ap[band][n], s->splits[band], q[idx], sample_rate);
 449         }
 450     }
 451
 452     switch (inlink->format) {
 453     case AV_SAMPLE_FMT_FLTP: s->filter_channels = filter_channels_fltp; break;
 454     case AV_SAMPLE_FMT_DBLP: s->filter_channels = filter_channels_dblp; break;
 455     }
 456
 457     s->xover = ff_get_audio_buffer(inlink, 2 * (ctx->nb_outputs * 10 + ctx->nb_outputs * 10 +
 458                                                 ctx->nb_outputs * ctx->nb_outputs * 10));
 459     if (!s->xover)
 460         return AVERROR(ENOMEM);
 461
 462     return 0;
 463 }
 464
 465 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
 466 {
 467     AVFilterContext *ctx = inlink->dst;
 468     AudioCrossoverContext *s = ctx->priv;
 469     AVFrame **frames = s->frames;
 470     int i, ret = 0;
 471
 472     for (i = 0; i < ctx->nb_outputs; i++) {
 473         frames[i] = ff_get_audio_buffer(ctx->outputs[i], in->nb_samples);
 474
 475         if (!frames[i]) {
 476             ret = AVERROR(ENOMEM);
 477             break;
 478         }
 479
 480         frames[i]->pts = in->pts;
 481     }
 482
 483     if (ret < 0)
 484         goto fail;
 485
 486     s->input_frame = in;
 487     ctx->internal->execute(ctx, s->filter_channels, NULL, NULL, FFMIN(inlink->channels,
 488                                                                       ff_filter_get_nb_threads(ctx)));
 489
 490     for (i = 0; i < ctx->nb_outputs; i++) {
 491         ret = ff_filter_frame(ctx->outputs[i], frames[i]);
 492         frames[i] = NULL;
 493         if (ret < 0)
 494             break;
 495     }
 496
 497 fail:
 498     for (i = 0; i < ctx->nb_outputs; i++)
 499         av_frame_free(&frames[i]);
 500     av_frame_free(&in);
 501     s->input_frame = NULL;
 502
 503     return ret;
 504 }
 505
 506 static av_cold void uninit(AVFilterContext *ctx)
 507 {
 508     AudioCrossoverContext *s = ctx->priv;
 509     int i;
 510
 511     av_freep(&s->fdsp);
 512     av_frame_free(&s->xover);
 513
 514     for (i = 0; i < ctx->nb_outputs; i++)
 515         av_freep(&ctx->output_pads[i].name);
 516 }
 517
 518 static const AVFilterPad inputs[] = {
 519     {
 520         .name         = "default",
 521         .type         = AVMEDIA_TYPE_AUDIO,
 522         .filter_frame = filter_frame,
 523         .config_props = config_input,
 524     },
 525     { NULL }
 526 };
 527
 528 AVFilter ff_af_acrossover = {
 529     .name           = "acrossover",
 530     .description    = NULL_IF_CONFIG_SMALL("Split audio into per-bands streams."),
 531     .priv_size      = sizeof(AudioCrossoverContext),
 532     .priv_class     = &acrossover_class,
 533     .init           = init,
 534     .uninit         = uninit,
 535     .query_formats  = query_formats,
 536     .inputs         = inputs,
 537     .outputs        = NULL,
 538     .flags          = AVFILTER_FLAG_DYNAMIC_OUTPUTS |
 539                       AVFILTER_FLAG_SLICE_THREADS,
 540 };