r3978 - trunk/bse
- From: timj svn gnome org
- To: svn-commits-list gnome org
- Subject: r3978 - trunk/bse
- Date: Mon, 16 Oct 2006 18:17:54 -0400 (EDT)
Author: timj
Date: 2006-10-16 18:17:52 -0400 (Mon, 16 Oct 2006)
New Revision: 3978
Modified:
trunk/bse/ChangeLog
trunk/bse/bseellipticfilter.c
Log:
Tue Oct 17 00:17:23 2006 Tim Janik <timj gtk org>
* bseellipticfilter.c: whitespace and paranthese fixes.
Modified: trunk/bse/ChangeLog
===================================================================
--- trunk/bse/ChangeLog 2006-10-16 22:06:39 UTC (rev 3977)
+++ trunk/bse/ChangeLog 2006-10-16 22:17:52 UTC (rev 3978)
@@ -1,3 +1,7 @@
+Tue Oct 17 00:17:23 2006 Tim Janik <timj gtk org>
+
+ * bseellipticfilter.c: whitespace and paranthese fixes.
+
Tue Oct 17 00:05:55 2006 Tim Janik <timj gtk org>
* bseellipticfilter.h, bseellipticfilter.c: white space fixups,
Modified: trunk/bse/bseellipticfilter.c
===================================================================
--- trunk/bse/bseellipticfilter.c 2006-10-16 22:06:39 UTC (rev 3977)
+++ trunk/bse/bseellipticfilter.c 2006-10-16 22:17:52 UTC (rev 3978)
@@ -41,7 +41,7 @@
}
#define EVERBOSE VERBOSE
-//#define EVERBOSE(...) do{}while(0)
+//#define EVERBOSE(...) do{}while (0)
#if 0 // FIXME: increase precision by using:
@@ -67,7 +67,7 @@
static void
init_constants (void)
{
- DECIBELL_FACTOR = 10.0 / log(10.0);
+ DECIBELL_FACTOR = 10.0 / log (10.0);
}
static const double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
static const double PI = 3.14159265358979323846; /* pi */
@@ -175,7 +175,7 @@
* int code;
* int math_set_error();
*
- * math_set_error( fctnam, code );
+ * math_set_error(fctnam, code);
*
* DESCRIPTION:
* This routine may be called to report one of the following
@@ -229,14 +229,14 @@
/* Display error message defined
* by the code argument.
*/
- if( (code <= 0) || (code >= 7) )
+ if ((code <= 0) || (code >= 7))
code = 0;
- printf( "%s error\n", ermsg[code] );
+ printf ("%s error\n", ermsg[code]);
/* Return to calling
* program
*/
- return( 0 );
+ return 0;
}
/* --- complex number arithmetic --- */
@@ -284,14 +284,14 @@
double p = b->r * a->r + b->i * a->i;
double q = b->i * a->r - b->r * a->i;
- if( y < 1.0 )
+ if (y < 1.0)
{
double w = MAXNUM * y;
- if( (fabs(p) > w) || (fabs(q) > w) || (y == 0.0) )
+ if ((fabs (p) > w) || (fabs (q) > w) || (y == 0.0))
{
c->r = MAXNUM;
c->i = MAXNUM;
- math_set_error( "Cdiv", MATH_ERROR_OVERFLOW );
+ math_set_error ("Cdiv", MATH_ERROR_OVERFLOW);
return;
}
}
@@ -329,12 +329,12 @@
* Complex z;
* double a;
*
- * a = Cabs( &z );
+ * a = Cabs (&z);
*
* DESCRIPTION:
* If z = x + iy
* then
- * a = sqrt( x**2 + y**2 ).
+ * a = sqrt (x**2 + y**2).
* Overflow and underflow are avoided by testing the magnitudes
* of x and y before squaring. If either is outside half of
* the floating point full scale range, both are rescaled.
@@ -356,61 +356,61 @@
double x, y, b, re, im;
int ex, ey, e;
- /* Note, Cabs(INFINITY,NAN) = INFINITY. */
- if( z->r == INFINITY || z->i == INFINITY
- || z->r == -INFINITY || z->i == -INFINITY )
- return( INFINITY );
+ /* Note, Cabs (INFINITY,NAN) = INFINITY. */
+ if (z->r == INFINITY || z->i == INFINITY
+ || z->r == -INFINITY || z->i == -INFINITY)
+ return INFINITY;
- if( isnan(z->r) )
- return(z->r);
- if( isnan(z->i) )
- return(z->i);
+ if (isnan (z->r))
+ return z->r;
+ if (isnan (z->i))
+ return z->i;
- re = fabs( z->r );
- im = fabs( z->i );
+ re = fabs (z->r);
+ im = fabs (z->i);
- if( re == 0.0 )
- return( im );
- if( im == 0.0 )
- return( re );
+ if (re == 0.0)
+ return im;
+ if (im == 0.0)
+ return re;
/* Get the exponents of the numbers */
- x = frexp( re, &ex );
- y = frexp( im, &ey );
+ x = frexp (re, &ex);
+ y = frexp (im, &ey);
/* Check if one number is tiny compared to the other */
e = ex - ey;
- if( e > PREC )
- return( re );
- if( e < -PREC )
- return( im );
+ if (e > PREC)
+ return re;
+ if (e < -PREC)
+ return im;
/* Find approximate exponent e of the geometric mean. */
e = (ex + ey) >> 1;
/* Rescale so mean is about 1 */
- x = ldexp( re, -e );
- y = ldexp( im, -e );
+ x = ldexp (re, -e);
+ y = ldexp (im, -e);
/* Hypotenuse of the right triangle */
- b = sqrt( x * x + y * y );
+ b = sqrt (x * x + y * y);
/* Compute the exponent of the answer. */
- y = frexp( b, &ey );
+ y = frexp (b, &ey);
ey = e + ey;
/* Check it for overflow and underflow. */
- if( ey > MAXEXP )
+ if (ey > MAXEXP)
{
- math_set_error( "Cabs", MATH_ERROR_OVERFLOW );
- return( INFINITY );
+ math_set_error ("Cabs", MATH_ERROR_OVERFLOW);
+ return INFINITY;
}
- if( ey < MINEXP )
- return(0.0);
+ if (ey < MINEXP)
+ return 0.0;
/* Undo the scaling */
- b = ldexp( b, e );
- return( b );
+ b = ldexp (b, e);
+ return b;
}
/* Csqrt() - Complex square root
@@ -418,7 +418,7 @@
* SYNOPSIS:
* void Csqrt();
* Complex z, w;
- * Csqrt( &z, &w );
+ * Csqrt (&z, &w);
*
* DESCRIPTION:
* If z = x + iy, r = |z|, then
@@ -438,7 +438,7 @@
* DEC -10,+10 25000 3.2e-17 9.6e-18
* IEEE -10,+10 100000 3.2e-16 7.7e-17
* 2
- * Also tested by Csqrt( z ) = z, and tested by arguments
+ * Also tested by Csqrt (z) = z, and tested by arguments
* close to the real axis.
*/
static void
@@ -450,28 +450,28 @@
x = z->r;
y = z->i;
- if( y == 0.0 )
+ if (y == 0.0)
{
- if( x < 0.0 )
+ if (x < 0.0)
{
w->r = 0.0;
- w->i = sqrt(-x);
+ w->i = sqrt (-x);
return;
}
else
{
- w->r = sqrt(x);
+ w->r = sqrt (x);
w->i = 0.0;
return;
}
}
- if( x == 0.0 )
+ if (x == 0.0)
{
- r = fabs(y);
- r = sqrt(0.5*r);
- if( y > 0 )
+ r = fabs (y);
+ r = sqrt (0.5*r);
+ if (y > 0)
w->r = r;
else
w->r = -r;
@@ -479,26 +479,26 @@
return;
}
- /* Approximate sqrt(x^2+y^2) - x = y^2/2x - y^4/24x^3 + ... .
+ /* Approximate sqrt (x^2+y^2) - x = y^2/2x - y^4/24x^3 + ... .
* The relative error in the first term is approximately y^2/12x^2 .
*/
- if( (fabs(y) < 2.e-4 * fabs(x))
- && (x > 0) )
+ if ((fabs (y) < 2.e-4 * fabs (x))
+ && (x > 0))
{
t = 0.25*y*(y/x);
}
else
{
- r = Cabs(z);
+ r = Cabs (z);
t = 0.5*(r - x);
}
- r = sqrt(t);
+ r = sqrt (t);
q.i = r;
q.r = y/(2.0*r);
/* Heron iteration in complex arithmetic */
- Cdiv( &q, z, &s );
- Cadd( &q, &s, w );
+ Cdiv (&q, z, &s);
+ Cadd (&q, &s, w);
w->r *= 0.5;
w->i *= 0.5;
}
@@ -508,7 +508,7 @@
*
* SYNOPSIS:
* double phi, m, y, ellik();
- * y = ellik( phi, m );
+ * y = ellik (phi, m);
*
* DESCRIPTION:
* Approximates the integral
@@ -518,7 +518,7 @@
* | dt
* F(phi_\m) = | ------------------
* | 2
- * | | sqrt( 1 - m sin t )
+ * | | sqrt(1 - m sin t)
* -
* 0
* of amplitude phi and modulus m, using the arithmetic -
@@ -536,76 +536,76 @@
double a, b, c, e, temp, t, K;
int d, mod, sign, npio2;
- if( m == 0.0 )
- return( phi );
+ if (m == 0.0)
+ return phi;
a = 1.0 - m;
- if( a == 0.0 )
+ if (a == 0.0)
{
- if( fabs(phi) >= PIO2 )
+ if (fabs (phi) >= PIO2)
{
- math_set_error( "ellik", MATH_ERROR_SING );
- return( MAXNUM );
+ math_set_error ("ellik", MATH_ERROR_SING);
+ return MAXNUM;
}
- return( log( tan( (PIO2 + phi)/2.0 ) ) );
+ return log ( tan ((PIO2 + phi) / 2.0) );
}
- npio2 = floor( phi/PIO2 );
- if( npio2 & 1 )
+ npio2 = floor (phi/PIO2);
+ if (npio2 & 1)
npio2 += 1;
- if( npio2 )
+ if (npio2)
{
- K = ellpk( a );
+ K = ellpk (a);
phi = phi - npio2 * PIO2;
}
else
K = 0.0;
- if( phi < 0.0 )
+ if (phi < 0.0)
{
phi = -phi;
sign = -1;
}
else
sign = 0;
- b = sqrt(a);
- t = tan( phi );
- if( fabs(t) > 10.0 )
+ b = sqrt (a);
+ t = tan (phi);
+ if (fabs (t) > 10.0)
{
/* Transform the amplitude */
e = 1.0/(b*t);
/* ... but avoid multiple recursions. */
- if( fabs(e) < 10.0 )
+ if (fabs (e) < 10.0)
{
- e = atan(e);
- if( npio2 == 0 )
- K = ellpk( a );
- temp = K - ellik( e, m );
+ e = atan (e);
+ if (npio2 == 0)
+ K = ellpk (a);
+ temp = K - ellik (e, m);
goto done;
}
}
a = 1.0;
- c = sqrt(m);
+ c = sqrt (m);
d = 1;
mod = 0;
- while( fabs(c/a) > MACHEP )
+ while (fabs (c/a) > MACHEP)
{
temp = b/a;
- phi = phi + atan(t*temp) + mod * PI;
+ phi = phi + atan (t*temp) + mod * PI;
mod = (phi + PIO2)/PI;
- t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
- c = ( a - b )/2.0;
- temp = sqrt( a * b );
- a = ( a + b )/2.0;
+ t = t * (1.0 + temp)/(1.0 - temp * t * t);
+ c = (a - b)/2.0;
+ temp = sqrt (a * b);
+ a = (a + b)/2.0;
b = temp;
d += d;
}
- temp = (atan(t) + mod * PI)/(d * a);
+ temp = (atan (t) + mod * PI)/(d * a);
done:
- if( sign < 0 )
+ if (sign < 0)
temp = -temp;
temp += npio2 * K;
- return( temp );
+ return temp;
}
/* ellpj - Jacobian Elliptic Functions
@@ -613,7 +613,7 @@
* SYNOPSIS:
* double u, m, sn, cn, dn, phi;
* int ellpj();
- * ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
+ * ellpj (u, m, _&sn, _&cn, _&dn, _&phi);
*
* DESCRIPTION:
* Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
@@ -657,60 +657,60 @@
double a[9], c[9];
int i;
/* Check for special cases */
- if( m < 0.0 || m > 1.0 )
+ if (m < 0.0 || m > 1.0)
{
- math_set_error( "ellpj", MATH_ERROR_DOMAIN );
+ math_set_error ("ellpj", MATH_ERROR_DOMAIN);
*sn = 0.0;
*cn = 0.0;
*ph = 0.0;
*dn = 0.0;
- return(-1);
+ return -1;
}
- if( m < 1.0e-9 )
+ if (m < 1.0e-9)
{
- t = sin(u);
- b = cos(u);
+ t = sin (u);
+ b = cos (u);
ai = 0.25 * m * (u - t*b);
*sn = t - ai*b;
*cn = b + ai*t;
*ph = u - ai;
*dn = 1.0 - 0.5*m*t*t;
- return(0);
+ return 0;
}
- if( m >= 0.9999999999 )
+ if (m >= 0.9999999999)
{
ai = 0.25 * (1.0-m);
- b = cosh(u);
- t = tanh(u);
+ b = cosh (u);
+ t = tanh (u);
phi = 1.0/b;
- twon = b * sinh(u);
+ twon = b * sinh (u);
*sn = t + ai * (twon - u)/(b*b);
- *ph = 2.0*atan(exp(u)) - PIO2 + ai*(twon - u)/b;
+ *ph = 2.0*atan (exp (u)) - PIO2 + ai*(twon - u)/b;
ai *= t * phi;
*cn = phi - ai * (twon - u);
*dn = phi + ai * (twon + u);
- return(0);
+ return 0;
}
/* A. G. M. scale */
a[0] = 1.0;
- b = sqrt(1.0 - m);
- c[0] = sqrt(m);
+ b = sqrt (1.0 - m);
+ c[0] = sqrt (m);
twon = 1.0;
i = 0;
- while( fabs(c[i]/a[i]) > MACHEP )
+ while (fabs (c[i]/a[i]) > MACHEP)
{
- if( i > 7 )
+ if (i > 7)
{
- math_set_error( "ellpj", MATH_ERROR_OVERFLOW );
+ math_set_error ("ellpj", MATH_ERROR_OVERFLOW);
goto done;
}
ai = a[i];
++i;
- c[i] = ( ai - b )/2.0;
- t = sqrt( ai * b );
- a[i] = ( ai + b )/2.0;
+ c[i] = (ai - b)/2.0;
+ t = sqrt (ai * b);
+ a[i] = (ai + b)/2.0;
b = t;
twon *= 2.0;
}
@@ -720,25 +720,25 @@
phi = twon * a[i] * u;
do
{
- t = c[i] * sin(phi) / a[i];
+ t = c[i] * sin (phi) / a[i];
b = phi;
- phi = (asin(t) + phi)/2.0;
+ phi = (asin (t) + phi)/2.0;
}
- while( --i );
+ while (--i);
- *sn = sin(phi);
- t = cos(phi);
+ *sn = sin (phi);
+ t = cos (phi);
*cn = t;
- *dn = t/cos(phi-b);
+ *dn = t/cos (phi-b);
*ph = phi;
- return(0);
+ return 0;
}
/* ellpk - Complete elliptic integral of the first kind
*
* SYNOPSIS:
* double m1, y, ellpk();
- * y = ellpk( m1 );
+ * y = ellpk (m1);
*
* DESCRIPTION:
* Approximates the integral
@@ -748,7 +748,7 @@
* | dt
* K(m) = | ------------------
* | 2
- * | | sqrt( 1 - m sin t )
+ * | | sqrt(1 - m sin t)
* -
* 0
* where m = 1 - m1, using the approximation
@@ -799,26 +799,26 @@
};
static const double C1_ellpk = 1.3862943611198906188E0; /* log(4) */
- if( (x < 0.0) || (x > 1.0) )
+ if ((x < 0.0) || (x > 1.0))
{
- math_set_error( "ellpk", MATH_ERROR_DOMAIN );
- return( 0.0 );
+ math_set_error ("ellpk", MATH_ERROR_DOMAIN);
+ return 0.0;
}
- if( x > MACHEP )
+ if (x > MACHEP)
{
- return( polevl(x,P_ellpk,10) - log(x) * polevl(x,Q_ellpk,10) );
+ return polevl (x,P_ellpk,10) - log (x) * polevl (x,Q_ellpk,10);
}
else
{
- if( x == 0.0 )
+ if (x == 0.0)
{
- math_set_error( "ellpk", MATH_ERROR_SING );
- return( MAXNUM );
+ math_set_error ("ellpk", MATH_ERROR_SING);
+ return MAXNUM;
}
else
{
- return( C1_ellpk - 0.5 * log(x) );
+ return C1_ellpk - 0.5 * log (x);
}
}
}
@@ -829,7 +829,7 @@
* SYNOPSIS:
* int N;
* double x, y, coef[N+1], polevl[];
- * y = polevl( x, coef, N );
+ * y = polevl(x, coef, N);
*
* DESCRIPTION:
* Evaluates polynomial of degree N:
@@ -860,9 +860,9 @@
do
ans = ans * x + *p++;
- while( --i );
+ while (--i);
- return( ans );
+ return ans;
}
#if 1
@@ -923,7 +923,7 @@
else
zgain = an / (pn * ds->gain_scale);
VERBOSE ("# constant mygain factor %23.13E\n", zgain); // BSE info
- VERBOSE ("# z plane Denominator Numerator\n" ); // BSE info
+ VERBOSE ("# 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, aa[j], pp[j] * zgain); // BSE info
@@ -951,20 +951,20 @@
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 );
+ phi = exp (0.5 * ifr->passband_ripple_db / DECIBELL_FACTOR);
- if( (ifr->order & 1) == 0 )
+ if ((ifr->order & 1) == 0)
ds->gain_scale = phi;
else
ds->gain_scale = 1.0;
}
else
{ /* elliptic */
- ds->ripple_epsilon = exp( ifr->passband_ripple_db / DECIBELL_FACTOR );
+ 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->order & 1) == 0)
+ ds->gain_scale = sqrt (ds->ripple_epsilon);
+ ds->ripple_epsilon = sqrt (ds->ripple_epsilon - 1.0);
}
}
@@ -978,7 +978,7 @@
if (passband_edge1 >= ds->nyquist_frequency)
return "passband_edge1 too high";
- if( (ifr->type & 1) == 0 )
+ if ((ifr->type & 1) == 0)
{
if (passband_edge0 <= 0.0)
return "passband_edge too small";
@@ -1008,7 +1008,7 @@
}
/* Frequency correspondence for bilinear transformation
*
- * Wanalog = tan( 2 pi Fdigital T / 2 )
+ * Wanalog = tan(2 pi Fdigital T / 2)
*
* where T = 1/ifr->sampling_frequency
*/
@@ -1016,16 +1016,16 @@
double sang;
double cang = cos (ang);
ds->tan_angle_frequency = sin (ang) / cang; /* Wanalog */
- if( ifr->kind != 3 )
+ if (ifr->kind != 3)
{
ds->wc = ds->tan_angle_frequency;
- /*printf( "cos( 1/2 (Whigh-Wlow) T ) = %.5e, wc = %.5e\n", cang, ds->wc );*/
+ /*printf("cos(1/2 (Whigh-Wlow) T) = %.5e, wc = %.5e\n", cang, ds->wc);*/
}
- if( ifr->kind == 3 )
+ if (ifr->kind == 3)
{ /* elliptic */
- double tmp_cgam = cos( (high_edge + passband_edge0) * PI / ifr->sampling_frequency ) / cang;
+ 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;
@@ -1034,20 +1034,20 @@
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 );
+ 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 Kk1 = ellpk (m1p);
+ double Kpk1 = ellpk (m1);
+ double q = exp (-PI * Kpk1 / (ifr->order * Kk1));
fixme2local_1 = jacobi_theta_by_nome (q);
- if( ifr->type >= 3 )
+ if (ifr->type >= 3)
ds->wr = fixme2local_1;
else
ds->wr = 1.0 / fixme2local_1;
- if( ifr->type & 1 )
+ if (ifr->type & 1)
{
- ds->stopband_edge = atan( ds->tan_angle_frequency * ds->wr ) * ifr->sampling_frequency / PI;
+ ds->stopband_edge = atan (ds->tan_angle_frequency * ds->wr) * ifr->sampling_frequency / PI;
}
else
{
@@ -1055,11 +1055,11 @@
fixme2local_a = ds->tan_angle_frequency * ds->wr;
fixme2local_a *= fixme2local_a;
double b = fixme2local_a * (1.0 - ds->cgam * ds->cgam) + fixme2local_a * fixme2local_a;
- b = (ds->cgam + sqrt(b))/(1.0 + fixme2local_a);
- ds->stopband_edge = (PI / 2.0 - asin(b)) * ifr->sampling_frequency / (2.0 * PI);
+ b = (ds->cgam + sqrt (b))/(1.0 + fixme2local_a);
+ ds->stopband_edge = (PI / 2.0 - asin (b)) * ifr->sampling_frequency / (2.0 * PI);
}
}
- switch( ifr->type )
+ switch (ifr->type)
{
case 1:
if (ds->stopband_edge <= passband_edge1)
@@ -1081,10 +1081,10 @@
break;
}
ang = ds->stopband_edge * PI / ifr->sampling_frequency;
- cang = cos(ang);
- sang = sin(ang);
+ cang = cos (ang);
+ sang = sin (ang);
- if( ifr->type & 1 )
+ if (ifr->type & 1)
{
ds->wr = sang/(cang * ds->tan_angle_frequency);
}
@@ -1096,26 +1096,26 @@
ds->wr = (ds->cgam - cang)/(sang * ds->tan_angle_frequency);
}
- if( ifr->type >= 3 )
+ if (ifr->type >= 3)
ds->wr = 1.0 / ds->wr;
- if( ds->wr < 0.0 )
+ if (ds->wr < 0.0)
ds->wr = -ds->wr;
y[0] = 1.0;
y[1] = ds->wr;
/* ds->chebyshev_band_cbp = ds->wr; */
- if( ifr->type >= 3 )
+ if (ifr->type >= 3)
y[1] = 1.0 / y[1];
- if( ifr->type & 1 )
+ if (ifr->type & 1)
{
int i;
for (i = 1; i <= 2; i++)
{
- aa[i] = atan( ds->tan_angle_frequency * y[i - 1] ) * ifr->sampling_frequency / PI ;
+ aa[i] = atan (ds->tan_angle_frequency * y[i - 1]) * ifr->sampling_frequency / PI ;
}
- printf( "pass band %.9E\n", aa[1] );
- printf( "stop band %.9E\n", aa[2] );
+ printf ("pass band %.9E\n", aa[1]);
+ printf ("stop band %.9E\n", aa[2]);
}
else
{
@@ -1123,46 +1123,46 @@
for (i = 1; i <= 2; i++)
{
double a = ds->tan_angle_frequency * y[i - 1];
- double b = atan(a);
- double q = sqrt( 1.0 + a * a - ds->cgam * ds->cgam );
- q = atan2( q, ds->cgam );
+ double b = atan (a);
+ double q = sqrt (1.0 + a * a - ds->cgam * ds->cgam);
+ q = atan2 (q, ds->cgam);
aa[i] = (q + b) * ds->nyquist_frequency / PI;
pp[i] = (q - b) * ds->nyquist_frequency / PI;
}
- printf( "pass band %.9E %.9E\n", pp[1], aa[1] );
- printf( "stop band %.9E %.9E\n", pp[2], aa[2] );
+ printf ("pass band %.9E %.9E\n", pp[1], aa[1]);
+ printf ("stop band %.9E %.9E\n", pp[2], aa[2]);
}
ds->wc = 1.0;
find_elliptic_locations_in_lambda_plane (ifr, ds); /* find locations in lambda plane */
- if( (2 * ifr->order + 2) > ARRSIZ )
+ if ((2 * ifr->order + 2) > ARRSIZ)
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 )
+ * cos(1/2 (Whigh+Wlow) T)
+ * cos(Wcenter T) = ----------------------
+ * cos(1/2 (Whigh-Wlow) T)
*
*
* Band edges
- * cos( Wcenter T) - cos( Wdigital T )
+ * cos(Wcenter T) - cos(Wdigital T)
* Wanalog = -----------------------------------
- * sin( Wdigital T )
+ * sin(Wdigital T)
*/
- if( ifr->kind == 2 )
+ if (ifr->kind == 2)
{ /* Chebyshev */
double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
- ds->cgam = cos(a) / cang;
+ 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 )
+ if (ifr->kind == 1)
{ /* Butterworth */
double a = PI * (high_edge + passband_edge0) / ifr->sampling_frequency ;
- ds->cgam = cos(a) / cang;
+ 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;
@@ -1175,7 +1175,7 @@
find_s_plane_poles_and_zeros (ifr, ds); /* find s plane poles and zeros */
- if( ((ifr->type & 1) == 0) && ((4 * ifr->order + 2) > ARRSIZ) )
+ if (((ifr->type & 1) == 0) && ((4 * ifr->order + 2) > ARRSIZ))
goto toosml;
convert_s_plane_to_z_plane (ifr, ds); /* convert s plane to z plane */
@@ -1201,29 +1201,29 @@
Kk = ellpk (1.0 - fixme2local_m);
Kpk = ellpk (fixme2local_m);
EVERBOSE ("check: k=%.20g m=%.20g Kk=%.20g Kpk=%.20g\n", fixme2local_k, fixme2local_m, Kk, Kpk); // BSE info
- double q = exp( -PI * ifr->order * Kpk / Kk ); /* the nome of k1 */
+ double q = exp (-PI * ifr->order * Kpk / Kk); /* the nome of k1 */
double m1 = jacobi_theta_by_nome (q); /* see below */
- /* Note m1 = ds->ripple_epsilon / sqrt( A*A - 1.0 ) */
+ /* 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 );
+ a = 10.0 * log (a) / log (10.0);
+ printf ("dbdown %.9E\n", a);
a = 180.0 * asin (fixme2local_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 );
+ b = sqrt (1.0 - 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 Kk1 = ellpk (m1p);
+ double Kpk1 = ellpk (m1);
double r = Kpk1 * Kk / (Kk1 * Kpk);
- printf( "consistency check: n= %.14E\n", r );
+ 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 )
+ * sn j/ds->ripple_epsilon\m = j ellik(atan(1/ds->ripple_epsilon), m)
*/
b = 1.0 / ds->ripple_epsilon;
- phi = atan( b );
+ phi = atan (b);
fixme2local_2 = ellik (phi, m1p);
EVERBOSE ("phi=%.20g m=%.20g u=%.20g\n", phi, m1p, fixme2local_2);
/* consistency check on inverse sn */
@@ -1257,8 +1257,8 @@
for (i = 0; i < ds->n_poles; i++)
{ /* poles */
int lr = i + i;
- zs[lr] = -cos(m);
- zs[lr + 1] = sin(m);
+ zs[lr] = -cos (m);
+ zs[lr + 1] = sin (m);
m += PI / ifr->order;
}
/* high pass or band reject
@@ -1277,7 +1277,7 @@
/* The zeros at infinity map to the origin.
*/
ds->n_zeros = ds->n_poles;
- if( ifr->type == 4 )
+ if (ifr->type == 4)
{
ds->n_zeros += ifr->order / 2;
}
@@ -1290,14 +1290,14 @@
}
}
}
- if( ifr->kind == 2 )
+ if (ifr->kind == 2)
{
/* 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 */
- 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
+ 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
@@ -1313,13 +1313,13 @@
for (i = 0; i < ds->n_poles; i++)
{ /* poles */
int lr = i + i;
- zs[lr] = -a * cos(m);
- zs[lr + 1] = b * sin(m);
+ zs[lr] = -a * cos (m);
+ zs[lr + 1] = b * sin (m);
m += PI / ifr->order;
}
/* high pass or band reject
*/
- if( ifr->type >= 3 )
+ if (ifr->type >= 3)
{
/* map s => 1/s */
for (j = 0; j < ds->n_poles; j++)
@@ -1333,7 +1333,7 @@
/* The zeros at infinity map to the origin.
*/
ds->n_zeros = ds->n_poles;
- if( ifr->type == 4 )
+ if (ifr->type == 4)
{
ds->n_zeros += ifr->order / 2;
}
@@ -1346,7 +1346,7 @@
}
}
}
- if( ifr->kind == 3 ) /* elliptic filter -- stw */
+ if (ifr->kind == 3) /* elliptic filter -- stw */
{
double m = ds->elliptic_m;
ds->n_zeros = ifr->order / 2;
@@ -1357,7 +1357,7 @@
{ /* zeros */
double a = ifr->order - 1 - i - i;
double b = (Kk * a) / ifr->order;
- ellpj( b, m, &sn, &cn, &dn, &phi );
+ ellpj (b, m, &sn, &cn, &dn, &phi);
int lr = 2 * ds->n_poles + 2 * i;
zs[ lr ] = 0.0;
a = ds->wc / (ds->elliptic_k * sn); /* elliptic_k = sqrt(m) */
@@ -1367,7 +1367,7 @@
{ /* poles */
double a = ifr->order - 1 - i - i;
double b = a * Kk / ifr->order;
- ellpj( b, m, &sn, &cn, &dn, &phi );
+ ellpj (b, m, &sn, &cn, &dn, &phi);
double r = ds->elliptic_k * sn * sn1;
b = cn1 * cn1 + r * r;
a = -ds->wc * cn * dn * sn1 * cn1 / b;
@@ -1376,7 +1376,7 @@
b = ds->wc * sn * dn1 / b;
zs[lr + 1] = b;
}
- if( ifr->type >= 3 )
+ if (ifr->type >= 3)
{
int nt = ds->n_poles + ds->n_zeros;
for (j = 0; j < nt; j++)
@@ -1397,7 +1397,7 @@
}
}
}
- printf( "s plane poles:\n" );
+ printf ("s plane poles:\n");
j = 0;
for (i = 0; i < ds->n_poles + ds->n_zeros; i++)
{
@@ -1405,9 +1405,9 @@
++j;
double b = zs[j];
++j;
- printf( "%.9E %.9E\n", a, b );
+ printf ("%.9E %.9E\n", a, b);
if (i == ds->n_poles - 1)
- printf( "s plane zeros:\n" );
+ printf ("s plane zeros:\n");
}
return 0;
}
@@ -1418,17 +1418,17 @@
* AMS55 #16.38.5, 16.38.7
*
* 1/2
- * ( 2K ) 4 9
- * ( -- ) = 1 + 2q + 2q + 2q + ... = Theta (0,q)
- * ( pi ) 3
+ * (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
+ * (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) ).
+ * The nome q(m) = exp(- pi K(1-m)/K(m)).
*
* 1/2
* Given q, this program returns m .
@@ -1490,15 +1490,15 @@
r.r = zs[ir];
r.i = zs[ii];
- switch( ifr->type )
+ switch (ifr->type)
{
case 1:
case 3:
/* Substitute s - r = s/wc - r = (1/wc)(z-1)/(z+1) - r
*
- * 1 1 - r wc ( 1 + r wc )
- * = --- -------- ( z - -------- )
- * z+1 wc ( 1 - r wc )
+ * 1 1 - r wc ( 1 + r wc)
+ * = --- -------- (z - --------)
+ * z+1 wc ( 1 - r wc)
*
* giving the root in the z plane.
*/
@@ -1507,8 +1507,8 @@
cden.r = 1 - C * r.r;
cden.i = -C * r.i;
jt += 1;
- Cdiv( &cden, &cnum, &z[jt] );
- if( r.i != 0.0 )
+ Cdiv (&cden, &cnum, &z[jt]);
+ if (r.i != 0.0)
{
/* fill in complex conjugate root */
jt += 1;
@@ -1531,43 +1531,43 @@
*
* and solve for the roots in the z plane.
*/
- if( ifr->kind == 2 )
+ if (ifr->kind == 2)
cwc.r = ds->chebyshev_band_cbp;
else
cwc.r = ds->tan_angle_frequency;
cwc.i = 0.0;
- Cmul( &r, &cwc, &cnum ); /* r wc */
- Csub( &cnum, &COMPLEX_ONE, &ca ); /* a = 1 - r wc */
- Cmul( &cnum, &cnum, &b4ac ); /* 1 - (r wc)^2 */
- Csub( &b4ac, &COMPLEX_ONE, &b4ac );
+ Cmul (&r, &cwc, &cnum); /* r wc */
+ Csub (&cnum, &COMPLEX_ONE, &ca); /* a = 1 - r wc */
+ Cmul (&cnum, &cnum, &b4ac); /* 1 - (r wc)^2 */
+ Csub (&b4ac, &COMPLEX_ONE, &b4ac);
b4ac.r *= 4.0; /* 4ac */
b4ac.i *= 4.0;
cb.r = -2.0 * ds->cgam; /* b */
cb.i = 0.0;
- Cmul( &cb, &cb, &cnum ); /* b^2 */
- Csub( &b4ac, &cnum, &b4ac ); /* b^2 - 4 ac */
- Csqrt( &b4ac, &b4ac );
+ Cmul (&cb, &cb, &cnum); /* b^2 */
+ Csub (&b4ac, &cnum, &b4ac); /* b^2 - 4 ac */
+ Csqrt (&b4ac, &b4ac);
cb.r = -cb.r; /* -b */
cb.i = -cb.i;
ca.r *= 2.0; /* 2a */
ca.i *= 2.0;
- Cadd( &b4ac, &cb, &cnum ); /* -b + sqrt( b^2 - 4ac) */
- Cdiv( &ca, &cnum, &cnum ); /* ... /2a */
+ Cadd (&b4ac, &cb, &cnum); /* -b + sqrt(b^2 - 4ac) */
+ Cdiv (&ca, &cnum, &cnum); /* ... /2a */
jt += 1;
- Cmov( &cnum, &z[jt] );
- if( cnum.i != 0.0 )
+ Cmov (&cnum, &z[jt]);
+ if (cnum.i != 0.0)
{
jt += 1;
z[jt].r = cnum.r;
z[jt].i = -cnum.i;
}
- if( (r.i != 0.0) || (cnum.i == 0) )
+ if ((r.i != 0.0) || (cnum.i == 0))
{
- Csub( &b4ac, &cb, &cnum ); /* -b - sqrt( b^2 - 4ac) */
- Cdiv( &ca, &cnum, &cnum ); /* ... /2a */
+ Csub (&b4ac, &cb, &cnum); /* -b - sqrt(b^2 - 4ac) */
+ Cdiv (&ca, &cnum, &cnum); /* ... /2a */
jt += 1;
- Cmov( &cnum, &z[jt] );
- if( cnum.i != 0.0 )
+ Cmov (&cnum, &z[jt]);
+ if (cnum.i != 0.0)
{
jt += 1;
z[jt].r = cnum.r;
@@ -1576,15 +1576,15 @@
}
} /* end switch */
}
- while( --nc > 0 );
+ while (--nc > 0);
- if( icnt == 0 )
+ if (icnt == 0)
{
ds->n_solved_poles = jt + 1;
- if( ds->n_zeros <= 0 )
+ if (ds->n_zeros <= 0)
{
- if( ifr->kind != 3 )
- return(0);
+ if (ifr->kind != 3)
+ return 0;
else
break;
}
@@ -1603,21 +1603,21 @@
lin[1].r = 1.0;
lin[1].i = 0.0;
- if( ifr->kind != 3 )
+ if (ifr->kind != 3)
{ /* Butterworth or Chebyshev */
/* generate the remaining zeros */
- while( 2 * ds->n_solved_poles - 1 > jt )
+ while (2 * ds->n_solved_poles - 1 > jt)
{
- if( ifr->type != 3 )
+ if (ifr->type != 3)
{
- printf( "adding zero at Nyquist frequency\n" );
+ printf ("adding zero at Nyquist frequency\n");
jt += 1;
z[jt].r = -1.0; /* zero at Nyquist frequency */
z[jt].i = 0.0;
}
- if( (ifr->type == 2) || (ifr->type == 3) )
+ if ((ifr->type == 2) || (ifr->type == 3))
{
- printf( "adding zero at 0 Hz\n" );
+ printf ("adding zero at 0 Hz\n");
jt += 1;
z[jt].r = 1.0; /* zero at 0 Hz */
z[jt].i = 0.0;
@@ -1626,12 +1626,12 @@
}
else
{ /* elliptic */
- while( 2 * ds->n_solved_poles - 1 > jt )
+ while (2 * ds->n_solved_poles - 1 > jt)
{
jt += 1;
z[jt].r = -1.0; /* zero at Nyquist frequency */
z[jt].i = 0.0;
- if( (ifr->type == 2) || (ifr->type == 4) )
+ if ((ifr->type == 2) || (ifr->type == 4))
{
jt += 1;
z[jt].r = 1.0; /* zero at 0 Hz */
@@ -1639,7 +1639,7 @@
}
}
}
- printf( "order = %d\n", ds->n_solved_poles );
+ printf ("order = %d\n", ds->n_solved_poles);
int j;
@@ -1657,7 +1657,7 @@
for (j=0; j < ds->n_solved_poles; j++)
{
int jj = j;
- if( icnt )
+ if (icnt)
jj += ds->n_solved_poles;
double a = z[jj].r;
double b = z[jj].i;
@@ -1669,7 +1669,7 @@
y[jh + 1] = y[jh + 1] - b * pp[jh] - a * y[jh];
}
}
- if( icnt == 0 )
+ if (icnt == 0)
{
for (j=0; j <= ds->n_solved_poles; j++)
aa[j] = pp[j];
@@ -1677,7 +1677,7 @@
}
/* Scale factors of the pole and zero polynomials */
double a = 1.0;
- switch( ifr->type )
+ switch (ifr->type)
{
case 3:
a = -1.0;
@@ -1695,12 +1695,12 @@
break;
case 2:
- gam = PI / 2.0 - asin( ds->cgam ); /* = acos( cgam ) */
+ gam = PI / 2.0 - asin (ds->cgam); /* = acos(cgam) */
int mh = ds->n_solved_poles / 2;
pn = pp[mh];
an = aa[mh];
double ai = 0.0;
- if( mh > ((ds->n_solved_poles / 4)*2) )
+ if (mh > ((ds->n_solved_poles / 4)*2))
{
ai = 1.0;
pn = 0.0;
@@ -1725,37 +1725,37 @@
{
int j;
gain = an/(pn * ds->gain_scale);
- if( (ifr->kind != 3) && (pn == 0) )
+ if ((ifr->kind != 3) && (pn == 0))
gain = 1.0;
- printf( "constant gain factor %23.13E\n", gain );
+ printf ("constant gain factor %23.13E\n", gain);
for (j = 0; j <= ds->n_solved_poles; j++)
pp[j] = gain * pp[j];
- printf( "z plane Denominator Numerator\n" );
+ printf ("z plane Denominator Numerator\n");
for (j=0; j <= ds->n_solved_poles; j++)
{
- printf( "%2d %17.9E %17.9E\n", j, aa[j], pp[j] );
+ printf ("%2d %17.9E %17.9E\n", j, aa[j], pp[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" );
+ printf ("poles and zeros with corresponding quadratic factors\n");
for (j=0; j < ds->n_solved_poles; j++)
{
double a = z[j].r;
double b = z[j].i;
- if( b >= 0.0 )
+ if (b >= 0.0)
{
- printf( "pole %23.13E %23.13E\n", a, b );
+ 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 = z[jj].r;
b = z[jj].i;
- if( b >= 0.0 )
+ if (b >= 0.0)
{
- printf( "zero %23.13E %23.13E\n", a, b );
+ printf ("zero %23.13E %23.13E\n", a, b);
print_quadratic_factors (ifr, ds, a, b, 0);
}
}
@@ -1771,7 +1771,7 @@
{
double a, b, r, f, g, g0;
- if( y > 1.0e-16 )
+ if (y > 1.0e-16)
{
a = -2.0 * x;
b = x * x + y * y;
@@ -1781,17 +1781,17 @@
a = -x;
b = 0.0;
}
- printf( "q. f.\nz**2 %23.13E\nz**1 %23.13E\n", b, a );
- if( b != 0.0 )
+ printf ("q. f.\nz**2 %23.13E\nz**1 %23.13E\n", b, a);
+ if (b != 0.0)
{
/* resonant frequency */
- r = sqrt(b);
- f = PI / 2.0 - asin( -a/(2.0 * r) );
- f = f * ifr->sampling_frequency / (2.0 * PI );
+ r = sqrt (b);
+ f = PI / 2.0 - asin (-a/(2.0 * r));
+ f = f * ifr->sampling_frequency / (2.0 * PI);
/* gain at resonance */
g = 1.0 + r;
g = g * g - (a * a / r);
- g = (1.0 - r) * sqrt(g);
+ g = (1.0 - r) * sqrt (g);
g0 = 1.0 + a + b; /* gain at d.c. */
}
else
@@ -1804,18 +1804,18 @@
g0 = 1.0 + a;
}
- if( pzflg )
+ if (pzflg)
{
- if( g != 0.0 )
+ if (g != 0.0)
g = 1.0 / g;
else
g = MAXNUM;
- if( g0 != 0.0 )
+ if (g0 != 0.0)
g0 = 1.0 / g0;
else
g = MAXNUM;
}
- printf( "f0 %16.8E gain %12.4E DC gain %12.4E\n\n", f, g, g0 );
+ printf ("f0 %16.8E gain %12.4E DC gain %12.4E\n\n", f, g, g0);
return 0;
}
@@ -1829,11 +1829,11 @@
for (f=0; f < limit; f += limit / 21.)
{
double r = response (ifr, ds, f, gain);
- if( r <= 0.0 )
+ if (r <= 0.0)
r = -999.99;
else
- r = 2.0 * DECIBELL_FACTOR * log( r );
- printf( "%10.1f %10.2f\n", f, r );
+ r = 2.0 * DECIBELL_FACTOR * log (r);
+ printf ("%10.1f %10.2f\n", f, r);
// f = f + 0.05 * ds->nyquist_frequency;
}
}
@@ -1848,10 +1848,10 @@
double u;
int j;
- /* exp( j omega T ) */
+ /* exp(j omega T) */
u = 2.0 * PI * f /ifr->sampling_frequency;
- x.r = cos(u);
- x.i = sin(u);
+ x.r = cos (u);
+ x.i = sin (u);
num.r = 1.0;
num.i = 0.0;
@@ -1859,16 +1859,16 @@
den.i = 0.0;
for (j=0; j < ds->n_solved_poles; j++)
{
- Csub( &z[j], &x, &w );
- Cmul( &w, &den, &den );
- Csub( &z[j + ds->n_solved_poles], &x, &w );
- Cmul( &w, &num, &num );
+ Csub (&z[j], &x, &w);
+ Cmul (&w, &den, &den);
+ Csub (&z[j + ds->n_solved_poles], &x, &w);
+ Cmul (&w, &num, &num);
}
- Cdiv( &den, &num, &w );
+ Cdiv (&den, &num, &w);
w.r *= amp;
w.i *= amp;
- u = Cabs( &w );
- return(u);
+ u = Cabs (&w);
+ return u;
}
static double
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
