1 /* 

     2  * GStreamer

     3  * Copyright (C) 2007-2009 Sebastian Dröge <sebastian.droege@collabora.co.uk>

     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  */

    30 

    31 /**

    32  * SECTION:element-audiocheblimit

    33  *

    34  * Attenuates all frequencies above the cutoff frequency (low-pass) or all frequencies below the

    35  * cutoff frequency (high-pass). The number of poles and the ripple parameter control the rolloff.

    36  *

    37  * This element has the advantage over the windowed sinc lowpass and highpass filter that it is

    38  * much faster and produces almost as good results. It's only disadvantages are the highly

    39  * non-linear phase and the slower rolloff compared to a windowed sinc filter with a large kernel.

    40  *

    41  * For type 1 the ripple parameter specifies how much ripple in dB is allowed in the passband, i.e.

    42  * some frequencies in the passband will be amplified by that value. A higher ripple value will allow

    43  * a faster rolloff.

    44  *

    45  * For type 2 the ripple parameter specifies the stopband attenuation. In the stopband the gain will

    46  * be at most this value. A lower ripple value will allow a faster rolloff.

    47  *

    48  * As a special case, a Chebyshev type 1 filter with no ripple is a Butterworth filter.

    49  * </para>

    50  * <note><para>

    51  * Be warned that a too large number of poles can produce noise. The most poles are possible with

    52  * a cutoff frequency at a quarter of the sampling rate.

    53  * </para></note>

    54  * <para>

    55  * <refsect2>

    56  * <title>Example launch line</title>

    57  * |[

    58  * gst-launch audiotestsrc freq=1500 ! audioconvert ! audiocheblimit mode=low-pass cutoff=1000 poles=4 ! audioconvert ! alsasink

    59  * gst-launch filesrc location="melo1.ogg" ! oggdemux ! vorbisdec ! audioconvert ! audiocheblimit mode=high-pass cutoff=400 ripple=0.2 ! audioconvert ! alsasink

    60  * gst-launch audiotestsrc wave=white-noise ! audioconvert ! audiocheblimit mode=low-pass cutoff=800 type=2 ! audioconvert ! alsasink

    61  * ]|

    62  * </refsect2>

    63  */

    64 

    65 #ifdef HAVE_CONFIG_H

    66 #include "config.h"

    67 #endif

    68 

    69 #include <gst/gst.h>

    70 #include <gst/base/gstbasetransform.h>

    71 #include <gst/audio/audio.h>

    72 #include <gst/audio/gstaudiofilter.h>

    73 #include <gst/controller/gstcontroller.h>

    74 

    75 #include <math.h>

    76 

    77 #include "math_compat.h"

    78 

    79 #include "audiocheblimit.h"

    80 

    81 #define GST_CAT_DEFAULT gst_audio_cheb_limit_debug

    82 GST_DEBUG_CATEGORY_STATIC (GST_CAT_DEFAULT);

    83 

    84 enum

    85 {

    86   PROP_0,

    87   PROP_MODE,

    88   PROP_TYPE,

    89   PROP_CUTOFF,

    90   PROP_RIPPLE,

    91   PROP_POLES

    92 };

    93 

    94 #define DEBUG_INIT(bla) \

    95   GST_DEBUG_CATEGORY_INIT (gst_audio_cheb_limit_debug, "audiocheblimit", 0, "audiocheblimit element");

    96 

    97 GST_BOILERPLATE_FULL (GstAudioChebLimit,

    98     gst_audio_cheb_limit, GstAudioFXBaseIIRFilter,

    99     GST_TYPE_AUDIO_FX_BASE_IIR_FILTER, DEBUG_INIT);

   100 

   101 static void gst_audio_cheb_limit_set_property (GObject * object,

   102     guint prop_id, const GValue * value, GParamSpec * pspec);

   103 static void gst_audio_cheb_limit_get_property (GObject * object,

   104     guint prop_id, GValue * value, GParamSpec * pspec);

   105 static void gst_audio_cheb_limit_finalize (GObject * object);

   106 

   107 static gboolean gst_audio_cheb_limit_setup (GstAudioFilter * filter,

   108     GstRingBufferSpec * format);

   109 

   110 enum

   111 {

   112   MODE_LOW_PASS = 0,

   113   MODE_HIGH_PASS

   114 };

   115 

   116 #define GST_TYPE_AUDIO_CHEBYSHEV_FREQ_LIMIT_MODE (gst_audio_cheb_limit_mode_get_type ())

   117 static GType

   118 gst_audio_cheb_limit_mode_get_type (void)

   119 {

   120   static GType gtype = 0;

   121 

   122   if (gtype == 0) {

   123     static const GEnumValue values[] = {

   124       {MODE_LOW_PASS, "Low pass (default)",

   125           "low-pass"},

   126       {MODE_HIGH_PASS, "High pass",

   127           "high-pass"},

   128       {0, NULL, NULL}

   129     };

   130 

   131     gtype = g_enum_register_static ("GstAudioChebLimitMode", values);

   132   }

   133   return gtype;

   134 }

   135 

   136 /* GObject vmethod implementations */

   137 

   138 static void

   139 gst_audio_cheb_limit_base_init (gpointer klass)

   140 {

   141   GstElementClass *element_class = GST_ELEMENT_CLASS (klass);

   142 

   143   gst_element_class_set_details_simple (element_class,

   144       "Low pass & high pass filter",

   145       "Filter/Effect/Audio",

   146       "Chebyshev low pass and high pass filter",

   147       "Sebastian Dröge <sebastian.droege@collabora.co.uk>");

   148 }

   149 

   150 static void

   151 gst_audio_cheb_limit_class_init (GstAudioChebLimitClass * klass)

   152 {

   153   GObjectClass *gobject_class = (GObjectClass *) klass;

   154   GstAudioFilterClass *filter_class = (GstAudioFilterClass *) klass;

   155 

   156   gobject_class->set_property = gst_audio_cheb_limit_set_property;

   157   gobject_class->get_property = gst_audio_cheb_limit_get_property;

   158   gobject_class->finalize = gst_audio_cheb_limit_finalize;

   159 

   160   g_object_class_install_property (gobject_class, PROP_MODE,

   161       g_param_spec_enum ("mode", "Mode",

   162           "Low pass or high pass mode",

   163           GST_TYPE_AUDIO_CHEBYSHEV_FREQ_LIMIT_MODE, MODE_LOW_PASS,

   164           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));

   165   g_object_class_install_property (gobject_class, PROP_TYPE,

   166       g_param_spec_int ("type", "Type", "Type of the chebychev filter", 1, 2, 1,

   167           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));

   168 

   169   /* FIXME: Don't use the complete possible range but restrict the upper boundary

   170    * so automatically generated UIs can use a slider without */

   171   g_object_class_install_property (gobject_class, PROP_CUTOFF,

   172       g_param_spec_float ("cutoff", "Cutoff", "Cut off frequency (Hz)", 0.0,

   173           100000.0, 0.0,

   174           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));

   175   g_object_class_install_property (gobject_class, PROP_RIPPLE,

   176       g_param_spec_float ("ripple", "Ripple", "Amount of ripple (dB)", 0.0,

   177           200.0, 0.25,

   178           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));

   179 

   180   /* FIXME: What to do about this upper boundary? With a cutoff frequency of

   181    * rate/4 32 poles are completely possible, with a cutoff frequency very low

   182    * or very high 16 poles already produces only noise */

   183   g_object_class_install_property (gobject_class, PROP_POLES,

   184       g_param_spec_int ("poles", "Poles",

   185           "Number of poles to use, will be rounded up to the next even number",

   186           2, 32, 4,

   187           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));

   188 

   189   filter_class->setup = GST_DEBUG_FUNCPTR (gst_audio_cheb_limit_setup);

   190 }

   191 

   192 static void

   193 gst_audio_cheb_limit_init (GstAudioChebLimit * filter,

   194     GstAudioChebLimitClass * klass)

   195 {

   196   filter->cutoff = 0.0;

   197   filter->mode = MODE_LOW_PASS;

   198   filter->type = 1;

   199   filter->poles = 4;

   200   filter->ripple = 0.25;

   201 

   202   filter->lock = g_mutex_new ();

   203 }

   204 

   205 static void

   206 generate_biquad_coefficients (GstAudioChebLimit * filter,

   207     gint p, gdouble * a0, gdouble * a1, gdouble * a2,

   208     gdouble * b1, gdouble * b2)

   209 {

   210   gint np = filter->poles;

   211   gdouble ripple = filter->ripple;

   212 

   213   /* pole location in s-plane */

   214   gdouble rp, ip;

   215 

   216   /* zero location in s-plane */

   217   gdouble iz = 0.0;

   218 

   219   /* transfer function coefficients for the z-plane */

   220   gdouble x0, x1, x2, y1, y2;

   221   gint type = filter->type;

   222 

   223   /* Calculate pole location for lowpass at frequency 1 */

   224   {

   225     gdouble angle = (M_PI / 2.0) * (2.0 * p - 1) / np;

   226 

   227     rp = -sin (angle);

   228     ip = cos (angle);

   229   }

   230 

   231   /* If we allow ripple, move the pole from the unit

   232    * circle to an ellipse and keep cutoff at frequency 1 */

   233   if (ripple > 0 && type == 1) {

   234     gdouble es, vx;

   235 

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

   237 

   238     vx = (1.0 / np) * asinh (1.0 / es);

   239     rp = rp * sinh (vx);

   240     ip = ip * cosh (vx);

   241   } else if (type == 2) {

   242     gdouble es, vx;

   243 

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

   245     vx = (1.0 / np) * asinh (es);

   246     rp = rp * sinh (vx);

   247     ip = ip * cosh (vx);

   248   }

   249 

   250   /* Calculate inverse of the pole location to convert from

   251    * type I to type II */

   252   if (type == 2) {

   253     gdouble mag2 = rp * rp + ip * ip;

   254 

   255     rp /= mag2;

   256     ip /= mag2;

   257   }

   258 

   259   /* Calculate zero location for frequency 1 on the

   260    * unit circle for type 2 */

   261   if (type == 2) {

   262     gdouble angle = M_PI / (np * 2.0) + ((p - 1) * M_PI) / (np);

   263     gdouble mag2;

   264 

   265     iz = cos (angle);

   266     mag2 = iz * iz;

   267     iz /= mag2;

   268   }

   269 

   270   /* Convert from s-domain to z-domain by

   271    * using the bilinear Z-transform, i.e.

   272    * substitute s by (2/t)*((z-1)/(z+1))

   273    * with t = 2 * tan(0.5).

   274    */

   275   if (type == 1) {

   276     gdouble t, m, d;

   277 

   278     t = 2.0 * tan (0.5);

   279     m = rp * rp + ip * ip;

   280     d = 4.0 - 4.0 * rp * t + m * t * t;

   281 

   282     x0 = (t * t) / d;

   283     x1 = 2.0 * x0;

   284     x2 = x0;

   285     y1 = (8.0 - 2.0 * m * t * t) / d;

   286     y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;

   287   } else {

   288     gdouble t, m, d;

   289 

   290     t = 2.0 * tan (0.5);

   291     m = rp * rp + ip * ip;

   292     d = 4.0 - 4.0 * rp * t + m * t * t;

   293 

   294     x0 = (t * t * iz * iz + 4.0) / d;

   295     x1 = (-8.0 + 2.0 * iz * iz * t * t) / d;

   296     x2 = x0;

   297     y1 = (8.0 - 2.0 * m * t * t) / d;

   298     y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;

   299   }

   300 

   301   /* Convert from lowpass at frequency 1 to either lowpass

   302    * or highpass.

   303    *

   304    * For lowpass substitute z^(-1) with:

   305    *  -1

   306    * z   - k

   307    * ------------

   308    *          -1

   309    * 1 - k * z

   310    *

   311    * k = sin((1-w)/2) / sin((1+w)/2)

   312    *

   313    * For highpass substitute z^(-1) with:

   314    *

   315    *   -1

   316    * -z   - k

   317    * ------------

   318    *          -1

   319    * 1 + k * z

   320    *

   321    * k = -cos((1+w)/2) / cos((1-w)/2)

   322    *

   323    */

   324   {

   325     gdouble k, d;

   326     gdouble omega =

   327         2.0 * M_PI * (filter->cutoff / GST_AUDIO_FILTER (filter)->format.rate);

   328 

   329     if (filter->mode == MODE_LOW_PASS)

   330       k = sin ((1.0 - omega) / 2.0) / sin ((1.0 + omega) / 2.0);

   331     else

   332       k = -cos ((omega + 1.0) / 2.0) / cos ((omega - 1.0) / 2.0);

   333 

   334     d = 1.0 + y1 * k - y2 * k * k;

   335     *a0 = (x0 + k * (-x1 + k * x2)) / d;

   336     *a1 = (x1 + k * k * x1 - 2.0 * k * (x0 + x2)) / d;

   337     *a2 = (x0 * k * k - x1 * k + x2) / d;

   338     *b1 = (2.0 * k + y1 + y1 * k * k - 2.0 * y2 * k) / d;

   339     *b2 = (-k * k - y1 * k + y2) / d;

   340 

   341     if (filter->mode == MODE_HIGH_PASS) {

   342       *a1 = -*a1;

   343       *b1 = -*b1;

   344     }

   345   }

   346 }

   347 

   348 static void

   349 generate_coefficients (GstAudioChebLimit * filter)

   350 {

   351   if (GST_AUDIO_FILTER (filter)->format.rate == 0) {

   352     gdouble *a = g_new0 (gdouble, 1);

   353 

   354     a[0] = 1.0;

   355     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER

   356         (filter), a, 1, NULL, 0);

   357 

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

   359     return;

   360   }

   361 

   362   if (filter->cutoff >= GST_AUDIO_FILTER (filter)->format.rate / 2.0) {

   363     gdouble *a = g_new0 (gdouble, 1);

   364 

   365     a[0] = (filter->mode == MODE_LOW_PASS) ? 1.0 : 0.0;

   366     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER

   367         (filter), a, 1, NULL, 0);

   368     GST_LOG_OBJECT (filter, "cutoff was higher than nyquist frequency");

   369     return;

   370   } else if (filter->cutoff <= 0.0) {

   371     gdouble *a = g_new0 (gdouble, 1);

   372 

   373     a[0] = (filter->mode == MODE_LOW_PASS) ? 0.0 : 1.0;

   374     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER

   375         (filter), a, 1, NULL, 0);

   376     GST_LOG_OBJECT (filter, "cutoff is lower than zero");

   377     return;

   378   }

   379 

   380   /* Calculate coefficients for the chebyshev filter */

   381   {

   382     gint np = filter->poles;

   383     gdouble *a, *b;

   384     gint i, p;

   385 

   386     a = g_new0 (gdouble, np + 3);

   387     b = g_new0 (gdouble, np + 3);

   388 

   389     /* Calculate transfer function coefficients */

   390     a[2] = 1.0;

   391     b[2] = 1.0;

   392 

   393     for (p = 1; p <= np / 2; p++) {

   394       gdouble a0, a1, a2, b1, b2;

   395       gdouble *ta = g_new0 (gdouble, np + 3);

   396       gdouble *tb = g_new0 (gdouble, np + 3);

   397 

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

   399 

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

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

   402 

   403       /* add the new coefficients for the new two poles

   404        * to the cascade by multiplication of the transfer

   405        * functions */

   406       for (i = 2; i < np + 3; i++) {

   407         a[i] = a0 * ta[i] + a1 * ta[i - 1] + a2 * ta[i - 2];

   408         b[i] = tb[i] - b1 * tb[i - 1] - b2 * tb[i - 2];

   409       }

   410       g_free (ta);

   411       g_free (tb);

   412     }

   413 

   414     /* Move coefficients to the beginning of the array

   415      * and multiply the b coefficients with -1 to move from

   416      * the transfer function's coefficients to the difference

   417      * equation's coefficients */

   418     b[2] = 0.0;

   419     for (i = 0; i <= np; i++) {

   420       a[i] = a[i + 2];

   421       b[i] = -b[i + 2];

   422     }

   423 

   424     /* Normalize to unity gain at frequency 0 for lowpass

   425      * and frequency 0.5 for highpass */

   426     {

   427       gdouble gain;

   428 

   429       if (filter->mode == MODE_LOW_PASS)

   430         gain =

   431             gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b, np + 1,

   432             1.0, 0.0);

   433       else

   434         gain =

   435             gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b, np + 1,

   436             -1.0, 0.0);

   437 

   438       for (i = 0; i <= np; i++) {

   439         a[i] /= gain;

   440       }

   441     }

   442 

   443     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER

   444         (filter), a, np + 1, b, np + 1);

   445 

   446     GST_LOG_OBJECT (filter,

   447         "Generated IIR coefficients for the Chebyshev filter");

   448     GST_LOG_OBJECT (filter,

   449         "mode: %s, type: %d, poles: %d, cutoff: %.2f Hz, ripple: %.2f dB",

   450         (filter->mode == MODE_LOW_PASS) ? "low-pass" : "high-pass",

   451         filter->type, filter->poles, filter->cutoff, filter->ripple);

   452     GST_LOG_OBJECT (filter, "%.2f dB gain @ 0 Hz",

   453         20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b,

   454                 np + 1, 1.0, 0.0)));

   455 

   456 #ifndef GST_DISABLE_GST_DEBUG

   457     {

   458       gdouble wc =

   459           2.0 * M_PI * (filter->cutoff /

   460           GST_AUDIO_FILTER (filter)->format.rate);

   461       gdouble zr = cos (wc), zi = sin (wc);

   462 

   463       GST_LOG_OBJECT (filter, "%.2f dB gain @ %d Hz",

   464           20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1,

   465                   b, np + 1, zr, zi)), (int) filter->cutoff);

   466     }

   467 #endif

   468 

   469     GST_LOG_OBJECT (filter, "%.2f dB gain @ %d Hz",

   470         20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b,

   471                 np + 1, -1.0, 0.0)),

   472         GST_AUDIO_FILTER (filter)->format.rate / 2);

   473   }

   474 }

   475 

   476 static void

   477 gst_audio_cheb_limit_finalize (GObject * object)

   478 {

   479   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);

   480 

   481   g_mutex_free (filter->lock);

   482   filter->lock = NULL;

   483 

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

   485 }

   486 

   487 static void

   488 gst_audio_cheb_limit_set_property (GObject * object, guint prop_id,

   489     const GValue * value, GParamSpec * pspec)

   490 {

   491   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);

   492 

   493   switch (prop_id) {

   494     case PROP_MODE:

   495       g_mutex_lock (filter->lock);

   496       filter->mode = g_value_get_enum (value);

   497       generate_coefficients (filter);

   498       g_mutex_unlock (filter->lock);

   499       break;

   500     case PROP_TYPE:

   501       g_mutex_lock (filter->lock);

   502       filter->type = g_value_get_int (value);

   503       generate_coefficients (filter);

   504       g_mutex_unlock (filter->lock);

   505       break;

   506     case PROP_CUTOFF:

   507       g_mutex_lock (filter->lock);

   508       filter->cutoff = g_value_get_float (value);

   509       generate_coefficients (filter);

   510       g_mutex_unlock (filter->lock);

   511       break;

   512     case PROP_RIPPLE:

   513       g_mutex_lock (filter->lock);

   514       filter->ripple = g_value_get_float (value);

   515       generate_coefficients (filter);

   516       g_mutex_unlock (filter->lock);

   517       break;

   518     case PROP_POLES:

   519       g_mutex_lock (filter->lock);

   520       filter->poles = GST_ROUND_UP_2 (g_value_get_int (value));

   521       generate_coefficients (filter);

   522       g_mutex_unlock (filter->lock);

   523       break;

   524     default:

   525       G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);

   526       break;

   527   }

   528 }

   529 

   530 static void

   531 gst_audio_cheb_limit_get_property (GObject * object, guint prop_id,

   532     GValue * value, GParamSpec * pspec)

   533 {

   534   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);

   535 

   536   switch (prop_id) {

   537     case PROP_MODE:

   538       g_value_set_enum (value, filter->mode);

   539       break;

   540     case PROP_TYPE:

   541       g_value_set_int (value, filter->type);

   542       break;

   543     case PROP_CUTOFF:

   544       g_value_set_float (value, filter->cutoff);

   545       break;

   546     case PROP_RIPPLE:

   547       g_value_set_float (value, filter->ripple);

   548       break;

   549     case PROP_POLES:

   550       g_value_set_int (value, filter->poles);

   551       break;

   552     default:

   553       G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);

   554       break;

   555   }

   556 }

   557 

   558 /* GstAudioFilter vmethod implementations */

   559 

   560 static gboolean

   561 gst_audio_cheb_limit_setup (GstAudioFilter * base, GstRingBufferSpec * format)

   562 {

   563   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (base);

   564 

   565   generate_coefficients (filter);

   566 

   567   return GST_AUDIO_FILTER_CLASS (parent_class)->setup (base, format);

   568 }