avfilter/af_aiir: add analog transfer function format

1 files changed, 65 insertions, 7 deletions

diff --git a/libavfilter/af_aiir.c b/libavfilter/af_aiir.c
index f6ca2438d7..c5df4b1202 100644
--- a/libavfilter/af_aiir.c
+++ b/libavfilter/af_aiir.c
@@ -428,7 +428,7 @@ static int read_channels(AVFilterContext *ctx, int channels, uint8_t *item_str,
             return AVERROR(ENOMEM);
         }
 
-        if (s->format) {
+        if (s->format > 0) {
             ret = read_zp_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab], format[s->format]);
         } else {
             ret = read_tf_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab]);
@@ -898,6 +898,60 @@ static void convert_sp2zp(AVFilterContext *ctx, int channels)
     }
 }
 
+static double fact(double i)
+{
+    if (i <= 0.)
+        return 1.;
+    return i * fact(i - 1.);
+}
+
+static double coef_sf2zf(double *a, int N, int n, double fs)
+{
+    double z = 0.;
+
+    for (int i = 0; i <= N; i++) {
+        double acc = 0.;
+
+        for (int k = FFMAX(n - N + i, 0); k <= FFMIN(i, n); k++) {
+            acc += ((fact(i) * fact(N - i)) /
+                    (fact(k) * fact(i - k) * fact(n - k) * fact(N - i - n + k))) *
+                   ((k & 1) ? -1. : 1.);;
+        }
+
+        z += a[i] * pow(2., i) * acc;
+    }
+
+    return z;
+}
+
+static void convert_sf2tf(AVFilterContext *ctx, int channels, int sample_rate)
+{
+    AudioIIRContext *s = ctx->priv;
+    int ch;
+
+    for (ch = 0; ch < channels; ch++) {
+        IIRChannel *iir = &s->iir[ch];
+        double *temp0 = av_calloc(iir->nb_ab[0], sizeof(*temp0));
+        double *temp1 = av_calloc(iir->nb_ab[1], sizeof(*temp1));
+
+        if (!temp0 || !temp1)
+            goto next;
+
+        memcpy(temp0, iir->ab[0], iir->nb_ab[0] * sizeof(*temp0));
+        memcpy(temp1, iir->ab[1], iir->nb_ab[1] * sizeof(*temp1));
+
+        for (int n = 0; n < iir->nb_ab[0]; n++)
+            iir->ab[0][n] = coef_sf2zf(temp0, iir->nb_ab[0] - 1, n, sample_rate);
+
+        for (int n = 0; n < iir->nb_ab[1]; n++)
+            iir->ab[1][n] = coef_sf2zf(temp1, iir->nb_ab[1] - 1, n, sample_rate);
+
+next:
+        av_free(temp0);
+        av_free(temp1);
+    }
+}
+
 static void convert_pd2zp(AVFilterContext *ctx, int channels)
 {
     AudioIIRContext *s = ctx->priv;
@@ -1186,7 +1240,10 @@ static int config_output(AVFilterLink *outlink)
     if (ret < 0)
         return ret;
 
-    if (s->format == 2) {
+    if (s->format == -1) {
+        convert_sf2tf(ctx, inlink->channels, inlink->sample_rate);
+        s->format = 0;
+    } else if (s->format == 2) {
         convert_pr2zp(ctx, inlink->channels);
     } else if (s->format == 3) {
         convert_pd2zp(ctx, inlink->channels);
@@ -1207,7 +1264,7 @@ static int config_output(AVFilterLink *outlink)
     }
 
     if (s->format == 0)
-        av_log(ctx, AV_LOG_WARNING, "tf coefficients format is not recommended for too high number of zeros/poles.\n");
+        av_log(ctx, AV_LOG_WARNING, "transfer function coefficients format is not recommended for too high number of zeros/poles.\n");
 
     if (s->format > 0 && s->process == 0) {
         av_log(ctx, AV_LOG_WARNING, "Direct processsing is not recommended for zp coefficients format.\n");
@@ -1215,10 +1272,10 @@ static int config_output(AVFilterLink *outlink)
         ret = convert_zp2tf(ctx, inlink->channels);
         if (ret < 0)
             return ret;
-    } else if (s->format == 0 && s->process == 1) {
+    } else if (s->format <= 0 && s->process == 1) {
         av_log(ctx, AV_LOG_ERROR, "Serial processing is not implemented for transfer function.\n");
         return AVERROR_PATCHWELCOME;
-    } else if (s->format == 0 && s->process == 2) {
+    } else if (s->format <= 0 && s->process == 2) {
         av_log(ctx, AV_LOG_ERROR, "Parallel processing is not implemented for transfer function.\n");
         return AVERROR_PATCHWELCOME;
     } else if (s->format > 0 && s->process == 1) {
@@ -1416,8 +1473,9 @@ static const AVOption aiir_options[] = {
     { "k", "set channels gains",                   OFFSET(g_str),    AV_OPT_TYPE_STRING, {.str="1|1"}, 0, 0, AF },
     { "dry", "set dry gain",                       OFFSET(dry_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1},     0, 1, AF },
     { "wet", "set wet gain",                       OFFSET(wet_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1},     0, 1, AF },
-    { "format", "set coefficients format",         OFFSET(format),   AV_OPT_TYPE_INT,    {.i64=1},     0, 4, AF, "format" },
-    { "f", "set coefficients format",              OFFSET(format),   AV_OPT_TYPE_INT,    {.i64=1},     0, 4, AF, "format" },
+    { "format", "set coefficients format",         OFFSET(format),   AV_OPT_TYPE_INT,    {.i64=1},    -1, 4, AF, "format" },
+    { "f", "set coefficients format",              OFFSET(format),   AV_OPT_TYPE_INT,    {.i64=1},    -1, 4, AF, "format" },
+    { "sf", "analog transfer function",            0,                AV_OPT_TYPE_CONST,  {.i64=-1},    0, 0, AF, "format" },
     { "tf", "digital transfer function",           0,                AV_OPT_TYPE_CONST,  {.i64=0},     0, 0, AF, "format" },
     { "zp", "Z-plane zeros/poles",                 0,                AV_OPT_TYPE_CONST,  {.i64=1},     0, 0, AF, "format" },
     { "pr", "Z-plane zeros/poles (polar radians)", 0,                AV_OPT_TYPE_CONST,  {.i64=2},     0, 0, AF, "format" },