From: Paul Brossier <piem@piem.org>
Date: Mon, 19 Oct 2009 12:45:13 +0000 (+0200)
Subject: src/temporal: derive biquad from filter, use in peakpicker
X-Git-Tag: 0.4.0-beta1~623
X-Git-Url: https://git.aubio.org/?a=commitdiff_plain;h=b01bd4a8331c5e24f88bb0659a1f723fe7bf12ca;p=aubio.git

src/temporal: derive biquad from filter, use in peakpicker
---

diff --git a/src/aubio.h b/src/aubio.h
index 379d6ef1..6e300082 100644
--- a/src/aubio.h
+++ b/src/aubio.h
@@ -71,8 +71,8 @@ extern "C"
 #if HAVE_SAMPLERATE
 #include "temporal/resampler.h"
 #endif /* HAVE_SAMPLERATE */
-#include "temporal/biquad.h"
 #include "temporal/filter.h"
+#include "temporal/biquad.h"
 #include "temporal/a_weighting.h"
 #include "temporal/c_weighting.h"
 #include "spectral/filterbank.h"
diff --git a/src/onset/peakpick.c b/src/onset/peakpick.c
index 17b9a6a6..fefe2359 100644
--- a/src/onset/peakpick.c
+++ b/src/onset/peakpick.c
@@ -20,6 +20,8 @@
 #include "aubio_priv.h"
 #include "fvec.h"
 #include "mathutils.h"
+#include "lvec.h"
+#include "temporal/filter.h"
 #include "temporal/biquad.h"
 #include "onset/peakpick.h"
 
@@ -42,7 +44,7 @@ struct _aubio_peakpicker_t {
 	aubio_pickerfn_t pickerfn;
 
 	/** biquad lowpass filter */
-	aubio_biquad_t * biquad;
+	aubio_filter_t * biquad;
 	/** original onsets */
 	fvec_t * onset_keep;
 	/** modified onsets */
@@ -94,7 +96,7 @@ aubio_peakpicker_do (aubio_peakpicker_t * p, fvec_t * onset, fvec_t * out)
 
   /* filter onset_proc */
   /** \bug filtfilt calculated post+pre times, should be only once !? */
-  //aubio_biquad_do_filtfilt(p->biquad,onset_proc,scratch);
+  aubio_filter_do_filtfilt (p->biquad, onset_proc, scratch);
 
   for (i = 0; i < p->channels; i++) {
     /* calculate mean and median for onset_proc */
@@ -159,14 +161,17 @@ aubio_peakpicker_t * new_aubio_peakpicker(uint_t channels) {
 	t->onset_proc = new_fvec(t->win_post+t->win_pre+1, channels);
 	t->onset_peek = new_fvec(3, channels);
 
-	/* cutoff: low-pass filter cutoff [0.34, 1] */
-	/* t->cutoff=0.34; */
-	t->biquad = new_aubio_biquad(0.1600,0.3200,0.1600,-0.5949,0.2348);
+	/* cutoff: low-pass filter with cutoff reduced frequency at 0.34
+   generated with octave butter function: [b,a] = butter(2, 0.34);
+   */
+  t->biquad = new_aubio_filter_biquad (0.15998789, 0.31997577, 0.15998789,
+    -0.59488894, 0.23484048, channels);
+
 	return t;
 }
 
 void del_aubio_peakpicker(aubio_peakpicker_t * p) {
-	del_aubio_biquad(p->biquad);
+	del_aubio_filter(p->biquad);
 	del_fvec(p->onset_keep);
 	del_fvec(p->onset_proc);
 	del_fvec(p->onset_peek);
diff --git a/src/temporal/biquad.c b/src/temporal/biquad.c
index 7674ccaa..d715b412 100644
--- a/src/temporal/biquad.c
+++ b/src/temporal/biquad.c
@@ -19,96 +19,34 @@
 
 #include "aubio_priv.h"
 #include "fvec.h"
-#include "mathutils.h"
-#include "temporal/biquad.h"
-
-/** \note this file needs to be in double or more less precision would lead to large
- * errors in the output 
- */
-
-struct _aubio_biquad_t {
-  lsmp_t a2;
-  lsmp_t a3;
-  lsmp_t b1;
-  lsmp_t b2;
-  lsmp_t b3;
-  lsmp_t o1;
-  lsmp_t o2;
-  lsmp_t i1;
-  lsmp_t i2;
-};
-
-void aubio_biquad_do(aubio_biquad_t * b, fvec_t * in) {
-  uint_t i,j;
-  lsmp_t i1 = b->i1;
-  lsmp_t i2 = b->i2;
-  lsmp_t o1 = b->o1;
-  lsmp_t o2 = b->o2;
-  lsmp_t a2 = b->a2;
-  lsmp_t a3 = b->a3;
-  lsmp_t b1 = b->b1;
-  lsmp_t b2 = b->b2;
-  lsmp_t b3 = b->b3;
-
-  i=0; // works in mono only !!!
-  //for (i=0;i<in->channels;i++) {
-  for (j = 0; j < in->length; j++) {
-    lsmp_t i0 = in->data[i][j];
-    lsmp_t o0 = b1 * i0 + b2 * i1 + b3 * i2
-      - a2 * o1 - a3 * o2;// + 1e-37;
-    in->data[i][j] = o0;
-    i2 = i1;
-    i1 = i0;
-    o2 = o1;
-    o1 = o0;
+#include "lvec.h"
+#include "temporal/filter.h"
+
+uint_t
+aubio_filter_set_biquad (aubio_filter_t * f, lsmp_t b0, lsmp_t b1, lsmp_t b2,
+    lsmp_t a1, lsmp_t a2) {
+  uint_t order = aubio_filter_get_order (f);
+  lvec_t *bs = aubio_filter_get_feedforward (f);
+  lvec_t *as = aubio_filter_get_feedback (f);
+
+  if (order != 3) {
+    AUBIO_ERROR ("order of biquad filter must be 3, not %d\n", order);
+    return AUBIO_FAIL;
   }
-  b->i2 = i2;
-  b->i1 = i1;
-  b->o2 = o2;
-  b->o1 = o1;
-  //}
-}
-
-void aubio_biquad_do_filtfilt(aubio_biquad_t * b, fvec_t * in, fvec_t * tmp) {
-  uint_t j,i=0;
-  uint_t length = in->length;
-  lsmp_t mir;
-  /* mirroring */
-  mir = 2*in->data[i][0];
-  b->i1 = mir - in->data[i][2];
-  b->i2 = mir - in->data[i][1];
-  /* apply filtering */
-  aubio_biquad_do(b,in);
-  /* invert  */
-  for (j = 0; j < length; j++)
-    tmp->data[i][length-j-1] = in->data[i][j];
-  /* mirror again */
-  mir = 2*tmp->data[i][0];
-  b->i1 = mir - tmp->data[i][2];
-  b->i2 = mir - tmp->data[i][1];
-  /* apply filtering */
-  aubio_biquad_do(b,tmp);
-  /* invert back */
-  for (j = 0; j < length; j++)
-    in->data[i][j] = tmp->data[i][length-j-1];
-}
-
-aubio_biquad_t * new_aubio_biquad(
-    lsmp_t b1, lsmp_t b2, lsmp_t b3, 
-    lsmp_t a2, lsmp_t a3) {
-  aubio_biquad_t * b = AUBIO_NEW(aubio_biquad_t);
-  b->a2 = a2;
-  b->a3 = a3;
-  b->b1 = b1;
-  b->b2 = b2;
-  b->b3 = b3;
-  b->i1 = 0.;
-  b->i2 = 0.;
-  b->o1 = 0.;
-  b->o2 = 0.;
-  return b;
+  bs->data[0][0] = b0;
+  bs->data[0][1] = b1;
+  bs->data[0][2] = b2;
+  as->data[0][0] = 1.;
+  as->data[0][1] = a1;
+  as->data[0][1] = a2;
+  return AUBIO_OK;
 }
 
-void del_aubio_biquad(aubio_biquad_t * b) {
-  AUBIO_FREE(b);
+aubio_filter_t *
+new_aubio_filter_biquad (lsmp_t b0, lsmp_t b1, lsmp_t b2, lsmp_t a1, lsmp_t a2,
+    uint_t channels)
+{
+  aubio_filter_t *f = new_aubio_filter (3, channels);
+  aubio_filter_set_biquad (f, b0, b1, b2, a1, a2);
+  return f;
 }
diff --git a/src/temporal/biquad.h b/src/temporal/biquad.h
index 2ee3a810..94f6fd9a 100644
--- a/src/temporal/biquad.h
+++ b/src/temporal/biquad.h
@@ -26,52 +26,45 @@
 
   This file implements a normalised biquad filter (second order IIR):
  
-  \f$ y[n] = b_1 x[n] + b_2 x[n-1] + b_3 x[n-2] - a_2 y[n-1] - a_3 y[n-2] \f$
+  \f$ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2] \f$
 
   The filtfilt version runs the filter twice, forward and backward, to
   compensate the phase shifting of the forward operation.
 
+  See also <a href="http://en.wikipedia.org/wiki/Digital_biquad_filter">Digital
+  biquad filter</a> on wikipedia.
+
 */
 
 #ifdef __cplusplus
 extern "C" {
 #endif
 
-/** biquad filter object */
-typedef struct _aubio_biquad_t aubio_biquad_t;
-
-/** filter input vector
+/** set coefficients of a biquad filter
 
-  \param b biquad object as returned by new_aubio_biquad
-  \param in input vector to filter
+  \param b0 forward filter coefficient
+  \param b1 forward filter coefficient
+  \param b2 forward filter coefficient
+  \param a1 feedback filter coefficient
+  \param a2 feedback filter coefficient
 
 */
-void aubio_biquad_do(aubio_biquad_t * b, fvec_t * in);
-/** filter input vector forward and backward
-
-  \param b biquad object as returned by new_aubio_biquad
-  \param in input vector to filter
-  \param tmp memory space to use for computation
+uint_t
+aubio_filter_set_biquad (aubio_filter_t * f, lsmp_t b0, lsmp_t b1, lsmp_t b2,
+    lsmp_t a1, lsmp_t a2);
 
-*/
-void aubio_biquad_do_filtfilt(aubio_biquad_t * b, fvec_t * in, fvec_t * tmp);
 /** create new biquad filter
 
+  \param b0 forward filter coefficient
   \param b1 forward filter coefficient
   \param b2 forward filter coefficient
-  \param b3 forward filter coefficient
+  \param a1 feedback filter coefficient
   \param a2 feedback filter coefficient
-  \param a3 feedback filter coefficient
-
-*/
-aubio_biquad_t * new_aubio_biquad(lsmp_t b1, lsmp_t b2, lsmp_t b3, lsmp_t a2, lsmp_t a3);
-
-/** delete biquad filter 
- 
-  \param b biquad object to delete 
 
 */
-void del_aubio_biquad(aubio_biquad_t * b);
+aubio_filter_t *
+new_aubio_filter_biquad (lsmp_t b0, lsmp_t b1, lsmp_t b2, lsmp_t a1, lsmp_t a2,
+    uint_t channels);
 
 #ifdef __cplusplus
 }
diff --git a/swig/aubio.i b/swig/aubio.i
index a314ae5c..3f87c8c2 100644
--- a/swig/aubio.i
+++ b/swig/aubio.i
@@ -84,10 +84,8 @@ aubio_filter_t * new_aubio_filter_c_weighting (uint_t channels, uint_t samplerat
 uint_t aubio_filter_set_c_weighting (aubio_filter_t * b, uint_t samplerate);
 
 /* biquad */
-aubio_biquad_t * new_aubio_biquad(lsmp_t b1, lsmp_t b2, lsmp_t b3, lsmp_t a2, lsmp_t a3);
-void aubio_biquad_do(aubio_biquad_t * b, fvec_t * in);
-void aubio_biquad_do_filtfilt(aubio_biquad_t * b, fvec_t * in, fvec_t * tmp);
-void del_aubio_biquad(aubio_biquad_t * b);
+aubio_filter_t * new_aubio_filter_biquad(lsmp_t b1, lsmp_t b2, lsmp_t b3, lsmp_t a2, lsmp_t a3, uint_t channels);
+uint_t aubio_filter_set_biquad (aubio_filter_t * b, lsmp_t b1, lsmp_t b2, lsmp_t b3, lsmp_t a2, lsmp_t a3);
 
 /* hist */
 aubio_hist_t * new_aubio_hist(smpl_t flow, smpl_t fhig, uint_t nelems, uint_t channels);
diff --git a/tests/src/test-biquad.c b/tests/src/test-biquad.c
index 319e2307..78d73188 100644
--- a/tests/src/test-biquad.c
+++ b/tests/src/test-biquad.c
@@ -5,12 +5,12 @@ int main(){
         uint_t win_s      = 1024;                       /* window size */
         uint_t channels   = 1;                          /* number of channel */
         fvec_t * in       = new_fvec (win_s, channels); /* input buffer */
-        aubio_biquad_t * o = new_aubio_biquad(0.3,0.2,0.1,0.2,0.3);
+        aubio_filter_t * o = new_aubio_filter_biquad(0.3,0.2,0.1,0.2,0.3, channels);
 
-        aubio_biquad_do_filtfilt(o,in,in);
-        aubio_biquad_do(o,in);
+        aubio_filter_do_filtfilt(o,in,in);
+        aubio_filter_do(o,in);
 
-        del_aubio_biquad(o);
+        del_aubio_filter(o);
         del_fvec(in);
         return 0;
 }