home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ellf / ellf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  14.6 KB  |  911 lines
  1. /* ellf.c
  2.  * 
  3.  * Read ellf.doc before attempting to compile this program.
  4.  */
  5.  
  6.  
  7. /* size of arrays: */
  8. #define ARRSIZ 50
  9.  
  10. double atan2();
  11.  
  12. /* 1 = yes, 0 = no: */
  13. #define DECPDP 0
  14. #define UNK 1
  15. #define CEPHES 0
  16.  
  17. #if DECPDP
  18. #define asin arcsin
  19. #define atan arctan
  20. #define artan2(a,b) (atan2(a,b))
  21. #endif
  22.  
  23.  
  24. #if CEPHES
  25. #define gtlin gets
  26. #define artan2(a,b) (atan2(a,b))
  27. #endif
  28.  
  29. #if UNK
  30. #define gtlin gets
  31. #define artan2(a,b) (atan2(b,a))
  32. double PI   = 3.14159265358979323846;
  33. double PIO2 = 1.57079632679489661923;
  34. /* MACHEP is the relative resolution of your computer's
  35.  * floating point arithmetic.  It will be 1.1E-16 for IEEE
  36.  * double precision, 1.39E-17 for DEC VAX or PDP-11. */
  37. double MACHEP = 1.1e-16;
  38. double MAXNUM = 1.5e38;
  39. #endif
  40.  
  41. extern double PI;
  42. extern double PIO2;
  43. extern double MACHEP;
  44.  
  45. static double aa[ARRSIZ] = {0.0};
  46. static double pp[ARRSIZ] = {0.0};
  47. static double y[ARRSIZ] = {0.0};
  48. static double zs[ARRSIZ] = {0.0};
  49. static double z[2][ARRSIZ] = {0.0};
  50. static double wr = 0.0;
  51. static double wc = 0.0;
  52. static double rn = 8.0;
  53. static double c = 0.0;
  54. static double cgam = 0.0;
  55. static double scale = 0.0;
  56. static double fs = 1.0e4;
  57. static double dbr = 0.5;
  58. static double dbd = -40.0;
  59. static double f1 = 1.5;
  60. static double f2 = 2.0e3;
  61. static double f3 = 2.4e3;
  62. static double dbfac = 0.0;
  63. static double a = 0.0;
  64. static double b = 0.0;
  65. static double p = 0.0;
  66. static double q = 0.0;
  67. static double r = 0.0;
  68. static double ri = 0.0;
  69.  
  70. static double t1 = 0.0;
  71. static double t2 = 0.0;
  72. static double u = 0.0;
  73. static double k = 0.0;
  74. static double kp = 0.0;
  75. static double m = 0.0;
  76. static double Kk = 0.0;
  77. static double Kk1 = 0.0;
  78. static double Kpk = 0.0;
  79. static double Kpk1 = 0.0;
  80. static double eps = 0.0;
  81. static double rho = 0.0;
  82. static double phi = 0.0;
  83. static double sn = 0.0;
  84. static double cn = 0.0;
  85. static double dn = 0.0;
  86. static double sn1 = 0.0;
  87. static double cn1 = 0.0;
  88. static double dn1 = 0.0;
  89. static double phi1 = 0.0;
  90. static double m1 = 0.0;
  91. static double m1p = 0.0;
  92. static double cang = 0.0;
  93. static double sang = 0.0;
  94. static double bw = 0.0;
  95. static double ang = 0.0;
  96. static double fnyq = 0.0;
  97. static double rr = 0.0;
  98. static double ar = 0.0;
  99. static double ai = 0.0;
  100. static double br = 0.0;
  101. static double bi = 0.0;
  102. static double qr = 0.0;
  103. static double qm = 0.0;
  104. static double qa = 0.0;
  105. static double sqr = 0.0;
  106. static double sqi = 0.0;
  107. static double qi = 0.0;
  108. static double pn = 0.0;
  109. static double an = 0.0;
  110. static double gam = 0.0;
  111. static double cng = 0.0;
  112. static int lr = 0;
  113. static int nt = 0;
  114. static int i = 0;
  115. static int j = 0;
  116. static int jt = 0;
  117. static int nc = 0;
  118. static int ii = 0;
  119. static int ir = 0;
  120. static int ni = 0;
  121. static int nr = 0;
  122.  
  123. static int zord = 0;
  124. static int icnt = 0;
  125. static int mh = 0;
  126. static int jj = 0;
  127. static int jh = 0;
  128. static int jl = 0;
  129. static int n = 8;
  130. static int np = 0;
  131. static int nz = 0;
  132. static int type = 1;
  133. static char salut[] =
  134.  
  135. {"Filter type:\n1 low pass\n2 band pass\n3 high pass\n4 band stop\n"};
  136.  
  137. double abs(), exp(), log(), cos(), sin(), sqrt();
  138. double ellpk(), ellik(), asin(), atan();
  139. double cay();
  140.  
  141.  
  142. main()
  143. {
  144. char str[80];
  145.  
  146. dbfac = 10.0/log(10.0);
  147.  
  148. top:
  149.  
  150. printf( "%s ? ", salut );    /* ask for filter type */
  151. gtlin( str );            /* gtlin() is like gets() */
  152. sscanf( str, "%d", &type );
  153. printf( "%d\n", type );
  154. if( (type <= 0) || (type > 4) )
  155.     exit(0);
  156.  
  157. getnum( "Order of filter", &rn ); /* see below for getnum() */
  158. n = rn;
  159. if( n <= 0 )
  160.     {
  161. specerr:
  162.     printf( "? Specification error\n" );
  163.     goto top;
  164.     }
  165. rn = n;    /* ensure it is an integer */
  166. getnum( "Passband ripple, db", &dbr );
  167. if( dbr <= 0.0 )
  168.     goto specerr;
  169. eps = exp( dbr/dbfac );
  170. scale = 1.0;
  171. if( (n & 1) == 0 )
  172.     scale = sqrt( eps );
  173. eps = sqrt( eps - 1.0 );
  174.  
  175. getnum( "Sampling frequency", &fs );
  176. if( fs <= 0.0 )
  177.     goto specerr;
  178.  
  179. fnyq = fs/2.0;
  180.  
  181. getnum( "Passband edge", &f2 );
  182. if( (f2 <= 0.0) || (f2 >= fnyq) )
  183.     goto specerr;
  184.  
  185. if( (type & 1) == 0 )
  186.     {
  187.     getnum( "Other passband edge", &f1 );
  188.     if( (f1 <= 0.0) || (f1 >= fnyq) )
  189.         goto specerr;
  190.     }
  191. else
  192.     f1 = 0.0;
  193.  
  194. if( f2 < f1 )
  195.     {
  196.     a = f2;
  197.     f2 = f1;
  198.     f1 = a;
  199.     }
  200. if( type == 3 )    /* high pass */
  201.     {
  202.     bw = f2;
  203.     a = fnyq;
  204.     }
  205. else
  206.     {
  207.     bw = f2 - f1;
  208.     a = f2;
  209.     }
  210. ang = bw * PI / fs;
  211. cang = cos( ang );
  212. c = sin(ang) / cang;
  213. cgam = cos( (a+f1) * PI / fs ) / cang;
  214.  
  215. getnum( "Stop band edge or -(db down)", &dbd );
  216. if( dbd > 0.0 )
  217.     f3 = dbd;
  218. else
  219.     { /* calculate band edge from db down */
  220.     a = exp( -dbd/dbfac );
  221.     m1 = eps/sqrt( a - 1.0 );
  222.     m1 *= m1;
  223.     m1p = 1.0 - m1;
  224.     Kk1 = ellpk( m1p );
  225.     Kpk1 = ellpk( m1 );
  226.     q = exp( -PI * Kpk1 / (rn * Kk1) );
  227.     k = cay(q);
  228.     if( type >= 3 )
  229.         wr = k;
  230.     else
  231.         wr = 1.0/k;
  232.     if( type & 1 )
  233.         {
  234.         f3 = atan( c * wr ) * fs / PI;
  235.         }
  236.     else
  237.         {
  238.         a = c * wr;
  239.         a *= a;
  240.         b = a * (1.0 - cgam * cgam) + a * a;
  241.         b = (cgam + sqrt(b))/(1.0 + a);
  242.         f3 = (PI/2.0 - asin(b)) * fs / (2.0*PI);
  243.         }
  244.     }
  245.  
  246. switch( type )
  247.     {
  248.     case 1:
  249.     if( f3 <= f2 )
  250.         goto specerr;
  251.     break;
  252.     case 2:
  253.     if( (f3 > f2) || (f3 < f1) )
  254.         break;
  255.     goto specerr;
  256.     case 3:
  257.     if( f3 >= f2 )
  258.         goto specerr;
  259.     break;
  260.     case 4:
  261.     if( (f3 <= f1) || (f3 >= f2) )
  262.         goto specerr;
  263.     break;
  264.     }
  265.  
  266. ang = f3 * PI / fs;
  267. cang = cos(ang);
  268. sang = sin(ang);
  269.  
  270. if( type & 1 )
  271.     {
  272.     wr = sang/(cang*c);
  273.     }
  274. else
  275.     {
  276.     q = cang * cang  -  sang * sang;
  277.     sang = 2.0 * cang * sang;
  278.     cang = q;
  279.     wr = (cgam - cang)/(sang * c);
  280.     }
  281.  
  282. if( type >= 3 )
  283.     wr = 1.0/wr;
  284. if( wr < 0.0 )
  285.     wr = -wr;
  286. y[0] = 1.0;
  287. y[1] = wr;
  288. if( type >= 3 )
  289.     y[1] = 1.0/y[1];
  290.  
  291. if( type & 1 )
  292.     {
  293.     for( i=1; i<=2; i++ )
  294.         {
  295.         aa[i] = atan( c * y[i-1] ) * fs / PI ;
  296.         }
  297.     printf( "pass band %.9E\n", aa[1] );
  298.     printf( "stop band %.9E\n", aa[2] );
  299.     }
  300. else
  301.     {
  302.     for( i=1; i<=2; i++ )
  303.         {
  304.         a = c * y[i-1];
  305.         b = atan(a);
  306.         q = sqrt( 1.0 + a * a  -  cgam * cgam );
  307.         q = artan2( cgam, q );
  308.         aa[i] = (q + b) * fnyq / PI;
  309.         pp[i] = (q - b) * fnyq / PI;
  310.         }
  311.     printf( "pass band %.9E %.9E\n", pp[1], aa[1] );
  312.     printf( "stop band %.9E %.9E\n", pp[2], aa[2] );
  313.     }
  314.  
  315. lampln();    /* find locations in lambda plane */
  316. if( (2*n+2) > ARRSIZ )
  317.     goto toosml;
  318.  
  319. spln();        /* find s plane poles and zeros */
  320.  
  321. if( ((type & 1) == 0) && ((4*n+2) > ARRSIZ) )
  322.     goto toosml;
  323.  
  324. zplna();    /* convert s plane to z plane */
  325. zplnb();
  326. zplnc();
  327. goto top;
  328.  
  329. toosml:
  330. printf( "Cannot continue, storage arrays too small\n" );
  331. goto top;
  332. }
  333.  
  334.  
  335. lampln()
  336. {
  337.  
  338. wc = 1.0;
  339. k = wc/wr;
  340.  
  341. m = k * k;
  342. Kk = ellpk( 1.0 - m );
  343. Kpk = ellpk( m );
  344.  
  345. q = exp( -PI * rn * Kpk / Kk );    /* the nome of k1 */
  346. m1 = cay(q); /* see below */
  347.  
  348. /* Note m1 = eps / sqrt( A*A - 1.0 ) */
  349. a = eps/m1;
  350. a =  a * a + 1;
  351. a = 10.0 * log(a) / log(10.0);
  352. printf( "dbdown %.9E\n", a );
  353.  
  354. a = 180.0 * asin( k ) / PI;
  355. b = 1.0/(1.0 + eps*eps);
  356. b = sqrt( 1.0 - b );
  357. printf( "theta %.9E, rho %.9E\n", a, b );
  358. m1 *= m1;
  359. m1p = 1.0 - m1;
  360. Kk1 = ellpk( m1p );
  361. Kpk1 = ellpk( m1 );
  362.  
  363. r = Kpk1 * Kk / (Kk1 * Kpk);
  364. printf( "consistency check: n= %.14E\n", r );
  365.  
  366. /*   -1
  367.  * sn   j/eps\m  =  j ellik( atan(1/eps), m )
  368.  */
  369. b = 1.0/eps;
  370. phi = atan( b );
  371. u = ellik( phi, m1p );
  372. printf( "phi %.7e m %.7e u %.7e\n", phi, m1p, u );
  373.  
  374. /* consistency check on inverse sn */
  375.  
  376. ellpj( u, m1p, &sn, &cn, &dn, &phi );
  377. a = sn/cn;
  378. printf( "consistency check: sn/cn = %.9E = %.9E = 1/eps\n", a, b );
  379.  
  380.  
  381. u = u * Kk / (rn * Kk1);    /* or, u = u * Kpk / Kpk1 */
  382. }
  383.  
  384.  
  385. /* calculate s plane poles and zeros, normalized to wc = 1 */
  386. spln()
  387. {
  388. np = (n+1)/2;
  389. nz = n/2;
  390. ellpj( u, 1.0-m, &sn1, &cn1, &dn1, &phi1 );
  391. for( i=0; i<ARRSIZ; i++ )
  392.     zs[i] = 0.0;
  393. for( i=0; i<nz; i++ )
  394.     {    /* zeros */
  395.     a = n - 1 - i - i;
  396.     b = (Kk * a) / rn;
  397.     ellpj( b, m, &sn, &cn, &dn, &phi );
  398.     lr = 2*np + 2*i;
  399.     zs[ lr ] = 0.0;
  400.     a = wc/(k*sn);    /* k = sqrt(m) */
  401.     zs[ lr + 1 ] = a;
  402.     }
  403. for( i=0; i<np; i++ )
  404.     {    /* poles */
  405.     a = n - 1 - i - i;
  406.     b = a * Kk / rn;        
  407.     ellpj( b, m, &sn, &cn, &dn, &phi );
  408.     r = k * sn * sn1;
  409.     b = cn1*cn1 + r*r;
  410.     a = -wc*cn*dn*sn1*cn1/b;
  411.     lr = i + i;
  412.     zs[lr] = a;
  413.     b = wc*sn*dn1/b;
  414.     zs[lr+1] = b;
  415.     }    
  416. if( type >= 3 )
  417. {
  418. nt = np + nz;
  419. for( j=0; j<nt; j++ )
  420.     {
  421.     ir = j + j;
  422.     ii = ir + 1;
  423.     b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
  424.     zs[ir] = zs[ir] / b;
  425.     zs[ii] = zs[ii] / b;
  426.     }
  427. while( np > nz )
  428.     {
  429.     ir = ii + 1;
  430.     ii = ir + 1;
  431.     nz += 1;
  432.     zs[ir] = 0.0;
  433.     zs[ii] = 0.0;
  434.     }
  435. }
  436.  
  437. printf( "s plane poles:\n" );
  438. j = 0;
  439. for( i=0; i<np+nz; i++ )
  440.     {
  441.     a = zs[j];
  442.     ++j;
  443.     b = zs[j];
  444.     ++j;
  445.     printf( "%.9E %.9E\n", a, b );
  446.     if( i == np-1 )
  447.         printf( "s plane zeros:\n" );
  448.     }
  449. }
  450.  
  451. getnum( line, val )
  452. char *line;
  453. double *val;
  454. {
  455. char s[40];
  456.  
  457. printf( "%s = %.9E ? ", line, *val );
  458. gtlin( s );
  459. if( s[0] != '\0' )
  460.     {
  461.     sscanf( s, "%lf", val );
  462.     printf( "%.9E\n", *val );
  463.     }
  464. }
  465.  
  466. /*        cay()
  467.  *
  468.  * Find parameter corresponding to given nome by expansion
  469.  * in theta functions:
  470.  * AMS55 #16.38.5, 16.38.7
  471.  *
  472.  *       1/2
  473.  * ( 2K )                   4     9
  474.  * ( -- )     =  1 + 2q + 2q  + 2q  + ...  =  Theta (0,q)
  475.  * ( pi )                                          3
  476.  *
  477.  *
  478.  *       1/2
  479.  * ( 2K )     1/4       1/4        2    6    12    20
  480.  * ( -- )    m     =  2q    ( 1 + q  + q  + q   + q   + ...) = Theta (0,q)
  481.  * ( pi )                                                           2
  482.  *
  483.  * The nome q(m) = exp( - pi K(1-m)/K(m) ).
  484.  *
  485.  *                                1/2
  486.  * Given q, this program returns m   .
  487.  */
  488.  
  489. double cay(q)
  490. double q;
  491. {
  492. double a, b, p, r;
  493. double t1, t2;
  494. double abs(), sqrt();
  495.  
  496. a = 1.0;
  497. b = 1.0;
  498. r = 1.0;
  499. p = q;
  500.  
  501. do
  502. {
  503. r *= p;
  504. a += 2.0 * r;
  505. t1 = abs( r/a );
  506.  
  507. r *= p;
  508. b += r;
  509. p *= q;
  510. t2 = abs( r/b );
  511. if( t2 > t1 )
  512.     t1 = t2;
  513. }
  514. while( t1 > MACHEP );
  515.  
  516. a = b/a;
  517. a = 4.0 * sqrt(q) * a * a;    /* see above formulas, solved for m */
  518. return(a);
  519. }
  520.  
  521. /*        zpln.c
  522.  * Program to convert s plane poles and zeros to the z plane.
  523.  */
  524.  
  525.  
  526. zplna()
  527. {
  528.  
  529. for( i=0; i<ARRSIZ; i++ )
  530.     {
  531.     z[0][i] = 0.0;
  532.     z[1][i] = 0.0;
  533.     }
  534.  
  535.  
  536. nc = np;
  537. jt = -1;
  538. ii = -1;
  539.  
  540. for( icnt=0; icnt<2; icnt++ )
  541. {
  542.     /* The maps from s plane to z plane */
  543. do
  544. {
  545. ir = ii + 1;
  546. ii = ir + 1;
  547. rr = zs[ir];
  548. ri = zs[ii];
  549. b = 1.0 - rr * c;
  550. b *= b;
  551. b += ri * ri * c * c;
  552.  
  553. switch( type )
  554. {
  555. case 1:
  556. case 3:
  557.  
  558. jt += 1;
  559. z[0][jt] = (1.0 - (rr * rr + ri * ri)*c*c)/b;
  560. z[1][jt] = -2.0 * ri * c / b;
  561. if( ri == 0.0 )
  562.     break;
  563.  
  564. jt += 1;
  565. z[0][jt] = z[0][jt-1 ];
  566. z[1][jt] = -z[1][jt-1 ];
  567. break;
  568.  
  569. case 2:
  570. case 4:
  571.  
  572. rr = c * rr;
  573. ri = c * ri;
  574. b = 1.0 - rr;
  575. b *= b;
  576. b += ri * ri;
  577. ar = (1.0 - rr) * cgam / b;
  578. ai = cgam * ri / b;
  579. br = (1.0 - rr * rr - ri * ri) / b;
  580. bi = 2.0 * ri / b;
  581. qr = ar * ar - ai * ai - br;
  582. qi = 2.0 * ar * ai - bi;
  583. qm = sqrt( qr * qr + qi * qi );
  584. qm = sqrt( qm );
  585. qa = artan2( qr, qi ) / 2.0;
  586. sqr = qm * cos( qa );
  587. sqi = qm * sin( qa );
  588. jt += 1;
  589. z[0][jt] = ar + sqr;
  590. z[1][jt] = ai + sqi;
  591. jt += 1;
  592. z[0][jt] = ar - sqr;
  593. z[1][jt] = ai - sqi;
  594.  
  595. if( ri > 0.0 )
  596.     {
  597.     jt += 2;
  598.     for( i=0; i<=1; i++ )
  599.         {
  600.         jj = jt - 1 + i;
  601.         z[0][jj] = z[0][jj-2];
  602.         z[1][jj] = -z[1][jj-2];
  603.         }
  604.     }
  605. } /* end switch */
  606.  
  607. }
  608. while( --nc > 0 );
  609.  
  610. if( icnt == 0 )
  611.     {
  612.     zord = jt+1;
  613.     if( nz <= 0 )
  614.         break;
  615.     }
  616. nc = nz;
  617. }    /* end for() loop */
  618.  
  619. }
  620.  
  621. zplnb()
  622. {
  623.  
  624. while( 2*zord - 1 > jt )
  625.     {
  626.     jt += 1;
  627.     z[0][jt] = -1.0;
  628.     z[1][jt] = 0.0;
  629.     if( (type == 2) || (type == 4) )
  630.         {
  631.         jt += 1;
  632.         z[0][jt] = 1.0;
  633.         z[1][jt] = 0.0;
  634.         }
  635.     }
  636. printf( "order = %d\n", zord );
  637.  
  638. /* Expand the poles and zeros into numerator and
  639.  * denominator polynomials
  640.  */
  641. for( icnt=0; icnt<2; icnt++ )
  642. {
  643. for( j=0; j<ARRSIZ; j++ )
  644.     {
  645.     pp[j] = 0.0;
  646.     y[j] = 0.0;
  647.     }
  648.  
  649. pp[0] = 1.0;
  650.  
  651. for( j=0; j<zord; j++ )
  652.     {
  653.     jj = j;
  654.     if( icnt )
  655.         jj += zord;
  656.     a = z[0][jj];
  657.     b = z[1][jj];
  658.     for( i=0; i<=j; i++ )
  659.         {
  660.         jh = j - i;
  661.         pp[jh+1] = pp[jh+1] - a * pp[jh] + b * y[jh];
  662.         y[jh+1] =  y[jh+1]  - b * pp[jh] - a * y[jh];
  663.         }
  664.     }
  665.  
  666. if( icnt == 0 )
  667.     {
  668.     for( j=0; j<=zord; j++ )
  669.         aa[j] = pp[j];
  670.     }
  671. }
  672.  
  673. /* Scale factors of the pole and zero polynomials */
  674.  
  675. a = 1.0;
  676.  
  677. switch( type )
  678. {
  679.  
  680. case 3:
  681.  
  682. a = -1.0;
  683.  
  684. case 1:
  685. case 4:
  686.  
  687. pn = 1.0;
  688. an = 1.0;
  689. for( j=1; j<=zord; j++ )
  690.     {
  691.     pn = a * pn + pp[j];
  692.     an = a * an + aa[j];
  693.     }
  694. break;
  695.  
  696. case 2:
  697.  
  698. gam = PI/2.0 - asin( cgam );
  699. mh = zord/2;
  700. pn = pp[mh];
  701. an = aa[mh];
  702. ai = 0.0;
  703.  
  704. if( mh > ((zord/4)*2) )
  705.     {
  706.     ai = 1.0;
  707.     pn = 0.0;
  708.     an = 0.0;
  709.     }
  710.  
  711. for( j=1; j<=mh; j++ )
  712.     {
  713.     a = gam * j - ai * PI / 2.0;
  714.     cng = cos(a);
  715.     jh = mh + j;
  716.     jl = mh - j;
  717.     pn = pn + cng * (pp[jh] + (1.0 - 2.0 * ai) * pp[jl]);
  718.     an = an + cng * (aa[jh] + (1.0 - 2.0 * ai) * aa[jl]);
  719.     }
  720. }
  721. }
  722.  
  723. zplnc()
  724. {
  725.  
  726. q = an/(pn*scale);
  727. printf( "constant gain factor %23.13E\n", q );
  728. for( j=0; j<=zord; j++ )
  729.     pp[j] = q * pp[j];
  730.  
  731. printf( "z plane Denominator      Numerator\n" );
  732. for( j=0; j<=zord; j++ )
  733.     {
  734.     printf( "%2d %17.9E %17.9E\n", j, aa[j], pp[j] );
  735.     }
  736. printf( "poles and zeros with corresponding quadratic factors\n" );
  737. for( j=0; j<zord; j++ )
  738.     {
  739.     a = z[0][j];
  740.     b = z[1][j];
  741.     if( b >= 0.0 )
  742.         {
  743.         printf( "pole  %23.13E %23.13E\n", a, b );
  744.         quadf( a, b, 1 );
  745.         }
  746.     jj = j + zord;
  747.     a = z[0][jj];
  748.     b = z[1][jj];
  749.     if( b >= 0.0 )
  750.         {
  751.         printf( "zero  %23.13E %23.13E\n", a, b );
  752.         quadf( a, b, 0 );
  753.         }
  754.     }
  755. }
  756.  
  757. /* display quadratic factors */
  758.  
  759. quadf( x, y, pzflg )
  760. double x, y;
  761. int pzflg;    /* 1 if poles, 0 if zeros */
  762. {
  763. double a, b, r, f, g, g0;
  764. double asin(), sqrt();
  765.  
  766. if( y > 1.0e-16 )
  767.     {
  768.     a = -2.0 * x;
  769.     b = x*x + y*y;
  770.     }
  771. else
  772.     {
  773.     a = -x;
  774.     b = 0.0;
  775.     }
  776. printf( "q. f.\nz**2 %23.13E\nz**1 %23.13E\n", b, a );
  777. if( b != 0.0 )
  778.     {
  779.     /* resonant frequency */
  780.     r = sqrt(b);
  781.     f = PI/2.0 - asin( -a/(2.0*r) );
  782.     f = f * fs / (2.0 * PI );
  783.     /* gain at resonance */
  784.     g = 1.0 + r;
  785.     g = g*g - (a*a/r);
  786.     g = (1.0 - r) * sqrt(g);
  787.     g0 = 1.0 + a + b;    /* gain at d.c. */
  788.     }
  789. else
  790.     {
  791.     /* It is really a first-order network.
  792.      * Give the gain at fnyq and D.C.
  793.      */
  794.  
  795.     f = fnyq;
  796.     g = 1.0 - a;
  797.     g0 = 1.0 + a;
  798.     }
  799.  
  800. if( pzflg )
  801.     {
  802.     g = 1.0/g;
  803.     g0 = 1.0/g0;
  804.     }
  805. printf( "f0 %16.8E  gain %12.4E  DC gain %12.4E\n\n", f, g, g0 );
  806. }
  807.  
  808. /* double precision absolute value */
  809.  
  810. double abs(x)
  811. double x;
  812. {
  813. if( x < 0.0 )
  814.     return( -x );
  815. else
  816.     return( x );
  817. }
  818. #if !CEPHES
  819. /*                            polevl.c
  820.  *                            p1evl.c
  821.  *
  822.  *    Evaluate polynomial
  823.  *
  824.  *
  825.  *
  826.  * SYNOPSIS:
  827.  *
  828.  * int N;
  829.  * double x, y, coef[N+1], polevl[];
  830.  *
  831.  * y = polevl( C, x, N );
  832.  *
  833.  *
  834.  *
  835.  * DESCRIPTION:
  836.  *
  837.  * Evaluates polynomial of degree N:
  838.  *
  839.  *                     2          N
  840.  * y  =  C  + C x + C x  +...+ C x
  841.  *        0    1     2          N
  842.  *
  843.  * Coefficients are stored in reverse order:
  844.  *
  845.  * coef[0] = C  , ..., coef[N] = C  .
  846.  *            N                   0
  847.  *
  848.  *  The function p1evl() assumes that coef[N] = 1.0 and is
  849.  * omitted from the array.  Its calling arguments are
  850.  * otherwise the same as polevl().
  851.  *
  852.  *
  853.  * SPEED:
  854.  *
  855.  * In the interest of speed, there are no checks for out
  856.  * of bounds arithmetic.  This routine is used by most of
  857.  * the functions in the library.  Depending on available
  858.  * equipment features, the user may wish to rewrite the
  859.  * program in microcode or assembly language.
  860.  *
  861.  */
  862. /*                            polevl.c */
  863.  
  864.  
  865. double polevl( x, coef, N )
  866. double x;
  867. double coef[];
  868. int N;
  869. {
  870. double ans;
  871. short i;
  872. register double *p;
  873.  
  874. p = coef;
  875. ans = 0.0;
  876. i = N+1;
  877.  
  878. do
  879.     ans = ans * x  +  *p++;
  880. while( --i );
  881.  
  882. return( ans );
  883. }
  884.  
  885. /*                            p1evl()    */
  886. /*                                          N
  887.  * Evaluate polynomial when coefficient of x  is 1.0.
  888.  * Otherwise same as polevl.
  889.  */
  890.  
  891. double p1evl( x, coef, N )
  892. double x;
  893. double coef[];
  894. int N;
  895. {
  896. double ans;
  897. register double *p;
  898. short i;
  899.  
  900. p = coef;
  901. ans = x + *p++;
  902. i = N-1;
  903.  
  904. do
  905.     ans = ans * x  + *p++;
  906. while( --i );
  907.  
  908. return( ans );
  909. }
  910. #endif    /*!CEPHES*/
  911.