audiocheblimit.c 16.9 KB
Newer Older
1 2
/* 
 * GStreamer
3
 * Copyright (C) 2007-2009 Sebastian Dröge <sebastian.droege@collabora.co.uk>
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/* 
 * Chebyshev type 1 filter design based on
 * "The Scientist and Engineer's Guide to DSP", Chapter 20.
 * http://www.dspguide.com/
 *
 * For type 2 and Chebyshev filters in general read
 * http://en.wikipedia.org/wiki/Chebyshev_filter
 *
 */

/**
32
 * SECTION:element-audiocheblimit
33 34 35
 *
 * Attenuates all frequencies above the cutoff frequency (low-pass) or all frequencies below the
 * cutoff frequency (high-pass). The number of poles and the ripple parameter control the rolloff.
36
 *
37 38 39
 * This element has the advantage over the windowed sinc lowpass and highpass filter that it is
 * much faster and produces almost as good results. It's only disadvantages are the highly
 * non-linear phase and the slower rolloff compared to a windowed sinc filter with a large kernel.
40
 *
41 42 43
 * For type 1 the ripple parameter specifies how much ripple in dB is allowed in the passband, i.e.
 * some frequencies in the passband will be amplified by that value. A higher ripple value will allow
 * a faster rolloff.
44
 *
45 46
 * For type 2 the ripple parameter specifies the stopband attenuation. In the stopband the gain will
 * be at most this value. A lower ripple value will allow a faster rolloff.
47
 *
48 49
 * As a special case, a Chebyshev type 1 filter with no ripple is a Butterworth filter.
 * </para>
50
 * <note><para>
51 52
 * Be warned that a too large number of poles can produce noise. The most poles are possible with
 * a cutoff frequency at a quarter of the sampling rate.
53
 * </para></note>
54
 * <para>
55 56 57
 * <refsect2>
 * <title>Example launch line</title>
 * |[
58 59 60
 * gst-launch audiotestsrc freq=1500 ! audioconvert ! audiocheblimit mode=low-pass cutoff=1000 poles=4 ! audioconvert ! alsasink
 * gst-launch filesrc location="melo1.ogg" ! oggdemux ! vorbisdec ! audioconvert ! audiocheblimit mode=high-pass cutoff=400 ripple=0.2 ! audioconvert ! alsasink
 * gst-launch audiotestsrc wave=white-noise ! audioconvert ! audiocheblimit mode=low-pass cutoff=800 type=2 ! audioconvert ! alsasink
61
 * ]|
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
 * </refsect2>
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <gst/gst.h>
#include <gst/base/gstbasetransform.h>
#include <gst/audio/audio.h>
#include <gst/audio/gstaudiofilter.h>
#include <gst/controller/gstcontroller.h>

#include <math.h>

77 78
#include "math_compat.h"

79
#include "audiocheblimit.h"
80

81 82
#include "gst/glib-compat-private.h"

83
#define GST_CAT_DEFAULT gst_audio_cheb_limit_debug
84 85 86 87 88 89 90 91 92 93 94 95 96
GST_DEBUG_CATEGORY_STATIC (GST_CAT_DEFAULT);

enum
{
  PROP_0,
  PROP_MODE,
  PROP_TYPE,
  PROP_CUTOFF,
  PROP_RIPPLE,
  PROP_POLES
};

#define DEBUG_INIT(bla) \
97
  GST_DEBUG_CATEGORY_INIT (gst_audio_cheb_limit_debug, "audiocheblimit", 0, "audiocheblimit element");
98

99
GST_BOILERPLATE_FULL (GstAudioChebLimit,
100 101
    gst_audio_cheb_limit, GstAudioFXBaseIIRFilter,
    GST_TYPE_AUDIO_FX_BASE_IIR_FILTER, DEBUG_INIT);
102

103
static void gst_audio_cheb_limit_set_property (GObject * object,
104
    guint prop_id, const GValue * value, GParamSpec * pspec);
105
static void gst_audio_cheb_limit_get_property (GObject * object,
106
    guint prop_id, GValue * value, GParamSpec * pspec);
107
static void gst_audio_cheb_limit_finalize (GObject * object);
108

109
static gboolean gst_audio_cheb_limit_setup (GstAudioFilter * filter,
110 111 112 113 114 115 116 117
    GstRingBufferSpec * format);

enum
{
  MODE_LOW_PASS = 0,
  MODE_HIGH_PASS
};

118
#define GST_TYPE_AUDIO_CHEBYSHEV_FREQ_LIMIT_MODE (gst_audio_cheb_limit_mode_get_type ())
119
static GType
120
gst_audio_cheb_limit_mode_get_type (void)
121 122 123 124 125 126 127 128 129 130 131 132
{
  static GType gtype = 0;

  if (gtype == 0) {
    static const GEnumValue values[] = {
      {MODE_LOW_PASS, "Low pass (default)",
          "low-pass"},
      {MODE_HIGH_PASS, "High pass",
          "high-pass"},
      {0, NULL, NULL}
    };

133
    gtype = g_enum_register_static ("GstAudioChebLimitMode", values);
134 135 136 137 138 139 140
  }
  return gtype;
}

/* GObject vmethod implementations */

static void
141
gst_audio_cheb_limit_base_init (gpointer klass)
142 143 144
{
  GstElementClass *element_class = GST_ELEMENT_CLASS (klass);

145 146 147 148 149
  gst_element_class_set_details_simple (element_class,
      "Low pass & high pass filter",
      "Filter/Effect/Audio",
      "Chebyshev low pass and high pass filter",
      "Sebastian Dröge <sebastian.droege@collabora.co.uk>");
150 151 152
}

static void
153
gst_audio_cheb_limit_class_init (GstAudioChebLimitClass * klass)
154
{
155 156
  GObjectClass *gobject_class = (GObjectClass *) klass;
  GstAudioFilterClass *filter_class = (GstAudioFilterClass *) klass;
157

158 159
  gobject_class->set_property = gst_audio_cheb_limit_set_property;
  gobject_class->get_property = gst_audio_cheb_limit_get_property;
160
  gobject_class->finalize = gst_audio_cheb_limit_finalize;
161 162 163 164 165

  g_object_class_install_property (gobject_class, PROP_MODE,
      g_param_spec_enum ("mode", "Mode",
          "Low pass or high pass mode",
          GST_TYPE_AUDIO_CHEBYSHEV_FREQ_LIMIT_MODE, MODE_LOW_PASS,
166
          G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
167 168
  g_object_class_install_property (gobject_class, PROP_TYPE,
      g_param_spec_int ("type", "Type", "Type of the chebychev filter", 1, 2, 1,
169
          G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
170 171 172

  /* FIXME: Don't use the complete possible range but restrict the upper boundary
   * so automatically generated UIs can use a slider without */
173 174
  g_object_class_install_property (gobject_class, PROP_CUTOFF,
      g_param_spec_float ("cutoff", "Cutoff", "Cut off frequency (Hz)", 0.0,
175 176
          100000.0, 0.0,
          G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
177 178
  g_object_class_install_property (gobject_class, PROP_RIPPLE,
      g_param_spec_float ("ripple", "Ripple", "Amount of ripple (dB)", 0.0,
179 180
          200.0, 0.25,
          G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
181 182 183 184

  /* FIXME: What to do about this upper boundary? With a cutoff frequency of
   * rate/4 32 poles are completely possible, with a cutoff frequency very low
   * or very high 16 poles already produces only noise */
185 186 187
  g_object_class_install_property (gobject_class, PROP_POLES,
      g_param_spec_int ("poles", "Poles",
          "Number of poles to use, will be rounded up to the next even number",
188 189
          2, 32, 4,
          G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
190

191
  filter_class->setup = GST_DEBUG_FUNCPTR (gst_audio_cheb_limit_setup);
192 193 194
}

static void
195 196
gst_audio_cheb_limit_init (GstAudioChebLimit * filter,
    GstAudioChebLimitClass * klass)
197 198 199 200 201 202
{
  filter->cutoff = 0.0;
  filter->mode = MODE_LOW_PASS;
  filter->type = 1;
  filter->poles = 4;
  filter->ripple = 0.25;
203 204

  filter->lock = g_mutex_new ();
205 206 207
}

static void
208
generate_biquad_coefficients (GstAudioChebLimit * filter,
209 210 211 212 213 214 215 216 217 218
    gint p, gdouble * a0, gdouble * a1, gdouble * a2,
    gdouble * b1, gdouble * b2)
{
  gint np = filter->poles;
  gdouble ripple = filter->ripple;

  /* pole location in s-plane */
  gdouble rp, ip;

  /* zero location in s-plane */
219
  gdouble iz = 0.0;
220 221 222 223 224 225 226

  /* transfer function coefficients for the z-plane */
  gdouble x0, x1, x2, y1, y2;
  gint type = filter->type;

  /* Calculate pole location for lowpass at frequency 1 */
  {
David Schleef's avatar
David Schleef committed
227
    gdouble angle = (G_PI / 2.0) * (2.0 * p - 1) / np;
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263

    rp = -sin (angle);
    ip = cos (angle);
  }

  /* If we allow ripple, move the pole from the unit
   * circle to an ellipse and keep cutoff at frequency 1 */
  if (ripple > 0 && type == 1) {
    gdouble es, vx;

    es = sqrt (pow (10.0, ripple / 10.0) - 1.0);

    vx = (1.0 / np) * asinh (1.0 / es);
    rp = rp * sinh (vx);
    ip = ip * cosh (vx);
  } else if (type == 2) {
    gdouble es, vx;

    es = sqrt (pow (10.0, ripple / 10.0) - 1.0);
    vx = (1.0 / np) * asinh (es);
    rp = rp * sinh (vx);
    ip = ip * cosh (vx);
  }

  /* Calculate inverse of the pole location to convert from
   * type I to type II */
  if (type == 2) {
    gdouble mag2 = rp * rp + ip * ip;

    rp /= mag2;
    ip /= mag2;
  }

  /* Calculate zero location for frequency 1 on the
   * unit circle for type 2 */
  if (type == 2) {
David Schleef's avatar
David Schleef committed
264
    gdouble angle = G_PI / (np * 2.0) + ((p - 1) * G_PI) / (np);
265 266 267
    gdouble mag2;

    iz = cos (angle);
268
    mag2 = iz * iz;
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
    iz /= mag2;
  }

  /* Convert from s-domain to z-domain by
   * using the bilinear Z-transform, i.e.
   * substitute s by (2/t)*((z-1)/(z+1))
   * with t = 2 * tan(0.5).
   */
  if (type == 1) {
    gdouble t, m, d;

    t = 2.0 * tan (0.5);
    m = rp * rp + ip * ip;
    d = 4.0 - 4.0 * rp * t + m * t * t;

    x0 = (t * t) / d;
    x1 = 2.0 * x0;
    x2 = x0;
    y1 = (8.0 - 2.0 * m * t * t) / d;
    y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;
  } else {
    gdouble t, m, d;

    t = 2.0 * tan (0.5);
    m = rp * rp + ip * ip;
    d = 4.0 - 4.0 * rp * t + m * t * t;

    x0 = (t * t * iz * iz + 4.0) / d;
    x1 = (-8.0 + 2.0 * iz * iz * t * t) / d;
    x2 = x0;
    y1 = (8.0 - 2.0 * m * t * t) / d;
    y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;
  }

  /* Convert from lowpass at frequency 1 to either lowpass
   * or highpass.
   *
   * For lowpass substitute z^(-1) with:
   *  -1
   * z   - k
   * ------------
   *          -1
   * 1 - k * z
   *
   * k = sin((1-w)/2) / sin((1+w)/2)
   *
   * For highpass substitute z^(-1) with:
   *
   *   -1
   * -z   - k
   * ------------
   *          -1
   * 1 + k * z
   *
   * k = -cos((1+w)/2) / cos((1-w)/2)
   *
   */
  {
    gdouble k, d;
    gdouble omega =
David Schleef's avatar
David Schleef committed
329
        2.0 * G_PI * (filter->cutoff / GST_AUDIO_FILTER (filter)->format.rate);
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350

    if (filter->mode == MODE_LOW_PASS)
      k = sin ((1.0 - omega) / 2.0) / sin ((1.0 + omega) / 2.0);
    else
      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;

    if (filter->mode == MODE_HIGH_PASS) {
      *a1 = -*a1;
      *b1 = -*b1;
    }
  }
}

static void
351
generate_coefficients (GstAudioChebLimit * filter)
352
{
353 354
  if (GST_AUDIO_FILTER (filter)->format.rate == 0) {
    gdouble *a = g_new0 (gdouble, 1);
355

356 357 358
    a[0] = 1.0;
    gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
        (filter), a, 1, NULL, 0);
359 360 361 362 363 364

    GST_LOG_OBJECT (filter, "rate was not set yet");
    return;
  }

  if (filter->cutoff >= GST_AUDIO_FILTER (filter)->format.rate / 2.0) {
365 366 367 368 369
    gdouble *a = g_new0 (gdouble, 1);

    a[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);
370 371 372
    GST_LOG_OBJECT (filter, "cutoff was higher than nyquist frequency");
    return;
  } else if (filter->cutoff <= 0.0) {
373 374 375 376 377
    gdouble *a = g_new0 (gdouble, 1);

    a[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);
378 379 380 381 382 383 384 385 386 387
    GST_LOG_OBJECT (filter, "cutoff is lower than zero");
    return;
  }

  /* Calculate coefficients for the chebyshev filter */
  {
    gint np = filter->poles;
    gdouble *a, *b;
    gint i, p;

388 389
    a = g_new0 (gdouble, np + 3);
    b = g_new0 (gdouble, np + 3);
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431

    /* Calculate transfer function coefficients */
    a[2] = 1.0;
    b[2] = 1.0;

    for (p = 1; p <= np / 2; p++) {
      gdouble a0, a1, a2, b1, b2;
      gdouble *ta = g_new0 (gdouble, np + 3);
      gdouble *tb = g_new0 (gdouble, np + 3);

      generate_biquad_coefficients (filter, p, &a0, &a1, &a2, &b1, &b2);

      memcpy (ta, a, sizeof (gdouble) * (np + 3));
      memcpy (tb, b, sizeof (gdouble) * (np + 3));

      /* add the new coefficients for the new two poles
       * 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];
      }
      g_free (ta);
      g_free (tb);
    }

    /* Move coefficients to the beginning of the array
     * and multiply the b coefficients with -1 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];
    }

    /* Normalize to unity gain at frequency 0 for lowpass
     * and frequency 0.5 for highpass */
    {
      gdouble gain;

      if (filter->mode == MODE_LOW_PASS)
432 433 434
        gain =
            gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b, np + 1,
            1.0, 0.0);
435
      else
436 437 438
        gain =
            gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b, np + 1,
            -1.0, 0.0);
439 440 441 442 443 444

      for (i = 0; i <= np; i++) {
        a[i] /= gain;
      }
    }

445 446 447
    gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
        (filter), a, np + 1, b, np + 1);

448 449 450 451 452 453 454
    GST_LOG_OBJECT (filter,
        "Generated IIR coefficients for the Chebyshev filter");
    GST_LOG_OBJECT (filter,
        "mode: %s, type: %d, poles: %d, cutoff: %.2f Hz, ripple: %.2f dB",
        (filter->mode == MODE_LOW_PASS) ? "low-pass" : "high-pass",
        filter->type, filter->poles, filter->cutoff, filter->ripple);
    GST_LOG_OBJECT (filter, "%.2f dB gain @ 0 Hz",
455 456
        20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b,
                np + 1, 1.0, 0.0)));
457 458

#ifndef GST_DISABLE_GST_DEBUG
459 460
    {
      gdouble wc =
David Schleef's avatar
David Schleef committed
461
          2.0 * G_PI * (filter->cutoff /
462 463 464 465
          GST_AUDIO_FILTER (filter)->format.rate);
      gdouble zr = cos (wc), zi = sin (wc);

      GST_LOG_OBJECT (filter, "%.2f dB gain @ %d Hz",
466 467
          20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1,
                  b, np + 1, zr, zi)), (int) filter->cutoff);
468
    }
469 470
#endif

471
    GST_LOG_OBJECT (filter, "%.2f dB gain @ %d Hz",
472 473
        20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b,
                np + 1, -1.0, 0.0)),
474 475 476 477
        GST_AUDIO_FILTER (filter)->format.rate / 2);
  }
}

478 479 480 481 482 483 484 485 486 487 488
static void
gst_audio_cheb_limit_finalize (GObject * object)
{
  GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);

  g_mutex_free (filter->lock);
  filter->lock = NULL;

  G_OBJECT_CLASS (parent_class)->finalize (object);
}

489
static void
490
gst_audio_cheb_limit_set_property (GObject * object, guint prop_id,
491 492
    const GValue * value, GParamSpec * pspec)
{
493
  GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);
494 495 496

  switch (prop_id) {
    case PROP_MODE:
497
      g_mutex_lock (filter->lock);
498 499
      filter->mode = g_value_get_enum (value);
      generate_coefficients (filter);
500
      g_mutex_unlock (filter->lock);
501 502
      break;
    case PROP_TYPE:
503
      g_mutex_lock (filter->lock);
504 505
      filter->type = g_value_get_int (value);
      generate_coefficients (filter);
506
      g_mutex_unlock (filter->lock);
507 508
      break;
    case PROP_CUTOFF:
509
      g_mutex_lock (filter->lock);
510 511
      filter->cutoff = g_value_get_float (value);
      generate_coefficients (filter);
512
      g_mutex_unlock (filter->lock);
513 514
      break;
    case PROP_RIPPLE:
515
      g_mutex_lock (filter->lock);
516 517
      filter->ripple = g_value_get_float (value);
      generate_coefficients (filter);
518
      g_mutex_unlock (filter->lock);
519 520
      break;
    case PROP_POLES:
521
      g_mutex_lock (filter->lock);
522 523
      filter->poles = GST_ROUND_UP_2 (g_value_get_int (value));
      generate_coefficients (filter);
524
      g_mutex_unlock (filter->lock);
525 526 527 528 529 530 531 532
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
533
gst_audio_cheb_limit_get_property (GObject * object, guint prop_id,
534 535
    GValue * value, GParamSpec * pspec)
{
536
  GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);
537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562

  switch (prop_id) {
    case PROP_MODE:
      g_value_set_enum (value, filter->mode);
      break;
    case PROP_TYPE:
      g_value_set_int (value, filter->type);
      break;
    case PROP_CUTOFF:
      g_value_set_float (value, filter->cutoff);
      break;
    case PROP_RIPPLE:
      g_value_set_float (value, filter->ripple);
      break;
    case PROP_POLES:
      g_value_set_int (value, filter->poles);
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

/* GstAudioFilter vmethod implementations */

static gboolean
563
gst_audio_cheb_limit_setup (GstAudioFilter * base, GstRingBufferSpec * format)
564
{
565
  GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (base);
566

567
  generate_coefficients (filter);
568

569
  return GST_AUDIO_FILTER_CLASS (parent_class)->setup (base, format);
570
}