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;
-}
-
/* display quadratic factors */
static int
print_quadratic_factors (const BseIIRFilterRequirements *ifr,
@@ -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)
*/
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
