home *** CD-ROM | disk | FTP | other *** search
/ Peanuts NeXT Software Archives / Peanuts-2.iso / Developer / hardware / dsp / drbub / analog / ellip.c next >
Encoding:
C/C++ Source or Header  |  1991-09-07  |  5.6 KB  |  176 lines
  1. /********************************************************************
  2. *                       ELLIPTICAL FILTER DESIGN
  3. * Given inputs of lamda (frequency transistion range), epsilon ( a measuer
  4. * of the pass band ripple, and the order of the sysetm, this program
  5. * will generate the necessary poles and zeros.  It is based on a algorithmn
  6. * by Sindy Darlington.
  7. *                                       steve punte 10-20-87
  8. */
  9. #include        <stdio.h>
  10. #include        <math.h>
  11. #define         TRUE            1
  12. #define         FALSE           0
  13. #define         MX              50      /* Max number of poles or zeros */
  14. #define         PIE             3.141592653
  15. #define         scanf           my_scanf
  16.  
  17. struct    cmplx_num{
  18.             double    real;
  19.             double    imag;
  20.          };
  21.  
  22. #define         sqr(X)          ( X )*( X )
  23. #define         cset(X, Y, Z)   X.real = Y; X.imag = Z
  24. #define         cmpl(X)         X.imag = - X.imag
  25. #define         csub(X, Y, Z)   Z.real = X.real - Y.real;\
  26.                                 Z.imag = X.imag - Y.imag
  27. #define         cmul(X, Y, Z)   { \
  28.                                 struct cmplx_num cmplx_tmp; \
  29.                                 cmplx_tmp.real = X.real * Y.real - X.imag * Y.imag ; \
  30.                                 cmplx_tmp.imag = X.real * Y.imag + X.imag * Y.real ; \
  31.                                 Z.real = cmplx_tmp.real; \
  32.                                 Z.imag = cmplx_tmp.imag; \
  33.                                 }
  34. #define         cdiv(X, Y, Z)   { \
  35.                                 struct cmplx_num cmplx_tmp; \
  36.                                 double den; \
  37.                                 den = sqr( Y.real ) + sqr( Y.imag ); \
  38.                                 cmplx_tmp.real = X.real * Y.real + X.imag * Y.imag ; \
  39.                                 cmplx_tmp.imag = X.imag * Y.real - X.real * Y.imag ; \
  40.                                 Z.imag = cmplx_tmp.imag / den; \
  41.                                 Z.real = cmplx_tmp.real / den; \
  42.                                 }
  43.  
  44. /* Common Recurrisions used in the algorithm */
  45. #define         recurs10(A)     A*A + sqrt( A*A*A*A - 1 )
  46. #define         recurs12(Y, A)  ( Y + 1/Y )/( 2.0 * A)
  47. #define         recurs16(J)     sqrt( 0.5 * ( J + 1/J ) )
  48. #define         recurs18(JS)    JS + sqrt( JS*JS + 1 )
  49. #define         recurs20(P, A, PP)      { \
  50.                                 struct cmplx_num one, A_tmp, this_tmp; \
  51.                                 cset( one, 1.0, 0.0); \
  52.                                 cdiv( one, P, this_tmp); \
  53.                                 csub( P, this_tmp, this_tmp); \
  54.                                 cset( A_tmp, ( 1 / ( 2.0 * A ) ), 0.0); \
  55.                                 cmul( this_tmp, A_tmp, PP ); \
  56.                                 }
  57. main()
  58. {
  59. double lamda, epsilon;
  60. int M, i;
  61. double A0, A1, A2, A3, A4;
  62. double Y4[MX], Y3[MX], Y2[MX], Y1[MX], Y0[MX];
  63. double S0, S1, S2, S3, S4, S5;
  64. double J0, J1, J2, J3, J4, J5;
  65. double SJ0, SJ1, SJ2, SJ3, SJ4, SJ5;
  66. struct cmplx_num num;
  67. double den, K_coeff;
  68. struct cmplx_num P5[MX], P4[MX], P3[MX], P2[MX], P1[MX], P0[MX];
  69.  
  70. /* Fetch necessary input data */
  71. fprintf( stderr, "Enter value of Lamda ->");
  72. scanf( "%f", &lamda );
  73. fprintf( stderr, "Enter value of Epsilon ->");
  74. scanf( "%f", &epsilon );
  75. do{
  76.         fprintf( stderr, "Enter Order (odd) of system ->");
  77.         scanf( "%d", &M );
  78.         } while( M == 2 * (int)(M/2));
  79.  
  80.  
  81. /* Calcualtes the zeros */
  82. A0 = sqrt( lamda );
  83. A1 = recurs10( A0 );
  84. A2 = recurs10( A1 );
  85. A3 = recurs10( A2 );
  86. A4 = recurs10( A3 );
  87.  
  88. for( i = 1; i <= M/2 ; i++ ) {
  89.         Y4[i] = A4 / cos( (double)( 2.0 * i - 1.0 ) * PIE * 0.5 / (double)M );
  90.         Y3[i] = recurs12( Y4[i], A3);
  91.         Y2[i] = recurs12( Y3[i], A2);
  92.         Y1[i] = recurs12( Y2[i], A1);
  93.         Y0[i] = recurs12( Y1[i], A0);
  94.  
  95.         printf( "ZC  0.0  %f ; \n", Y0[i] );
  96.         }
  97.  
  98. /* Calculates the Poles */
  99. J4 = exp( (double)(M - 1) * log( 2.0 ) + (double)( M ) * log( A4 ) );
  100. J3 = recurs16( J4 );
  101. J2 = recurs16( J3 );
  102. J1 = recurs16( J2 );
  103. J0 = recurs16( J1 );
  104.  
  105. fprintf( stderr, "Stop Band Atten = %f dB \n", 10 * log10( 1 + sqr( epsilon * sqr( J0 ))));
  106. SJ0 = 1 / epsilon ;
  107. SJ1 = J1 * (recurs18( SJ0 ));
  108. SJ2 = J2 * (recurs18( SJ1 ));
  109. SJ3 = J3 * (recurs18( SJ2 ));
  110. S4  = 2 * SJ3 ;
  111.  
  112. S5 = exp( log( J4 / S4 + sqrt( ( J4 / S4 ) * ( J4 / S4 ) + 1 )) / (double)( M ) );
  113.  
  114.  
  115. for( i = 0; i <= M/2 ; i++ ) {
  116.         double theta;
  117.  
  118.         theta = (double)i * PIE / (double)(M);
  119.         cset( P5[i], (S5 * cos( theta) ), (S5 * sin( theta) ) );
  120.  
  121.         recurs20( P5[i], A4, P4[i]);
  122.         recurs20( P4[i], A3, P3[i]);
  123.         recurs20( P3[i], A2, P2[i]);
  124.         recurs20( P2[i], A1, P1[i]);
  125.         recurs20( P1[i], A0, P0[i]);
  126.  
  127.         if( fabs( P0[i].imag / P0[i].real ) < 0.0000001 ) {
  128.                 printf( "P  %f  ; \n", P0[i].real);
  129.                 }
  130.         else printf( "PC  %f  %f ; \n", P0[i].real, P0[i].imag );
  131.         }
  132.  
  133.  
  134. /* Calcualte the scaling constant */
  135. cset( num, 1.0, 0.0);
  136. for( i = 0; i <= M/2 ; i++ ) {
  137.         if( fabs( P0[i].imag / P0[i].real ) < 0.0000001 ) { cmul( num, P0[i], num); }
  138.         else {
  139.                 cmul( num, P0[i], num);
  140.                 cmpl( P0[i]);
  141.                 cmul( num, P0[i], num);
  142.                 }
  143.         }
  144.  
  145. den = 1.0 ;
  146. for( i = 1; i <= M/2 ; i++ ) {
  147.         den *= sqr( Y0[i] );
  148.         }
  149. K_coeff = fabs(num.real / den);
  150. printf( "K  %f ; \n", K_coeff );
  151. }
  152.  
  153.  
  154.  
  155. #undef          scanf
  156. /* A lousy fucking patch for micro-soft C.  It doesn't
  157. *  seem to read floating point numbers properly.      
  158. */
  159. my_scanf( str, reg)
  160. char *str;
  161. union {
  162.         int x;
  163.         double y;
  164.         char t;
  165.         } *reg;
  166.         
  167. {
  168. double atof();
  169. if( !strcmp( str, "%f") ) {
  170.         char temp[40];
  171.         scanf( "%s", temp);
  172.         reg->y = atof( temp);
  173.         }
  174. else scanf( str, reg);
  175. }
  176.