avcodec/vc1: correct ff_vc1_mbmode_intfrp
[ffmpeg.git] / libavfilter / af_aiir.c
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 static int query_formats(AVFilterContext *ctx)
70 {
71     AudioIIRContext *s = ctx->priv;
72     AVFilterFormats *formats;
73     AVFilterChannelLayouts *layouts;
74     enum AVSampleFormat sample_fmts[] = {
75         AV_SAMPLE_FMT_DBLP,
76         AV_SAMPLE_FMT_NONE
77     };
78     int ret;
79
80     layouts = ff_all_channel_counts();
81     if (!layouts)
82         return AVERROR(ENOMEM);
83     ret = ff_set_common_channel_layouts(ctx, layouts);
84     if (ret < 0)
85         return ret;
86
87     sample_fmts[0] = s->sample_format;
88     formats = ff_make_format_list(sample_fmts);
89     if (!formats)
90         return AVERROR(ENOMEM);
91     ret = ff_set_common_formats(ctx, formats);
92     if (ret < 0)
93         return ret;
94
95     formats = ff_all_samplerates();
96     if (!formats)
97         return AVERROR(ENOMEM);
98     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 IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
150 IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
151 IIR_CH(fltp, float,         -1.,        1., 0)
152 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 SERIAL_IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
210 SERIAL_IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
211 SERIAL_IIR_CH(fltp, float,         -1.,        1., 0)
212 SERIAL_IIR_CH(dblp, double,        -1.,        1., 0)
213
214 static void count_coefficients(char *item_str, int *nb_items)
215 {
216     char *p;
217
218     if (!item_str)
219         return;
220
221     *nb_items = 1;
222     for (p = item_str; *p && *p != '|'; p++) {
223         if (*p == ' ')
224             (*nb_items)++;
225     }
226 }
227
228 static int read_gains(AVFilterContext *ctx, char *item_str, int nb_items)
229 {
230     AudioIIRContext *s = ctx->priv;
231     char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
232     int i;
233
234     p = old_str = av_strdup(item_str);
235     if (!p)
236         return AVERROR(ENOMEM);
237     for (i = 0; i < nb_items; i++) {
238         if (!(arg = av_strtok(p, "|", &saveptr)))
239             arg = prev_arg;
240
241         if (!arg) {
242             av_freep(&old_str);
243             return AVERROR(EINVAL);
244         }
245
246         p = NULL;
247         if (sscanf(arg, "%lf", &s->iir[i].g) != 1) {
248             av_log(ctx, AV_LOG_ERROR, "Invalid gains supplied: %s\n", arg);
249             av_freep(&old_str);
250             return AVERROR(EINVAL);
251         }
252
253         prev_arg = arg;
254     }
255
256     av_freep(&old_str);
257
258     return 0;
259 }
260
261 static int read_tf_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst)
262 {
263     char *p, *arg, *old_str, *saveptr = NULL;
264     int i;
265
266     p = old_str = av_strdup(item_str);
267     if (!p)
268         return AVERROR(ENOMEM);
269     for (i = 0; i < nb_items; i++) {
270         if (!(arg = av_strtok(p, " ", &saveptr)))
271             break;
272
273         p = NULL;
274         if (sscanf(arg, "%lf", &dst[i]) != 1) {
275             av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
276             av_freep(&old_str);
277             return AVERROR(EINVAL);
278         }
279     }
280
281     av_freep(&old_str);
282
283     return 0;
284 }
285
286 static int read_zp_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst, const char *format)
287 {
288     char *p, *arg, *old_str, *saveptr = NULL;
289     int i;
290
291     p = old_str = av_strdup(item_str);
292     if (!p)
293         return AVERROR(ENOMEM);
294     for (i = 0; i < nb_items; i++) {
295         if (!(arg = av_strtok(p, " ", &saveptr)))
296             break;
297
298         p = NULL;
299         if (sscanf(arg, format, &dst[i*2], &dst[i*2+1]) != 2) {
300             av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
301             av_freep(&old_str);
302             return AVERROR(EINVAL);
303         }
304     }
305
306     av_freep(&old_str);
307
308     return 0;
309 }
310
311 static const char *format[] = { "%lf", "%lf %lfi", "%lf %lfr", "%lf %lfd" };
312
313 static int read_channels(AVFilterContext *ctx, int channels, uint8_t *item_str, int ab)
314 {
315     AudioIIRContext *s = ctx->priv;
316     char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
317     int i, ret;
318
319     p = old_str = av_strdup(item_str);
320     if (!p)
321         return AVERROR(ENOMEM);
322     for (i = 0; i < channels; i++) {
323         IIRChannel *iir = &s->iir[i];
324
325         if (!(arg = av_strtok(p, "|", &saveptr)))
326             arg = prev_arg;
327
328         if (!arg) {
329             av_freep(&old_str);
330             return AVERROR(EINVAL);
331         }
332
333         count_coefficients(arg, &iir->nb_ab[ab]);
334
335         p = NULL;
336         iir->cache[ab] = av_calloc(iir->nb_ab[ab] + 1, sizeof(double));
337         iir->ab[ab] = av_calloc(iir->nb_ab[ab] * (!!s->format + 1), sizeof(double));
338         if (!iir->ab[ab] || !iir->cache[ab]) {
339             av_freep(&old_str);
340             return AVERROR(ENOMEM);
341         }
342
343         if (s->format) {
344             ret = read_zp_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab], format[s->format]);
345         } else {
346             ret = read_tf_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab]);
347         }
348         if (ret < 0) {
349             av_freep(&old_str);
350             return ret;
351         }
352         prev_arg = arg;
353     }
354
355     av_freep(&old_str);
356
357     return 0;
358 }
359
360 static void multiply(double wre, double wim, int npz, double *coeffs)
361 {
362     double nwre = -wre, nwim = -wim;
363     double cre, cim;
364     int i;
365
366     for (i = npz; i >= 1; i--) {
367         cre = coeffs[2 * i + 0];
368         cim = coeffs[2 * i + 1];
369
370         coeffs[2 * i + 0] = (nwre * cre - nwim * cim) + coeffs[2 * (i - 1) + 0];
371         coeffs[2 * i + 1] = (nwre * cim + nwim * cre) + coeffs[2 * (i - 1) + 1];
372     }
373
374     cre = coeffs[0];
375     cim = coeffs[1];
376     coeffs[0] = nwre * cre - nwim * cim;
377     coeffs[1] = nwre * cim + nwim * cre;
378 }
379
380 static int expand(AVFilterContext *ctx, double *pz, int nb, double *coeffs)
381 {
382     int i;
383
384     coeffs[0] = 1.0;
385     coeffs[1] = 0.0;
386
387     for (i = 0; i < nb; i++) {
388         coeffs[2 * (i + 1)    ] = 0.0;
389         coeffs[2 * (i + 1) + 1] = 0.0;
390     }
391
392     for (i = 0; i < nb; i++)
393         multiply(pz[2 * i], pz[2 * i + 1], nb, coeffs);
394
395     for (i = 0; i < nb + 1; i++) {
396         if (fabs(coeffs[2 * i + 1]) > FLT_EPSILON) {
397             av_log(ctx, AV_LOG_ERROR, "coeff: %lf of z^%d is not real; poles/zeros are not complex conjugates.\n",
398                    coeffs[2 * i + 1], i);
399             return AVERROR(EINVAL);
400         }
401     }
402
403     return 0;
404 }
405
406 static int convert_zp2tf(AVFilterContext *ctx, int channels)
407 {
408     AudioIIRContext *s = ctx->priv;
409     int ch, i, j, ret = 0;
410
411     for (ch = 0; ch < channels; ch++) {
412         IIRChannel *iir = &s->iir[ch];
413         double *topc, *botc;
414
415         topc = av_calloc((iir->nb_ab[0] + 1) * 2, sizeof(*topc));
416         botc = av_calloc((iir->nb_ab[1] + 1) * 2, sizeof(*botc));
417         if (!topc || !botc) {
418             ret = AVERROR(ENOMEM);
419             goto fail;
420         }
421
422         ret = expand(ctx, iir->ab[0], iir->nb_ab[0], botc);
423         if (ret < 0) {
424             goto fail;
425         }
426
427         ret = expand(ctx, iir->ab[1], iir->nb_ab[1], topc);
428         if (ret < 0) {
429             goto fail;
430         }
431
432         for (j = 0, i = iir->nb_ab[1]; i >= 0; j++, i--) {
433             iir->ab[1][j] = topc[2 * i];
434         }
435         iir->nb_ab[1]++;
436
437         for (j = 0, i = iir->nb_ab[0]; i >= 0; j++, i--) {
438             iir->ab[0][j] = botc[2 * i];
439         }
440         iir->nb_ab[0]++;
441
442 fail:
443         av_free(topc);
444         av_free(botc);
445         if (ret < 0)
446             break;
447     }
448
449     return ret;
450 }
451
452 static int decompose_zp2biquads(AVFilterContext *ctx, int channels)
453 {
454     AudioIIRContext *s = ctx->priv;
455     int ch, ret;
456
457     for (ch = 0; ch < channels; ch++) {
458         IIRChannel *iir = &s->iir[ch];
459         int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2;
460         int current_biquad = 0;
461
462         iir->biquads = av_calloc(nb_biquads, sizeof(BiquadContext));
463         if (!iir->biquads)
464             return AVERROR(ENOMEM);
465
466         while (nb_biquads--) {
467             Pair outmost_pole = { -1, -1 };
468             Pair nearest_zero = { -1, -1 };
469             double zeros[4] = { 0 };
470             double poles[4] = { 0 };
471             double b[6] = { 0 };
472             double a[6] = { 0 };
473             double min_distance = DBL_MAX;
474             double max_mag = 0;
475             int i;
476
477             for (i = 0; i < iir->nb_ab[0]; i++) {
478                 double mag;
479
480                 if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
481                     continue;
482                 mag = hypot(iir->ab[0][2 * i], iir->ab[0][2 * i + 1]);
483
484                 if (mag > max_mag) {
485                     max_mag = mag;
486                     outmost_pole.a = i;
487                 }
488             }
489
490             for (i = 0; i < iir->nb_ab[1]; i++) {
491                 if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
492                     continue;
493
494                 if (iir->ab[0][2 * i    ] ==  iir->ab[0][2 * outmost_pole.a    ] &&
495                     iir->ab[0][2 * i + 1] == -iir->ab[0][2 * outmost_pole.a + 1]) {
496                     outmost_pole.b = i;
497                     break;
498                 }
499             }
500
501             av_log(ctx, AV_LOG_VERBOSE, "outmost_pole is %d.%d\n", outmost_pole.a, outmost_pole.b);
502
503             if (outmost_pole.a < 0 || outmost_pole.b < 0)
504                 return AVERROR(EINVAL);
505
506             for (i = 0; i < iir->nb_ab[1]; i++) {
507                 double distance;
508
509                 if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
510                     continue;
511                 distance = hypot(iir->ab[0][2 * outmost_pole.a    ] - iir->ab[1][2 * i    ],
512                                  iir->ab[0][2 * outmost_pole.a + 1] - iir->ab[1][2 * i + 1]);
513
514                 if (distance < min_distance) {
515                     min_distance = distance;
516                     nearest_zero.a = i;
517                 }
518             }
519
520             for (i = 0; i < iir->nb_ab[1]; i++) {
521                 if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
522                     continue;
523
524                 if (iir->ab[1][2 * i    ] ==  iir->ab[1][2 * nearest_zero.a    ] &&
525                     iir->ab[1][2 * i + 1] == -iir->ab[1][2 * nearest_zero.a + 1]) {
526                     nearest_zero.b = i;
527                     break;
528                 }
529             }
530
531             av_log(ctx, AV_LOG_VERBOSE, "nearest_zero is %d.%d\n", nearest_zero.a, nearest_zero.b);
532
533             if (nearest_zero.a < 0 || nearest_zero.b < 0)
534                 return AVERROR(EINVAL);
535
536             poles[0] = iir->ab[0][2 * outmost_pole.a    ];
537             poles[1] = iir->ab[0][2 * outmost_pole.a + 1];
538
539             zeros[0] = iir->ab[1][2 * nearest_zero.a    ];
540             zeros[1] = iir->ab[1][2 * nearest_zero.a + 1];
541
542             if (nearest_zero.a == nearest_zero.b && outmost_pole.a == outmost_pole.b) {
543                 zeros[2] = 0;
544                 zeros[3] = 0;
545
546                 poles[2] = 0;
547                 poles[3] = 0;
548             } else {
549                 poles[2] = iir->ab[0][2 * outmost_pole.b    ];
550                 poles[3] = iir->ab[0][2 * outmost_pole.b + 1];
551
552                 zeros[2] = iir->ab[1][2 * nearest_zero.b    ];
553                 zeros[3] = iir->ab[1][2 * nearest_zero.b + 1];
554             }
555
556             ret = expand(ctx, zeros, 2, b);
557             if (ret < 0)
558                 return ret;
559
560             ret = expand(ctx, poles, 2, a);
561             if (ret < 0)
562                 return ret;
563
564             iir->ab[0][2 * outmost_pole.a] = iir->ab[0][2 * outmost_pole.a + 1] = NAN;
565             iir->ab[0][2 * outmost_pole.b] = iir->ab[0][2 * outmost_pole.b + 1] = NAN;
566             iir->ab[1][2 * nearest_zero.a] = iir->ab[1][2 * nearest_zero.a + 1] = NAN;
567             iir->ab[1][2 * nearest_zero.b] = iir->ab[1][2 * nearest_zero.b + 1] = NAN;
568
569             iir->biquads[current_biquad].a0 = 1.0;
570             iir->biquads[current_biquad].a1 = a[2] / a[4];
571             iir->biquads[current_biquad].a2 = a[0] / a[4];
572             iir->biquads[current_biquad].b0 = b[4] / a[4] * (current_biquad ? 1.0 : iir->g);
573             iir->biquads[current_biquad].b1 = b[2] / a[4] * (current_biquad ? 1.0 : iir->g);
574             iir->biquads[current_biquad].b2 = b[0] / a[4] * (current_biquad ? 1.0 : iir->g);
575
576             av_log(ctx, AV_LOG_VERBOSE, "a=%lf %lf %lf:b=%lf %lf %lf\n",
577                    iir->biquads[current_biquad].a0,
578                    iir->biquads[current_biquad].a1,
579                    iir->biquads[current_biquad].a2,
580                    iir->biquads[current_biquad].b0,
581                    iir->biquads[current_biquad].b1,
582                    iir->biquads[current_biquad].b2);
583
584             current_biquad++;
585         }
586     }
587
588     return 0;
589 }
590
591 static void convert_pr2zp(AVFilterContext *ctx, int channels)
592 {
593     AudioIIRContext *s = ctx->priv;
594     int ch;
595
596     for (ch = 0; ch < channels; ch++) {
597         IIRChannel *iir = &s->iir[ch];
598         int n;
599
600         for (n = 0; n < iir->nb_ab[0]; n++) {
601             double r = iir->ab[0][2*n];
602             double angle = iir->ab[0][2*n+1];
603
604             iir->ab[0][2*n]   = r * cos(angle);
605             iir->ab[0][2*n+1] = r * sin(angle);
606         }
607
608         for (n = 0; n < iir->nb_ab[1]; n++) {
609             double r = iir->ab[1][2*n];
610             double angle = iir->ab[1][2*n+1];
611
612             iir->ab[1][2*n]   = r * cos(angle);
613             iir->ab[1][2*n+1] = r * sin(angle);
614         }
615     }
616 }
617
618 static void convert_pd2zp(AVFilterContext *ctx, int channels)
619 {
620     AudioIIRContext *s = ctx->priv;
621     int ch;
622
623     for (ch = 0; ch < channels; ch++) {
624         IIRChannel *iir = &s->iir[ch];
625         int n;
626
627         for (n = 0; n < iir->nb_ab[0]; n++) {
628             double r = iir->ab[0][2*n];
629             double angle = M_PI*iir->ab[0][2*n+1]/180.;
630
631             iir->ab[0][2*n]   = r * cos(angle);
632             iir->ab[0][2*n+1] = r * sin(angle);
633         }
634
635         for (n = 0; n < iir->nb_ab[1]; n++) {
636             double r = iir->ab[1][2*n];
637             double angle = M_PI*iir->ab[1][2*n+1]/180.;
638
639             iir->ab[1][2*n]   = r * cos(angle);
640             iir->ab[1][2*n+1] = r * sin(angle);
641         }
642     }
643 }
644
645 static int config_output(AVFilterLink *outlink)
646 {
647     AVFilterContext *ctx = outlink->src;
648     AudioIIRContext *s = ctx->priv;
649     AVFilterLink *inlink = ctx->inputs[0];
650     int ch, ret, i;
651
652     s->channels = inlink->channels;
653     s->iir = av_calloc(s->channels, sizeof(*s->iir));
654     if (!s->iir)
655         return AVERROR(ENOMEM);
656
657     ret = read_gains(ctx, s->g_str, inlink->channels);
658     if (ret < 0)
659         return ret;
660
661     ret = read_channels(ctx, inlink->channels, s->a_str, 0);
662     if (ret < 0)
663         return ret;
664
665     ret = read_channels(ctx, inlink->channels, s->b_str, 1);
666     if (ret < 0)
667         return ret;
668
669     if (s->format == 2) {
670         convert_pr2zp(ctx, inlink->channels);
671     } else if (s->format == 3) {
672         convert_pd2zp(ctx, inlink->channels);
673     }
674
675     if (s->format == 0)
676         av_log(ctx, AV_LOG_WARNING, "tf coefficients format is not recommended for too high number of zeros/poles.\n");
677
678     if (s->format > 0 && s->process == 0) {
679         av_log(ctx, AV_LOG_WARNING, "Direct processsing is not recommended for zp coefficients format.\n");
680
681         ret = convert_zp2tf(ctx, inlink->channels);
682         if (ret < 0)
683             return ret;
684     } else if (s->format == 0 && s->process == 1) {
685         av_log(ctx, AV_LOG_ERROR, "Serial cascading is not implemented for transfer function.\n");
686         return AVERROR_PATCHWELCOME;
687     } else if (s->format > 0 && s->process == 1) {
688         if (inlink->format == AV_SAMPLE_FMT_S16P)
689             av_log(ctx, AV_LOG_WARNING, "Serial cascading is not recommended for i16 precision.\n");
690
691         ret = decompose_zp2biquads(ctx, inlink->channels);
692         if (ret < 0)
693             return ret;
694     }
695
696     for (ch = 0; ch < inlink->channels; ch++) {
697         IIRChannel *iir = &s->iir[ch];
698
699         for (i = 1; i < iir->nb_ab[0]; i++) {
700             iir->ab[0][i] /= iir->ab[0][0];
701         }
702
703         for (i = 0; i < iir->nb_ab[1]; i++) {
704             iir->ab[1][i] *= iir->g / iir->ab[0][0];
705         }
706     }
707
708     switch (inlink->format) {
709     case AV_SAMPLE_FMT_DBLP: s->iir_channel = s->process == 1 ? iir_ch_serial_dblp : iir_ch_dblp; break;
710     case AV_SAMPLE_FMT_FLTP: s->iir_channel = s->process == 1 ? iir_ch_serial_fltp : iir_ch_fltp; break;
711     case AV_SAMPLE_FMT_S32P: s->iir_channel = s->process == 1 ? iir_ch_serial_s32p : iir_ch_s32p; break;
712     case AV_SAMPLE_FMT_S16P: s->iir_channel = s->process == 1 ? iir_ch_serial_s16p : iir_ch_s16p; break;
713     }
714
715     return 0;
716 }
717
718 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
719 {
720     AVFilterContext *ctx = inlink->dst;
721     AudioIIRContext *s = ctx->priv;
722     AVFilterLink *outlink = ctx->outputs[0];
723     ThreadData td;
724     AVFrame *out;
725     int ch;
726
727     if (av_frame_is_writable(in)) {
728         out = in;
729     } else {
730         out = ff_get_audio_buffer(outlink, in->nb_samples);
731         if (!out) {
732             av_frame_free(&in);
733             return AVERROR(ENOMEM);
734         }
735         av_frame_copy_props(out, in);
736     }
737
738     td.in  = in;
739     td.out = out;
740     ctx->internal->execute(ctx, s->iir_channel, &td, NULL, outlink->channels);
741
742     for (ch = 0; ch < outlink->channels; ch++) {
743         if (s->iir[ch].clippings > 0)
744             av_log(ctx, AV_LOG_WARNING, "Channel %d clipping %d times. Please reduce gain.\n",
745                    ch, s->iir[ch].clippings);
746         s->iir[ch].clippings = 0;
747     }
748
749     if (in != out)
750         av_frame_free(&in);
751
752     return ff_filter_frame(outlink, out);
753 }
754
755 static av_cold int init(AVFilterContext *ctx)
756 {
757     AudioIIRContext *s = ctx->priv;
758
759     if (!s->a_str || !s->b_str || !s->g_str) {
760         av_log(ctx, AV_LOG_ERROR, "Valid coefficients are mandatory.\n");
761         return AVERROR(EINVAL);
762     }
763
764     switch (s->precision) {
765     case 0: s->sample_format = AV_SAMPLE_FMT_DBLP; break;
766     case 1: s->sample_format = AV_SAMPLE_FMT_FLTP; break;
767     case 2: s->sample_format = AV_SAMPLE_FMT_S32P; break;
768     case 3: s->sample_format = AV_SAMPLE_FMT_S16P; break;
769     default: return AVERROR_BUG;
770     }
771
772     return 0;
773 }
774
775 static av_cold void uninit(AVFilterContext *ctx)
776 {
777     AudioIIRContext *s = ctx->priv;
778     int ch;
779
780     if (s->iir) {
781         for (ch = 0; ch < s->channels; ch++) {
782             IIRChannel *iir = &s->iir[ch];
783             av_freep(&iir->ab[0]);
784             av_freep(&iir->ab[1]);
785             av_freep(&iir->cache[0]);
786             av_freep(&iir->cache[1]);
787             av_freep(&iir->biquads);
788         }
789     }
790     av_freep(&s->iir);
791 }
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 };