r4053 - in trunk/bse: . tests
- From: timj svn gnome org
- To: svn-commits-list gnome org
- Subject: r4053 - in trunk/bse: . tests
- Date: Thu, 2 Nov 2006 16:45:28 -0500 (EST)
Author: timj
Date: 2006-11-02 16:45:13 -0500 (Thu, 02 Nov 2006)
New Revision: 4053
Modified:
trunk/bse/ChangeLog
trunk/bse/bsefilter-ellf.c
trunk/bse/bsefilter.cc
trunk/bse/bsemath.c
trunk/bse/tests/filtertest.cc
Log:
Thu Nov 2 22:44:16 2006 Tim Janik <timj gtk org>
* bsemath.c: temporary fix against ring buffer overflow.
* bsefilter-ellf.c: replaced overly verbose debugging function
implementation by simple printf like macros. moved ellf compat
constant definitions.
* bsefilter.cc: cosmetic changes to bse_iir_filter_request_string().
* tests/filtertest.cc: compare filter zero/pole coefficients
instead of filter transfer function coefficients which are useless
for high order filters. zero/pole pairs for comparsion are
determined through nearest-neighbour searches.
use TABORT_HANDLER() to hook up print_filter_on_abort(), to print
out mismatching filters when tests fail. use T*_CMP() test macros
for more informative printouts. extended generic filter test to
cover more filter properties.
added predesigned filters where ellf miscalculates gain.
Modified: trunk/bse/ChangeLog
===================================================================
--- trunk/bse/ChangeLog 2006-11-02 21:31:58 UTC (rev 4052)
+++ trunk/bse/ChangeLog 2006-11-02 21:45:13 UTC (rev 4053)
@@ -1,3 +1,23 @@
+Thu Nov 2 22:44:16 2006 Tim Janik <timj gtk org>
+
+ * bsemath.c: temporary fix against ring buffer overflow.
+
+ * bsefilter-ellf.c: replaced overly verbose debugging function
+ implementation by simple printf like macros. moved ellf compat
+ constant definitions.
+
+ * bsefilter.cc: cosmetic changes to bse_iir_filter_request_string().
+
+ * tests/filtertest.cc: compare filter zero/pole coefficients
+ instead of filter transfer function coefficients which are useless
+ for high order filters. zero/pole pairs for comparsion are
+ determined through nearest-neighbour searches.
+ use TABORT_HANDLER() to hook up print_filter_on_abort(), to print
+ out mismatching filters when tests fail. use T*_CMP() test macros
+ for more informative printouts. extended generic filter test to
+ cover more filter properties.
+ added predesigned filters where ellf miscalculates gain.
+
Tue Oct 31 13:52:17 2006 Stefan Westerfeld <stefan space twc de>
* tests/filtertest.cc: Raise coefficient comparision epsilon ceps
Modified: trunk/bse/bsefilter-ellf.c
===================================================================
--- trunk/bse/bsefilter-ellf.c 2006-11-02 21:31:58 UTC (rev 4052)
+++ trunk/bse/bsefilter-ellf.c 2006-11-02 21:45:13 UTC (rev 4053)
@@ -118,28 +118,7 @@
double zd[BSE_IIR_CARRAY_SIZE]; /* denominator coefficients [order+1] */
} EllfDesignState;
-static void __attribute__ ((__format__ (__printf__, 1, 2)))
-VERBOSE (const char *format,
- ...)
-{
- char buffer[4096];
- va_list args;
- va_start (args, format);
- vsnprintf (buffer, sizeof (buffer), format, args);
- va_end (args);
- printf ("# %s", buffer);
- fflush (stdout);
-}
-
-#define EVERBOSE VERBOSE
-//#define EVERBOSE(...) do{}while (0)
-
-static const char* ellf_filter_design (const BseIIRFilterRequest *ifr,
- EllfDesignState *ds);
-
-
-#ifndef ELLF_BEHAVIOR
-
+#if 0 // precision
//#include "bseieee754.h"
#define PI (3.141592653589793238462643383279502884197) // pi
#define BSE_PI_DIV_2 (1.570796326794896619231321691639751442099) // pi/2
@@ -155,15 +134,51 @@
#define PIO2 (BSE_PI_DIV_2) /* pi/2 */
#define MAXNUM (BSE_DOUBLE_MAX_NORMAL) /* 2**1024*(1-MACHEP) */
+#else
+
+static double DECIBELL_FACTOR = -1;
+static void
+init_constants (void)
+{
+ DECIBELL_FACTOR = 10.0 / log (10.0);
+}
+static const double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
+#undef PI
+static const double PI = 3.14159265358979323846; /* pi */
+static const double PIO2 = 1.57079632679489661923; /* pi/2 */
+static const double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
+#endif
+
+
+
+#if 0
+#define error_printf(...) fprintf (stderr, __VA_ARGS__)
+#define ellf_debugf(...) fprintf (stderr, __VA_ARGS__)
+#define ellf_outputf(...) fprintf (stdout, __VA_ARGS__)
+#define ellf_inputf(...) fprintf (stdout, __VA_ARGS__)
+#else
+#define error_printf(...) while (0) { printf (__VA_ARGS__); }
+#define ellf_debugf(...) while (0) { printf (__VA_ARGS__); }
+#define ellf_outputf(...) while (0) { printf (__VA_ARGS__); }
+#define ellf_inputf(...) while (0) { printf (__VA_ARGS__); }
+#endif
+
+static const char* ellf_filter_design (const BseIIRFilterRequest *ifr,
+ EllfDesignState *ds);
+
+
+#ifndef ELLF_BEHAVIOR
+
bool
_bse_filter_design_ellf (const BseIIRFilterRequest *ifr,
BseIIRFilterDesign *fid)
{
+ init_constants();
EllfDesignState ds = { 0, };
const char *errmsg = ellf_filter_design (ifr, &ds);
fflush (stdout);
fflush (stderr);
- // VERBOSE ("DEBUG: %.20g %.20g %.20g %.20g %.20g\n", a, cos(a), cang, ds->cgam, 0.);
+ // ellf_debugf ("DEBUG: %.20g %.20g %.20g %.20g %.20g\n", a, cos(a), cang, ds->cgam, 0.);
if (errmsg)
return false;
fid->order = ds.n_solved_poles;
@@ -191,22 +206,10 @@
#else
-static double DECIBELL_FACTOR = -1;
-static void
-init_constants (void)
-{
- DECIBELL_FACTOR = 10.0 / log (10.0);
-}
-static const double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
-static const double PI = 3.14159265358979323846; /* pi */
-static const double PIO2 = 1.57079632679489661923; /* pi/2 */
-static const double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
-
-
static double
my_getnum (const char *text)
{
- printf ("%s ? ", text);
+ ellf_inputf ("%s ? ", text);
char s[4096];
if (!fgets (s, sizeof (s), stdin))
exit (0);
@@ -240,14 +243,14 @@
ifr.passband_edge2 = my_getnum ("passband_edge2");
if (ifr.kind == BSE_IIR_FILTER_ELLIPTIC)
ifr.stopband_db = ifr.stopband_edge = my_getnum ("stopband_edge or stopband_db");
- printf ("\n");
+ ellf_inputf ("\n");
const char *errmsg = ellf_filter_design (&ifr, &ds);
fflush (stdout);
fflush (stderr);
- // VERBOSE ("DEBUG: %.20g %.20g %.20g %.20g %.20g\n", a, cos(a), cang, ds->cgam, 0.);
+ // ellf_debugf ("DEBUG: %.20g %.20g %.20g %.20g %.20g\n", a, cos(a), cang, ds->cgam, 0.);
if (errmsg)
{
- fprintf (stderr, "Invalid specification: %s\n", errmsg);
+ error_printf (stderr, "Invalid specification: %s\n", errmsg);
fflush (stderr);
return 1;
}
@@ -397,7 +400,7 @@
* which is supposed to be the name of the
* function in which the error occurred:
*/
- printf ("\n%s ", name); // FIXME
+ error_printf ("\n%s ", name); // FIXME
/* Set global error message word */
math_global_error = code;
@@ -407,7 +410,7 @@
*/
if ((code <= 0) || (code >= 7))
code = 0;
- printf ("%s error\n", ermsg[code]);
+ error_printf ("%s error\n", ermsg[code]);
/* Return to calling
* program
@@ -1094,11 +1097,11 @@
zgain = 1.0;
else
zgain = ds->denominator_accu / (ds->numerator_accu * ds->gain_scale);
- VERBOSE ("# constant mygain factor %23.13E\n", zgain); // BSE info
- VERBOSE ("# z plane Denominator Numerator\n"); // BSE info
+ ellf_outputf ("# constant mygain factor %23.13E\n", zgain); // BSE info
+ ellf_outputf ("# z plane Denominator Numerator\n"); // BSE info
int j;
for (j = 0; j <= ds->n_solved_poles; j++)
- VERBOSE ("%2d %17.9E %17.9E\n", j, ds->zd[j], ds->zn[j] * zgain); // BSE info
+ ellf_outputf ("# %2d %17.9E %17.9E\n", j, ds->zd[j], ds->zn[j] * zgain); // BSE info
}
@@ -1110,37 +1113,37 @@
double m = k * k;
ds->elliptic_Kk = ellpk (1.0 - m);
ds->elliptic_Kpk = ellpk (m);
- EVERBOSE ("check: k=%.20g m=%.20g Kk=%.20g Kpk=%.20g\n", k, m, ds->elliptic_Kk, ds->elliptic_Kpk); // BSE info
+ ellf_debugf ("check: k=%.20g m=%.20g Kk=%.20g Kpk=%.20g\n", k, m, ds->elliptic_Kk, ds->elliptic_Kpk); // BSE info
double q = exp (-PI * ifr->order * ds->elliptic_Kpk / ds->elliptic_Kk); /* the nome of k1 */
double m1 = jacobi_theta_by_nome (q); /* see below */
/* Note m1 = ds->ripple_epsilon / sqrt(A*A - 1.0) */
double a = ds->ripple_epsilon / m1;
a = a * a + 1;
a = 10.0 * log (a) / log (10.0);
- printf ("dbdown %.9E\n", a);
+ ellf_debugf ("dbdown %.9E\n", a);
a = 180.0 * asin (k) / PI;
double b = 1.0 / (1.0 + ds->ripple_epsilon * ds->ripple_epsilon);
b = sqrt (1.0 - b);
- printf ("theta=%.9E, rho=%.9E\n", a, b);
+ ellf_debugf ("theta=%.9E, rho=%.9E\n", a, b);
m1 *= m1;
double m1p = 1.0 - m1;
double Kk1 = ellpk (m1p);
double Kpk1 = ellpk (m1);
double r = Kpk1 * ds->elliptic_Kk / (Kk1 * ds->elliptic_Kpk);
- printf ("consistency check: n= %.14E\n", r);
- EVERBOSE ("consistency check: r=%.20g Kpk1=%.20g Kk1=%.20g m1=%.20g m1p=%.20g\n", r, Kpk1, Kk1, m1, m1p); // BSE info
+ ellf_debugf ("consistency check: n= %.14E\n", r);
+ ellf_debugf ("consistency check: r=%.20g Kpk1=%.20g Kk1=%.20g m1=%.20g m1p=%.20g\n", r, Kpk1, Kk1, m1, m1p); // BSE info
/* -1
* sn j/ds->ripple_epsilon\m = j ellik(atan(1/ds->ripple_epsilon), m)
*/
b = 1.0 / ds->ripple_epsilon;
ds->elliptic_phi = atan (b);
double u = ellik (ds->elliptic_phi, m1p);
- EVERBOSE ("phi=%.20g m=%.20g u=%.20g\n", ds->elliptic_phi, m1p, u);
+ ellf_debugf ("phi=%.20g m=%.20g u=%.20g\n", ds->elliptic_phi, m1p, u);
/* consistency check on inverse sn */
double sn, cn, dn;
ellpj (u, m1p, &sn, &cn, &dn, &ds->elliptic_phi);
a = sn / cn;
- EVERBOSE ("consistency check: sn/cn = %.20g = %.20g = 1/ripple\n", a, b);
+ ellf_debugf ("consistency check: sn/cn = %.20g = %.20g = 1/ripple\n", a, b);
ds->elliptic_k = k;
ds->elliptic_u = u * ds->elliptic_Kk / (ifr->order * Kk1); /* or, u = u * Kpk / Kpk1 */
ds->elliptic_m = m;
@@ -1213,9 +1216,9 @@
/* sqrt(1 + 1/ds->ripple_epsilon^2) + 1/ds->ripple_epsilon = {sqrt(1 + ds->ripple_epsilon^2) + 1} / ds->ripple_epsilon
*/
ds->chebyshev_phi = (ds->chebyshev_phi + 1.0) / ds->ripple_epsilon;
- EVERBOSE ("Chebychev: phi-before=%.20g ripple=%.20g\n", ds->chebyshev_phi, ds->ripple_epsilon); // BSE info
+ ellf_debugf ("Chebychev: phi-before=%.20g ripple=%.20g\n", ds->chebyshev_phi, ds->ripple_epsilon); // BSE info
ds->chebyshev_phi = pow (ds->chebyshev_phi, 1.0 / ifr->order); /* raise to the 1/n power */
- EVERBOSE ("Chebychev: phi-raised=%.20g rn=%.20g\n", ds->chebyshev_phi, ifr->order * 1.0); // BSE info
+ ellf_debugf ("Chebychev: phi-raised=%.20g rn=%.20g\n", ds->chebyshev_phi, ifr->order * 1.0); // BSE info
double b = 0.5 * (ds->chebyshev_phi + 1.0 / ds->chebyshev_phi); /* y coordinates are on this circle */
double a = 0.5 * (ds->chebyshev_phi - 1.0 / ds->chebyshev_phi); /* x coordinates are on this circle */
double m;
@@ -1312,7 +1315,7 @@
}
}
}
- printf ("s plane poles:\n");
+ ellf_outputf ("s plane poles:\n");
j = 0;
for (i = 0; i < ds->n_poles + ds->n_zeros; i++)
{
@@ -1320,9 +1323,9 @@
++j;
double b = spz[j];
++j;
- printf ("%.9E %.9E\n", a, b);
+ ellf_outputf ("%.9E %.9E\n", a, b);
if (i == ds->n_poles - 1)
- printf ("s plane zeros:\n");
+ ellf_outputf ("s plane zeros:\n");
}
return 0;
}
@@ -1481,14 +1484,14 @@
{
if (ifr->type != BSE_IIR_FILTER_HIGH_PASS)
{
- printf ("adding zero at Nyquist frequency\n");
+ ellf_debugf ("adding zero at Nyquist frequency\n");
ds->z_counter += 1;
ds->zcpz[ds->z_counter].r = -1.0; /* zero at Nyquist frequency */
ds->zcpz[ds->z_counter].i = 0.0;
}
if (ifr->type == BSE_IIR_FILTER_BAND_PASS || ifr->type == BSE_IIR_FILTER_HIGH_PASS)
{
- printf ("adding zero at 0 Hz\n");
+ ellf_debugf ("adding zero at 0 Hz\n");
ds->z_counter += 1;
ds->zcpz[ds->z_counter].r = 1.0; /* zero at 0 Hz */
ds->zcpz[ds->z_counter].i = 0.0;
@@ -1510,7 +1513,7 @@
}
}
}
- printf ("order = %d\n", ds->n_solved_poles);
+ ellf_outputf ("order = %d\n", ds->n_solved_poles);
/* Expand the poles and zeros into numerator and
* denominator polynomials
@@ -1605,7 +1608,7 @@
a = -x;
b = 0.0;
}
- printf ("q. f.\nz**2 %23.13E\nz**1 %23.13E\n", b, a);
+ ellf_outputf ("q. f.\nz**2 %23.13E\nz**1 %23.13E\n", b, a);
if (b != 0.0)
{
/* resonant frequency */
@@ -1639,7 +1642,7 @@
else
g = MAXNUM;
}
- printf ("f0 %16.8E gain %12.4E DC gain %12.4E\n\n", f, g, g0);
+ ellf_outputf ("f0 %16.8E gain %12.4E DC gain %12.4E\n\n", f, g, g0);
return 0;
}
@@ -1651,34 +1654,34 @@
ds->gain = ds->denominator_accu / (ds->numerator_accu * ds->gain_scale);
if ((ifr->kind == BSE_IIR_FILTER_BUTTERWORTH || ifr->kind == BSE_IIR_FILTER_CHEBYSHEV1) && ds->numerator_accu == 0)
ds->gain = 1.0;
- printf ("constant gain factor %23.13E\n", ds->gain);
+ ellf_outputf ("constant gain factor %23.13E\n", ds->gain);
for (j = 0; j <= ds->n_solved_poles; j++)
ds->zn[j] = ds->gain * ds->zn[j];
- printf ("z plane Denominator Numerator\n");
+ ellf_outputf ("z plane Denominator Numerator\n");
for (j = 0; j <= ds->n_solved_poles; j++)
{
- printf ("%2d %17.9E %17.9E\n", j, ds->zd[j], ds->zn[j]);
+ ellf_outputf ("%2d %17.9E %17.9E\n", j, ds->zd[j], ds->zn[j]);
}
/* I /think/ at this point the polynomial is factorized in 2nd order filters,
* so that it can be implemented without stability problems -- stw
*/
- printf ("poles and zeros with corresponding quadratic factors\n");
+ ellf_outputf ("poles and zeros with corresponding quadratic factors\n");
for (j = 0; j < ds->n_solved_poles; j++)
{
double a = ds->zcpz[j].r;
double b = ds->zcpz[j].i;
if (b >= 0.0)
{
- printf ("pole %23.13E %23.13E\n", a, b);
+ ellf_outputf ("pole %23.13E %23.13E\n", a, b);
print_quadratic_factors (ifr, ds, a, b, true);
}
a = ds->zcpz[ds->n_solved_poles + j].r;
b = ds->zcpz[ds->n_solved_poles + j].i;
if (b >= 0.0)
{
- printf ("zero %23.13E %23.13E\n", a, b);
+ ellf_outputf ("zero %23.13E %23.13E\n", a, b);
print_quadratic_factors (ifr, ds, a, b, false);
}
}
@@ -1732,7 +1735,7 @@
r = -999.99;
else
r = 2.0 * DECIBELL_FACTOR * log (r);
- printf ("%10.1f %10.2f\n", f, r);
+ ellf_debugf ("%10.1f %10.2f\n", f, r);
// f = f + 0.05 * ds->nyquist_frequency;
}
}
@@ -1832,7 +1835,7 @@
if (ifr->kind == BSE_IIR_FILTER_BUTTERWORTH || ifr->kind == BSE_IIR_FILTER_CHEBYSHEV1)
{
ds->wc = ds->tan_angle_frequency;
- /*printf("cos(1/2 (Whigh-Wlow) T) = %.5e, wc = %.5e\n", cang, ds->wc);*/
+ /*ellf_debugf("cos(1/2 (Whigh-Wlow) T) = %.5e, wc = %.5e\n", cang, ds->wc);*/
}
if (ifr->kind == BSE_IIR_FILTER_ELLIPTIC)
@@ -1926,8 +1929,8 @@
double tmp_y = i == 1 ? tmp_y0 : tmp_y1;
ds->zd[i] = atan (ds->tan_angle_frequency * tmp_y) * ifr->sampling_frequency / PI ;
}
- printf ("pass band %.9E\n", ds->zd[1]);
- printf ("stop band %.9E\n", ds->zd[2]);
+ ellf_debugf ("pass band %.9E\n", ds->zd[1]);
+ ellf_debugf ("stop band %.9E\n", ds->zd[2]);
}
else
{
@@ -1942,8 +1945,8 @@
ds->zd[i] = (q + b) * ds->nyquist_frequency / PI;
ds->zn[i] = (q - b) * ds->nyquist_frequency / PI;
}
- printf ("pass band %.9E %.9E\n", ds->zn[1], ds->zd[1]);
- printf ("stop band %.9E %.9E\n", ds->zn[2], ds->zd[2]);
+ ellf_debugf ("pass band %.9E %.9E\n", ds->zn[1], ds->zd[1]);
+ ellf_debugf ("stop band %.9E %.9E\n", ds->zn[2], ds->zd[2]);
}
ds->wc = 1.0;
find_elliptic_locations_in_lambda_plane (ifr, ds); /* find locations in lambda plane */
@@ -1980,12 +1983,12 @@
/* ds->chebyshev_band_cbp = (ds->cgam - cos (a)) / sin (a); */
ds->gain_scale = 1.0;
}
-
- EVERBOSE ("State: gain_scale=%.20g ripple_epsilon=%.20g nyquist_frequency=%.20g " // BSE info
- "tan_angle_frequency=%.20g stopband_edge=%.20g wc=%.20g wr=%.20g cgam=%.20g\n",
- ds->gain_scale, ds->ripple_epsilon, ds->nyquist_frequency,
- ds->tan_angle_frequency, ds->stopband_edge, ds->wc, ds->wr, ds->cgam);
-
+
+ ellf_debugf ("State: gain_scale=%.20g ripple_epsilon=%.20g nyquist_frequency=%.20g " // BSE info
+ "tan_angle_frequency=%.20g stopband_edge=%.20g wc=%.20g wr=%.20g cgam=%.20g\n",
+ ds->gain_scale, ds->ripple_epsilon, ds->nyquist_frequency,
+ ds->tan_angle_frequency, ds->stopband_edge, ds->wc, ds->wr, ds->cgam);
+
find_s_plane_poles_and_zeros (ifr, ds); /* find s plane poles and zeros */
if ((ifr->type == BSE_IIR_FILTER_BAND_PASS || ifr->type == BSE_IIR_FILTER_BAND_STOP) && 4 * ifr->order + 2 > BSE_IIR_CARRAY_SIZE)
@@ -1994,7 +1997,7 @@
convert_s_plane_to_z_plane (ifr, ds); /* convert s plane to z plane */
// volatile_sink ("x");
z_plane_zeros_poles_to_numerator_denomerator (ifr, ds);
- EVERBOSE ("an=%.20g pn=%.20g scale=%.20g\n", ds->denominator_accu, ds->numerator_accu, ds->gain_scale); // BSE info
+ ellf_debugf ("an=%.20g pn=%.20g scale=%.20g\n", ds->denominator_accu, ds->numerator_accu, ds->gain_scale); // BSE info
print_z_fraction_before_zplnc (ifr, ds);
gainscale_and_print_deno_nume_zeros2_poles2 (ifr, ds);
print_filter_table (ifr, ds); /* tabulate transfer function */
Modified: trunk/bse/bsefilter.cc
===================================================================
--- trunk/bse/bsefilter.cc 2006-11-02 21:31:58 UTC (rev 4052)
+++ trunk/bse/bsefilter.cc 2006-11-02 21:45:13 UTC (rev 4053)
@@ -59,15 +59,15 @@
s += bse_iir_filter_type_string (ifr->type);
s += " order=" + string_from_int (ifr->order);
s += " sample-rate=" + string_from_float (ifr->sampling_frequency);
- s += " passband-edge=" + string_from_float (ifr->passband_edge);
if (ifr->kind == BSE_IIR_FILTER_CHEBYSHEV1 || ifr->kind == BSE_IIR_FILTER_ELLIPTIC)
s += " passband-ripple-db=" + string_from_float (ifr->passband_ripple_db);
+ s += " passband-edge=" + string_from_float (ifr->passband_edge);
if (ifr->type == BSE_IIR_FILTER_BAND_PASS || ifr->type == BSE_IIR_FILTER_BAND_STOP)
s += " passband-edge2=" + string_from_float (ifr->passband_edge2);
+ if (ifr->kind == BSE_IIR_FILTER_ELLIPTIC && ifr->stopband_db < 0)
+ s += " stopband-db=" + string_from_float (ifr->stopband_db);
if (ifr->kind == BSE_IIR_FILTER_ELLIPTIC && ifr->stopband_edge > 0)
s += " stopband-edge=" + string_from_float (ifr->stopband_edge);
- if (ifr->kind == BSE_IIR_FILTER_ELLIPTIC && ifr->stopband_db < 0)
- s += " stopband-db=" + string_from_float (ifr->stopband_db);
return g_strdup (s.c_str());
}
Modified: trunk/bse/bsemath.c
===================================================================
--- trunk/bse/bsemath.c 2006-11-02 21:31:58 UTC (rev 4052)
+++ trunk/bse/bsemath.c 2006-11-02 21:45:13 UTC (rev 4053)
@@ -21,7 +21,7 @@
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
-#define RING_BUFFER_LENGTH (16)
+#define RING_BUFFER_LENGTH (256) // FIXME: simlpy dup strings in the API
#define PRINTF_DIGITS "1270"
#define FLOAT_STRING_SIZE (2048)
Modified: trunk/bse/tests/filtertest.cc
===================================================================
--- trunk/bse/tests/filtertest.cc 2006-11-02 21:31:58 UTC (rev 4052)
+++ trunk/bse/tests/filtertest.cc 2006-11-02 21:45:13 UTC (rev 4053)
@@ -21,6 +21,8 @@
#include <sfi/sfitests.h>
#include <bse/bsefilter.h>
#include <bse/bsemain.h>
+#include <bse/gslfilter.h> // FIXME
+#include <bse/bseglobals.h> // FIXME
#include "topconfig.h"
#include <math.h>
@@ -28,7 +30,57 @@
using std::max;
using std::min;
+static inline double
+sqr (register double a)
+{
+ return a * a;
+}
+
+static inline uint
+complex_find_nearest (const Complex *zp,
+ uint n_zps,
+ const Complex *zps)
+{
+ const Complex z = *zp;
+ uint j = 0;
+ double last = sqr (real (z) - real (zps[j])) + sqr (imag (z) - imag (zps[j]));
+ if (last == 0.0)
+ return j;
+ for (uint i = 1; i < n_zps; i++)
+ {
+ double dist2 = sqr (real (z) - real (zps[i])) + sqr (imag (z) - imag (zps[i]));
+ if (dist2 < last)
+ {
+ last = dist2;
+ j = i;
+ if (last == 0.0)
+ return j;
+ }
+ }
+ return j;
+}
+
static double
+compare_zeros (uint n_zeros,
+ const Complex *czeros,
+ const double *rizeros)
+{
+ Complex *z2 = (Complex*) alloca (sizeof (Complex) * n_zeros);
+ for (uint i = 0; i < n_zeros; i++)
+ z2[i] = Complex (rizeros[i * 2], rizeros[i * 2 + 1]);
+ double max_eps = 0;
+ for (uint i = 0; i < n_zeros; i++)
+ {
+ uint j = complex_find_nearest (&czeros[i], n_zeros - i, z2);
+ double reps = fabs (real (z2[j]) - real (czeros[i]));
+ double ieps = fabs (imag (z2[j]) - imag (czeros[i]));
+ z2[j] = z2[n_zeros - i - 1]; // optimization of: swap (z2[last], z2[j]);
+ max_eps = max (max_eps, max (reps, ieps));
+ }
+ return max_eps;
+}
+
+static double
compare_coefficients (const BseIIRFilterDesign *fdes,
uint n_dncoeffs,
const double *dncoeffs)
@@ -154,15 +206,15 @@
}
static void
-exit_with_iir_filter_gnuplot (const BseIIRFilterRequest *fireq,
- const BseIIRFilterDesign *fdes,
- const char *fname,
- double passband_ripple_db = NAN,
- double passband_edge = NAN,
- double passband_edge2 = NAN,
- double stopband_db = NAN,
- double stopband_edge = NAN,
- double stopband_edge2 = NAN)
+noexit_dump_iir_filter_gnuplot (const BseIIRFilterRequest *fireq,
+ const BseIIRFilterDesign *fdes,
+ const char *fname,
+ double passband_ripple_db = NAN,
+ double passband_edge = NAN,
+ double passband_edge2 = NAN,
+ double stopband_db = NAN,
+ double stopband_edge = NAN,
+ double stopband_edge2 = NAN)
{
double consts[64];
double arrows[64];
@@ -185,6 +237,20 @@
g_printerr ("Filter: %s\n", bse_iir_filter_request_string (fireq));
g_printerr ("Design: %s\n", bse_iir_filter_design_string (fdes));
g_printerr ("GnuplotDump: wrote %s.gp and %s.dat use: gnuplot %s.gp\n", fname, fname, fname);
+}
+
+static void
+exit_with_iir_filter_gnuplot (const BseIIRFilterRequest *fireq,
+ const BseIIRFilterDesign *fdes,
+ const char *fname,
+ double passband_ripple_db = NAN,
+ double passband_edge = NAN,
+ double passband_edge2 = NAN,
+ double stopband_db = NAN,
+ double stopband_edge = NAN,
+ double stopband_edge2 = NAN)
+{
+ noexit_dump_iir_filter_gnuplot (fireq, fdes, fname, passband_ripple_db, passband_edge, passband_edge2, stopband_db, stopband_edge, stopband_edge2);
exit (0);
}
@@ -270,9 +336,8 @@
double start_freq,
double end_freq)
{
- double res1 = max_band_damping_ztrans (fdes, start_freq, end_freq);
- double res2 = max_band_damping_zp (fdes, start_freq, end_freq);
- return min (res1, res2);
+ // double res1 = max_band_damping_ztrans (fdes, start_freq, end_freq);
+ return max_band_damping_zp (fdes, MIN (start_freq, end_freq), MAX (start_freq, end_freq));
}
static double
@@ -280,12 +345,22 @@
double start_freq,
double end_freq)
{
- double res1 = min_band_damping_ztrans (fdes, start_freq, end_freq);
- double res2 = min_band_damping_zp (fdes, start_freq, end_freq);
- return max (res1, res2);
+ // double res1 = min_band_damping_ztrans (fdes, start_freq, end_freq);
+ return min_band_damping_zp (fdes, MIN (start_freq, end_freq), MAX (start_freq, end_freq));
}
static void
+print_filter_on_abort (void *data)
+{
+ void **adata = (void**) data;
+ const BseIIRFilterRequest *req = (BseIIRFilterRequest*) adata[0];
+ const BseIIRFilterDesign *fdes = (BseIIRFilterDesign*) adata[1];
+ noexit_dump_iir_filter_gnuplot (req, fdes, "tmpfilter",
+ -fabs(req->passband_ripple_db), req->passband_edge, req->passband_edge2,
+ req->stopband_db != 0 ? req->stopband_db : NAN, req->stopband_edge, NAN);
+}
+
+static void
butterwoth_tests ()
{
TSTART ("Butterworth");
@@ -295,6 +370,10 @@
BseIIRFilterDesign fdes;
BseIIRFilterRequest req = { BseIIRFilterKind (0), };
req.kind = BSE_IIR_FILTER_BUTTERWORTH;
+ void *abort_data[2];
+ abort_data[0] = (void*) &req;
+ abort_data[1] = &fdes;
+ TABORT_HANDLER (print_filter_on_abort, abort_data);
TOK();
{
@@ -322,9 +401,9 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 0, 2000) < gaineps);
- TASSERT (max_band_damping (&fdes, 0, 2000) > -3.0103);
- TASSERT (min_band_damping (&fdes, 3500, 5000) < -68);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 2000), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 0, 2000), >, -3.0103);
+ TASSERT_CMP (min_band_damping (&fdes, 3500, 5000), <, -68);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -3.0103, 2000, NAN, -68, 3500);
}
@@ -348,9 +427,9 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 2000, 5000) < gaineps);
- TASSERT (max_band_damping (&fdes, 2000, 5000) > -3.0103);
- TASSERT (min_band_damping (&fdes, 0, 600) < -80);
+ TASSERT_CMP (min_band_damping (&fdes, 2000, 5000), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 2000, 5000), >, -3.0103);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 600), <, -80);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -3.0103, 2000, NAN, -80, 600);
}
@@ -386,10 +465,10 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 1500, 3500) < gaineps);
- TASSERT (max_band_damping (&fdes, 1500, 3500) > -3.0103);
- TASSERT (min_band_damping (&fdes, 0, 1000) < -49.5);
- TASSERT (min_band_damping (&fdes, 4000, 5000) < -49.5);
+ TASSERT_CMP (min_band_damping (&fdes, 1500, 3500), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 1500, 3500), >, -3.0103);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 1000), <, -49.5);
+ TASSERT_CMP (min_band_damping (&fdes, 4000, 5000), <, -49.5);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -3.0103, 1500, 3500, -49.5, 1000, 4000);
}
@@ -435,10 +514,10 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 0, 1000) < gaineps);
- TASSERT (max_band_damping (&fdes, 0, 1000) > -3.0103);
- TASSERT (min_band_damping (&fdes, 1500, 3500) < -77);
- TASSERT (max_band_damping (&fdes, 4000, 5000) > -3.0103);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 1000), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 0, 1000), >, -3.0103);
+ TASSERT_CMP (min_band_damping (&fdes, 1500, 3500), <, -77);
+ TASSERT_CMP (max_band_damping (&fdes, 4000, 5000), >, -3.0103);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -3.0103, 1000, 4000, -77, 1500, 3500);
}
@@ -450,20 +529,24 @@
static void
chebychev1_tests ()
{
- TSTART ("Chebyshev1");
bool success;
double eps;
- const double ceps = 8e-13, gaineps = 1e-7;
+ TSTART ("Chebyshev1");
BseIIRFilterDesign fdes;
BseIIRFilterRequest req = { BseIIRFilterKind (0), };
req.kind = BSE_IIR_FILTER_CHEBYSHEV1;
+ void *abort_data[2];
+ abort_data[0] = (void*) &req;
+ abort_data[1] = &fdes;
+ TABORT_HANDLER (print_filter_on_abort, abort_data);
+ const double ceps = 8e-13, gaineps = 1e-7;
TOK();
{
req.type = BSE_IIR_FILTER_LOW_PASS;
req.order = 8;
req.sampling_frequency = 20000;
- req.passband_ripple_db = 1.55;
+ req.passband_ripple_db = 1.5500;
req.passband_edge = 3000;
success = bse_iir_filter_design (&req, &fdes);
TASSERT (success == true);
@@ -480,9 +563,9 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 0, 3000) < gaineps);
- TASSERT (max_band_damping (&fdes, 0, 3000) > -1.55);
- TASSERT (min_band_damping (&fdes, 5000, 10000) < -80);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 3000), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 0, 3000), >, -1.5501);
+ TASSERT_CMP (min_band_damping (&fdes, 5000, 10000), <, -80);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -1.55, 3000, NAN, -80, 5000);
}
@@ -507,9 +590,9 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 600, 5000) < gaineps);
- TASSERT (max_band_damping (&fdes, 600, 5000) >= -0.101);
- TASSERT (min_band_damping (&fdes, 0, 250) < -70);
+ TASSERT_CMP (min_band_damping (&fdes, 600, 5000), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 600, 5000), >, -0.101);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 250), <, -70);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -0.1, 600, NAN, -70, 250);
}
@@ -548,10 +631,10 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 3500, 9500) < gaineps);
- TASSERT (max_band_damping (&fdes, 3500, 9500) > -1.801);
- TASSERT (min_band_damping (&fdes, 0, 3000) < -55);
- TASSERT (min_band_damping (&fdes, 10200, 15000) < -55);
+ TASSERT_CMP (min_band_damping (&fdes, 3500, 9500), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 3500, 9500), >, -1.801);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 3000), <, -55);
+ TASSERT_CMP (min_band_damping (&fdes, 10200, 15000), <, -55);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -1.801, 3500, 9500, -55, 3000, 10200);
}
@@ -596,10 +679,10 @@
};
eps = compare_coefficients (&fdes, BIRNET_ARRAY_SIZE (dncoeffs), dncoeffs);
TASSERT (eps <= ceps);
- TASSERT (min_band_damping (&fdes, 0, 8000) < gaineps);
- TASSERT (max_band_damping (&fdes, 0, 8000) > -1.001);
- TASSERT (min_band_damping (&fdes, 8500, 11500) < -78);
- TASSERT (max_band_damping (&fdes, 12000, 20000) > -1.001);
+ TASSERT_CMP (min_band_damping (&fdes, 0, 8000), <, gaineps);
+ TASSERT_CMP (max_band_damping (&fdes, 0, 8000), >, -1.001);
+ TASSERT_CMP (min_band_damping (&fdes, 8500, 11500), <, -78);
+ TASSERT_CMP (max_band_damping (&fdes, 12000, 20000), >, -1.001);
if (0)
exit_with_iir_filter_gnuplot (&req, &fdes, "tmpfilter", -1.001, 8000, 12000, -78, 8500, 11500);
}
@@ -615,32 +698,131 @@
const BseIIRFilterRequest *filter_request;
const double *filter_coefficients;
uint n_filter_coefficients;
- } filters[1000000];
+ double gain;
+ uint n_zeros;
+ const double *zeros;
+ uint n_poles;
+ const double *poles;
+ } filters[100000] = { { 0, } };
uint index = 0;
- // #include "../../r+d-files/tmp.c"
+
+ { // ellf gain is too high
+ static const BseIIRFilterRequest filter_request = {
+ /* kind = */ BSE_IIR_FILTER_ELLIPTIC,
+ /* type = */ BSE_IIR_FILTER_HIGH_PASS,
+ /* order = */ 20,
+ /* sampling_frequency = */ 1073741824.0000,
+ /* passband_ripple_db = */ 1.3719436899999999,
+ /* passband_edge = */ 483183820.79999995,
+ /* passband_edge2 = */ 0,
+ /* stopband_edge = */ 0,
+ /* stopband_db = */ -26,
+ };
+ filters[index].filter_request = &filter_request;
+ filters[index].gain = +6.50971316315766546e-02;
+ index++;
+ }
+ { // ellf gain is too high
+ static const BseIIRFilterRequest filter_request = {
+ /* kind = */ BSE_IIR_FILTER_ELLIPTIC,
+ /* type = */ BSE_IIR_FILTER_BAND_PASS,
+ /* order = */ 19,
+ /* sampling_frequency = */ 3333.0000,
+ /* passband_ripple_db = */ 0.01,
+ /* passband_edge = */ 1333.2,
+ /* passband_edge2 = */ 999.89999999999998,
+ /* stopband_edge = */ 0,
+ /* stopband_db = */ -120,
+ };
+ filters[index].filter_request = &filter_request;
+ filters[index].gain = +1.21552945102360275e-05;
+ index++;
+ }
+ { // ellf gain is too high
+ static const BseIIRFilterRequest filter_request = {
+ /* kind = */ BSE_IIR_FILTER_ELLIPTIC,
+ /* type = */ BSE_IIR_FILTER_BAND_PASS,
+ /* order = */ 20,
+ /* sampling_frequency = */ 8000.0000,
+ /* passband_ripple_db = */ 0.01,
+ /* passband_edge = */ 3200,
+ /* passband_edge2 = */ 2400,
+ /* stopband_edge = */ 0,
+ /* stopband_db = */ -26,
+ };
+ filters[index].filter_request = &filter_request;
+ filters[index].gain = +1.12649838606861938e-01;
+ index++;
+ }
+
+ if (0) { // FIXME: bse gain = 4.3037825362077964
+ static const BseIIRFilterRequest filter_request = {
+ /* kind = */ BSE_IIR_FILTER_ELLIPTIC,
+ /* type = */ BSE_IIR_FILTER_BAND_PASS,
+ /* order = */ 36,
+ /* sampling_frequency = */ 3333.0000,
+ /* passband_ripple_db = */ 1.3719436899999999,
+ /* passband_edge = */ 943.23900000000003,
+ /* passband_edge2 = */ 499.94999999999999,
+ /* stopband_edge = */ 0,
+ /* stopband_db = */ -26,
+ };
+ filters[index].filter_request = &filter_request;
+ filters[index].gain = -3.83263619423376167e-01;
+ index++;
+ }
+
+#include "../../r+d-files/tmp.c"
+
uint i;
- const double coefficients_epsilon = 1e-9;
+ const double coefficients_epsilon = 1e-7;
+ const double max_gain = 0.01; // maximum gain in passband
TSTART ("Bruteforce filter checks");
+ void *abort_data[2];
+ TABORT_HANDLER (print_filter_on_abort, abort_data);
+ const uint n_filters_tok = 13;
for (i = 0; i < index; i++)
{
const BseIIRFilterRequest *req = filters[i].filter_request;
- const uint n_coefficients = filters[i].n_filter_coefficients;
- const double *coefficients = filters[i].filter_coefficients;
BseIIRFilterDesign fdes;
+ abort_data[0] = (void*) req;
+ abort_data[1] = &fdes;
bool success = bse_iir_filter_design (req, &fdes);
- TASSERT (success == true);
- double eps = compare_coefficients (&fdes, n_coefficients, coefficients);
- if (eps > coefficients_epsilon)
+ TCHECK (success == true);
+ if (filters[i].zeros)
+ TCHECK_CMP (filters[i].n_zeros, ==, fdes.n_zeros);
+ if (filters[i].poles)
+ TCHECK_CMP (filters[i].n_poles, ==, fdes.n_poles);
+ double zeps = filters[i].zeros ? compare_zeros (fdes.n_zeros, fdes.zz, filters[i].zeros) : 0;
+ double peps = filters[i].poles ? compare_zeros (fdes.n_poles, fdes.zp, filters[i].poles) : 0;
+ TCHECK_CMP (zeps, <, coefficients_epsilon);
+ TCHECK_CMP (peps, <, coefficients_epsilon);
+ /* broad gain check */
+ TCHECK_CMP (min_band_damping (&fdes, 0, 0.5 * req->sampling_frequency), <, max_gain);
+ /* finer gain checks */
+ double passband_ripple_db = req->kind == BSE_IIR_FILTER_BUTTERWORTH ? -3.010299956639811952 : -req->passband_ripple_db;
+ passband_ripple_db -= 0.01; // add a bit of fudge for evaluation artefacts
+ if (req->type == BSE_IIR_FILTER_LOW_PASS || req->type == BSE_IIR_FILTER_BAND_STOP)
{
- g_printerr ("EPSILON FAILED: %.16g %.16g\n", eps, coefficients_epsilon);
- g_printerr ("EPSILON FAILED: %.16g %.16g\n", eps, coefficients_epsilon);
- g_printerr ("EPSILON FAILED: %.16g %.16g\n", eps, coefficients_epsilon);
- exit_with_iir_filter_gnuplot (req, &fdes, "tmpfilter",
- req->passband_ripple_db, req->passband_edge, req->passband_edge2,
- req->stopband_db != 0 ? req->stopband_db : NAN, req->stopband_edge, NAN);
- TASSERT (eps <= coefficients_epsilon);
+ TCHECK_CMP (min_band_damping (&fdes, 0, req->passband_edge), <, max_gain);
+ TCHECK_CMP (max_band_damping (&fdes, 0, req->passband_edge), >, passband_ripple_db);
}
- else
+ if (req->type == BSE_IIR_FILTER_HIGH_PASS)
+ {
+ TCHECK_CMP (min_band_damping (&fdes, req->passband_edge, 0.5 * req->sampling_frequency), <, max_gain);
+ TCHECK_CMP (max_band_damping (&fdes, req->passband_edge, 0.5 * req->sampling_frequency), >, passband_ripple_db);
+ }
+ if (req->type == BSE_IIR_FILTER_BAND_PASS)
+ {
+ TCHECK_CMP (min_band_damping (&fdes, req->passband_edge, req->passband_edge2), <, max_gain);
+ TCHECK_CMP (max_band_damping (&fdes, req->passband_edge, req->passband_edge2), >, passband_ripple_db);
+ }
+ if (req->type == BSE_IIR_FILTER_BAND_STOP)
+ {
+ TCHECK_CMP (min_band_damping (&fdes, req->passband_edge2, 0.5 * req->sampling_frequency), <, max_gain);
+ TCHECK_CMP (max_band_damping (&fdes, req->passband_edge2, 0.5 * req->sampling_frequency), >, passband_ripple_db);
+ }
+ if (i % n_filters_tok == 0)
TOK();
}
TDONE();
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
