gst/audiofx/audiochebband.c

   1 /*
   2  * GStreamer
   3  * Copyright (C) 2007 Sebastian Dröge <slomo@circular-chaos.org>
   4  *
   5  * This library is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU Library General Public
   7  * License as published by the Free Software Foundation; either
   8  * version 2 of the License, or (at your option) any later version.
   9  *
  10  * This library is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * Library General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU Library General Public
  16  * License along with this library; if not, write to the
  17  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  18  * Boston, MA 02111-1307, USA.
  19  */
  20
  21 /*
  22  * Chebyshev type 1 filter design based on
  23  * "The Scientist and Engineer's Guide to DSP", Chapter 20.
  24  * http://www.dspguide.com/
  25  *
  26  * For type 2 and Chebyshev filters in general read
  27  * http://en.wikipedia.org/wiki/Chebyshev_filter
  28  *
  29  * Transformation from lowpass to bandpass/bandreject:
  30  * http://docs.dewresearch.com/DspHelp/html/IDH_LinearSystems_LowpassToBandPassZ.htm
  31  * http://docs.dewresearch.com/DspHelp/html/IDH_LinearSystems_LowpassToBandStopZ.htm
  32  *
  33  */
  34
  35 /**
  36  * SECTION:element-audiochebyshevfreqband
  37  * @short_description: Chebyshev band pass and band reject filter
  38  *
  39  * <refsect2>
  40  * <para>
  41  * Attenuates all frequencies outside (bandpass) or inside (bandreject) of a frequency
  42  * band. The number of poles and the ripple parameter control the rolloff.
  43  * </para>
  44  * <para>
  45  * This element has the advantage over the windowed sinc bandpass and bandreject filter that it is
  46  * much faster and produces almost as good results. It's only disadvantages are the highly
  47  * non-linear phase and the slower rolloff compared to a windowed sinc filter with a large kernel.
  48  * </para>
  49  * <para>
  50  * For type 1 the ripple parameter specifies how much ripple in dB is allowed in the passband, i.e.
  51  * some frequencies in the passband will be amplified by that value. A higher ripple value will allow
  52  * a faster rolloff.
  53  * </para>
  54  * <para>
  55  * For type 2 the ripple parameter specifies the stopband attenuation. In the stopband the gain will
  56  * be at most this value. A lower ripple value will allow a faster rolloff.
  57  * </para>
  58  * <para>
  59  * As a special case, a Chebyshev type 1 filter with no ripple is a Butterworth filter.
  60  * </para>
  61  * <para><note>
  62  * Be warned that a too large number of poles can produce noise. The most poles are possible with
  63  * a cutoff frequency at a quarter of the sampling rate.
  64  * </note></para>
  65  * <title>Example launch line</title>
  66  * <para>
  67  * <programlisting>
  68  * gst-launch audiotestsrc freq=1500 ! audioconvert ! audiochebyshevfreqband mode=band-pass lower-frequency=1000 upper-frequenc=6000 poles=4 ! audioconvert ! alsasink
  69  * gst-launch filesrc location="melo1.ogg" ! oggdemux ! vorbisdec ! audioconvert ! audiochebyshevfreqband mode=band-reject lower-frequency=1000 upper-frequency=4000 ripple=0.2 ! audioconvert ! alsasink
  70  * gst-launch audiotestsrc wave=white-noise ! audioconvert ! audiochebyshevfreqband mode=band-pass lower-frequency=1000 upper-frequency=4000 type=2 ! audioconvert ! alsasink
  71  * </programlisting>
  72  * </para>
  73  * </refsect2>
  74  */
  75
  76 #ifdef HAVE_CONFIG_H
  77 #include "config.h"
  78 #endif
  79
  80 #include <gst/gst.h>
  81 #include <gst/base/gstbasetransform.h>
  82 #include <gst/audio/audio.h>
  83 #include <gst/audio/gstaudiofilter.h>
  84 #include <gst/controller/gstcontroller.h>
  85
  86 #include <math.h>
  87
  88 #include "audiochebyshevfreqband.h"
  89
  90 #define GST_CAT_DEFAULT gst_audio_chebyshev_freq_band_debug
  91 GST_DEBUG_CATEGORY_STATIC (GST_CAT_DEFAULT);
  92
  93 static const GstElementDetails element_details =
  94 GST_ELEMENT_DETAILS ("AudioChebyshevFreqBand",
  95     "Filter/Effect/Audio",
  96     "Chebyshev band pass and band reject filter",
  97     "Sebastian Dröge <slomo@circular-chaos.org>");
  98
  99 /* Filter signals and args */
 100 enum
 101 {
 102   /* FILL ME */
 103   LAST_SIGNAL
 104 };
 105
 106 enum
 107 {
 108   PROP_0,
 109   PROP_MODE,
 110   PROP_TYPE,
 111   PROP_LOWER_FREQUENCY,
 112   PROP_UPPER_FREQUENCY,
 113   PROP_RIPPLE,
 114   PROP_POLES
 115 };
 116
 117 #define ALLOWED_CAPS \
 118     "audio/x-raw-float,"                                              \
 119     " width = (int) { 32, 64 }, "                                     \
 120     " endianness = (int) BYTE_ORDER,"                                 \
 121     " rate = (int) [ 1, MAX ],"                                       \
 122     " channels = (int) [ 1, MAX ]"
 123
 124 #define DEBUG_INIT(bla) \
 125   GST_DEBUG_CATEGORY_INIT (gst_audio_chebyshev_freq_band_debug, "audiochebyshevfreqband", 0, "audiochebyshevfreqband element");
 126
 127 GST_BOILERPLATE_FULL (GstAudioChebyshevFreqBand, gst_audio_chebyshev_freq_band,
 128     GstAudioFilter, GST_TYPE_AUDIO_FILTER, DEBUG_INIT);
 129
 130 static void gst_audio_chebyshev_freq_band_set_property (GObject * object,
 131     guint prop_id, const GValue * value, GParamSpec * pspec);
 132 static void gst_audio_chebyshev_freq_band_get_property (GObject * object,
 133     guint prop_id, GValue * value, GParamSpec * pspec);
 134
 135 static gboolean gst_audio_chebyshev_freq_band_setup (GstAudioFilter * filter,
 136     GstRingBufferSpec * format);
 137 static GstFlowReturn
 138 gst_audio_chebyshev_freq_band_transform_ip (GstBaseTransform * base,
 139     GstBuffer * buf);
 140 static gboolean gst_audio_chebyshev_freq_band_start (GstBaseTransform * base);
 141
 142 static void process_64 (GstAudioChebyshevFreqBand * filter,
 143     gdouble * data, guint num_samples);
 144 static void process_32 (GstAudioChebyshevFreqBand * filter,
 145     gfloat * data, guint num_samples);
 146
 147 enum
 148 {
 149   MODE_BAND_PASS = 0,
 150   MODE_BAND_REJECT
 151 };
 152
 153 #define GST_TYPE_AUDIO_CHEBYSHEV_FREQ_BAND_MODE (gst_audio_chebyshev_freq_band_mode_get_type ())
 154 static GType
 155 gst_audio_chebyshev_freq_band_mode_get_type (void)
 156 {
 157   static GType gtype = 0;
 158
 159   if (gtype == 0) {
 160     static const GEnumValue values[] = {
 161       {MODE_BAND_PASS, "Band pass (default)",
 162           "band-pass"},
 163       {MODE_BAND_REJECT, "Band reject",
 164           "band-reject"},
 165       {0, NULL, NULL}
 166     };
 167
 168     gtype = g_enum_register_static ("GstAudioChebyshevFreqBandMode", values);
 169   }
 170   return gtype;
 171 }
 172
 173 /* GObject vmethod implementations */
 174
 175 static void
 176 gst_audio_chebyshev_freq_band_base_init (gpointer klass)
 177 {
 178   GstElementClass *element_class = GST_ELEMENT_CLASS (klass);
 179   GstCaps *caps;
 180
 181   gst_element_class_set_details (element_class, &element_details);
 182
 183   caps = gst_caps_from_string (ALLOWED_CAPS);
 184   gst_audio_filter_class_add_pad_templates (GST_AUDIO_FILTER_CLASS (klass),
 185       caps);
 186   gst_caps_unref (caps);
 187 }
 188
 189 static void
 190 gst_audio_chebyshev_freq_band_dispose (GObject * object)
 191 {
 192   GstAudioChebyshevFreqBand *filter = GST_AUDIO_CHEBYSHEV_FREQ_BAND (object);
 193
 194   if (filter->a) {
 195     g_free (filter->a);
 196     filter->a = NULL;
 197   }
 198
 199   if (filter->b) {
 200     g_free (filter->b);
 201     filter->b = NULL;
 202   }
 203
 204   if (filter->channels) {
 205     GstAudioChebyshevFreqBandChannelCtx *ctx;
 206     gint i, channels = GST_AUDIO_FILTER (filter)->format.channels;
 207
 208     for (i = 0; i < channels; i++) {
 209       ctx = &filter->channels[i];
 210       g_free (ctx->x);
 211       g_free (ctx->y);
 212     }
 213
 214     g_free (filter->channels);
 215     filter->channels = NULL;
 216   }
 217
 218   G_OBJECT_CLASS (parent_class)->dispose (object);
 219 }
 220
 221 static void
 222 gst_audio_chebyshev_freq_band_class_init (GstAudioChebyshevFreqBandClass *
 223     klass)
 224 {
 225   GObjectClass *gobject_class;
 226   GstBaseTransformClass *trans_class;
 227   GstAudioFilterClass *filter_class;
 228
 229   gobject_class = (GObjectClass *) klass;
 230   trans_class = (GstBaseTransformClass *) klass;
 231   filter_class = (GstAudioFilterClass *) klass;
 232
 233   gobject_class->set_property = gst_audio_chebyshev_freq_band_set_property;
 234   gobject_class->get_property = gst_audio_chebyshev_freq_band_get_property;
 235   gobject_class->dispose = gst_audio_chebyshev_freq_band_dispose;
 236
 237   g_object_class_install_property (gobject_class, PROP_MODE,
 238       g_param_spec_enum ("mode", "Mode",
 239           "Low pass or high pass mode", GST_TYPE_AUDIO_CHEBYSHEV_FREQ_BAND_MODE,
 240           MODE_BAND_PASS, G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE));
 241   g_object_class_install_property (gobject_class, PROP_TYPE,
 242       g_param_spec_int ("type", "Type",
 243           "Type of the chebychev filter", 1, 2,
 244           1, G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE));
 245
 246   /* FIXME: Don't use the complete possible range but restrict the upper boundary
 247    * so automatically generated UIs can use a slider without */
 248   g_object_class_install_property (gobject_class, PROP_LOWER_FREQUENCY,
 249       g_param_spec_float ("lower-frequency", "Lower frequency",
 250           "Start frequency of the band (Hz)", 0.0, 100000.0,
 251           0.0, G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE));
 252   g_object_class_install_property (gobject_class, PROP_UPPER_FREQUENCY,
 253       g_param_spec_float ("upper-frequency", "Upper frequency",
 254           "Stop frequency of the band (Hz)", 0.0, 100000.0,
 255           0.0, G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE));
 256   g_object_class_install_property (gobject_class, PROP_RIPPLE,
 257       g_param_spec_float ("ripple", "Ripple",
 258           "Amount of ripple (dB)", 0.0, 200.0,
 259           0.25, G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE));
 260   /* FIXME: What to do about this upper boundary? With a frequencies near
 261    * rate/4 32 poles are completely possible, with frequencies very low
 262    * or very high 16 poles already produces only noise */
 263   g_object_class_install_property (gobject_class, PROP_POLES,
 264       g_param_spec_int ("poles", "Poles",
 265           "Number of poles to use, will be rounded up to the next multiply of four",
 266           4, 32, 4, G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE));
 267
 268   filter_class->setup = GST_DEBUG_FUNCPTR (gst_audio_chebyshev_freq_band_setup);
 269   trans_class->transform_ip =
 270       GST_DEBUG_FUNCPTR (gst_audio_chebyshev_freq_band_transform_ip);
 271   trans_class->start = GST_DEBUG_FUNCPTR (gst_audio_chebyshev_freq_band_start);
 272 }
 273
 274 static void
 275 gst_audio_chebyshev_freq_band_init (GstAudioChebyshevFreqBand * filter,
 276     GstAudioChebyshevFreqBandClass * klass)
 277 {
 278   filter->lower_frequency = filter->upper_frequency = 0.0;
 279   filter->mode = MODE_BAND_PASS;
 280   filter->type = 1;
 281   filter->poles = 4;
 282   filter->ripple = 0.25;
 283   gst_base_transform_set_in_place (GST_BASE_TRANSFORM (filter), TRUE);
 284
 285   filter->have_coeffs = FALSE;
 286   filter->num_a = 0;
 287   filter->num_b = 0;
 288   filter->channels = NULL;
 289 }
 290
 291 static void
 292 generate_biquad_coefficients (GstAudioChebyshevFreqBand * filter,
 293     gint p, gdouble * a0, gdouble * a1, gdouble * a2, gdouble * a3,
 294     gdouble * a4, gdouble * b1, gdouble * b2, gdouble * b3, gdouble * b4)
 295 {
 296   gint np = filter->poles / 2;
 297   gdouble ripple = filter->ripple;
 298
 299   /* pole location in s-plane */
 300   gdouble rp, ip;
 301
 302   /* zero location in s-plane */
 303   gdouble rz = 0.0, iz = 0.0;
 304
 305   /* transfer function coefficients for the z-plane */
 306   gdouble x0, x1, x2, y1, y2;
 307   gint type = filter->type;
 308
 309   /* Calculate pole location for lowpass at frequency 1 */
 310   {
 311     gdouble angle = (M_PI / 2.0) * (2.0 * p - 1) / np;
 312
 313     rp = -sin (angle);
 314     ip = cos (angle);
 315   }
 316
 317   /* If we allow ripple, move the pole from the unit
 318    * circle to an ellipse and keep cutoff at frequency 1 */
 319   if (ripple > 0 && type == 1) {
 320     gdouble es, vx;
 321
 322     es = sqrt (pow (10.0, ripple / 10.0) - 1.0);
 323
 324     vx = (1.0 / np) * asinh (1.0 / es);
 325     rp = rp * sinh (vx);
 326     ip = ip * cosh (vx);
 327   } else if (type == 2) {
 328     gdouble es, vx;
 329
 330     es = sqrt (pow (10.0, ripple / 10.0) - 1.0);
 331     vx = (1.0 / np) * asinh (es);
 332     rp = rp * sinh (vx);
 333     ip = ip * cosh (vx);
 334   }
 335
 336   /* Calculate inverse of the pole location to move from
 337    * type I to type II */
 338   if (type == 2) {
 339     gdouble mag2 = rp * rp + ip * ip;
 340
 341     rp /= mag2;
 342     ip /= mag2;
 343   }
 344
 345   /* Calculate zero location for frequency 1 on the
 346    * unit circle for type 2 */
 347   if (type == 2) {
 348     gdouble angle = M_PI / (np * 2.0) + ((p - 1) * M_PI) / (np);
 349     gdouble mag2;
 350
 351     rz = 0.0;
 352     iz = cos (angle);
 353     mag2 = rz * rz + iz * iz;
 354     rz /= mag2;
 355     iz /= mag2;
 356   }
 357
 358   /* Convert from s-domain to z-domain by
 359    * using the bilinear Z-transform, i.e.
 360    * substitute s by (2/t)*((z-1)/(z+1))
 361    * with t = 2 * tan(0.5).
 362    */
 363   if (type == 1) {
 364     gdouble t, m, d;
 365
 366     t = 2.0 * tan (0.5);
 367     m = rp * rp + ip * ip;
 368     d = 4.0 - 4.0 * rp * t + m * t * t;
 369
 370     x0 = (t * t) / d;
 371     x1 = 2.0 * x0;
 372     x2 = x0;
 373     y1 = (8.0 - 2.0 * m * t * t) / d;
 374     y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;
 375   } else {
 376     gdouble t, m, d;
 377
 378     t = 2.0 * tan (0.5);
 379     m = rp * rp + ip * ip;
 380     d = 4.0 - 4.0 * rp * t + m * t * t;
 381
 382     x0 = (t * t * iz * iz + 4.0) / d;
 383     x1 = (-8.0 + 2.0 * iz * iz * t * t) / d;
 384     x2 = x0;
 385     y1 = (8.0 - 2.0 * m * t * t) / d;
 386     y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;
 387   }
 388
 389   /* Convert from lowpass at frequency 1 to either bandpass
 390    * or band reject.
 391    *
 392    * For bandpass substitute z^(-1) with:
 393    *
 394    *   -2            -1
 395    * -z   + alpha * z   - beta
 396    * ----------------------------
 397    *         -2            -1
 398    * beta * z   - alpha * z   + 1
 399    *
 400    * alpha = (2*a*b)/(1+b)
 401    * beta = (b-1)/(b+1)
 402    * a = cos((w1 + w0)/2) / cos((w1 - w0)/2)
 403    * b = tan(1/2) * cot((w1 - w0)/2)
 404    *
 405    * For bandreject substitute z^(-1) with:
 406    *
 407    *  -2            -1
 408    * z   - alpha * z   + beta
 409    * ----------------------------
 410    *         -2            -1
 411    * beta * z   - alpha * z   + 1
 412    *
 413    * alpha = (2*a)/(1+b)
 414    * beta = (1-b)/(1+b)
 415    * a = cos((w1 + w0)/2) / cos((w1 - w0)/2)
 416    * b = tan(1/2) * tan((w1 - w0)/2)
 417    *
 418    */
 419   {
 420     gdouble a, b, d;
 421     gdouble alpha, beta;
 422     gdouble w0 =
 423         2.0 * M_PI * (filter->lower_frequency /
 424         GST_AUDIO_FILTER (filter)->format.rate);
 425     gdouble w1 =
 426         2.0 * M_PI * (filter->upper_frequency /
 427         GST_AUDIO_FILTER (filter)->format.rate);
 428
 429     if (filter->mode == MODE_BAND_PASS) {
 430       a = cos ((w1 + w0) / 2.0) / cos ((w1 - w0) / 2.0);
 431       b = tan (1.0 / 2.0) / tan ((w1 - w0) / 2.0);
 432
 433       alpha = (2.0 * a * b) / (1.0 + b);
 434       beta = (b - 1.0) / (b + 1.0);
 435
 436       d = 1.0 + beta * (y1 - beta * y2);
 437
 438       *a0 = (x0 + beta * (-x1 + beta * x2)) / d;
 439       *a1 = (alpha * (-2.0 * x0 + x1 + beta * x1 - 2.0 * beta * x2)) / d;
 440       *a2 =
 441           (-x1 - beta * beta * x1 + 2.0 * beta * (x0 + x2) +
 442           alpha * alpha * (x0 - x1 + x2)) / d;
 443       *a3 = (alpha * (x1 + beta * (-2.0 * x0 + x1) - 2.0 * x2)) / d;
 444       *a4 = (beta * (beta * x0 - x1) + x2) / d;
 445       *b1 = (alpha * (2.0 + y1 + beta * y1 - 2.0 * beta * y2)) / d;
 446       *b2 =
 447           (-y1 - beta * beta * y1 - alpha * alpha * (1.0 + y1 - y2) +
 448           2.0 * beta * (-1.0 + y2)) / d;
 449       *b3 = (alpha * (y1 + beta * (2.0 + y1) - 2.0 * y2)) / d;
 450       *b4 = (-beta * beta - beta * y1 + y2) / d;
 451     } else {
 452       a = cos ((w1 + w0) / 2.0) / cos ((w1 - w0) / 2.0);
 453       b = tan (1.0 / 2.0) * tan ((w1 - w0) / 2.0);
 454
 455       alpha = (2.0 * a) / (1.0 + b);
 456       beta = (1.0 - b) / (1.0 + b);
 457
 458       d = -1.0 + beta * (beta * y2 + y1);
 459
 460       *a0 = (-x0 - beta * x1 - beta * beta * x2) / d;
 461       *a1 = (alpha * (2.0 * x0 + x1 + beta * x1 + 2.0 * beta * x2)) / d;
 462       *a2 =
 463           (-x1 - beta * beta * x1 - 2.0 * beta * (x0 + x2) -
 464           alpha * alpha * (x0 + x1 + x2)) / d;
 465       *a3 = (alpha * (x1 + beta * (2.0 * x0 + x1) + 2.0 * x2)) / d;
 466       *a4 = (-beta * beta * x0 - beta * x1 - x2) / d;
 467       *b1 = (alpha * (-2.0 + y1 + beta * y1 + 2.0 * beta * y2)) / d;
 468       *b2 =
 469           -(y1 + beta * beta * y1 + 2.0 * beta * (-1.0 + y2) +
 470           alpha * alpha * (-1.0 + y1 + y2)) / d;
 471       *b3 = (alpha * (beta * (-2.0 + y1) + y1 + 2.0 * y2)) / d;
 472       *b4 = -(-beta * beta + beta * y1 + y2) / d;
 473     }
 474   }
 475 }
 476
 477 /* Evaluate the transfer function that corresponds to the IIR
 478  * coefficients at zr + zi*I and return the magnitude */
 479 static gdouble
 480 calculate_gain (gdouble * a, gdouble * b, gint num_a, gint num_b, gdouble zr,
 481     gdouble zi)
 482 {
 483   gdouble sum_ar, sum_ai;
 484   gdouble sum_br, sum_bi;
 485   gdouble gain_r, gain_i;
 486
 487   gdouble sum_r_old;
 488   gdouble sum_i_old;
 489
 490   gint i;
 491
 492   sum_ar = 0.0;
 493   sum_ai = 0.0;
 494   for (i = num_a; i >= 0; i--) {
 495     sum_r_old = sum_ar;
 496     sum_i_old = sum_ai;
 497
 498     sum_ar = (sum_r_old * zr - sum_i_old * zi) + a[i];
 499     sum_ai = (sum_r_old * zi + sum_i_old * zr) + 0.0;
 500   }
 501
 502   sum_br = 0.0;
 503   sum_bi = 0.0;
 504   for (i = num_b; i >= 0; i--) {
 505     sum_r_old = sum_br;
 506     sum_i_old = sum_bi;
 507
 508     sum_br = (sum_r_old * zr - sum_i_old * zi) - b[i];
 509     sum_bi = (sum_r_old * zi + sum_i_old * zr) - 0.0;
 510   }
 511   sum_br += 1.0;
 512   sum_bi += 0.0;
 513
 514   gain_r =
 515       (sum_ar * sum_br + sum_ai * sum_bi) / (sum_br * sum_br + sum_bi * sum_bi);
 516   gain_i =
 517       (sum_ai * sum_br - sum_ar * sum_bi) / (sum_br * sum_br + sum_bi * sum_bi);
 518
 519   return (sqrt (gain_r * gain_r + gain_i * gain_i));
 520 }
 521
 522 static void
 523 generate_coefficients (GstAudioChebyshevFreqBand * filter)
 524 {
 525   gint channels = GST_AUDIO_FILTER (filter)->format.channels;
 526
 527   if (filter->a) {
 528     g_free (filter->a);
 529     filter->a = NULL;
 530   }
 531
 532   if (filter->b) {
 533     g_free (filter->b);
 534     filter->b = NULL;
 535   }
 536
 537   if (filter->channels) {
 538     GstAudioChebyshevFreqBandChannelCtx *ctx;
 539     gint i;
 540
 541     for (i = 0; i < channels; i++) {
 542       ctx = &filter->channels[i];
 543       g_free (ctx->x);
 544       g_free (ctx->y);
 545     }
 546
 547     g_free (filter->channels);
 548     filter->channels = NULL;
 549   }
 550
 551   if (GST_AUDIO_FILTER (filter)->format.rate == 0) {
 552     filter->num_a = 1;
 553     filter->a = g_new0 (gdouble, 1);
 554     filter->a[0] = 1.0;
 555     filter->num_b = 0;
 556     filter->channels = g_new0 (GstAudioChebyshevFreqBandChannelCtx, channels);
 557     GST_LOG_OBJECT (filter, "rate was not set yet");
 558     return;
 559   }
 560
 561   filter->have_coeffs = TRUE;
 562
 563   if (filter->upper_frequency <= filter->lower_frequency) {
 564     filter->num_a = 1;
 565     filter->a = g_new0 (gdouble, 1);
 566     filter->a[0] = (filter->mode == MODE_BAND_PASS) ? 0.0 : 1.0;
 567     filter->num_b = 0;
 568     filter->channels = g_new0 (GstAudioChebyshevFreqBandChannelCtx, channels);
 569     GST_LOG_OBJECT (filter, "frequency band had no or negative dimension");
 570     return;
 571   }
 572
 573   if (filter->upper_frequency > GST_AUDIO_FILTER (filter)->format.rate / 2) {
 574     filter->upper_frequency = GST_AUDIO_FILTER (filter)->format.rate / 2;
 575     GST_LOG_OBJECT (filter, "clipped upper frequency to nyquist frequency");
 576   }
 577
 578   if (filter->lower_frequency < 0.0) {
 579     filter->lower_frequency = 0.0;
 580     GST_LOG_OBJECT (filter, "clipped lower frequency to 0.0");
 581   }
 582
 583   /* Calculate coefficients for the chebyshev filter */
 584   {
 585     gint np = filter->poles;
 586     gdouble *a, *b;
 587     gint i, p;
 588
 589     filter->num_a = np + 1;
 590     filter->a = a = g_new0 (gdouble, np + 5);
 591     filter->num_b = np + 1;
 592     filter->b = b = g_new0 (gdouble, np + 5);
 593
 594     filter->channels = g_new0 (GstAudioChebyshevFreqBandChannelCtx, channels);
 595     for (i = 0; i < channels; i++) {
 596       GstAudioChebyshevFreqBandChannelCtx *ctx = &filter->channels[i];
 597
 598       ctx->x = g_new0 (gdouble, np + 1);
 599       ctx->y = g_new0 (gdouble, np + 1);
 600     }
 601
 602     /* Calculate transfer function coefficients */
 603     a[4] = 1.0;
 604     b[4] = 1.0;
 605
 606     for (p = 1; p <= np / 4; p++) {
 607       gdouble a0, a1, a2, a3, a4, b1, b2, b3, b4;
 608       gdouble *ta = g_new0 (gdouble, np + 5);
 609       gdouble *tb = g_new0 (gdouble, np + 5);
 610
 611       generate_biquad_coefficients (filter, p, &a0, &a1, &a2, &a3, &a4, &b1,
 612           &b2, &b3, &b4);
 613
 614       memcpy (ta, a, sizeof (gdouble) * (np + 5));
 615       memcpy (tb, b, sizeof (gdouble) * (np + 5));
 616
 617       /* add the new coefficients for the new two poles
 618        * to the cascade by multiplication of the transfer
 619        * functions */
 620       for (i = 4; i < np + 5; i++) {
 621         a[i] =
 622             a0 * ta[i] + a1 * ta[i - 1] + a2 * ta[i - 2] + a3 * ta[i - 3] +
 623             a4 * ta[i - 4];
 624         b[i] =
 625             tb[i] - b1 * tb[i - 1] - b2 * tb[i - 2] - b3 * tb[i - 3] -
 626             b4 * tb[i - 4];
 627       }
 628       g_free (ta);
 629       g_free (tb);
 630     }
 631
 632     /* Move coefficients to the beginning of the array
 633      * and multiply the b coefficients with -1 to move from
 634      * the transfer function's coefficients to the difference
 635      * equation's coefficients */
 636     b[4] = 0.0;
 637     for (i = 0; i <= np; i++) {
 638       a[i] = a[i + 4];
 639       b[i] = -b[i + 4];
 640     }
 641
 642     /* Normalize to unity gain at frequency 0 and frequency
 643      * 0.5 for bandreject and unity gain at band center frequency
 644      * for bandpass */
 645     if (filter->mode == MODE_BAND_REJECT) {
 646       /* gain is sqrt(H(0)*H(0.5)) */
 647
 648       gdouble gain1 = calculate_gain (a, b, np, np, 1.0, 0.0);
 649       gdouble gain2 = calculate_gain (a, b, np, np, -1.0, 0.0);
 650
 651       gain1 = sqrt (gain1 * gain2);
 652
 653       for (i = 0; i <= np; i++) {
 654         a[i] /= gain1;
 655       }
 656     } else {
 657       /* gain is H(wc), wc = center frequency */
 658
 659       gdouble w1 =
 660           2.0 * M_PI * (filter->lower_frequency /
 661           GST_AUDIO_FILTER (filter)->format.rate);
 662       gdouble w2 =
 663           2.0 * M_PI * (filter->upper_frequency /
 664           GST_AUDIO_FILTER (filter)->format.rate);
 665       gdouble w0 = (w2 + w1) / 2.0;
 666       gdouble zr = cos (w0), zi = sin (w0);
 667       gdouble gain = calculate_gain (a, b, np, np, zr, zi);
 668
 669       for (i = 0; i <= np; i++) {
 670         a[i] /= gain;
 671       }
 672     }
 673
 674     GST_LOG_OBJECT (filter,
 675         "Generated IIR coefficients for the Chebyshev filter");
 676     GST_LOG_OBJECT (filter,
 677         "mode: %s, type: %d, poles: %d, lower-frequency: %.2f Hz, upper-frequency: %.2f Hz, ripple: %.2f dB",
 678         (filter->mode == MODE_BAND_PASS) ? "band-pass" : "band-reject",
 679         filter->type, filter->poles, filter->lower_frequency,
 680         filter->upper_frequency, filter->ripple);
 681
 682     GST_LOG_OBJECT (filter, "%.2f dB gain @ 0Hz",
 683         20.0 * log10 (calculate_gain (a, b, np, np, 1.0, 0.0)));
 684     {
 685       gdouble w1 =
 686           2.0 * M_PI * (filter->lower_frequency /
 687           GST_AUDIO_FILTER (filter)->format.rate);
 688       gdouble w2 =
 689           2.0 * M_PI * (filter->upper_frequency /
 690           GST_AUDIO_FILTER (filter)->format.rate);
 691       gdouble w0 = (w2 + w1) / 2.0;
 692       gdouble zr, zi;
 693
 694       zr = cos (w1);
 695       zi = sin (w1);
 696       GST_LOG_OBJECT (filter, "%.2f dB gain @ %dHz",
 697           20.0 * log10 (calculate_gain (a, b, np, np, zr, zi)),
 698           (int) filter->lower_frequency);
 699       zr = cos (w0);
 700       zi = sin (w0);
 701       GST_LOG_OBJECT (filter, "%.2f dB gain @ %dHz",
 702           20.0 * log10 (calculate_gain (a, b, np, np, zr, zi)),
 703           (int) ((filter->lower_frequency + filter->upper_frequency) / 2.0));
 704       zr = cos (w2);
 705       zi = sin (w2);
 706       GST_LOG_OBJECT (filter, "%.2f dB gain @ %dHz",
 707           20.0 * log10 (calculate_gain (a, b, np, np, zr, zi)),
 708           (int) filter->upper_frequency);
 709     }
 710     GST_LOG_OBJECT (filter, "%.2f dB gain @ %dHz",
 711         20.0 * log10 (calculate_gain (a, b, np, np, -1.0, 0.0)),
 712         GST_AUDIO_FILTER (filter)->format.rate / 2);
 713   }
 714 }
 715
 716 static void
 717 gst_audio_chebyshev_freq_band_set_property (GObject * object, guint prop_id,
 718     const GValue * value, GParamSpec * pspec)
 719 {
 720   GstAudioChebyshevFreqBand *filter = GST_AUDIO_CHEBYSHEV_FREQ_BAND (object);
 721
 722   switch (prop_id) {
 723     case PROP_MODE:
 724       GST_BASE_TRANSFORM_LOCK (filter);
 725       filter->mode = g_value_get_enum (value);
 726       generate_coefficients (filter);
 727       GST_BASE_TRANSFORM_UNLOCK (filter);
 728       break;
 729     case PROP_TYPE:
 730       GST_BASE_TRANSFORM_LOCK (filter);
 731       filter->type = g_value_get_int (value);
 732       generate_coefficients (filter);
 733       GST_BASE_TRANSFORM_UNLOCK (filter);
 734       break;
 735     case PROP_LOWER_FREQUENCY:
 736       GST_BASE_TRANSFORM_LOCK (filter);
 737       filter->lower_frequency = g_value_get_float (value);
 738       generate_coefficients (filter);
 739       GST_BASE_TRANSFORM_UNLOCK (filter);
 740       break;
 741     case PROP_UPPER_FREQUENCY:
 742       GST_BASE_TRANSFORM_LOCK (filter);
 743       filter->upper_frequency = g_value_get_float (value);
 744       generate_coefficients (filter);
 745       GST_BASE_TRANSFORM_UNLOCK (filter);
 746       break;
 747     case PROP_RIPPLE:
 748       GST_BASE_TRANSFORM_LOCK (filter);
 749       filter->ripple = g_value_get_float (value);
 750       generate_coefficients (filter);
 751       GST_BASE_TRANSFORM_UNLOCK (filter);
 752       break;
 753     case PROP_POLES:
 754       GST_BASE_TRANSFORM_LOCK (filter);
 755       filter->poles = GST_ROUND_UP_4 (g_value_get_int (value));
 756       generate_coefficients (filter);
 757       GST_BASE_TRANSFORM_UNLOCK (filter);
 758       break;
 759     default:
 760       G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
 761       break;
 762   }
 763 }
 764
 765 static void
 766 gst_audio_chebyshev_freq_band_get_property (GObject * object, guint prop_id,
 767     GValue * value, GParamSpec * pspec)
 768 {
 769   GstAudioChebyshevFreqBand *filter = GST_AUDIO_CHEBYSHEV_FREQ_BAND (object);
 770
 771   switch (prop_id) {
 772     case PROP_MODE:
 773       g_value_set_enum (value, filter->mode);
 774       break;
 775     case PROP_TYPE:
 776       g_value_set_int (value, filter->type);
 777       break;
 778     case PROP_LOWER_FREQUENCY:
 779       g_value_set_float (value, filter->lower_frequency);
 780       break;
 781     case PROP_UPPER_FREQUENCY:
 782       g_value_set_float (value, filter->upper_frequency);
 783       break;
 784     case PROP_RIPPLE:
 785       g_value_set_float (value, filter->ripple);
 786       break;
 787     case PROP_POLES:
 788       g_value_set_int (value, filter->poles);
 789       break;
 790     default:
 791       G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
 792       break;
 793   }
 794 }
 795
 796 /* GstAudioFilter vmethod implementations */
 797
 798 static gboolean
 799 gst_audio_chebyshev_freq_band_setup (GstAudioFilter * base,
 800     GstRingBufferSpec * format)
 801 {
 802   GstAudioChebyshevFreqBand *filter = GST_AUDIO_CHEBYSHEV_FREQ_BAND (base);
 803   gboolean ret = TRUE;
 804
 805   if (format->width == 32)
 806     filter->process = (GstAudioChebyshevFreqBandProcessFunc)
 807         process_32;
 808   else if (format->width == 64)
 809     filter->process = (GstAudioChebyshevFreqBandProcessFunc)
 810         process_64;
 811   else
 812     ret = FALSE;
 813
 814   filter->have_coeffs = FALSE;
 815
 816   return ret;
 817 }
 818
 819 static inline gdouble
 820 process (GstAudioChebyshevFreqBand * filter,
 821     GstAudioChebyshevFreqBandChannelCtx * ctx, gdouble x0)
 822 {
 823   gdouble val = filter->a[0] * x0;
 824   gint i, j;
 825
 826   for (i = 1, j = ctx->x_pos; i < filter->num_a; i++) {
 827     val += filter->a[i] * ctx->x[j];
 828     j--;
 829     if (j < 0)
 830       j = filter->num_a - 1;
 831   }
 832
 833   for (i = 1, j = ctx->y_pos; i < filter->num_b; i++) {
 834     val += filter->b[i] * ctx->y[j];
 835     j--;
 836     if (j < 0)
 837       j = filter->num_b - 1;
 838   }
 839
 840   if (ctx->x) {
 841     ctx->x_pos++;
 842     if (ctx->x_pos > filter->num_a - 1)
 843       ctx->x_pos = 0;
 844     ctx->x[ctx->x_pos] = x0;
 845   }
 846
 847   if (ctx->y) {
 848     ctx->y_pos++;
 849     if (ctx->y_pos > filter->num_b - 1)
 850       ctx->y_pos = 0;
 851
 852     ctx->y[ctx->y_pos] = val;
 853   }
 854
 855   return val;
 856 }
 857
 858 #define DEFINE_PROCESS_FUNC(width,ctype) \
 859 static void \
 860 process_##width (GstAudioChebyshevFreqBand * filter, \
 861     g##ctype * data, guint num_samples) \
 862 { \
 863   gint i, j, channels = GST_AUDIO_FILTER (filter)->format.channels; \
 864   gdouble val; \
 865   \
 866   for (i = 0; i < num_samples / channels; i++) { \
 867     for (j = 0; j < channels; j++) { \
 868       val = process (filter, &filter->channels[j], *data); \
 869       *data++ = val; \
 870     } \
 871   } \
 872 }
 873
 874 DEFINE_PROCESS_FUNC (32, float);
 875 DEFINE_PROCESS_FUNC (64, double);
 876
 877 #undef DEFINE_PROCESS_FUNC
 878
 879 /* GstBaseTransform vmethod implementations */
 880 static GstFlowReturn
 881 gst_audio_chebyshev_freq_band_transform_ip (GstBaseTransform * base,
 882     GstBuffer * buf)
 883 {
 884   GstAudioChebyshevFreqBand *filter = GST_AUDIO_CHEBYSHEV_FREQ_BAND (base);
 885   guint num_samples =
 886       GST_BUFFER_SIZE (buf) / (GST_AUDIO_FILTER (filter)->format.width / 8);
 887
 888   if (!gst_buffer_is_writable (buf))
 889     return GST_FLOW_OK;
 890
 891   if (GST_CLOCK_TIME_IS_VALID (GST_BUFFER_TIMESTAMP (buf)))
 892     gst_object_sync_values (G_OBJECT (filter), GST_BUFFER_TIMESTAMP (buf));
 893
 894   if (!filter->have_coeffs)
 895     generate_coefficients (filter);
 896
 897   filter->process (filter, GST_BUFFER_DATA (buf), num_samples);
 898
 899   return GST_FLOW_OK;
 900 }
 901
 902 static gboolean
 903 gst_audio_chebyshev_freq_band_start (GstBaseTransform * base)
 904 {
 905   GstAudioChebyshevFreqBand *filter = GST_AUDIO_CHEBYSHEV_FREQ_BAND (base);
 906   gint channels = GST_AUDIO_FILTER (filter)->format.channels;
 907   GstAudioChebyshevFreqBandChannelCtx *ctx;
 908   gint i;
 909
 910   /* Reset the history of input and output values if
 911    * already existing */
 912   if (channels && filter->channels) {
 913     for (i = 0; i < channels; i++) {
 914       ctx = &filter->channels[i];
 915       if (ctx->x)
 916         memset (ctx->x, 0, (filter->poles + 1) * sizeof (gdouble));
 917       if (ctx->y)
 918         memset (ctx->y, 0, (filter->poles + 1) * sizeof (gdouble));
 919     }
 920   }
 921   return TRUE;
 922 }