GNOME.org
 

r3982 - trunk/bse

Author: timj
Date: 2006-10-16 20:14:42 -0400 (Mon, 16 Oct 2006)
New Revision: 3982

Modified:
   trunk/bse/ChangeLog
   trunk/bse/bseellipticfilter.c
   trunk/bse/bseellipticfilter.h
Log:
Tue Oct 17 02:12:24 2006  Tim Janik  <timj gtk org>                                                                                                           
                                                                                                                                                              
        * bseellipticfilter.h, bseellipticfilter.c: moved elliptic integral                                                                                   
        values, chebyshev/elliptic phi and gain into the filter state structure.                                                                              
        eliminated other global variables by declaring them in local scopes.                                                                                  
                                                                                                                                                              


Modified: trunk/bse/ChangeLog
===================================================================
--- trunk/bse/ChangeLog	2006-10-16 23:50:28 UTC (rev 3981)
+++ trunk/bse/ChangeLog	2006-10-17 00:14:42 UTC (rev 3982)
@@ -1,3 +1,9 @@
+Tue Oct 17 02:12:24 2006  Tim Janik  <timj gtk org>
+
+	* bseellipticfilter.h, bseellipticfilter.c: moved elliptic integral
+	values, chebyshev/elliptic phi and gain into the filter state structure.
+	eliminated other global variables by declaring them in local scopes.
+
 Tue Oct 17 01:43:26 2006  Tim Janik  <timj gtk org>
 
 	* bseellipticfilter.h: define MAX_ORDER and MAX_COEFFICIENT_ARRAY_SIZE.

Modified: trunk/bse/bseellipticfilter.c
===================================================================
--- trunk/bse/bseellipticfilter.c	2006-10-16 23:50:28 UTC (rev 3981)
+++ trunk/bse/bseellipticfilter.c	2006-10-17 00:14:42 UTC (rev 3982)
@@ -862,21 +862,14 @@
 static double aa[MAX_COEFFICIENT_ARRAY_SIZE];
 static double pp[MAX_COEFFICIENT_ARRAY_SIZE];
 static Complex z[MAX_COEFFICIENT_ARRAY_SIZE];
-static double Kk = 0.0;
-static double Kpk = 0.0;
-static double rho = 0.0;
-static double phi = 0.0;
 static double sn = 0.0;
 static double cn = 0.0;
 static double dn = 0.0;
 static double sn1 = 0.0;
 static double cn1 = 0.0;
 static double dn1 = 0.0;
-static double phi1 = 0.0;
 static double pn = 0.0;
 static double an = 0.0;
-static double gam = 0.0;
-static double gain = 0.0;
 #endif
 
 /* --- prototypes --- */
@@ -938,10 +931,10 @@
       if (ifr->kind == 2)
         {
           /* For Chebyshev filter, ripples go from 1.0 to 1/sqrt(1+ds->ripple_epsilon^2) */
-          phi = exp (0.5 * ifr->passband_ripple_db / DECIBELL_FACTOR);
+          ds->chebyshev_phi = exp (0.5 * ifr->passband_ripple_db / DECIBELL_FACTOR);
 
           if ((ifr->order & 1) == 0)
-            ds->gain_scale = phi;
+            ds->gain_scale = ds->chebyshev_phi;
           else
             ds->gain_scale = 1.0;
         }
@@ -1186,10 +1179,10 @@
 {
   double k = ds->wc / ds->wr;
   double m = k * k;
-  Kk = ellpk (1.0 - m);
-  Kpk = ellpk (m);
-  EVERBOSE ("check: k=%.20g m=%.20g Kk=%.20g Kpk=%.20g\n", k, m, Kk, Kpk); // BSE info
-  double q = exp (-PI * ifr->order * Kpk / Kk);	/* the nome of k1 */
+  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
+  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;
@@ -1199,27 +1192,27 @@
   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);
+  printf ("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 * Kk / (Kk1 * Kpk);
+  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
   /*   -1
    * sn   j/ds->ripple_epsilon\m  =  j ellik(atan(1/ds->ripple_epsilon), m)
    */
   b = 1.0 / ds->ripple_epsilon;
-  phi = atan (b);
-  double u = ellik (phi, m1p);
-  EVERBOSE ("phi=%.20g m=%.20g u=%.20g\n", phi, m1p, u);
+  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);
   /* consistency check on inverse sn */
-  ellpj (u, m1p, &sn, &cn, &dn, &phi);
+  ellpj (u, m1p, &sn, &cn, &dn, &ds->elliptic_phi);
   a = sn / cn;
   EVERBOSE ("consistency check: sn/cn = %.20g = %.20g = 1/ripple\n", a, b);
   ds->elliptic_k = k;
-  ds->elliptic_u = u * Kk / (ifr->order * Kk1);	/* or, u = u * Kpk / Kpk1 */
+  ds->elliptic_u = u * ds->elliptic_Kk / (ifr->order * Kk1);	/* or, u = u * Kpk / Kpk1 */
   ds->elliptic_m = m;
   return 0;
 }
@@ -1285,16 +1278,16 @@
       /* For Chebyshev, find radii of two Butterworth circles
        * See Gold & Rader, page 60
        */
-      rho = (phi - 1.0)*(phi + 1);  /* rho = ds->ripple_epsilon^2 = {sqrt(1+ds->ripple_epsilon^2)}^2 - 1 */
+      double rho = (ds->chebyshev_phi - 1.0) * (ds->chebyshev_phi + 1);  /* rho = ds->ripple_epsilon^2 = {sqrt(1+ds->ripple_epsilon^2)}^2 - 1 */
       ds->ripple_epsilon = sqrt (rho);
       /* sqrt(1 + 1/ds->ripple_epsilon^2) + 1/ds->ripple_epsilon  = {sqrt(1 + ds->ripple_epsilon^2)  +  1} / ds->ripple_epsilon
        */
-      phi = (phi + 1.0) / ds->ripple_epsilon;
-      EVERBOSE ("Chebychev: phi-before=%.20g ripple=%.20g\n", phi, ds->ripple_epsilon); // BSE info
-      phi = pow (phi, 1.0 / ifr->order);  /* raise to the 1/n power */
-      EVERBOSE ("Chebychev: phi-raised=%.20g rn=%.20g\n", phi, ifr->order * 1.0); // BSE info
-      double b = 0.5 * (phi + 1.0 / phi); /* y coordinates are on this circle */
-      double a = 0.5 * (phi - 1.0 / phi); /* x coordinates are on this circle */
+      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
+      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
+      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;
       if (ifr->order & 1)
         m = 0.0;
@@ -1339,14 +1332,14 @@
     }
   if (ifr->kind == 3)   /* elliptic filter -- stw */
     {
-      double m = ds->elliptic_m;
+      double phi1 = 0.0, m = ds->elliptic_m;
       ds->n_zeros = ifr->order / 2;
       ellpj (ds->elliptic_u, 1.0 - m, &sn1, &cn1, &dn1, &phi1);
       for (i = 0; i < ds->n_zeros; i++)
         {	/* zeros */
           double a = ifr->order - 1 - i - i;
-          double b = (Kk * a) / ifr->order;
-          ellpj (b, m, &sn, &cn, &dn, &phi);
+          double b = (ds->elliptic_Kk * a) / ifr->order;
+          ellpj (b, m, &sn, &cn, &dn, &ds->elliptic_phi);
           int lr = 2 * ds->n_poles + 2 * i;
           zs[lr] = 0.0;
           a = ds->wc / (ds->elliptic_k * sn);	/* elliptic_k = sqrt(m) */
@@ -1355,8 +1348,8 @@
       for (i = 0; i < ds->n_poles; i++)
         {	/* poles */
           double a = ifr->order - 1 - i - i;
-          double b = a * Kk / ifr->order;
-          ellpj (b, m, &sn, &cn, &dn, &phi);
+          double b = a * ds->elliptic_Kk / ifr->order;
+          ellpj (b, m, &sn, &cn, &dn, &ds->elliptic_phi);
           double r = ds->elliptic_k * sn * sn1;
           b = cn1 * cn1 + r * r;
           a = -ds->wc * cn * dn * sn1 * cn1 / b;
@@ -1664,6 +1657,8 @@
   double a = 1.0;
   switch (ifr->type)
     {
+      double gam;
+
     case 3:
       a = -1.0;
       
@@ -1709,12 +1704,12 @@
                                              DesignState                    *ds)
 {
   int j;
-  gain = an/(pn * ds->gain_scale);
+  ds->gain = an / (pn * ds->gain_scale);
   if ((ifr->kind != 3) && (pn == 0))
-    gain = 1.0;
-  printf ("constant gain factor %23.13E\n", gain);
+    ds->gain = 1.0;
+  printf ("constant gain factor %23.13E\n", ds->gain);
   for (j = 0; j <= ds->n_solved_poles; j++)
-    pp[j] = gain * pp[j];
+    pp[j] = ds->gain * pp[j];
   
   printf ("z plane Denominator      Numerator\n");
   for (j = 0; j <= ds->n_solved_poles; j++)
@@ -1813,7 +1808,7 @@
   
   for (f=0; f < limit; f += limit / 21.)
     {
-      double r = response (ifr, ds, f, gain);
+      double r = response (ifr, ds, f, ds->gain);
       if (r <= 0.0)
         r = -999.99;
       else

Modified: trunk/bse/bseellipticfilter.h
===================================================================
--- trunk/bse/bseellipticfilter.h	2006-10-16 23:50:28 UTC (rev 3981)
+++ trunk/bse/bseellipticfilter.h	2006-10-17 00:14:42 UTC (rev 3982)
@@ -56,6 +56,7 @@
   int    n_zeros;
   int    z_counter;	/* incremented as z^N coefficients are found, indexes poles and zeros */
   int    n_solved_poles;
+  /* common state */
   double gain_scale;
   double ripple_epsilon;
   double nyquist_frequency;
@@ -64,10 +65,18 @@
   double cgam; /* angle frequency temporary */
   double stopband_edge; /* derived from ifr->stopband_edge or ifr->stopband_db */
   double wr;
+  /* chebyshev state */
+  double chebyshev_phi;
+  double chebyshev_band_cbp;
+  /* elliptic state */
+  double elliptic_phi;
   double elliptic_k;
   double elliptic_u;
   double elliptic_m;
-  double chebyshev_band_cbp;
+  double elliptic_Kk;  /* complete elliptic integral of the first kind of 1-elliptic_m */
+  double elliptic_Kpk; /* complete elliptic integral of the first kind of elliptic_m */
+  /* common output */
+  double gain;
   double zs[MAX_COEFFICIENT_ARRAY_SIZE];	/* s-plane poles and zeros */
 } DesignState;
 
@@ -84,10 +93,15 @@
   .cgam = 0.0,
   .stopband_edge = 2400,
   .wr = 0.0,
+  .chebyshev_phi = 0.0,
+  .chebyshev_band_cbp = 0.0,
+  .elliptic_phi = 0.0,
   .elliptic_k = 0.0,
   .elliptic_u = 0.0,
   .elliptic_m = 0.0,
-  .chebyshev_band_cbp = 0.0,
+  .elliptic_Kk = 0.0,
+  .elliptic_Kpk = 0.0,
+  .gain = 0.0,
   .zs = { 0, },
 };
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]