audiofx: Use most common convention for definitions of IIR filter coefficients.

Most signal processing texts, including MATLAB, use the following convention for IIR filter coefficients:

a_0 y[n] + a_1 y[n-1] + ... + a_M y[n-M] = b_0 x[n] + b_1 x[n-1] + ... + b[N] x[n-N]

Usually, a_0 is set to 1 because the coefficients can always be rescaled, giving

y[n] = b_0 x[n] + b_1 x[n-1] + ... + b[N] x[n-N] - a_1 y[n-1] - ... - a_M y[n-M]

The convention that was previously used by audiofxbaseiirfilter and derived class had the a and b coefficients swapped, and did not have the minus signs.

This change makes the audiofx plugin use the more common convention described above.

5 files changed, 93 insertions, 88 deletions

diff --git a/gst/audiofx/audiochebband.c b/gst/audiofx/audiochebband.c
index cf431618e..c37cdbb37 100644
--- a/gst/audiofx/audiochebband.c
+++ b/gst/audiofx/audiochebband.c
@@ -209,8 +209,8 @@ gst_audio_cheb_band_init (GstAudioChebBand * filter)
 
 static void
 generate_biquad_coefficients (GstAudioChebBand * filter,
-    gint p, gdouble * a0, gdouble * a1, gdouble * a2, gdouble * a3,
-    gdouble * a4, gdouble * b1, gdouble * b2, gdouble * b3, gdouble * b4)
+    gint p, gdouble * b0, gdouble * b1, gdouble * b2, gdouble * b3,
+    gdouble * b4, gdouble * a1, gdouble * a2, gdouble * a3, gdouble * a4)
 {
   gint np = filter->poles / 2;
   gdouble ripple = filter->ripple;
@@ -349,19 +349,19 @@ generate_biquad_coefficients (GstAudioChebBand * filter,
 
       d = 1.0 + beta * (y1 - beta * y2);
 
-      *a0 = (x0 + beta * (-x1 + beta * x2)) / d;
-      *a1 = (alpha * (-2.0 * x0 + x1 + beta * x1 - 2.0 * beta * x2)) / d;
-      *a2 =
+      *b0 = (x0 + beta * (-x1 + beta * x2)) / d;
+      *b1 = (alpha * (-2.0 * x0 + x1 + beta * x1 - 2.0 * beta * x2)) / d;
+      *b2 =
           (-x1 - beta * beta * x1 + 2.0 * beta * (x0 + x2) +
           alpha * alpha * (x0 - x1 + x2)) / d;
-      *a3 = (alpha * (x1 + beta * (-2.0 * x0 + x1) - 2.0 * x2)) / d;
-      *a4 = (beta * (beta * x0 - x1) + x2) / d;
-      *b1 = (alpha * (2.0 + y1 + beta * y1 - 2.0 * beta * y2)) / d;
-      *b2 =
+      *b3 = (alpha * (x1 + beta * (-2.0 * x0 + x1) - 2.0 * x2)) / d;
+      *b4 = (beta * (beta * x0 - x1) + x2) / d;
+      *a1 = (alpha * (2.0 + y1 + beta * y1 - 2.0 * beta * y2)) / d;
+      *a2 =
           (-y1 - beta * beta * y1 - alpha * alpha * (1.0 + y1 - y2) +
           2.0 * beta * (-1.0 + y2)) / d;
-      *b3 = (alpha * (y1 + beta * (2.0 + y1) - 2.0 * y2)) / d;
-      *b4 = (-beta * beta - beta * y1 + y2) / d;
+      *a3 = (alpha * (y1 + beta * (2.0 + y1) - 2.0 * y2)) / d;
+      *a4 = (-beta * beta - beta * y1 + y2) / d;
     } else {
       a = cos ((w1 + w0) / 2.0) / cos ((w1 - w0) / 2.0);
       b = tan (1.0 / 2.0) * tan ((w1 - w0) / 2.0);
@@ -371,19 +371,19 @@ generate_biquad_coefficients (GstAudioChebBand * filter,
 
       d = -1.0 + beta * (beta * y2 + y1);
 
-      *a0 = (-x0 - beta * x1 - beta * beta * x2) / d;
-      *a1 = (alpha * (2.0 * x0 + x1 + beta * x1 + 2.0 * beta * x2)) / d;
-      *a2 =
+      *b0 = (-x0 - beta * x1 - beta * beta * x2) / d;
+      *b1 = (alpha * (2.0 * x0 + x1 + beta * x1 + 2.0 * beta * x2)) / d;
+      *b2 =
           (-x1 - beta * beta * x1 - 2.0 * beta * (x0 + x2) -
           alpha * alpha * (x0 + x1 + x2)) / d;
-      *a3 = (alpha * (x1 + beta * (2.0 * x0 + x1) + 2.0 * x2)) / d;
-      *a4 = (-beta * beta * x0 - beta * x1 - x2) / d;
-      *b1 = (alpha * (-2.0 + y1 + beta * y1 + 2.0 * beta * y2)) / d;
-      *b2 =
+      *b3 = (alpha * (x1 + beta * (2.0 * x0 + x1) + 2.0 * x2)) / d;
+      *b4 = (-beta * beta * x0 - beta * x1 - x2) / d;
+      *a1 = (alpha * (-2.0 + y1 + beta * y1 + 2.0 * beta * y2)) / d;
+      *a2 =
           -(y1 + beta * beta * y1 + 2.0 * beta * (-1.0 + y2) +
           alpha * alpha * (-1.0 + y1 + y2)) / d;
-      *b3 = (alpha * (beta * (-2.0 + y1) + y1 + 2.0 * y2)) / d;
-      *b4 = -(-beta * beta + beta * y1 + y2) / d;
+      *a3 = (alpha * (beta * (-2.0 + y1) + y1 + 2.0 * y2)) / d;
+      *a4 = -(-beta * beta + beta * y1 + y2) / d;
     }
   }
 }
@@ -395,20 +395,24 @@ generate_coefficients (GstAudioChebBand * filter)
 
   if (rate == 0) {
     gdouble *a = g_new0 (gdouble, 1);
+    gdouble *b = g_new0 (gdouble, 1);
 
     a[0] = 1.0;
+    b[0] = 1.0;
     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
-        (filter), a, 1, NULL, 0);
+        (filter), a, 1, b, 1);
     GST_LOG_OBJECT (filter, "rate was not set yet");
     return;
   }
 
   if (filter->upper_frequency <= filter->lower_frequency) {
     gdouble *a = g_new0 (gdouble, 1);
+    gdouble *b = g_new0 (gdouble, 1);
 
-    a[0] = (filter->mode == MODE_BAND_PASS) ? 0.0 : 1.0;
+    a[0] = 1.0;
+    b[0] = (filter->mode == MODE_BAND_PASS) ? 0.0 : 1.0;
     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
-        (filter), a, 1, NULL, 0);
+        (filter), a, 1, b, 1);
 
     GST_LOG_OBJECT (filter, "frequency band had no or negative dimension");
     return;
@@ -438,12 +442,12 @@ generate_coefficients (GstAudioChebBand * filter)
     b[4] = 1.0;
 
     for (p = 1; p <= np / 4; p++) {
-      gdouble a0, a1, a2, a3, a4, b1, b2, b3, b4;
+      gdouble b0, b1, b2, b3, b4, a1, a2, a3, a4;
       gdouble *ta = g_new0 (gdouble, np + 5);
       gdouble *tb = g_new0 (gdouble, np + 5);
 
-      generate_biquad_coefficients (filter, p, &a0, &a1, &a2, &a3, &a4, &b1,
-          &b2, &b3, &b4);
+      generate_biquad_coefficients (filter, p, &b0, &b1, &b2, &b3, &b4, &a1,
+          &a2, &a3, &a4);
 
       memcpy (ta, a, sizeof (gdouble) * (np + 5));
       memcpy (tb, b, sizeof (gdouble) * (np + 5));
@@ -452,25 +456,23 @@ generate_coefficients (GstAudioChebBand * filter)
        * to the cascade by multiplication of the transfer
        * functions */
       for (i = 4; i < np + 5; i++) {
-        a[i] =
-            a0 * ta[i] + a1 * ta[i - 1] + a2 * ta[i - 2] + a3 * ta[i - 3] +
-            a4 * ta[i - 4];
         b[i] =
-            tb[i] - b1 * tb[i - 1] - b2 * tb[i - 2] - b3 * tb[i - 3] -
+            b0 * tb[i] + b1 * tb[i - 1] + b2 * tb[i - 2] + b3 * tb[i - 3] +
             b4 * tb[i - 4];
+        a[i] =
+            ta[i] - a1 * ta[i - 1] - a2 * ta[i - 2] - a3 * ta[i - 3] -
+            a4 * ta[i - 4];
       }
       g_free (ta);
       g_free (tb);
     }
 
-    /* Move coefficients to the beginning of the array
-     * and multiply the b coefficients with -1 to move from
+    /* Move coefficients to the beginning of the array to move from
      * the transfer function's coefficients to the difference
      * equation's coefficients */
-    b[4] = 0.0;
     for (i = 0; i <= np; i++) {
       a[i] = a[i + 4];
-      b[i] = -b[i + 4];
+      b[i] = b[i + 4];
     }
 
     /* Normalize to unity gain at frequency 0 and frequency
@@ -489,7 +491,7 @@ generate_coefficients (GstAudioChebBand * filter)
       gain1 = sqrt (gain1 * gain2);
 
       for (i = 0; i <= np; i++) {
-        a[i] /= gain1;
+        b[i] /= gain1;
       }
     } else {
       /* gain is H(wc), wc = center frequency */
@@ -503,7 +505,7 @@ generate_coefficients (GstAudioChebBand * filter)
           zi);
 
       for (i = 0; i <= np; i++) {
-        a[i] /= gain;
+        b[i] /= gain;
       }
     }
 
diff --git a/gst/audiofx/audiocheblimit.c b/gst/audiofx/audiocheblimit.c
index 9654a0564..894152fe9 100644
--- a/gst/audiofx/audiocheblimit.c
+++ b/gst/audiofx/audiocheblimit.c
@@ -201,8 +201,8 @@ gst_audio_cheb_limit_init (GstAudioChebLimit * filter)
 
 static void
 generate_biquad_coefficients (GstAudioChebLimit * filter,
-    gint p, gdouble * a0, gdouble * a1, gdouble * a2,
-    gdouble * b1, gdouble * b2)
+    gint p, gdouble * b0, gdouble * b1, gdouble * b2,
+    gdouble * a1, gdouble * a2)
 {
   gint np = filter->poles;
   gdouble ripple = filter->ripple;
@@ -329,11 +329,11 @@ generate_biquad_coefficients (GstAudioChebLimit * filter,
       k = -cos ((omega + 1.0) / 2.0) / cos ((omega - 1.0) / 2.0);
 
     d = 1.0 + y1 * k - y2 * k * k;
-    *a0 = (x0 + k * (-x1 + k * x2)) / d;
-    *a1 = (x1 + k * k * x1 - 2.0 * k * (x0 + x2)) / d;
-    *a2 = (x0 * k * k - x1 * k + x2) / d;
-    *b1 = (2.0 * k + y1 + y1 * k * k - 2.0 * y2 * k) / d;
-    *b2 = (-k * k - y1 * k + y2) / d;
+    *b0 = (x0 + k * (-x1 + k * x2)) / d;
+    *b1 = (x1 + k * k * x1 - 2.0 * k * (x0 + x2)) / d;
+    *b2 = (x0 * k * k - x1 * k + x2) / d;
+    *a1 = (2.0 * k + y1 + y1 * k * k - 2.0 * y2 * k) / d;
+    *a2 = (-k * k - y1 * k + y2) / d;
 
     if (filter->mode == MODE_HIGH_PASS) {
       *a1 = -*a1;
@@ -347,10 +347,12 @@ generate_coefficients (GstAudioChebLimit * filter)
 {
   if (GST_AUDIO_FILTER_RATE (filter) == 0) {
     gdouble *a = g_new0 (gdouble, 1);
+    gdouble *b = g_new0 (gdouble, 1);
 
     a[0] = 1.0;
+    b[0] = 1.0;
     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
-        (filter), a, 1, NULL, 0);
+        (filter), a, 1, b, 1);
 
     GST_LOG_OBJECT (filter, "rate was not set yet");
     return;
@@ -358,18 +360,22 @@ generate_coefficients (GstAudioChebLimit * filter)
 
   if (filter->cutoff >= GST_AUDIO_FILTER_RATE (filter) / 2.0) {
     gdouble *a = g_new0 (gdouble, 1);
+    gdouble *b = g_new0 (gdouble, 1);
 
-    a[0] = (filter->mode == MODE_LOW_PASS) ? 1.0 : 0.0;
+    a[0] = 1.0;
+    b[0] = (filter->mode == MODE_LOW_PASS) ? 1.0 : 0.0;
     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
-        (filter), a, 1, NULL, 0);
+        (filter), a, 1, b, 1);
     GST_LOG_OBJECT (filter, "cutoff was higher than nyquist frequency");
     return;
   } else if (filter->cutoff <= 0.0) {
     gdouble *a = g_new0 (gdouble, 1);
+    gdouble *b = g_new0 (gdouble, 1);
 
-    a[0] = (filter->mode == MODE_LOW_PASS) ? 0.0 : 1.0;
+    a[0] = 1.0;
+    b[0] = (filter->mode == MODE_LOW_PASS) ? 0.0 : 1.0;
     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
-        (filter), a, 1, NULL, 0);
+        (filter), a, 1, b, 1);
     GST_LOG_OBJECT (filter, "cutoff is lower than zero");
     return;
   }
@@ -388,11 +394,11 @@ generate_coefficients (GstAudioChebLimit * filter)
     b[2] = 1.0;
 
     for (p = 1; p <= np / 2; p++) {
-      gdouble a0, a1, a2, b1, b2;
+      gdouble b0, b1, b2, a1, a2;
       gdouble *ta = g_new0 (gdouble, np + 3);
       gdouble *tb = g_new0 (gdouble, np + 3);
 
-      generate_biquad_coefficients (filter, p, &a0, &a1, &a2, &b1, &b2);
+      generate_biquad_coefficients (filter, p, &b0, &b1, &b2, &a1, &a2);
 
       memcpy (ta, a, sizeof (gdouble) * (np + 3));
       memcpy (tb, b, sizeof (gdouble) * (np + 3));
@@ -401,21 +407,19 @@ generate_coefficients (GstAudioChebLimit * filter)
        * to the cascade by multiplication of the transfer
        * functions */
       for (i = 2; i < np + 3; i++) {
-        a[i] = a0 * ta[i] + a1 * ta[i - 1] + a2 * ta[i - 2];
-        b[i] = tb[i] - b1 * tb[i - 1] - b2 * tb[i - 2];
+        b[i] = b0 * tb[i] + b1 * tb[i - 1] + b2 * tb[i - 2];
+        a[i] = ta[i] - a1 * ta[i - 1] - a2 * ta[i - 2];
       }
       g_free (ta);
       g_free (tb);
     }
 
-    /* Move coefficients to the beginning of the array
-     * and multiply the b coefficients with -1 to move from
+    /* Move coefficients to the beginning of the array to move from
      * the transfer function's coefficients to the difference
      * equation's coefficients */
-    b[2] = 0.0;
     for (i = 0; i <= np; i++) {
       a[i] = a[i + 2];
-      b[i] = -b[i + 2];
+      b[i] = b[i + 2];
     }
 
     /* Normalize to unity gain at frequency 0 for lowpass
@@ -433,7 +437,7 @@ generate_coefficients (GstAudioChebLimit * filter)
             -1.0, 0.0);
 
       for (i = 0; i <= np; i++) {
-        a[i] /= gain;
+        b[i] /= gain;
       }
     }
 
diff --git a/gst/audiofx/audiofxbaseiirfilter.c b/gst/audiofx/audiofxbaseiirfilter.c
index 2d318829e..e30c586e2 100644
--- a/gst/audiofx/audiofxbaseiirfilter.c
+++ b/gst/audiofx/audiofxbaseiirfilter.c
@@ -132,7 +132,7 @@ gst_audio_fx_base_iir_filter_init (GstAudioFXBaseIIRFilter * filter)
 }
 
 /* Evaluate the transfer function that corresponds to the IIR
- * coefficients at zr + zi*I and return the magnitude */
+ * coefficients at (zr + zi*I)^-1 and return the magnitude */
 gdouble
 gst_audio_fx_base_iir_filter_calculate_gain (gdouble * a, guint na, gdouble * b,
     guint nb, gdouble zr, gdouble zi)
@@ -146,9 +146,9 @@ gst_audio_fx_base_iir_filter_calculate_gain (gdouble * a, guint na, gdouble * b,
 
   gint i;
 
-  sum_ar = 0.0;
+  sum_ar = a[na - 1];
   sum_ai = 0.0;
-  for (i = na - 1; i >= 0; i--) {
+  for (i = na - 2; i >= 0; i--) {
     sum_r_old = sum_ar;
     sum_i_old = sum_ai;
 
@@ -156,22 +156,20 @@ gst_audio_fx_base_iir_filter_calculate_gain (gdouble * a, guint na, gdouble * b,
     sum_ai = (sum_r_old * zi + sum_i_old * zr) + 0.0;
   }
 
-  sum_br = 0.0;
+  sum_br = b[nb - 1];
   sum_bi = 0.0;
-  for (i = nb - 1; i >= 0; i--) {
+  for (i = nb - 2; i >= 0; i--) {
     sum_r_old = sum_br;
     sum_i_old = sum_bi;
 
-    sum_br = (sum_r_old * zr - sum_i_old * zi) - b[i];
-    sum_bi = (sum_r_old * zi + sum_i_old * zr) - 0.0;
+    sum_br = (sum_r_old * zr - sum_i_old * zi) + b[i];
+    sum_bi = (sum_r_old * zi + sum_i_old * zr) + 0.0;
   }
-  sum_br += 1.0;
-  sum_bi += 0.0;
 
   gain_r =
-      (sum_ar * sum_br + sum_ai * sum_bi) / (sum_br * sum_br + sum_bi * sum_bi);
+      (sum_br * sum_ar + sum_bi * sum_ai) / (sum_ar * sum_ar + sum_ai * sum_ai);
   gain_i =
-      (sum_ai * sum_br - sum_ar * sum_bi) / (sum_br * sum_br + sum_bi * sum_bi);
+      (sum_bi * sum_ar - sum_br * sum_ai) / (sum_ar * sum_ar + sum_ai * sum_ai);
 
   return (sqrt (gain_r * gain_r + gain_i * gain_i));
 }
@@ -201,12 +199,12 @@ gst_audio_fx_base_iir_filter_set_coefficients (GstAudioFXBaseIIRFilter * filter,
       if (free)
         g_free (ctx->x);
       else
-        memset (ctx->x, 0, filter->na * sizeof (gdouble));
+        memset (ctx->x, 0, filter->nb * sizeof (gdouble));
 
       if (free)
         g_free (ctx->y);
       else
-        memset (ctx->y, 0, filter->nb * sizeof (gdouble));
+        memset (ctx->y, 0, filter->na * sizeof (gdouble));
     }
 
     g_free (filter->channels);
@@ -227,8 +225,8 @@ gst_audio_fx_base_iir_filter_set_coefficients (GstAudioFXBaseIIRFilter * filter,
     for (i = 0; i < filter->nchannels; i++) {
       ctx = &filter->channels[i];
 
-      ctx->x = g_new0 (gdouble, filter->na);
-      ctx->y = g_new0 (gdouble, filter->nb);
+      ctx->x = g_new0 (gdouble, filter->nb);
+      ctx->y = g_new0 (gdouble, filter->na);
     }
   }
 
@@ -279,8 +277,8 @@ gst_audio_fx_base_iir_filter_setup (GstAudioFilter * base,
     for (i = 0; i < channels; i++) {
       ctx = &filter->channels[i];
 
-      ctx->x = g_new0 (gdouble, filter->na);
-      ctx->y = g_new0 (gdouble, filter->nb);
+      ctx->x = g_new0 (gdouble, filter->nb);
+      ctx->y = g_new0 (gdouble, filter->na);
     }
     filter->nchannels = channels;
   }
@@ -292,32 +290,33 @@ static inline gdouble
 process (GstAudioFXBaseIIRFilter * filter,
     GstAudioFXBaseIIRFilterChannelCtx * ctx, gdouble x0)
 {
-  gdouble val = filter->a[0] * x0;
+  gdouble val = filter->b[0] * x0;
   gint i, j;
 
-  for (i = 1, j = ctx->x_pos; i < filter->na; i++) {
-    val += filter->a[i] * ctx->x[j];
+  for (i = 1, j = ctx->x_pos; i < filter->nb; i++) {
+    val += filter->b[i] * ctx->x[j];
     j--;
     if (j < 0)
-      j = filter->na - 1;
+      j = filter->nb - 1;
   }
 
-  for (i = 1, j = ctx->y_pos; i < filter->nb; i++) {
-    val += filter->b[i] * ctx->y[j];
+  for (i = 1, j = ctx->y_pos; i < filter->na; i++) {
+    val -= filter->a[i] * ctx->y[j];
     j--;
     if (j < 0)
-      j = filter->nb - 1;
+      j = filter->na - 1;
   }
+  val /= filter->a[0];
 
   if (ctx->x) {
     ctx->x_pos++;
-    if (ctx->x_pos >= filter->na)
+    if (ctx->x_pos >= filter->nb)
       ctx->x_pos = 0;
     ctx->x[ctx->x_pos] = x0;
   }
   if (ctx->y) {
     ctx->y_pos++;
-    if (ctx->y_pos >= filter->nb)
+    if (ctx->y_pos >= filter->na)
       ctx->y_pos = 0;
 
     ctx->y[ctx->y_pos] = val;
diff --git a/gst/audiofx/audioiirfilter.c b/gst/audiofx/audioiirfilter.c
index dc3c17efb..09d02829f 100644
--- a/gst/audiofx/audioiirfilter.c
+++ b/gst/audiofx/audioiirfilter.c
@@ -101,14 +101,14 @@ gst_audio_iir_filter_class_init (GstAudioIIRFilterClass * klass)
 
   g_object_class_install_property (gobject_class, PROP_A,
       g_param_spec_value_array ("a", "A",
-          "Filter coefficients (numerator of transfer function)",
+          "Filter coefficients (denominator of transfer function)",
           g_param_spec_double ("Coefficient", "Filter Coefficient",
               "Filter coefficient", -G_MAXDOUBLE, G_MAXDOUBLE, 0.0,
               G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS),
           G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS));
   g_object_class_install_property (gobject_class, PROP_B,
       g_param_spec_value_array ("b", "B",
-          "Filter coefficients (denominator of transfer function)",
+          "Filter coefficients (numerator of transfer function)",
           g_param_spec_double ("Coefficient", "Filter Coefficient",
               "Filter coefficient", -G_MAXDOUBLE, G_MAXDOUBLE, 0.0,
               G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS),
diff --git a/tests/check/elements/audioiirfilter.c b/tests/check/elements/audioiirfilter.c
index f9d7d1b14..7965b2099 100644
--- a/tests/check/elements/audioiirfilter.c
+++ b/tests/check/elements/audioiirfilter.c
@@ -76,17 +76,17 @@ on_rate_changed (GstElement * element, gint rate, gpointer user_data)
   g_value_array_append (va, &v);
   g_value_reset (&v);
 
-  g_object_set (G_OBJECT (element), "a", va, NULL);
+  g_object_set (G_OBJECT (element), "b", va, NULL);
 
   g_value_array_free (va);
 
   va = g_value_array_new (6);
 
-  g_value_set_double (&v, 0.0);
+  g_value_set_double (&v, 1.0);
   g_value_array_append (va, &v);
   g_value_reset (&v);
 
-  g_object_set (G_OBJECT (element), "b", va, NULL);
+  g_object_set (G_OBJECT (element), "a", va, NULL);
 
   g_value_array_free (va);
 }