GNOME.org
 

r3973 - in trunk: . tests/bse

Author: stw
Date: 2006-10-16 07:46:51 -0400 (Mon, 16 Oct 2006)
New Revision: 3973

Added:
   trunk/tests/bse/filtertest.cc
Modified:
   trunk/ChangeLog
   trunk/tests/bse/Makefile.am
Log:
Mon Oct 16 13:40:47 2006  Stefan Westerfeld  <stefan space twc de>

	* tests/bse/Makefile.am:
	* tests/bse/filtertest.cc: Added test program for testing IIR filters
	against specifications. It supports both: scanning the actual filter
	transfer function by filtering sine waves, and computing the filter
	transfer function by evaluating H(z).


Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog	2006-10-16 01:05:27 UTC (rev 3972)
+++ trunk/ChangeLog	2006-10-16 11:46:51 UTC (rev 3973)
@@ -1,3 +1,11 @@
+Mon Oct 16 13:40:47 2006  Stefan Westerfeld  <stefan space twc de>
+
+	* tests/bse/Makefile.am:
+	* tests/bse/filtertest.cc: Added test program for testing IIR filters
+	against specifications. It supports both: scanning the actual filter
+	transfer function by filtering sine waves, and computing the filter
+	transfer function by evaluating H(z).
+
 Thu Oct 12 22:11:49 2006  Stefan Westerfeld  <stefan space twc de>
 
 	* slowtests/testresampler.cc: Documentation improvements. Fixed

Modified: trunk/tests/bse/Makefile.am
===================================================================
--- trunk/tests/bse/Makefile.am	2006-10-16 01:05:27 UTC (rev 3972)
+++ trunk/tests/bse/Makefile.am	2006-10-16 11:46:51 UTC (rev 3973)
@@ -50,10 +50,12 @@
 # test programs
 #
 TESTS = cxxbinding
-noinst_PROGRAMS = $(TESTS)
+noinst_PROGRAMS = $(TESTS) filtertest
 progs_ldadd = $(top_builddir)/bse/libbse.la
 cxxbinding_SOURCES = cxxbinding.cc bsecxxapi.cc
 cxxbinding_LDADD = $(progs_ldadd)
+filtertest_SOURCES = filtertest.cc
+filtertest_LDADD = $(progs_ldadd)
 EXTRA_DIST += empty.ogg
 
 #

Added: trunk/tests/bse/filtertest.cc
===================================================================
--- trunk/tests/bse/filtertest.cc	2006-10-16 01:05:27 UTC (rev 3972)
+++ trunk/tests/bse/filtertest.cc	2006-10-16 11:46:51 UTC (rev 3973)
@@ -0,0 +1,262 @@
+/* Tool for testing IIR filters generated by BSE
+ * Copyright (C) 2006 Stefan Westerfeld
+ * 
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser 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.
+ */
+#include <string>
+#include <vector>
+#include <set>
+#include <complex>
+#include <bse/bsemain.h>
+#include <bse/bsemath.h>
+#include <bse/bsemathsignal.h>
+#include <bse/gslfilter.h>
+//#define TEST_VERBOSE
+#include <birnet/birnettests.h>
+
+using std::string;
+using std::vector;
+using std::set;
+using std::min;
+
+class FilterTest
+{
+public:
+  enum TestMode
+  {
+    TEST_NOTHING = 0,
+    TEST_COMPUTED_RESPONSE = 1,
+    TEST_SCANNED_RESPONSE = 2,
+    TEST_COMPUTED_AND_SCANNED_RESPONSE = TEST_COMPUTED_RESPONSE | TEST_SCANNED_RESPONSE
+  };
+
+private:
+  string         m_name;
+  guint          m_order;
+  TestMode       m_test_mode;
+  vector<double> m_a, m_b;
+  set<double>    m_gp_arrows;
+  set<double>    m_gp_lines;
+
+  static const double FS = 10000.0;
+  static const double delta_f = 10000.0 / 1000.0; //997; /* use a prime number to decrease chance of aliasing with mixfreq */
+  static const double MIN_DB = -1000;
+  static const double DB_EPSILON = 0.01;  /* for comparisions */
+
+  double
+  response (double f)
+  {
+    double u = 2.0 * PI * f / FS;
+    std::complex<double> z (cos (u), sin (u)); /* exp( j omega T ) */
+
+    guint o = m_order;
+    std::complex<double> num = m_a[o], den = m_b[o]; 
+    while (o--)
+      {
+	num = num * z + m_a[o];
+	den = den * z + m_b[o];
+      }
+    std::complex<double> w = num / den;
+    return abs (w);
+  }
+
+  double
+  scan_response (double f)
+  {
+    const double MAX_SCAN_FREQ = FS / 2 - 0.01;
+    f = min (f, MAX_SCAN_FREQ);
+    return gsl_filter_sine_scan (m_order, &m_a[0], &m_b[0], f / FS * 2 * PI, 10000);
+  }
+
+public:
+  FilterTest (const char   *name,
+	      guint         order,
+              const double *coefficients,
+	      TestMode      test_mode = TEST_COMPUTED_AND_SCANNED_RESPONSE) :
+    m_name (name),
+    m_order (order),
+    m_test_mode (test_mode)
+  {
+    for (guint i = 0; i < m_order * 2 + 2; i += 2)
+      {
+	m_b.push_back (coefficients[i]);
+	m_a.push_back (coefficients[i+1]);
+      }
+  }
+  void
+  check_response_db (double freq,
+                     double min_resp_db,
+		     double max_resp_db)
+  {
+    if (m_test_mode & TEST_COMPUTED_RESPONSE)
+      {
+	double resp = bse_db_from_factor (response (freq), MIN_DB);
+	if (!(resp > min_resp_db - DB_EPSILON) || !(resp < max_resp_db + DB_EPSILON))
+	  {
+	    g_printerr ("\n*** check_response_db: computed response at frequency %f is %f\n", freq, resp);
+	    g_printerr ("*** check_response_db: but should be in interval [%f..%f]\n", min_resp_db, max_resp_db);
+	  }
+	TASSERT (resp > min_resp_db - DB_EPSILON);
+	TASSERT (resp < max_resp_db + DB_EPSILON);
+      }
+
+    if (m_test_mode & TEST_SCANNED_RESPONSE)
+      {
+	double scan_resp = bse_db_from_factor (scan_response (freq), MIN_DB);
+	if (!(scan_resp > min_resp_db - DB_EPSILON) || !(scan_resp < max_resp_db + DB_EPSILON))
+	  {
+	    g_printerr ("\n*** check_response_db: scanned response at frequency %f is %f\n", freq, scan_resp);
+	    g_printerr ("*** check_response_db: but should be in interval [%f..%f]\n", min_resp_db, max_resp_db);
+	  }
+	TASSERT (scan_resp > min_resp_db - DB_EPSILON);
+	TASSERT (scan_resp < max_resp_db + DB_EPSILON);
+      }
+  }
+  void
+  check_band (double	freq_start,
+	      double	freq_end,
+	      double    min_resp_db,
+	      double	max_resp_db)
+  {
+    g_return_if_fail (freq_start <= freq_end);
+    g_return_if_fail (freq_end <= FS/2);
+
+    TPRINT ("checking band: response in interval [%f..%f] should be in interval [%f..%f] dB\n",
+	freq_start, freq_end, min_resp_db, max_resp_db);
+
+    for (double f = freq_start; f < freq_end; f += delta_f)
+      check_response_db (f, min_resp_db, max_resp_db);
+
+    if (freq_start != freq_end)
+      check_response_db (freq_end, min_resp_db, max_resp_db);
+  }
+  void
+  check_passband (double freq_start,
+		  double freq_end,
+		  double ripple_db)
+  {
+    m_gp_arrows.insert (freq_start);
+    m_gp_arrows.insert (freq_end);
+    m_gp_lines.insert (ripple_db);
+
+    check_band (freq_start, freq_end, ripple_db, 0);
+  }
+  void
+  check_stopband (double freq_start,
+		  double freq_end,
+		  double ripple_db)
+  {
+    m_gp_arrows.insert (freq_start);
+    m_gp_arrows.insert (freq_end);
+    m_gp_lines.insert (ripple_db);
+
+    check_band (freq_start, freq_end, MIN_DB, ripple_db);
+  }
+  /**
+   * creates two files, a datafile named filename_prefix ".data" and a gnuplot
+   * script called filename_prefix ".gp", which can be used to plot the filter,
+   * including the specification checks
+   */
+  bool
+  dump_gp (const string& filename_prefix)
+  {
+    FILE *data_file = fopen ((filename_prefix + ".data").c_str(), "w");
+    if (!data_file)
+      return false;
+    FILE *gp_file = fopen ((filename_prefix + ".gp").c_str(), "w");
+    if (!gp_file)
+      {
+	fclose (data_file);
+	return false;
+      }
+
+    for (double f = 0; f < FS/2; f += delta_f)
+      {
+	fprintf (data_file, "%f %f %f\n", f, bse_db_from_factor (response (f), -1000),
+					     bse_db_from_factor (scan_response (f), -1000));
+      }
+    fprintf (gp_file, "# Test order=%u  norm=%f:\n", m_order,
+	bse_poly_eval (m_order, &m_a[0], 1) / bse_poly_eval (m_order, &m_b[0], 1));
+    /* we don't need H(z) that because we do it ourselves, but we can print it
+     * in case somebody wants to play around with gnuplot
+     */
+    fprintf (gp_file, "H(z)=%s/%s\n", bse_poly_str (m_order, &m_a[0], "z"),
+				      bse_poly_str (m_order, &m_b[0], "z"));
+    fprintf (gp_file, "load '../../bse/tests/filter-defs.gp'\n");
+    fprintf (gp_file, "call '../../bse/tests/arrows.gp' %d", m_gp_arrows.size());
+    for (set<double>::iterator ai = m_gp_arrows.begin(); ai != m_gp_arrows.end(); ai++)
+      fprintf (gp_file, " %f", *ai);
+    fprintf (gp_file, "\n");
+    fprintf (gp_file, "plot 0, '%s.data' using ($1):($3), '%s.data' using ($1):($2) with lines",
+		      filename_prefix.c_str(), filename_prefix.c_str());
+    for (set<double>::iterator li = m_gp_lines.begin(); li != m_gp_lines.end(); li++)
+      fprintf (gp_file, ", %f", *li);
+    fprintf (gp_file, "\n");
+    fprintf (gp_file, "pause -1\n");
+
+    fclose (gp_file);
+    fclose (data_file);
+    return 0;
+  }
+
+};
+
+int
+main (int     argc,
+      char  **argv)
+{
+  bse_init_test (&argc, &argv, NULL);
+
+  TSTART ("Butterworth Lowpass at 2000 Hz, Order 8");
+
+  const double bw8_coefficients[] =
+  {
+    1.000000000E+00, 2.271840012E-03,
+   -1.590566496E+00, 1.817472010E-02,
+    2.083813300E+00, 6.361152035E-02,
+   -1.532625563E+00, 1.272230407E-01,
+    8.694409155E-01, 1.590288009E-01,
+   -3.191759433E-01, 1.272230407E-01,
+    8.209013157E-02, 6.361152035E-02,
+   -1.224667019E-02, 1.817472010E-02,
+    8.613683812E-04, 2.271840012E-03
+  };
+
+  FilterTest bw8 ("Butterworth Lowpass", 8, bw8_coefficients);
+  bw8.check_passband (0, 2000, bse_db_from_factor (1/sqrt(2), -30));
+  bw8.check_stopband (3500, 5000, -68);
+  bw8.dump_gp ("filtertest_bw8");
+
+  TDONE();
+/*
+  const double bw2_coefficients[] =
+  {
+    1.000000000E+00, 2.065720838E-01,
+   -3.695273774E-01, 4.131441677E-01,
+    1.958157127E-01, 2.065720838E-01
+  };
+  FilterTest bw2 ("Butterworth Lowpass", 2, bw2_coefficients, 1);
+  const double bw4_coefficients[] =
+  {
+    1.000000000E+00, 4.658290664E-02,
+   -7.820951980E-01, 1.863316265E-01,
+    6.799785269E-01, 2.794974398E-01,
+   -1.826756978E-01, 1.863316265E-01,
+    3.011887504E-02, 4.658290664E-02
+  };
+  FilterTest bw4 ("Butterworth Lowpass", 4, bw4_coefficients, 1);
+*/
+}
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]