LCOV - code coverage report
Current view: top level - libavfilter - af_aiir.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 0 384 0.0 %
Date: 2018-05-20 11:54:08 Functions: 0 24 0.0 %

          Line data    Source code
       1             : /*
       2             :  * Copyright (c) 2018 Paul B Mahol
       3             :  *
       4             :  * This file is part of FFmpeg.
       5             :  *
       6             :  * FFmpeg is free software; you can redistribute it and/or
       7             :  * modify it under the terms of the GNU Lesser General Public
       8             :  * License as published by the Free Software Foundation; either
       9             :  * version 2.1 of the License, or (at your option) any later version.
      10             :  *
      11             :  * FFmpeg is distributed in the hope that it will be useful,
      12             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      14             :  * Lesser General Public License for more details.
      15             :  *
      16             :  * You should have received a copy of the GNU Lesser General Public
      17             :  * License along with FFmpeg; if not, write to the Free Software
      18             :  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
      19             :  */
      20             : 
      21             : #include <float.h>
      22             : 
      23             : #include "libavutil/avassert.h"
      24             : #include "libavutil/avstring.h"
      25             : #include "libavutil/opt.h"
      26             : #include "audio.h"
      27             : #include "avfilter.h"
      28             : #include "internal.h"
      29             : 
      30             : typedef struct ThreadData {
      31             :     AVFrame *in, *out;
      32             : } ThreadData;
      33             : 
      34             : typedef struct Pair {
      35             :     int a, b;
      36             : } Pair;
      37             : 
      38             : typedef struct BiquadContext {
      39             :     double a0, a1, a2;
      40             :     double b0, b1, b2;
      41             :     double i1, i2;
      42             :     double o1, o2;
      43             : } BiquadContext;
      44             : 
      45             : typedef struct IIRChannel {
      46             :     int nb_ab[2];
      47             :     double *ab[2];
      48             :     double g;
      49             :     double *cache[2];
      50             :     BiquadContext *biquads;
      51             :     int clippings;
      52             : } IIRChannel;
      53             : 
      54             : typedef struct AudioIIRContext {
      55             :     const AVClass *class;
      56             :     char *a_str, *b_str, *g_str;
      57             :     double dry_gain, wet_gain;
      58             :     int format;
      59             :     int process;
      60             :     int precision;
      61             : 
      62             :     IIRChannel *iir;
      63             :     int channels;
      64             :     enum AVSampleFormat sample_format;
      65             : 
      66             :     int (*iir_channel)(AVFilterContext *ctx, void *arg, int ch, int nb_jobs);
      67             : } AudioIIRContext;
      68             : 
      69           0 : static int query_formats(AVFilterContext *ctx)
      70             : {
      71           0 :     AudioIIRContext *s = ctx->priv;
      72             :     AVFilterFormats *formats;
      73             :     AVFilterChannelLayouts *layouts;
      74           0 :     enum AVSampleFormat sample_fmts[] = {
      75             :         AV_SAMPLE_FMT_DBLP,
      76             :         AV_SAMPLE_FMT_NONE
      77             :     };
      78             :     int ret;
      79             : 
      80           0 :     layouts = ff_all_channel_counts();
      81           0 :     if (!layouts)
      82           0 :         return AVERROR(ENOMEM);
      83           0 :     ret = ff_set_common_channel_layouts(ctx, layouts);
      84           0 :     if (ret < 0)
      85           0 :         return ret;
      86             : 
      87           0 :     sample_fmts[0] = s->sample_format;
      88           0 :     formats = ff_make_format_list(sample_fmts);
      89           0 :     if (!formats)
      90           0 :         return AVERROR(ENOMEM);
      91           0 :     ret = ff_set_common_formats(ctx, formats);
      92           0 :     if (ret < 0)
      93           0 :         return ret;
      94             : 
      95           0 :     formats = ff_all_samplerates();
      96           0 :     if (!formats)
      97           0 :         return AVERROR(ENOMEM);
      98           0 :     return ff_set_common_samplerates(ctx, formats);
      99             : }
     100             : 
     101             : #define IIR_CH(name, type, min, max, need_clipping)                     \
     102             : static int iir_ch_## name(AVFilterContext *ctx, void *arg, int ch, int nb_jobs)  \
     103             : {                                                                       \
     104             :     AudioIIRContext *s = ctx->priv;                                     \
     105             :     const double ig = s->dry_gain;                                      \
     106             :     const double og = s->wet_gain;                                      \
     107             :     ThreadData *td = arg;                                               \
     108             :     AVFrame *in = td->in, *out = td->out;                               \
     109             :     const type *src = (const type *)in->extended_data[ch];              \
     110             :     double *ic = (double *)s->iir[ch].cache[0];                         \
     111             :     double *oc = (double *)s->iir[ch].cache[1];                         \
     112             :     const int nb_a = s->iir[ch].nb_ab[0];                               \
     113             :     const int nb_b = s->iir[ch].nb_ab[1];                               \
     114             :     const double *a = s->iir[ch].ab[0];                                 \
     115             :     const double *b = s->iir[ch].ab[1];                                 \
     116             :     int *clippings = &s->iir[ch].clippings;                             \
     117             :     type *dst = (type *)out->extended_data[ch];                         \
     118             :     int n;                                                              \
     119             :                                                                         \
     120             :     for (n = 0; n < in->nb_samples; n++) {                              \
     121             :         double sample = 0.;                                             \
     122             :         int x;                                                          \
     123             :                                                                         \
     124             :         memmove(&ic[1], &ic[0], (nb_b - 1) * sizeof(*ic));              \
     125             :         memmove(&oc[1], &oc[0], (nb_a - 1) * sizeof(*oc));              \
     126             :         ic[0] = src[n] * ig;                                            \
     127             :         for (x = 0; x < nb_b; x++)                                      \
     128             :             sample += b[x] * ic[x];                                     \
     129             :                                                                         \
     130             :         for (x = 1; x < nb_a; x++)                                      \
     131             :             sample -= a[x] * oc[x];                                     \
     132             :                                                                         \
     133             :         oc[0] = sample;                                                 \
     134             :         sample *= og;                                                   \
     135             :         if (need_clipping && sample < min) {                            \
     136             :             (*clippings)++;                                             \
     137             :             dst[n] = min;                                               \
     138             :         } else if (need_clipping && sample > max) {                     \
     139             :             (*clippings)++;                                             \
     140             :             dst[n] = max;                                               \
     141             :         } else {                                                        \
     142             :             dst[n] = sample;                                            \
     143             :         }                                                               \
     144             :     }                                                                   \
     145             :                                                                         \
     146             :     return 0;                                                           \
     147             : }
     148             : 
     149           0 : IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
     150           0 : IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
     151           0 : IIR_CH(fltp, float,         -1.,        1., 0)
     152           0 : IIR_CH(dblp, double,        -1.,        1., 0)
     153             : 
     154             : #define SERIAL_IIR_CH(name, type, min, max, need_clipping)                  \
     155             : static int iir_ch_serial_## name(AVFilterContext *ctx, void *arg, int ch, int nb_jobs)  \
     156             : {                                                                       \
     157             :     AudioIIRContext *s = ctx->priv;                                     \
     158             :     const double ig = s->dry_gain;                                      \
     159             :     const double og = s->wet_gain;                                      \
     160             :     ThreadData *td = arg;                                               \
     161             :     AVFrame *in = td->in, *out = td->out;                               \
     162             :     const type *src = (const type *)in->extended_data[ch];              \
     163             :     type *dst = (type *)out->extended_data[ch];                         \
     164             :     IIRChannel *iir = &s->iir[ch];                                      \
     165             :     int *clippings = &iir->clippings;                                   \
     166             :     int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2;     \
     167             :     int n, i;                                                           \
     168             :                                                                         \
     169             :     for (i = 0; i < nb_biquads; i++) {                                  \
     170             :         const double a1 = -iir->biquads[i].a1;                          \
     171             :         const double a2 = -iir->biquads[i].a2;                          \
     172             :         const double b0 = iir->biquads[i].b0;                           \
     173             :         const double b1 = iir->biquads[i].b1;                           \
     174             :         const double b2 = iir->biquads[i].b2;                           \
     175             :         double i1 = iir->biquads[i].i1;                                 \
     176             :         double i2 = iir->biquads[i].i2;                                 \
     177             :         double o1 = iir->biquads[i].o1;                                 \
     178             :         double o2 = iir->biquads[i].o2;                                 \
     179             :                                                                         \
     180             :         for (n = 0; n < in->nb_samples; n++) {                          \
     181             :             double sample = ig * (i ? dst[n] : src[n]);                 \
     182             :             double o0 = sample * b0 + i1 * b1 + i2 * b2 + o1 * a1 + o2 * a2; \
     183             :                                                                         \
     184             :             i2 = i1;                                                    \
     185             :             i1 = src[n];                                                \
     186             :             o2 = o1;                                                    \
     187             :             o1 = o0;                                                    \
     188             :             o0 *= og;                                                   \
     189             :                                                                         \
     190             :             if (need_clipping && o0 < min) {                            \
     191             :                 (*clippings)++;                                         \
     192             :                 dst[n] = min;                                           \
     193             :             } else if (need_clipping && o0 > max) {                     \
     194             :                 (*clippings)++;                                         \
     195             :                 dst[n] = max;                                           \
     196             :             } else {                                                    \
     197             :                 dst[n] = o0;                                            \
     198             :             }                                                           \
     199             :         }                                                               \
     200             :         iir->biquads[i].i1 = i1;                                        \
     201             :         iir->biquads[i].i2 = i2;                                        \
     202             :         iir->biquads[i].o1 = o1;                                        \
     203             :         iir->biquads[i].o2 = o2;                                        \
     204             :     }                                                                   \
     205             :                                                                         \
     206             :     return 0;                                                           \
     207             : }
     208             : 
     209           0 : SERIAL_IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
     210           0 : SERIAL_IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
     211           0 : SERIAL_IIR_CH(fltp, float,         -1.,        1., 0)
     212           0 : SERIAL_IIR_CH(dblp, double,        -1.,        1., 0)
     213             : 
     214           0 : static void count_coefficients(char *item_str, int *nb_items)
     215             : {
     216             :     char *p;
     217             : 
     218           0 :     if (!item_str)
     219           0 :         return;
     220             : 
     221           0 :     *nb_items = 1;
     222           0 :     for (p = item_str; *p && *p != '|'; p++) {
     223           0 :         if (*p == ' ')
     224           0 :             (*nb_items)++;
     225             :     }
     226             : }
     227             : 
     228           0 : static int read_gains(AVFilterContext *ctx, char *item_str, int nb_items)
     229             : {
     230           0 :     AudioIIRContext *s = ctx->priv;
     231           0 :     char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
     232             :     int i;
     233             : 
     234           0 :     p = old_str = av_strdup(item_str);
     235           0 :     if (!p)
     236           0 :         return AVERROR(ENOMEM);
     237           0 :     for (i = 0; i < nb_items; i++) {
     238           0 :         if (!(arg = av_strtok(p, "|", &saveptr)))
     239           0 :             arg = prev_arg;
     240             : 
     241           0 :         if (!arg) {
     242           0 :             av_freep(&old_str);
     243           0 :             return AVERROR(EINVAL);
     244             :         }
     245             : 
     246           0 :         p = NULL;
     247           0 :         if (sscanf(arg, "%lf", &s->iir[i].g) != 1) {
     248           0 :             av_log(ctx, AV_LOG_ERROR, "Invalid gains supplied: %s\n", arg);
     249           0 :             av_freep(&old_str);
     250           0 :             return AVERROR(EINVAL);
     251             :         }
     252             : 
     253           0 :         prev_arg = arg;
     254             :     }
     255             : 
     256           0 :     av_freep(&old_str);
     257             : 
     258           0 :     return 0;
     259             : }
     260             : 
     261           0 : static int read_tf_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst)
     262             : {
     263           0 :     char *p, *arg, *old_str, *saveptr = NULL;
     264             :     int i;
     265             : 
     266           0 :     p = old_str = av_strdup(item_str);
     267           0 :     if (!p)
     268           0 :         return AVERROR(ENOMEM);
     269           0 :     for (i = 0; i < nb_items; i++) {
     270           0 :         if (!(arg = av_strtok(p, " ", &saveptr)))
     271           0 :             break;
     272             : 
     273           0 :         p = NULL;
     274           0 :         if (sscanf(arg, "%lf", &dst[i]) != 1) {
     275           0 :             av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
     276           0 :             av_freep(&old_str);
     277           0 :             return AVERROR(EINVAL);
     278             :         }
     279             :     }
     280             : 
     281           0 :     av_freep(&old_str);
     282             : 
     283           0 :     return 0;
     284             : }
     285             : 
     286           0 : static int read_zp_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst, const char *format)
     287             : {
     288           0 :     char *p, *arg, *old_str, *saveptr = NULL;
     289             :     int i;
     290             : 
     291           0 :     p = old_str = av_strdup(item_str);
     292           0 :     if (!p)
     293           0 :         return AVERROR(ENOMEM);
     294           0 :     for (i = 0; i < nb_items; i++) {
     295           0 :         if (!(arg = av_strtok(p, " ", &saveptr)))
     296           0 :             break;
     297             : 
     298           0 :         p = NULL;
     299           0 :         if (sscanf(arg, format, &dst[i*2], &dst[i*2+1]) != 2) {
     300           0 :             av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
     301           0 :             av_freep(&old_str);
     302           0 :             return AVERROR(EINVAL);
     303             :         }
     304             :     }
     305             : 
     306           0 :     av_freep(&old_str);
     307             : 
     308           0 :     return 0;
     309             : }
     310             : 
     311             : static const char *format[] = { "%lf", "%lf %lfi", "%lf %lfr", "%lf %lfd" };
     312             : 
     313           0 : static int read_channels(AVFilterContext *ctx, int channels, uint8_t *item_str, int ab)
     314             : {
     315           0 :     AudioIIRContext *s = ctx->priv;
     316           0 :     char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
     317             :     int i, ret;
     318             : 
     319           0 :     p = old_str = av_strdup(item_str);
     320           0 :     if (!p)
     321           0 :         return AVERROR(ENOMEM);
     322           0 :     for (i = 0; i < channels; i++) {
     323           0 :         IIRChannel *iir = &s->iir[i];
     324             : 
     325           0 :         if (!(arg = av_strtok(p, "|", &saveptr)))
     326           0 :             arg = prev_arg;
     327             : 
     328           0 :         if (!arg) {
     329           0 :             av_freep(&old_str);
     330           0 :             return AVERROR(EINVAL);
     331             :         }
     332             : 
     333           0 :         count_coefficients(arg, &iir->nb_ab[ab]);
     334             : 
     335           0 :         p = NULL;
     336           0 :         iir->cache[ab] = av_calloc(iir->nb_ab[ab] + 1, sizeof(double));
     337           0 :         iir->ab[ab] = av_calloc(iir->nb_ab[ab] * (!!s->format + 1), sizeof(double));
     338           0 :         if (!iir->ab[ab] || !iir->cache[ab]) {
     339           0 :             av_freep(&old_str);
     340           0 :             return AVERROR(ENOMEM);
     341             :         }
     342             : 
     343           0 :         if (s->format) {
     344           0 :             ret = read_zp_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab], format[s->format]);
     345             :         } else {
     346           0 :             ret = read_tf_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab]);
     347             :         }
     348           0 :         if (ret < 0) {
     349           0 :             av_freep(&old_str);
     350           0 :             return ret;
     351             :         }
     352           0 :         prev_arg = arg;
     353             :     }
     354             : 
     355           0 :     av_freep(&old_str);
     356             : 
     357           0 :     return 0;
     358             : }
     359             : 
     360           0 : static void multiply(double wre, double wim, int npz, double *coeffs)
     361             : {
     362           0 :     double nwre = -wre, nwim = -wim;
     363             :     double cre, cim;
     364             :     int i;
     365             : 
     366           0 :     for (i = npz; i >= 1; i--) {
     367           0 :         cre = coeffs[2 * i + 0];
     368           0 :         cim = coeffs[2 * i + 1];
     369             : 
     370           0 :         coeffs[2 * i + 0] = (nwre * cre - nwim * cim) + coeffs[2 * (i - 1) + 0];
     371           0 :         coeffs[2 * i + 1] = (nwre * cim + nwim * cre) + coeffs[2 * (i - 1) + 1];
     372             :     }
     373             : 
     374           0 :     cre = coeffs[0];
     375           0 :     cim = coeffs[1];
     376           0 :     coeffs[0] = nwre * cre - nwim * cim;
     377           0 :     coeffs[1] = nwre * cim + nwim * cre;
     378           0 : }
     379             : 
     380           0 : static int expand(AVFilterContext *ctx, double *pz, int nb, double *coeffs)
     381             : {
     382             :     int i;
     383             : 
     384           0 :     coeffs[0] = 1.0;
     385           0 :     coeffs[1] = 0.0;
     386             : 
     387           0 :     for (i = 0; i < nb; i++) {
     388           0 :         coeffs[2 * (i + 1)    ] = 0.0;
     389           0 :         coeffs[2 * (i + 1) + 1] = 0.0;
     390             :     }
     391             : 
     392           0 :     for (i = 0; i < nb; i++)
     393           0 :         multiply(pz[2 * i], pz[2 * i + 1], nb, coeffs);
     394             : 
     395           0 :     for (i = 0; i < nb + 1; i++) {
     396           0 :         if (fabs(coeffs[2 * i + 1]) > FLT_EPSILON) {
     397           0 :             av_log(ctx, AV_LOG_ERROR, "coeff: %lf of z^%d is not real; poles/zeros are not complex conjugates.\n",
     398           0 :                    coeffs[2 * i + 1], i);
     399           0 :             return AVERROR(EINVAL);
     400             :         }
     401             :     }
     402             : 
     403           0 :     return 0;
     404             : }
     405             : 
     406           0 : static int convert_zp2tf(AVFilterContext *ctx, int channels)
     407             : {
     408           0 :     AudioIIRContext *s = ctx->priv;
     409           0 :     int ch, i, j, ret = 0;
     410             : 
     411           0 :     for (ch = 0; ch < channels; ch++) {
     412           0 :         IIRChannel *iir = &s->iir[ch];
     413             :         double *topc, *botc;
     414             : 
     415           0 :         topc = av_calloc((iir->nb_ab[0] + 1) * 2, sizeof(*topc));
     416           0 :         botc = av_calloc((iir->nb_ab[1] + 1) * 2, sizeof(*botc));
     417           0 :         if (!topc || !botc) {
     418           0 :             ret = AVERROR(ENOMEM);
     419           0 :             goto fail;
     420             :         }
     421             : 
     422           0 :         ret = expand(ctx, iir->ab[0], iir->nb_ab[0], botc);
     423           0 :         if (ret < 0) {
     424           0 :             goto fail;
     425             :         }
     426             : 
     427           0 :         ret = expand(ctx, iir->ab[1], iir->nb_ab[1], topc);
     428           0 :         if (ret < 0) {
     429           0 :             goto fail;
     430             :         }
     431             : 
     432           0 :         for (j = 0, i = iir->nb_ab[1]; i >= 0; j++, i--) {
     433           0 :             iir->ab[1][j] = topc[2 * i];
     434             :         }
     435           0 :         iir->nb_ab[1]++;
     436             : 
     437           0 :         for (j = 0, i = iir->nb_ab[0]; i >= 0; j++, i--) {
     438           0 :             iir->ab[0][j] = botc[2 * i];
     439             :         }
     440           0 :         iir->nb_ab[0]++;
     441             : 
     442           0 : fail:
     443           0 :         av_free(topc);
     444           0 :         av_free(botc);
     445           0 :         if (ret < 0)
     446           0 :             break;
     447             :     }
     448             : 
     449           0 :     return ret;
     450             : }
     451             : 
     452           0 : static int decompose_zp2biquads(AVFilterContext *ctx, int channels)
     453             : {
     454           0 :     AudioIIRContext *s = ctx->priv;
     455             :     int ch, ret;
     456             : 
     457           0 :     for (ch = 0; ch < channels; ch++) {
     458           0 :         IIRChannel *iir = &s->iir[ch];
     459           0 :         int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2;
     460           0 :         int current_biquad = 0;
     461             : 
     462           0 :         iir->biquads = av_calloc(nb_biquads, sizeof(BiquadContext));
     463           0 :         if (!iir->biquads)
     464           0 :             return AVERROR(ENOMEM);
     465             : 
     466           0 :         while (nb_biquads--) {
     467           0 :             Pair outmost_pole = { -1, -1 };
     468           0 :             Pair nearest_zero = { -1, -1 };
     469           0 :             double zeros[4] = { 0 };
     470           0 :             double poles[4] = { 0 };
     471           0 :             double b[6] = { 0 };
     472           0 :             double a[6] = { 0 };
     473           0 :             double min_distance = DBL_MAX;
     474           0 :             double max_mag = 0;
     475             :             int i;
     476             : 
     477           0 :             for (i = 0; i < iir->nb_ab[0]; i++) {
     478             :                 double mag;
     479             : 
     480           0 :                 if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
     481           0 :                     continue;
     482           0 :                 mag = hypot(iir->ab[0][2 * i], iir->ab[0][2 * i + 1]);
     483             : 
     484           0 :                 if (mag > max_mag) {
     485           0 :                     max_mag = mag;
     486           0 :                     outmost_pole.a = i;
     487             :                 }
     488             :             }
     489             : 
     490           0 :             for (i = 0; i < iir->nb_ab[1]; i++) {
     491           0 :                 if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
     492           0 :                     continue;
     493             : 
     494           0 :                 if (iir->ab[0][2 * i    ] ==  iir->ab[0][2 * outmost_pole.a    ] &&
     495           0 :                     iir->ab[0][2 * i + 1] == -iir->ab[0][2 * outmost_pole.a + 1]) {
     496           0 :                     outmost_pole.b = i;
     497           0 :                     break;
     498             :                 }
     499             :             }
     500             : 
     501           0 :             av_log(ctx, AV_LOG_VERBOSE, "outmost_pole is %d.%d\n", outmost_pole.a, outmost_pole.b);
     502             : 
     503           0 :             if (outmost_pole.a < 0 || outmost_pole.b < 0)
     504           0 :                 return AVERROR(EINVAL);
     505             : 
     506           0 :             for (i = 0; i < iir->nb_ab[1]; i++) {
     507             :                 double distance;
     508             : 
     509           0 :                 if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
     510           0 :                     continue;
     511           0 :                 distance = hypot(iir->ab[0][2 * outmost_pole.a    ] - iir->ab[1][2 * i    ],
     512           0 :                                  iir->ab[0][2 * outmost_pole.a + 1] - iir->ab[1][2 * i + 1]);
     513             : 
     514           0 :                 if (distance < min_distance) {
     515           0 :                     min_distance = distance;
     516           0 :                     nearest_zero.a = i;
     517             :                 }
     518             :             }
     519             : 
     520           0 :             for (i = 0; i < iir->nb_ab[1]; i++) {
     521           0 :                 if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
     522           0 :                     continue;
     523             : 
     524           0 :                 if (iir->ab[1][2 * i    ] ==  iir->ab[1][2 * nearest_zero.a    ] &&
     525           0 :                     iir->ab[1][2 * i + 1] == -iir->ab[1][2 * nearest_zero.a + 1]) {
     526           0 :                     nearest_zero.b = i;
     527           0 :                     break;
     528             :                 }
     529             :             }
     530             : 
     531           0 :             av_log(ctx, AV_LOG_VERBOSE, "nearest_zero is %d.%d\n", nearest_zero.a, nearest_zero.b);
     532             : 
     533           0 :             if (nearest_zero.a < 0 || nearest_zero.b < 0)
     534           0 :                 return AVERROR(EINVAL);
     535             : 
     536           0 :             poles[0] = iir->ab[0][2 * outmost_pole.a    ];
     537           0 :             poles[1] = iir->ab[0][2 * outmost_pole.a + 1];
     538             : 
     539           0 :             zeros[0] = iir->ab[1][2 * nearest_zero.a    ];
     540           0 :             zeros[1] = iir->ab[1][2 * nearest_zero.a + 1];
     541             : 
     542           0 :             if (nearest_zero.a == nearest_zero.b && outmost_pole.a == outmost_pole.b) {
     543           0 :                 zeros[2] = 0;
     544           0 :                 zeros[3] = 0;
     545             : 
     546           0 :                 poles[2] = 0;
     547           0 :                 poles[3] = 0;
     548             :             } else {
     549           0 :                 poles[2] = iir->ab[0][2 * outmost_pole.b    ];
     550           0 :                 poles[3] = iir->ab[0][2 * outmost_pole.b + 1];
     551             : 
     552           0 :                 zeros[2] = iir->ab[1][2 * nearest_zero.b    ];
     553           0 :                 zeros[3] = iir->ab[1][2 * nearest_zero.b + 1];
     554             :             }
     555             : 
     556           0 :             ret = expand(ctx, zeros, 2, b);
     557           0 :             if (ret < 0)
     558           0 :                 return ret;
     559             : 
     560           0 :             ret = expand(ctx, poles, 2, a);
     561           0 :             if (ret < 0)
     562           0 :                 return ret;
     563             : 
     564           0 :             iir->ab[0][2 * outmost_pole.a] = iir->ab[0][2 * outmost_pole.a + 1] = NAN;
     565           0 :             iir->ab[0][2 * outmost_pole.b] = iir->ab[0][2 * outmost_pole.b + 1] = NAN;
     566           0 :             iir->ab[1][2 * nearest_zero.a] = iir->ab[1][2 * nearest_zero.a + 1] = NAN;
     567           0 :             iir->ab[1][2 * nearest_zero.b] = iir->ab[1][2 * nearest_zero.b + 1] = NAN;
     568             : 
     569           0 :             iir->biquads[current_biquad].a0 = 1.0;
     570           0 :             iir->biquads[current_biquad].a1 = a[2] / a[4];
     571           0 :             iir->biquads[current_biquad].a2 = a[0] / a[4];
     572           0 :             iir->biquads[current_biquad].b0 = b[4] / a[4] * (current_biquad ? 1.0 : iir->g);
     573           0 :             iir->biquads[current_biquad].b1 = b[2] / a[4] * (current_biquad ? 1.0 : iir->g);
     574           0 :             iir->biquads[current_biquad].b2 = b[0] / a[4] * (current_biquad ? 1.0 : iir->g);
     575             : 
     576           0 :             av_log(ctx, AV_LOG_VERBOSE, "a=%lf %lf %lf:b=%lf %lf %lf\n",
     577           0 :                    iir->biquads[current_biquad].a0,
     578           0 :                    iir->biquads[current_biquad].a1,
     579           0 :                    iir->biquads[current_biquad].a2,
     580           0 :                    iir->biquads[current_biquad].b0,
     581           0 :                    iir->biquads[current_biquad].b1,
     582           0 :                    iir->biquads[current_biquad].b2);
     583             : 
     584           0 :             current_biquad++;
     585             :         }
     586             :     }
     587             : 
     588           0 :     return 0;
     589             : }
     590             : 
     591           0 : static void convert_pr2zp(AVFilterContext *ctx, int channels)
     592             : {
     593           0 :     AudioIIRContext *s = ctx->priv;
     594             :     int ch;
     595             : 
     596           0 :     for (ch = 0; ch < channels; ch++) {
     597           0 :         IIRChannel *iir = &s->iir[ch];
     598             :         int n;
     599             : 
     600           0 :         for (n = 0; n < iir->nb_ab[0]; n++) {
     601           0 :             double r = iir->ab[0][2*n];
     602           0 :             double angle = iir->ab[0][2*n+1];
     603             : 
     604           0 :             iir->ab[0][2*n]   = r * cos(angle);
     605           0 :             iir->ab[0][2*n+1] = r * sin(angle);
     606             :         }
     607             : 
     608           0 :         for (n = 0; n < iir->nb_ab[1]; n++) {
     609           0 :             double r = iir->ab[1][2*n];
     610           0 :             double angle = iir->ab[1][2*n+1];
     611             : 
     612           0 :             iir->ab[1][2*n]   = r * cos(angle);
     613           0 :             iir->ab[1][2*n+1] = r * sin(angle);
     614             :         }
     615             :     }
     616           0 : }
     617             : 
     618           0 : static void convert_pd2zp(AVFilterContext *ctx, int channels)
     619             : {
     620           0 :     AudioIIRContext *s = ctx->priv;
     621             :     int ch;
     622             : 
     623           0 :     for (ch = 0; ch < channels; ch++) {
     624           0 :         IIRChannel *iir = &s->iir[ch];
     625             :         int n;
     626             : 
     627           0 :         for (n = 0; n < iir->nb_ab[0]; n++) {
     628           0 :             double r = iir->ab[0][2*n];
     629           0 :             double angle = M_PI*iir->ab[0][2*n+1]/180.;
     630             : 
     631           0 :             iir->ab[0][2*n]   = r * cos(angle);
     632           0 :             iir->ab[0][2*n+1] = r * sin(angle);
     633             :         }
     634             : 
     635           0 :         for (n = 0; n < iir->nb_ab[1]; n++) {
     636           0 :             double r = iir->ab[1][2*n];
     637           0 :             double angle = M_PI*iir->ab[1][2*n+1]/180.;
     638             : 
     639           0 :             iir->ab[1][2*n]   = r * cos(angle);
     640           0 :             iir->ab[1][2*n+1] = r * sin(angle);
     641             :         }
     642             :     }
     643           0 : }
     644             : 
     645           0 : static int config_output(AVFilterLink *outlink)
     646             : {
     647           0 :     AVFilterContext *ctx = outlink->src;
     648           0 :     AudioIIRContext *s = ctx->priv;
     649           0 :     AVFilterLink *inlink = ctx->inputs[0];
     650             :     int ch, ret, i;
     651             : 
     652           0 :     s->channels = inlink->channels;
     653           0 :     s->iir = av_calloc(s->channels, sizeof(*s->iir));
     654           0 :     if (!s->iir)
     655           0 :         return AVERROR(ENOMEM);
     656             : 
     657           0 :     ret = read_gains(ctx, s->g_str, inlink->channels);
     658           0 :     if (ret < 0)
     659           0 :         return ret;
     660             : 
     661           0 :     ret = read_channels(ctx, inlink->channels, s->a_str, 0);
     662           0 :     if (ret < 0)
     663           0 :         return ret;
     664             : 
     665           0 :     ret = read_channels(ctx, inlink->channels, s->b_str, 1);
     666           0 :     if (ret < 0)
     667           0 :         return ret;
     668             : 
     669           0 :     if (s->format == 2) {
     670           0 :         convert_pr2zp(ctx, inlink->channels);
     671           0 :     } else if (s->format == 3) {
     672           0 :         convert_pd2zp(ctx, inlink->channels);
     673             :     }
     674             : 
     675           0 :     if (s->format == 0)
     676           0 :         av_log(ctx, AV_LOG_WARNING, "tf coefficients format is not recommended for too high number of zeros/poles.\n");
     677             : 
     678           0 :     if (s->format > 0 && s->process == 0) {
     679           0 :         av_log(ctx, AV_LOG_WARNING, "Direct processsing is not recommended for zp coefficients format.\n");
     680             : 
     681           0 :         ret = convert_zp2tf(ctx, inlink->channels);
     682           0 :         if (ret < 0)
     683           0 :             return ret;
     684           0 :     } else if (s->format == 0 && s->process == 1) {
     685           0 :         av_log(ctx, AV_LOG_ERROR, "Serial cascading is not implemented for transfer function.\n");
     686           0 :         return AVERROR_PATCHWELCOME;
     687           0 :     } else if (s->format > 0 && s->process == 1) {
     688           0 :         if (inlink->format == AV_SAMPLE_FMT_S16P)
     689           0 :             av_log(ctx, AV_LOG_WARNING, "Serial cascading is not recommended for i16 precision.\n");
     690             : 
     691           0 :         ret = decompose_zp2biquads(ctx, inlink->channels);
     692           0 :         if (ret < 0)
     693           0 :             return ret;
     694             :     }
     695             : 
     696           0 :     for (ch = 0; ch < inlink->channels; ch++) {
     697           0 :         IIRChannel *iir = &s->iir[ch];
     698             : 
     699           0 :         for (i = 1; i < iir->nb_ab[0]; i++) {
     700           0 :             iir->ab[0][i] /= iir->ab[0][0];
     701             :         }
     702             : 
     703           0 :         for (i = 0; i < iir->nb_ab[1]; i++) {
     704           0 :             iir->ab[1][i] *= iir->g / iir->ab[0][0];
     705             :         }
     706             :     }
     707             : 
     708           0 :     switch (inlink->format) {
     709           0 :     case AV_SAMPLE_FMT_DBLP: s->iir_channel = s->process == 1 ? iir_ch_serial_dblp : iir_ch_dblp; break;
     710           0 :     case AV_SAMPLE_FMT_FLTP: s->iir_channel = s->process == 1 ? iir_ch_serial_fltp : iir_ch_fltp; break;
     711           0 :     case AV_SAMPLE_FMT_S32P: s->iir_channel = s->process == 1 ? iir_ch_serial_s32p : iir_ch_s32p; break;
     712           0 :     case AV_SAMPLE_FMT_S16P: s->iir_channel = s->process == 1 ? iir_ch_serial_s16p : iir_ch_s16p; break;
     713             :     }
     714             : 
     715           0 :     return 0;
     716             : }
     717             : 
     718           0 : static int filter_frame(AVFilterLink *inlink, AVFrame *in)
     719             : {
     720           0 :     AVFilterContext *ctx = inlink->dst;
     721           0 :     AudioIIRContext *s = ctx->priv;
     722           0 :     AVFilterLink *outlink = ctx->outputs[0];
     723             :     ThreadData td;
     724             :     AVFrame *out;
     725             :     int ch;
     726             : 
     727           0 :     if (av_frame_is_writable(in)) {
     728           0 :         out = in;
     729             :     } else {
     730           0 :         out = ff_get_audio_buffer(outlink, in->nb_samples);
     731           0 :         if (!out) {
     732           0 :             av_frame_free(&in);
     733           0 :             return AVERROR(ENOMEM);
     734             :         }
     735           0 :         av_frame_copy_props(out, in);
     736             :     }
     737             : 
     738           0 :     td.in  = in;
     739           0 :     td.out = out;
     740           0 :     ctx->internal->execute(ctx, s->iir_channel, &td, NULL, outlink->channels);
     741             : 
     742           0 :     for (ch = 0; ch < outlink->channels; ch++) {
     743           0 :         if (s->iir[ch].clippings > 0)
     744           0 :             av_log(ctx, AV_LOG_WARNING, "Channel %d clipping %d times. Please reduce gain.\n",
     745           0 :                    ch, s->iir[ch].clippings);
     746           0 :         s->iir[ch].clippings = 0;
     747             :     }
     748             : 
     749           0 :     if (in != out)
     750           0 :         av_frame_free(&in);
     751             : 
     752           0 :     return ff_filter_frame(outlink, out);
     753             : }
     754             : 
     755           0 : static av_cold int init(AVFilterContext *ctx)
     756             : {
     757           0 :     AudioIIRContext *s = ctx->priv;
     758             : 
     759           0 :     if (!s->a_str || !s->b_str || !s->g_str) {
     760           0 :         av_log(ctx, AV_LOG_ERROR, "Valid coefficients are mandatory.\n");
     761           0 :         return AVERROR(EINVAL);
     762             :     }
     763             : 
     764           0 :     switch (s->precision) {
     765           0 :     case 0: s->sample_format = AV_SAMPLE_FMT_DBLP; break;
     766           0 :     case 1: s->sample_format = AV_SAMPLE_FMT_FLTP; break;
     767           0 :     case 2: s->sample_format = AV_SAMPLE_FMT_S32P; break;
     768           0 :     case 3: s->sample_format = AV_SAMPLE_FMT_S16P; break;
     769           0 :     default: return AVERROR_BUG;
     770             :     }
     771             : 
     772           0 :     return 0;
     773             : }
     774             : 
     775           0 : static av_cold void uninit(AVFilterContext *ctx)
     776             : {
     777           0 :     AudioIIRContext *s = ctx->priv;
     778             :     int ch;
     779             : 
     780           0 :     if (s->iir) {
     781           0 :         for (ch = 0; ch < s->channels; ch++) {
     782           0 :             IIRChannel *iir = &s->iir[ch];
     783           0 :             av_freep(&iir->ab[0]);
     784           0 :             av_freep(&iir->ab[1]);
     785           0 :             av_freep(&iir->cache[0]);
     786           0 :             av_freep(&iir->cache[1]);
     787           0 :             av_freep(&iir->biquads);
     788             :         }
     789             :     }
     790           0 :     av_freep(&s->iir);
     791           0 : }
     792             : 
     793             : static const AVFilterPad inputs[] = {
     794             :     {
     795             :         .name         = "default",
     796             :         .type         = AVMEDIA_TYPE_AUDIO,
     797             :         .filter_frame = filter_frame,
     798             :     },
     799             :     { NULL }
     800             : };
     801             : 
     802             : static const AVFilterPad outputs[] = {
     803             :     {
     804             :         .name         = "default",
     805             :         .type         = AVMEDIA_TYPE_AUDIO,
     806             :         .config_props = config_output,
     807             :     },
     808             :     { NULL }
     809             : };
     810             : 
     811             : #define OFFSET(x) offsetof(AudioIIRContext, x)
     812             : #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
     813             : 
     814             : static const AVOption aiir_options[] = {
     815             :     { "z", "set B/numerator/zeros coefficients",   OFFSET(b_str),    AV_OPT_TYPE_STRING, {.str="1+0i 1-0i"}, 0, 0, AF },
     816             :     { "p", "set A/denominator/poles coefficients", OFFSET(a_str),    AV_OPT_TYPE_STRING, {.str="1+0i 1-0i"}, 0, 0, AF },
     817             :     { "k", "set channels gains",                   OFFSET(g_str),    AV_OPT_TYPE_STRING, {.str="1|1"}, 0, 0, AF },
     818             :     { "dry", "set dry gain",                       OFFSET(dry_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1},     0, 1, AF },
     819             :     { "wet", "set wet gain",                       OFFSET(wet_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1},     0, 1, AF },
     820             :     { "f", "set coefficients format",              OFFSET(format),   AV_OPT_TYPE_INT,    {.i64=1},     0, 3, AF, "format" },
     821             :     { "tf", "transfer function",                   0,                AV_OPT_TYPE_CONST,  {.i64=0},     0, 0, AF, "format" },
     822             :     { "zp", "Z-plane zeros/poles",                 0,                AV_OPT_TYPE_CONST,  {.i64=1},     0, 0, AF, "format" },
     823             :     { "pr", "Z-plane zeros/poles (polar radians)", 0,                AV_OPT_TYPE_CONST,  {.i64=2},     0, 0, AF, "format" },
     824             :     { "pd", "Z-plane zeros/poles (polar degrees)", 0,                AV_OPT_TYPE_CONST,  {.i64=3},     0, 0, AF, "format" },
     825             :     { "r", "set kind of processing",               OFFSET(process),  AV_OPT_TYPE_INT,    {.i64=1},     0, 1, AF, "process" },
     826             :     { "d", "direct",                               0,                AV_OPT_TYPE_CONST,  {.i64=0},     0, 0, AF, "process" },
     827             :     { "s", "serial cascading",                     0,                AV_OPT_TYPE_CONST,  {.i64=1},     0, 0, AF, "process" },
     828             :     { "e", "set precision",                        OFFSET(precision),AV_OPT_TYPE_INT,    {.i64=0},     0, 3, AF, "precision" },
     829             :     { "dbl", "double-precision floating-point",    0,                AV_OPT_TYPE_CONST,  {.i64=0},     0, 0, AF, "precision" },
     830             :     { "flt", "single-precision floating-point",    0,                AV_OPT_TYPE_CONST,  {.i64=1},     0, 0, AF, "precision" },
     831             :     { "i32", "32-bit integers",                    0,                AV_OPT_TYPE_CONST,  {.i64=2},     0, 0, AF, "precision" },
     832             :     { "i16", "16-bit integers",                    0,                AV_OPT_TYPE_CONST,  {.i64=3},     0, 0, AF, "precision" },
     833             :     { NULL },
     834             : };
     835             : 
     836             : AVFILTER_DEFINE_CLASS(aiir);
     837             : 
     838             : AVFilter ff_af_aiir = {
     839             :     .name          = "aiir",
     840             :     .description   = NULL_IF_CONFIG_SMALL("Apply Infinite Impulse Response filter with supplied coefficients."),
     841             :     .priv_size     = sizeof(AudioIIRContext),
     842             :     .priv_class    = &aiir_class,
     843             :     .init          = init,
     844             :     .uninit        = uninit,
     845             :     .query_formats = query_formats,
     846             :     .inputs        = inputs,
     847             :     .outputs       = outputs,
     848             :     .flags         = AVFILTER_FLAG_SLICE_THREADS,
     849             : };

Generated by: LCOV version 1.13