# r3985 - trunk/bse

• From: timj svn gnome org
• To: svn-commits-list gnome org
• Subject: r3985 - trunk/bse
• Date: Tue, 17 Oct 2006 17:07:12 -0400 (EDT)

```Author: timj
Date: 2006-10-17 17:07:11 -0400 (Tue, 17 Oct 2006)
New Revision: 3985

Modified:
trunk/bse/bseellipticfilter.c
Log:
rearranged functions

Modified: trunk/bse/bseellipticfilter.c
===================================================================
--- trunk/bse/bseellipticfilter.c	2006-10-17 21:05:15 UTC (rev 3984)
+++ trunk/bse/bseellipticfilter.c	2006-10-17 21:07:11 UTC (rev 3985)
@@ -809,6 +809,50 @@
}
}

+/* jacobi_theta_by_nome():
+ * Find parameter corresponding to given nome by expansion
+ * in theta functions:
+ * AMS55 #16.38.5, 16.38.7
+ *
+ *       1/2
+ * (2K)                   4     9
+ * (--)     =  1 + 2q + 2q  + 2q  + ...  =  Theta (0,q)
+ * (pi)                                          3
+ *
+ *
+ *       1/2
+ * (2K)     1/4       1/4        2    6    12    20
+ * (--)    m     =  2q    (1 + q  + q  + q   + q   + ...) = Theta (0,q)
+ * (pi)                                                           2
+ *
+ * The nome q(m) = exp(- pi K(1-m)/K(m)).
+ *
+ *                                1/2
+ * Given q, this program returns m   .
+ */
+static double
+jacobi_theta_by_nome (double q)
+{
+  double t1, a = 1.0, b = 1.0, r = 1.0, p = q;
+  do
+    {
+      r *= p;
+      a += 2.0 * r;
+      t1 = fabs (r / a);
+
+      r *= p;
+      b += r;
+      p *= q;
+      double t2 = fabs (r / b);
+      if (t2 > t1)
+	t1 = t2;
+    }
+  while (t1 > MACHEP);
+  a = b / a;
+  a = 4.0 * sqrt (q) * a * a;	/* see above formulas, solved for m */
+  return a;
+}
+
/* --- misc utilities --- */
/* polevl - Evaluate polynomial
*
@@ -852,6 +896,7 @@
}

/* --- prototypes --- */
+#if 0
static int    find_elliptic_locations_in_lambda_plane      (const BseIIRFilterRequirements *ifr,
DesignState                    *ds);
static int    find_s_plane_poles_and_zeros                 (const BseIIRFilterRequirements *ifr,
@@ -871,6 +916,7 @@
static double response                                     (const BseIIRFilterRequirements *ifr,
DesignState                    *ds,
double f, double amp);
+#endif

/* --- functions --- */
static void
@@ -889,269 +935,7 @@
VERBOSE ("%2d %17.9E %17.9E\n", j, ds->cofd[j], ds->cofn[j] * zgain); // BSE info
}

-/* --- main IIR filter design function --- */
-static const char*
-iir_filter_design (const BseIIRFilterRequirements *ifr,
-                   DesignState                    *ds)
-{
-  double passband_edge1 = ifr->passband_edge;
-  double passband_edge0 = ifr->passband_edge2;
-
-  if (ifr->kind <= 0 || ifr->kind > 3)
-    return "unknown kind";
-  if (ifr->type <= 0 || ifr->type > 4)
-    return "unknown type";
-  if (ifr->order <= 0)
-    return "order too small";

-  if (ifr->kind > 1) /* not Butterworth */
-    {
-      if (ifr->passband_ripple_db <= 0.0)
-        return "passband_ripple_db too small";
-      if (ifr->kind == 2)
-        {
-          /* For Chebyshev filter, ripples go from 1.0 to 1/sqrt(1+ds->ripple_epsilon^2) */
-          ds->chebyshev_phi = exp (0.5 * ifr->passband_ripple_db / DECIBELL_FACTOR);
-
-          if ((ifr->order & 1) == 0)
-            ds->gain_scale = ds->chebyshev_phi;
-          else
-            ds->gain_scale = 1.0;
-        }
-      else
-        { /* elliptic */
-          ds->ripple_epsilon = exp (ifr->passband_ripple_db / DECIBELL_FACTOR);
-          ds->gain_scale = 1.0;
-          if ((ifr->order & 1) == 0)
-            ds->gain_scale = sqrt (ds->ripple_epsilon);
-          ds->ripple_epsilon = sqrt (ds->ripple_epsilon - 1.0);
-        }
-    }
-
-  if (ifr->sampling_frequency <= 0.0)
-    return "sampling_frequency too small";
-
-  ds->nyquist_frequency = 0.5 * ifr->sampling_frequency;
-
-  if (passband_edge1 <= 0.0)
-    return "passband_edge1 too small";
-  if (passband_edge1 >= ds->nyquist_frequency)
-    return "passband_edge1 too high";
-
-  if ((ifr->type & 1) == 0)
-    {
-      if (passband_edge0 <= 0.0)
-        return "passband_edge too small";
-      if (passband_edge0 >= ds->nyquist_frequency)
-        return "passband_edge too high";
-    }
-  else
-    {
-      passband_edge0 = 0.0;
-    }
-  if (passband_edge1 < passband_edge0)
-    {
-      double tmp = passband_edge1;
-      passband_edge1 = passband_edge0;
-      passband_edge0 = tmp;
-    }
-  double high_edge, band_width;
-  if (ifr->type == 3)	/* high pass */
-    {
-      band_width = passband_edge1;
-      high_edge = ds->nyquist_frequency;
-    }
-  else
-    {
-      band_width = passband_edge1 - passband_edge0;
-      high_edge = passband_edge1;
-    }
-  /* Frequency correspondence for bilinear transformation
-   *
-   *  Wanalog = tan(2 pi Fdigital T / 2)
-   *
-   * where T = 1/ifr->sampling_frequency
-   */
-  double ang = band_width * PI / ifr->sampling_frequency; /* angle frequency */
-  double sang;
-  double cang = cos (ang);
-  ds->tan_angle_frequency = sin (ang) / cang; /* Wanalog */
-  if (ifr->kind != 3)
-    {
-      ds->wc = ds->tan_angle_frequency;
-      /*printf("cos(1/2 (Whigh-Wlow) T) = %.5e, wc = %.5e\n", cang, ds->wc);*/
-    }
-
-
-  if (ifr->kind == 3)
-    { /* elliptic */
-      double tmp_cgam = cos ((high_edge + passband_edge0) * PI / ifr->sampling_frequency) / cang;
-      ds->cgam = tmp_cgam;
-      if (ifr->stopband_edge > 0.0)
-        ds->stopband_edge = ifr->stopband_edge;
-      else if (ifr->stopband_db >= 0.0)
-        return "need stopband_db or stopband_edge";
-      else /* stopband_db < 0.0 */
-        { /* calculate band edge from db down */
-          double a = exp (-ifr->stopband_db / DECIBELL_FACTOR);
-          double m1 = ds->ripple_epsilon / sqrt (a - 1.0);
-          m1 *= m1;
-          double m1p = 1.0 - m1;
-          double Kk1 = ellpk (m1p);
-          double Kpk1 = ellpk (m1);
-          double q = exp (-PI * Kpk1 / (ifr->order * Kk1));
-          double k = jacobi_theta_by_nome (q);
-          if (ifr->type >= 3)
-            ds->wr = k;
-          else
-            ds->wr = 1.0 / k;
-          if (ifr->type & 1)
-            {
-              ds->stopband_edge = atan (ds->tan_angle_frequency * ds->wr) * ifr->sampling_frequency / PI;
-            }
-          else
-            {
-              // FIXME: using tmp_cgam here increases precision
-              double a = ds->tan_angle_frequency * ds->wr;
-              a *= a;
-              double b = a * (1.0 - ds->cgam * ds->cgam) + a * a;
-              b = (ds->cgam + sqrt (b)) / (1.0 + a);
-              ds->stopband_edge = (PI / 2.0 - asin (b)) * ifr->sampling_frequency / (2.0 * PI);
-            }
-        }
-      switch (ifr->type)
-	{
-	case 1:
-          if (ds->stopband_edge <= passband_edge1)
-            return "need stopband_edge > passband_edge";
-          break;
-	case 2:
-          if (ds->stopband_edge >= passband_edge0 && ds->stopband_edge <= passband_edge1)
-            return "need stopband_edge < passband_edge or stopband_edge > passband_edge2";
-          break;
-	case 3:
-          if (ds->stopband_edge >= passband_edge1)
-            return "need stopband_edge < passband_edge";
-          break;
-	case 4:
-          if (ds->stopband_edge <= passband_edge0)
-            return "need stopband_edge > passband_edge2";
-          if (ds->stopband_edge >= passband_edge1)
-            return "need stopband_edge < passband_edge";
-          break;
-	}
-      ang = ds->stopband_edge * PI / ifr->sampling_frequency;
-      cang = cos (ang);
-      sang = sin (ang);
-
-      if (ifr->type & 1)
-	{
-          ds->wr = sang / (cang * ds->tan_angle_frequency);
-	}
-      else
-	{
-          double q = cang * cang  -  sang * sang;
-          sang = 2.0 * cang * sang;
-          cang = q;
-          ds->wr = (ds->cgam - cang) / (sang * ds->tan_angle_frequency);
-	}
-
-      if (ifr->type >= 3)
-	ds->wr = 1.0 / ds->wr;
-      if (ds->wr < 0.0)
-	ds->wr = -ds->wr;
-
-      const double tmp_y0 = 1.0;
-      double tmp_y1 = ds->wr;
-      /* ds->chebyshev_band_cbp = ds->wr; */
-      if (ifr->type >= 3)
-	tmp_y1 = 1.0 / tmp_y1;
-      if (ifr->type & 1)
-	{
-          int i;
-          for (i = 1; i <= 2; i++)
-            {
-              double tmp_y = i == 1 ? tmp_y0 : tmp_y1;
-              ds->cofd[i] = atan (ds->tan_angle_frequency * tmp_y) * ifr->sampling_frequency / PI ;
-            }
-          printf ("pass band %.9E\n", ds->cofd[1]);
-          printf ("stop band %.9E\n", ds->cofd[2]);
-	}
-      else
-	{
-          int i;
-          for (i = 1; i <= 2; i++)
-            {
-              double tmp_y = i == 1 ? tmp_y0 : tmp_y1;
-              double a = ds->tan_angle_frequency * tmp_y;
-              double b = atan (a);
-              double q = sqrt (1.0 + a * a  -  ds->cgam * ds->cgam);
-              q = atan2 (q, ds->cgam);
-              ds->cofd[i] = (q + b) * ds->nyquist_frequency / PI;
-              ds->cofn[i] = (q - b) * ds->nyquist_frequency / PI;
-            }
-          printf ("pass band %.9E %.9E\n", ds->cofn[1], ds->cofd[1]);
-          printf ("stop band %.9E %.9E\n", ds->cofn[2], ds->cofd[2]);
-	}
-      ds->wc = 1.0;
-      find_elliptic_locations_in_lambda_plane (ifr, ds);	/* find locations in lambda plane */
-      if ((2 * ifr->order + 2) > MAX_COEFFICIENT_ARRAY_SIZE)
-	goto toosml;
-    } /* elliptic */
-
-  /* Transformation from low-pass to band-pass critical frequencies
-   *
-   * Center frequency
-   *                     cos(1/2 (Whigh+Wlow) T)
-   *  cos(Wcenter T) = ----------------------
-   *                     cos(1/2 (Whigh-Wlow) T)
-   *
-   *
-   * Band edges
-   *            cos(Wcenter T) - cos(Wdigital T)
-   *  Wanalog = -----------------------------------
-   *                        sin(Wdigital T)
-   */
-
-  if (ifr->kind == 2)
-    { /* Chebyshev */
-      double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
-      ds->cgam = cos (a) / cang;
-      a = 2.0 * PI * passband_edge1 / ifr->sampling_frequency;
-      ds->chebyshev_band_cbp = (ds->cgam - cos (a)) / sin (a);
-    }
-  if (ifr->kind == 1)
-    { /* Butterworth */
-      double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
-      ds->cgam = cos (a) / cang;
-      a = 2.0 * PI * passband_edge1 / ifr->sampling_frequency;
-      /* 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);
-
-  find_s_plane_poles_and_zeros (ifr, ds);		/* find s plane poles and zeros */
-
-  if (((ifr->type & 1) == 0) && ((4 * ifr->order + 2) > MAX_COEFFICIENT_ARRAY_SIZE))
-    goto toosml;
-
-  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
-  print_z_fraction_before_zplnc (ifr, ds);
-  gainscale_and_print_deno_nume_zeros2_poles2 (ifr, ds);
-  print_filter_table (ifr, ds); /* tabulate transfer function */
-  return NULL;
-
- toosml:
-  return "storage arrays too small";
-}
-
static int
find_elliptic_locations_in_lambda_plane (const BseIIRFilterRequirements *ifr,
DesignState                    *ds)
@@ -1377,50 +1161,6 @@
return 0;
}

-/* jacobi_theta_by_nome():
- * Find parameter corresponding to given nome by expansion
- * in theta functions:
- * AMS55 #16.38.5, 16.38.7
- *
- *       1/2
- * (2K)                   4     9
- * (--)     =  1 + 2q + 2q  + 2q  + ...  =  Theta (0,q)
- * (pi)                                          3
- *
- *
- *       1/2
- * (2K)     1/4       1/4        2    6    12    20
- * (--)    m     =  2q    (1 + q  + q  + q   + q   + ...) = Theta (0,q)
- * (pi)                                                           2
- *
- * The nome q(m) = exp(- pi K(1-m)/K(m)).
- *
- *                                1/2
- * Given q, this program returns m   .
- */
-static double
-jacobi_theta_by_nome (double q)
-{
-  double t1, a = 1.0, b = 1.0, r = 1.0, p = q;
-  do
-    {
-      r *= p;
-      a += 2.0 * r;
-      t1 = fabs (r / a);
-
-      r *= p;
-      b += r;
-      p *= q;
-      double t2 = fabs (r / b);
-      if (t2 > t1)
-	t1 = t2;
-    }
-  while (t1 > MACHEP);
-  a = b / a;
-  a = 4.0 * sqrt (q) * a * a;	/* see above formulas, solved for m */
-  return a;
-}
-
/* convert s plane poles and zeros to the z plane. */
static int
convert_s_plane_to_z_plane (const BseIIRFilterRequirements *ifr,
@@ -1682,49 +1422,6 @@
return 0;
}

-static int
-gainscale_and_print_deno_nume_zeros2_poles2 (const BseIIRFilterRequirements *ifr, /* zplnc */
-                                             DesignState                    *ds)
-{
-  int j;
-  ds->gain = ds->denominator_accu / (ds->numerator_accu * ds->gain_scale);
-  if ((ifr->kind != 3) && (ds->numerator_accu == 0))
-    ds->gain = 1.0;
-  printf ("constant gain factor %23.13E\n", ds->gain);
-  for (j = 0; j <= ds->n_solved_poles; j++)
-    ds->cofn[j] = ds->gain * ds->cofn[j];
-
-  printf ("z plane Denominator      Numerator\n");
-  for (j = 0; j <= ds->n_solved_poles; j++)
-    {
-      printf ("%2d %17.9E %17.9E\n", j, ds->cofd[j], ds->cofn[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");
-  for (j = 0; j < ds->n_solved_poles; j++)
-    {
-      double a = ds->zz[j].r;
-      double b = ds->zz[j].i;
-      if (b >= 0.0)
-        {
-          printf ("pole  %23.13E %23.13E\n", a, b);
-          print_quadratic_factors (ifr, ds, a, b, 1);
-        }
-      int jj = j + ds->n_solved_poles;
-      a = ds->zz[jj].r;
-      b = ds->zz[jj].i;
-      if (b >= 0.0)
-        {
-          printf ("zero  %23.13E %23.13E\n", a, b);
-          print_quadratic_factors (ifr, ds, a, b, 0);
-        }
-    }
-  return 0;
-}
-
static int
@@ -1782,23 +1479,47 @@
return 0;
}

-/* Print table of filter frequency response */
-static void
-print_filter_table (const BseIIRFilterRequirements *ifr,
-                    DesignState                    *ds)
+static int
+gainscale_and_print_deno_nume_zeros2_poles2 (const BseIIRFilterRequirements *ifr, /* zplnc */
+                                             DesignState                    *ds)
{
-  double f, limit = 0.05 * ds->nyquist_frequency * 21;
+  int j;
+  ds->gain = ds->denominator_accu / (ds->numerator_accu * ds->gain_scale);
+  if ((ifr->kind != 3) && (ds->numerator_accu == 0))
+    ds->gain = 1.0;
+  printf ("constant gain factor %23.13E\n", ds->gain);
+  for (j = 0; j <= ds->n_solved_poles; j++)
+    ds->cofn[j] = ds->gain * ds->cofn[j];

-  for (f=0; f < limit; f += limit / 21.)
+  printf ("z plane Denominator      Numerator\n");
+  for (j = 0; j <= ds->n_solved_poles; j++)
{
-      double r = response (ifr, ds, f, ds->gain);
-      if (r <= 0.0)
-        r = -999.99;
-      else
-        r = 2.0 * DECIBELL_FACTOR * log (r);
-      printf ("%10.1f  %10.2f\n", f, r);
-      // f = f + 0.05 * ds->nyquist_frequency;
+      printf ("%2d %17.9E %17.9E\n", j, ds->cofd[j], ds->cofn[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");
+  for (j = 0; j < ds->n_solved_poles; j++)
+    {
+      double a = ds->zz[j].r;
+      double b = ds->zz[j].i;
+      if (b >= 0.0)
+        {
+          printf ("pole  %23.13E %23.13E\n", a, b);
+          print_quadratic_factors (ifr, ds, a, b, 1);
+        }
+      int jj = j + ds->n_solved_poles;
+      a = ds->zz[jj].r;
+      b = ds->zz[jj].i;
+      if (b >= 0.0)
+        {
+          printf ("zero  %23.13E %23.13E\n", a, b);
+          print_quadratic_factors (ifr, ds, a, b, 0);
+        }
+    }
+  return 0;
}

/* Calculate frequency response at f Hz mulitplied by amp */
@@ -1834,7 +1555,288 @@
return u;
}

+/* Print table of filter frequency response */
+static void
+print_filter_table (const BseIIRFilterRequirements *ifr,
+                    DesignState                    *ds)
+{
+  double f, limit = 0.05 * ds->nyquist_frequency * 21;
+
+  for (f=0; f < limit; f += limit / 21.)
+    {
+      double r = response (ifr, ds, f, ds->gain);
+      if (r <= 0.0)
+        r = -999.99;
+      else
+        r = 2.0 * DECIBELL_FACTOR * log (r);
+      printf ("%10.1f  %10.2f\n", f, r);
+      // f = f + 0.05 * ds->nyquist_frequency;
+    }
+}

+/* --- main IIR filter design function --- */
+static const char*
+iir_filter_design (const BseIIRFilterRequirements *ifr,
+                   DesignState                    *ds)
+{
+  double passband_edge1 = ifr->passband_edge;
+  double passband_edge0 = ifr->passband_edge2;
+
+  if (ifr->kind <= 0 || ifr->kind > 3)
+    return "unknown kind";
+  if (ifr->type <= 0 || ifr->type > 4)
+    return "unknown type";
+  if (ifr->order <= 0)
+    return "order too small";
+
+  if (ifr->kind > 1) /* not Butterworth */
+    {
+      if (ifr->passband_ripple_db <= 0.0)
+        return "passband_ripple_db too small";
+      if (ifr->kind == 2)
+        {
+          /* For Chebyshev filter, ripples go from 1.0 to 1/sqrt(1+ds->ripple_epsilon^2) */
+          ds->chebyshev_phi = exp (0.5 * ifr->passband_ripple_db / DECIBELL_FACTOR);
+
+          if ((ifr->order & 1) == 0)
+            ds->gain_scale = ds->chebyshev_phi;
+          else
+            ds->gain_scale = 1.0;
+        }
+      else
+        { /* elliptic */
+          ds->ripple_epsilon = exp (ifr->passband_ripple_db / DECIBELL_FACTOR);
+          ds->gain_scale = 1.0;
+          if ((ifr->order & 1) == 0)
+            ds->gain_scale = sqrt (ds->ripple_epsilon);
+          ds->ripple_epsilon = sqrt (ds->ripple_epsilon - 1.0);
+        }
+    }
+
+  if (ifr->sampling_frequency <= 0.0)
+    return "sampling_frequency too small";
+
+  ds->nyquist_frequency = 0.5 * ifr->sampling_frequency;
+
+  if (passband_edge1 <= 0.0)
+    return "passband_edge1 too small";
+  if (passband_edge1 >= ds->nyquist_frequency)
+    return "passband_edge1 too high";
+
+  if ((ifr->type & 1) == 0)
+    {
+      if (passband_edge0 <= 0.0)
+        return "passband_edge too small";
+      if (passband_edge0 >= ds->nyquist_frequency)
+        return "passband_edge too high";
+    }
+  else
+    {
+      passband_edge0 = 0.0;
+    }
+  if (passband_edge1 < passband_edge0)
+    {
+      double tmp = passband_edge1;
+      passband_edge1 = passband_edge0;
+      passband_edge0 = tmp;
+    }
+  double high_edge, band_width;
+  if (ifr->type == 3)	/* high pass */
+    {
+      band_width = passband_edge1;
+      high_edge = ds->nyquist_frequency;
+    }
+  else
+    {
+      band_width = passband_edge1 - passband_edge0;
+      high_edge = passband_edge1;
+    }
+  /* Frequency correspondence for bilinear transformation
+   *
+   *  Wanalog = tan(2 pi Fdigital T / 2)
+   *
+   * where T = 1/ifr->sampling_frequency
+   */
+  double ang = band_width * PI / ifr->sampling_frequency; /* angle frequency */
+  double sang;
+  double cang = cos (ang);
+  ds->tan_angle_frequency = sin (ang) / cang; /* Wanalog */
+  if (ifr->kind != 3)
+    {
+      ds->wc = ds->tan_angle_frequency;
+      /*printf("cos(1/2 (Whigh-Wlow) T) = %.5e, wc = %.5e\n", cang, ds->wc);*/
+    }
+
+
+  if (ifr->kind == 3)
+    { /* elliptic */
+      double tmp_cgam = cos ((high_edge + passband_edge0) * PI / ifr->sampling_frequency) / cang;
+      ds->cgam = tmp_cgam;
+      if (ifr->stopband_edge > 0.0)
+        ds->stopband_edge = ifr->stopband_edge;
+      else if (ifr->stopband_db >= 0.0)
+        return "need stopband_db or stopband_edge";
+      else /* stopband_db < 0.0 */
+        { /* calculate band edge from db down */
+          double a = exp (-ifr->stopband_db / DECIBELL_FACTOR);
+          double m1 = ds->ripple_epsilon / sqrt (a - 1.0);
+          m1 *= m1;
+          double m1p = 1.0 - m1;
+          double Kk1 = ellpk (m1p);
+          double Kpk1 = ellpk (m1);
+          double q = exp (-PI * Kpk1 / (ifr->order * Kk1));
+          double k = jacobi_theta_by_nome (q);
+          if (ifr->type >= 3)
+            ds->wr = k;
+          else
+            ds->wr = 1.0 / k;
+          if (ifr->type & 1)
+            {
+              ds->stopband_edge = atan (ds->tan_angle_frequency * ds->wr) * ifr->sampling_frequency / PI;
+            }
+          else
+            {
+              // FIXME: using tmp_cgam here increases precision
+              double a = ds->tan_angle_frequency * ds->wr;
+              a *= a;
+              double b = a * (1.0 - ds->cgam * ds->cgam) + a * a;
+              b = (ds->cgam + sqrt (b)) / (1.0 + a);
+              ds->stopband_edge = (PI / 2.0 - asin (b)) * ifr->sampling_frequency / (2.0 * PI);
+            }
+        }
+      switch (ifr->type)
+	{
+	case 1:
+          if (ds->stopband_edge <= passband_edge1)
+            return "need stopband_edge > passband_edge";
+          break;
+	case 2:
+          if (ds->stopband_edge >= passband_edge0 && ds->stopband_edge <= passband_edge1)
+            return "need stopband_edge < passband_edge or stopband_edge > passband_edge2";
+          break;
+	case 3:
+          if (ds->stopband_edge >= passband_edge1)
+            return "need stopband_edge < passband_edge";
+          break;
+	case 4:
+          if (ds->stopband_edge <= passband_edge0)
+            return "need stopband_edge > passband_edge2";
+          if (ds->stopband_edge >= passband_edge1)
+            return "need stopband_edge < passband_edge";
+          break;
+	}
+      ang = ds->stopband_edge * PI / ifr->sampling_frequency;
+      cang = cos (ang);
+      sang = sin (ang);
+
+      if (ifr->type & 1)
+	{
+          ds->wr = sang / (cang * ds->tan_angle_frequency);
+	}
+      else
+	{
+          double q = cang * cang  -  sang * sang;
+          sang = 2.0 * cang * sang;
+          cang = q;
+          ds->wr = (ds->cgam - cang) / (sang * ds->tan_angle_frequency);
+	}
+
+      if (ifr->type >= 3)
+	ds->wr = 1.0 / ds->wr;
+      if (ds->wr < 0.0)
+	ds->wr = -ds->wr;
+
+      const double tmp_y0 = 1.0;
+      double tmp_y1 = ds->wr;
+      /* ds->chebyshev_band_cbp = ds->wr; */
+      if (ifr->type >= 3)
+	tmp_y1 = 1.0 / tmp_y1;
+      if (ifr->type & 1)
+	{
+          int i;
+          for (i = 1; i <= 2; i++)
+            {
+              double tmp_y = i == 1 ? tmp_y0 : tmp_y1;
+              ds->cofd[i] = atan (ds->tan_angle_frequency * tmp_y) * ifr->sampling_frequency / PI ;
+            }
+          printf ("pass band %.9E\n", ds->cofd[1]);
+          printf ("stop band %.9E\n", ds->cofd[2]);
+	}
+      else
+	{
+          int i;
+          for (i = 1; i <= 2; i++)
+            {
+              double tmp_y = i == 1 ? tmp_y0 : tmp_y1;
+              double a = ds->tan_angle_frequency * tmp_y;
+              double b = atan (a);
+              double q = sqrt (1.0 + a * a  -  ds->cgam * ds->cgam);
+              q = atan2 (q, ds->cgam);
+              ds->cofd[i] = (q + b) * ds->nyquist_frequency / PI;
+              ds->cofn[i] = (q - b) * ds->nyquist_frequency / PI;
+            }
+          printf ("pass band %.9E %.9E\n", ds->cofn[1], ds->cofd[1]);
+          printf ("stop band %.9E %.9E\n", ds->cofn[2], ds->cofd[2]);
+	}
+      ds->wc = 1.0;
+      find_elliptic_locations_in_lambda_plane (ifr, ds);	/* find locations in lambda plane */
+      if ((2 * ifr->order + 2) > MAX_COEFFICIENT_ARRAY_SIZE)
+	goto toosml;
+    } /* elliptic */
+
+  /* Transformation from low-pass to band-pass critical frequencies
+   *
+   * Center frequency
+   *                     cos(1/2 (Whigh+Wlow) T)
+   *  cos(Wcenter T) = ----------------------
+   *                     cos(1/2 (Whigh-Wlow) T)
+   *
+   *
+   * Band edges
+   *            cos(Wcenter T) - cos(Wdigital T)
+   *  Wanalog = -----------------------------------
+   *                        sin(Wdigital T)
+   */
+
+  if (ifr->kind == 2)
+    { /* Chebyshev */
+      double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
+      ds->cgam = cos (a) / cang;
+      a = 2.0 * PI * passband_edge1 / ifr->sampling_frequency;
+      ds->chebyshev_band_cbp = (ds->cgam - cos (a)) / sin (a);
+    }
+  if (ifr->kind == 1)
+    { /* Butterworth */
+      double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
+      ds->cgam = cos (a) / cang;
+      a = 2.0 * PI * passband_edge1 / ifr->sampling_frequency;
+      /* 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);
+
+  find_s_plane_poles_and_zeros (ifr, ds);		/* find s plane poles and zeros */
+
+  if (((ifr->type & 1) == 0) && ((4 * ifr->order + 2) > MAX_COEFFICIENT_ARRAY_SIZE))
+    goto toosml;
+
+  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
+  print_z_fraction_before_zplnc (ifr, ds);
+  gainscale_and_print_deno_nume_zeros2_poles2 (ifr, ds);
+  print_filter_table (ifr, ds); /* tabulate transfer function */
+  return NULL;
+
+ toosml:
+  return "storage arrays too small";
+}
+
static double
my_getnum (const char *text)
{
@@ -1882,7 +1884,6 @@
return 0;
}

-
/* compile with: gcc -Wall -O2 -g bseellipticfilter.c -lm -o bseellipticfilter
* (use -ffloat-store for ellf.c compatibility)
*/

```