avfilter/af_aiir: add analog transfer function format

doc/filters.texi

@@ -1408,6 +1408,8 @@ Set output gain.
 Set coefficients format.
 
 @table @samp
+@item sf
+analog transfer function
 @item tf
 digital transfer function
 @item zp
@@ -1466,7 +1468,7 @@ displayed. This option is used only when @var{response} is enabled.
 Set video stream size. This option is used only when @var{response} is enabled.
 @end table
 
-Coefficients in @code{tf} format are separated by spaces and are in ascending
+Coefficients in @code{tf} and @code{sf} format are separated by spaces and are in ascending
 order.
 
 Coefficients in @code{zp} format are separated by spaces and order of coefficients
@@ -1491,6 +1493,12 @@ Same as above but in @code{zp} format:
 @example
 aiir=k=0.79575848078096756:z=0.80918701+0.58773007i 0.80918701-0.58773007i 0.80884700+0.58784055i 0.80884700-0.58784055i:p=0.63892345+0.59951235i 0.63892345-0.59951235i 0.79582691+0.44198673i 0.79582691-0.44198673i:f=zp:r=s
 @end example
+
+@item
+Apply 3-rd order analog normalized Butterworth low-pass filter, using analog transfer function format:
+@example
+aiir=z=1.3057 0 0 0:p=1.3057 2.3892 2.1860 1:f=sf:r=d
+@end example
 @end itemize
 
 @section alimiter

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" },
+    { "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},     0, 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" },