Merge branch 'master' into 0.11
[platform/upstream/gst-plugins-good.git] / gst / audiofx / audiocheblimit.c
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 gst_audio_cheb_limit_parent_class parent_class
95 G_DEFINE_TYPE (GstAudioChebLimit,
96     gst_audio_cheb_limit, GST_TYPE_AUDIO_FX_BASE_IIR_FILTER);
97
98 static void gst_audio_cheb_limit_set_property (GObject * object,
99     guint prop_id, const GValue * value, GParamSpec * pspec);
100 static void gst_audio_cheb_limit_get_property (GObject * object,
101     guint prop_id, GValue * value, GParamSpec * pspec);
102 static void gst_audio_cheb_limit_finalize (GObject * object);
103
104 static gboolean gst_audio_cheb_limit_setup (GstAudioFilter * filter,
105     const GstAudioInfo * info);
106
107 enum
108 {
109   MODE_LOW_PASS = 0,
110   MODE_HIGH_PASS
111 };
112
113 #define GST_TYPE_AUDIO_CHEBYSHEV_FREQ_LIMIT_MODE (gst_audio_cheb_limit_mode_get_type ())
114 static GType
115 gst_audio_cheb_limit_mode_get_type (void)
116 {
117   static GType gtype = 0;
118
119   if (gtype == 0) {
120     static const GEnumValue values[] = {
121       {MODE_LOW_PASS, "Low pass (default)",
122           "low-pass"},
123       {MODE_HIGH_PASS, "High pass",
124           "high-pass"},
125       {0, NULL, NULL}
126     };
127
128     gtype = g_enum_register_static ("GstAudioChebLimitMode", values);
129   }
130   return gtype;
131 }
132
133 /* GObject vmethod implementations */
134
135 static void
136 gst_audio_cheb_limit_class_init (GstAudioChebLimitClass * klass)
137 {
138   GObjectClass *gobject_class = (GObjectClass *) klass;
139   GstElementClass *gstelement_class = (GstElementClass *) klass;
140   GstAudioFilterClass *filter_class = (GstAudioFilterClass *) klass;
141
142   GST_DEBUG_CATEGORY_INIT (gst_audio_cheb_limit_debug, "audiocheblimit", 0,
143       "audiocheblimit element");
144
145   gobject_class->set_property = gst_audio_cheb_limit_set_property;
146   gobject_class->get_property = gst_audio_cheb_limit_get_property;
147   gobject_class->finalize = gst_audio_cheb_limit_finalize;
148
149   g_object_class_install_property (gobject_class, PROP_MODE,
150       g_param_spec_enum ("mode", "Mode",
151           "Low pass or high pass mode",
152           GST_TYPE_AUDIO_CHEBYSHEV_FREQ_LIMIT_MODE, MODE_LOW_PASS,
153           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
154   g_object_class_install_property (gobject_class, PROP_TYPE,
155       g_param_spec_int ("type", "Type", "Type of the chebychev filter", 1, 2, 1,
156           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
157
158   /* FIXME: Don't use the complete possible range but restrict the upper boundary
159    * so automatically generated UIs can use a slider without */
160   g_object_class_install_property (gobject_class, PROP_CUTOFF,
161       g_param_spec_float ("cutoff", "Cutoff", "Cut off frequency (Hz)", 0.0,
162           100000.0, 0.0,
163           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
164   g_object_class_install_property (gobject_class, PROP_RIPPLE,
165       g_param_spec_float ("ripple", "Ripple", "Amount of ripple (dB)", 0.0,
166           200.0, 0.25,
167           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
168
169   /* FIXME: What to do about this upper boundary? With a cutoff frequency of
170    * rate/4 32 poles are completely possible, with a cutoff frequency very low
171    * or very high 16 poles already produces only noise */
172   g_object_class_install_property (gobject_class, PROP_POLES,
173       g_param_spec_int ("poles", "Poles",
174           "Number of poles to use, will be rounded up to the next even number",
175           2, 32, 4,
176           G_PARAM_READWRITE | GST_PARAM_CONTROLLABLE | G_PARAM_STATIC_STRINGS));
177
178   gst_element_class_set_details_simple (gstelement_class,
179       "Low pass & high pass filter",
180       "Filter/Effect/Audio",
181       "Chebyshev low pass and high pass filter",
182       "Sebastian Dröge <sebastian.droege@collabora.co.uk>");
183
184   filter_class->setup = GST_DEBUG_FUNCPTR (gst_audio_cheb_limit_setup);
185 }
186
187 static void
188 gst_audio_cheb_limit_init (GstAudioChebLimit * filter)
189 {
190   filter->cutoff = 0.0;
191   filter->mode = MODE_LOW_PASS;
192   filter->type = 1;
193   filter->poles = 4;
194   filter->ripple = 0.25;
195
196   filter->lock = g_mutex_new ();
197 }
198
199 static void
200 generate_biquad_coefficients (GstAudioChebLimit * filter,
201     gint p, gdouble * a0, gdouble * a1, gdouble * a2,
202     gdouble * b1, gdouble * b2)
203 {
204   gint np = filter->poles;
205   gdouble ripple = filter->ripple;
206
207   /* pole location in s-plane */
208   gdouble rp, ip;
209
210   /* zero location in s-plane */
211   gdouble iz = 0.0;
212
213   /* transfer function coefficients for the z-plane */
214   gdouble x0, x1, x2, y1, y2;
215   gint type = filter->type;
216
217   /* Calculate pole location for lowpass at frequency 1 */
218   {
219     gdouble angle = (G_PI / 2.0) * (2.0 * p - 1) / np;
220
221     rp = -sin (angle);
222     ip = cos (angle);
223   }
224
225   /* If we allow ripple, move the pole from the unit
226    * circle to an ellipse and keep cutoff at frequency 1 */
227   if (ripple > 0 && type == 1) {
228     gdouble es, vx;
229
230     es = sqrt (pow (10.0, ripple / 10.0) - 1.0);
231
232     vx = (1.0 / np) * asinh (1.0 / es);
233     rp = rp * sinh (vx);
234     ip = ip * cosh (vx);
235   } else if (type == 2) {
236     gdouble es, vx;
237
238     es = sqrt (pow (10.0, ripple / 10.0) - 1.0);
239     vx = (1.0 / np) * asinh (es);
240     rp = rp * sinh (vx);
241     ip = ip * cosh (vx);
242   }
243
244   /* Calculate inverse of the pole location to convert from
245    * type I to type II */
246   if (type == 2) {
247     gdouble mag2 = rp * rp + ip * ip;
248
249     rp /= mag2;
250     ip /= mag2;
251   }
252
253   /* Calculate zero location for frequency 1 on the
254    * unit circle for type 2 */
255   if (type == 2) {
256     gdouble angle = G_PI / (np * 2.0) + ((p - 1) * G_PI) / (np);
257     gdouble mag2;
258
259     iz = cos (angle);
260     mag2 = iz * iz;
261     iz /= mag2;
262   }
263
264   /* Convert from s-domain to z-domain by
265    * using the bilinear Z-transform, i.e.
266    * substitute s by (2/t)*((z-1)/(z+1))
267    * with t = 2 * tan(0.5).
268    */
269   if (type == 1) {
270     gdouble t, m, d;
271
272     t = 2.0 * tan (0.5);
273     m = rp * rp + ip * ip;
274     d = 4.0 - 4.0 * rp * t + m * t * t;
275
276     x0 = (t * t) / d;
277     x1 = 2.0 * x0;
278     x2 = x0;
279     y1 = (8.0 - 2.0 * m * t * t) / d;
280     y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;
281   } else {
282     gdouble t, m, d;
283
284     t = 2.0 * tan (0.5);
285     m = rp * rp + ip * ip;
286     d = 4.0 - 4.0 * rp * t + m * t * t;
287
288     x0 = (t * t * iz * iz + 4.0) / d;
289     x1 = (-8.0 + 2.0 * iz * iz * t * t) / d;
290     x2 = x0;
291     y1 = (8.0 - 2.0 * m * t * t) / d;
292     y2 = (-4.0 - 4.0 * rp * t - m * t * t) / d;
293   }
294
295   /* Convert from lowpass at frequency 1 to either lowpass
296    * or highpass.
297    *
298    * For lowpass substitute z^(-1) with:
299    *  -1
300    * z   - k
301    * ------------
302    *          -1
303    * 1 - k * z
304    *
305    * k = sin((1-w)/2) / sin((1+w)/2)
306    *
307    * For highpass substitute z^(-1) with:
308    *
309    *   -1
310    * -z   - k
311    * ------------
312    *          -1
313    * 1 + k * z
314    *
315    * k = -cos((1+w)/2) / cos((1-w)/2)
316    *
317    */
318   {
319     gdouble k, d;
320     gdouble omega =
321         2.0 * G_PI * (filter->cutoff / GST_AUDIO_FILTER_RATE (filter));
322
323     if (filter->mode == MODE_LOW_PASS)
324       k = sin ((1.0 - omega) / 2.0) / sin ((1.0 + omega) / 2.0);
325     else
326       k = -cos ((omega + 1.0) / 2.0) / cos ((omega - 1.0) / 2.0);
327
328     d = 1.0 + y1 * k - y2 * k * k;
329     *a0 = (x0 + k * (-x1 + k * x2)) / d;
330     *a1 = (x1 + k * k * x1 - 2.0 * k * (x0 + x2)) / d;
331     *a2 = (x0 * k * k - x1 * k + x2) / d;
332     *b1 = (2.0 * k + y1 + y1 * k * k - 2.0 * y2 * k) / d;
333     *b2 = (-k * k - y1 * k + y2) / d;
334
335     if (filter->mode == MODE_HIGH_PASS) {
336       *a1 = -*a1;
337       *b1 = -*b1;
338     }
339   }
340 }
341
342 static void
343 generate_coefficients (GstAudioChebLimit * filter)
344 {
345   if (GST_AUDIO_FILTER_RATE (filter) == 0) {
346     gdouble *a = g_new0 (gdouble, 1);
347
348     a[0] = 1.0;
349     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
350         (filter), a, 1, NULL, 0);
351
352     GST_LOG_OBJECT (filter, "rate was not set yet");
353     return;
354   }
355
356   if (filter->cutoff >= GST_AUDIO_FILTER_RATE (filter) / 2.0) {
357     gdouble *a = g_new0 (gdouble, 1);
358
359     a[0] = (filter->mode == MODE_LOW_PASS) ? 1.0 : 0.0;
360     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
361         (filter), a, 1, NULL, 0);
362     GST_LOG_OBJECT (filter, "cutoff was higher than nyquist frequency");
363     return;
364   } else if (filter->cutoff <= 0.0) {
365     gdouble *a = g_new0 (gdouble, 1);
366
367     a[0] = (filter->mode == MODE_LOW_PASS) ? 0.0 : 1.0;
368     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
369         (filter), a, 1, NULL, 0);
370     GST_LOG_OBJECT (filter, "cutoff is lower than zero");
371     return;
372   }
373
374   /* Calculate coefficients for the chebyshev filter */
375   {
376     gint np = filter->poles;
377     gdouble *a, *b;
378     gint i, p;
379
380     a = g_new0 (gdouble, np + 3);
381     b = g_new0 (gdouble, np + 3);
382
383     /* Calculate transfer function coefficients */
384     a[2] = 1.0;
385     b[2] = 1.0;
386
387     for (p = 1; p <= np / 2; p++) {
388       gdouble a0, a1, a2, b1, b2;
389       gdouble *ta = g_new0 (gdouble, np + 3);
390       gdouble *tb = g_new0 (gdouble, np + 3);
391
392       generate_biquad_coefficients (filter, p, &a0, &a1, &a2, &b1, &b2);
393
394       memcpy (ta, a, sizeof (gdouble) * (np + 3));
395       memcpy (tb, b, sizeof (gdouble) * (np + 3));
396
397       /* add the new coefficients for the new two poles
398        * to the cascade by multiplication of the transfer
399        * functions */
400       for (i = 2; i < np + 3; i++) {
401         a[i] = a0 * ta[i] + a1 * ta[i - 1] + a2 * ta[i - 2];
402         b[i] = tb[i] - b1 * tb[i - 1] - b2 * tb[i - 2];
403       }
404       g_free (ta);
405       g_free (tb);
406     }
407
408     /* Move coefficients to the beginning of the array
409      * and multiply the b coefficients with -1 to move from
410      * the transfer function's coefficients to the difference
411      * equation's coefficients */
412     b[2] = 0.0;
413     for (i = 0; i <= np; i++) {
414       a[i] = a[i + 2];
415       b[i] = -b[i + 2];
416     }
417
418     /* Normalize to unity gain at frequency 0 for lowpass
419      * and frequency 0.5 for highpass */
420     {
421       gdouble gain;
422
423       if (filter->mode == MODE_LOW_PASS)
424         gain =
425             gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b, np + 1,
426             1.0, 0.0);
427       else
428         gain =
429             gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b, np + 1,
430             -1.0, 0.0);
431
432       for (i = 0; i <= np; i++) {
433         a[i] /= gain;
434       }
435     }
436
437     gst_audio_fx_base_iir_filter_set_coefficients (GST_AUDIO_FX_BASE_IIR_FILTER
438         (filter), a, np + 1, b, np + 1);
439
440     GST_LOG_OBJECT (filter,
441         "Generated IIR coefficients for the Chebyshev filter");
442     GST_LOG_OBJECT (filter,
443         "mode: %s, type: %d, poles: %d, cutoff: %.2f Hz, ripple: %.2f dB",
444         (filter->mode == MODE_LOW_PASS) ? "low-pass" : "high-pass",
445         filter->type, filter->poles, filter->cutoff, filter->ripple);
446     GST_LOG_OBJECT (filter, "%.2f dB gain @ 0 Hz",
447         20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b,
448                 np + 1, 1.0, 0.0)));
449
450 #ifndef GST_DISABLE_GST_DEBUG
451     {
452       gdouble wc =
453           2.0 * G_PI * (filter->cutoff / GST_AUDIO_FILTER_RATE (filter));
454       gdouble zr = cos (wc), zi = sin (wc);
455
456       GST_LOG_OBJECT (filter, "%.2f dB gain @ %d Hz",
457           20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1,
458                   b, np + 1, zr, zi)), (int) filter->cutoff);
459     }
460 #endif
461
462     GST_LOG_OBJECT (filter, "%.2f dB gain @ %d Hz",
463         20.0 * log10 (gst_audio_fx_base_iir_filter_calculate_gain (a, np + 1, b,
464                 np + 1, -1.0, 0.0)), GST_AUDIO_FILTER_RATE (filter) / 2);
465   }
466 }
467
468 static void
469 gst_audio_cheb_limit_finalize (GObject * object)
470 {
471   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);
472
473   g_mutex_free (filter->lock);
474   filter->lock = NULL;
475
476   G_OBJECT_CLASS (parent_class)->finalize (object);
477 }
478
479 static void
480 gst_audio_cheb_limit_set_property (GObject * object, guint prop_id,
481     const GValue * value, GParamSpec * pspec)
482 {
483   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);
484
485   switch (prop_id) {
486     case PROP_MODE:
487       g_mutex_lock (filter->lock);
488       filter->mode = g_value_get_enum (value);
489       generate_coefficients (filter);
490       g_mutex_unlock (filter->lock);
491       break;
492     case PROP_TYPE:
493       g_mutex_lock (filter->lock);
494       filter->type = g_value_get_int (value);
495       generate_coefficients (filter);
496       g_mutex_unlock (filter->lock);
497       break;
498     case PROP_CUTOFF:
499       g_mutex_lock (filter->lock);
500       filter->cutoff = g_value_get_float (value);
501       generate_coefficients (filter);
502       g_mutex_unlock (filter->lock);
503       break;
504     case PROP_RIPPLE:
505       g_mutex_lock (filter->lock);
506       filter->ripple = g_value_get_float (value);
507       generate_coefficients (filter);
508       g_mutex_unlock (filter->lock);
509       break;
510     case PROP_POLES:
511       g_mutex_lock (filter->lock);
512       filter->poles = GST_ROUND_UP_2 (g_value_get_int (value));
513       generate_coefficients (filter);
514       g_mutex_unlock (filter->lock);
515       break;
516     default:
517       G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
518       break;
519   }
520 }
521
522 static void
523 gst_audio_cheb_limit_get_property (GObject * object, guint prop_id,
524     GValue * value, GParamSpec * pspec)
525 {
526   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (object);
527
528   switch (prop_id) {
529     case PROP_MODE:
530       g_value_set_enum (value, filter->mode);
531       break;
532     case PROP_TYPE:
533       g_value_set_int (value, filter->type);
534       break;
535     case PROP_CUTOFF:
536       g_value_set_float (value, filter->cutoff);
537       break;
538     case PROP_RIPPLE:
539       g_value_set_float (value, filter->ripple);
540       break;
541     case PROP_POLES:
542       g_value_set_int (value, filter->poles);
543       break;
544     default:
545       G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
546       break;
547   }
548 }
549
550 /* GstAudioFilter vmethod implementations */
551
552 static gboolean
553 gst_audio_cheb_limit_setup (GstAudioFilter * base, const GstAudioInfo * info)
554 {
555   GstAudioChebLimit *filter = GST_AUDIO_CHEB_LIMIT (base);
556
557   generate_coefficients (filter);
558
559   return GST_AUDIO_FILTER_CLASS (parent_class)->setup (base, info);
560 }