This is ellip.c in view mode; [Download] [Up]
/********************************************************************
* ELLIPTICAL FILTER DESIGN
* Given inputs of lamda (frequency transistion range), epsilon ( a measuer
* of the pass band ripple, and the order of the sysetm, this program
* will generate the necessary poles and zeros. It is based on a algorithmn
* by Sindy Darlington.
* steve punte 10-20-87
*/
#include <stdio.h>
#include <math.h>
#define TRUE 1
#define FALSE 0
#define MX 50 /* Max number of poles or zeros */
#define PIE 3.141592653
#define scanf my_scanf
struct cmplx_num{
double real;
double imag;
};
#define sqr(X) ( X )*( X )
#define cset(X, Y, Z) X.real = Y; X.imag = Z
#define cmpl(X) X.imag = - X.imag
#define csub(X, Y, Z) Z.real = X.real - Y.real;\
Z.imag = X.imag - Y.imag
#define cmul(X, Y, Z) { \
struct cmplx_num cmplx_tmp; \
cmplx_tmp.real = X.real * Y.real - X.imag * Y.imag ; \
cmplx_tmp.imag = X.real * Y.imag + X.imag * Y.real ; \
Z.real = cmplx_tmp.real; \
Z.imag = cmplx_tmp.imag; \
}
#define cdiv(X, Y, Z) { \
struct cmplx_num cmplx_tmp; \
double den; \
den = sqr( Y.real ) + sqr( Y.imag ); \
cmplx_tmp.real = X.real * Y.real + X.imag * Y.imag ; \
cmplx_tmp.imag = X.imag * Y.real - X.real * Y.imag ; \
Z.imag = cmplx_tmp.imag / den; \
Z.real = cmplx_tmp.real / den; \
}
/* Common Recurrisions used in the algorithm */
#define recurs10(A) A*A + sqrt( A*A*A*A - 1 )
#define recurs12(Y, A) ( Y + 1/Y )/( 2.0 * A)
#define recurs16(J) sqrt( 0.5 * ( J + 1/J ) )
#define recurs18(JS) JS + sqrt( JS*JS + 1 )
#define recurs20(P, A, PP) { \
struct cmplx_num one, A_tmp, this_tmp; \
cset( one, 1.0, 0.0); \
cdiv( one, P, this_tmp); \
csub( P, this_tmp, this_tmp); \
cset( A_tmp, ( 1 / ( 2.0 * A ) ), 0.0); \
cmul( this_tmp, A_tmp, PP ); \
}
main()
{
double lamda, epsilon;
int M, i;
double A0, A1, A2, A3, A4;
double Y4[MX], Y3[MX], Y2[MX], Y1[MX], Y0[MX];
double S0, S1, S2, S3, S4, S5;
double J0, J1, J2, J3, J4, J5;
double SJ0, SJ1, SJ2, SJ3, SJ4, SJ5;
struct cmplx_num num;
double den, K_coeff;
struct cmplx_num P5[MX], P4[MX], P3[MX], P2[MX], P1[MX], P0[MX];
/* Fetch necessary input data */
fprintf( stderr, "Enter value of Lamda ->");
scanf( "%f", &lamda );
fprintf( stderr, "Enter value of Epsilon ->");
scanf( "%f", &epsilon );
do{
fprintf( stderr, "Enter Order (odd) of system ->");
scanf( "%d", &M );
} while( M == 2 * (int)(M/2));
/* Calcualtes the zeros */
A0 = sqrt( lamda );
A1 = recurs10( A0 );
A2 = recurs10( A1 );
A3 = recurs10( A2 );
A4 = recurs10( A3 );
for( i = 1; i <= M/2 ; i++ ) {
Y4[i] = A4 / cos( (double)( 2.0 * i - 1.0 ) * PIE * 0.5 / (double)M );
Y3[i] = recurs12( Y4[i], A3);
Y2[i] = recurs12( Y3[i], A2);
Y1[i] = recurs12( Y2[i], A1);
Y0[i] = recurs12( Y1[i], A0);
printf( "ZC 0.0 %f ; \n", Y0[i] );
}
/* Calculates the Poles */
J4 = exp( (double)(M - 1) * log( 2.0 ) + (double)( M ) * log( A4 ) );
J3 = recurs16( J4 );
J2 = recurs16( J3 );
J1 = recurs16( J2 );
J0 = recurs16( J1 );
fprintf( stderr, "Stop Band Atten = %f dB \n", 10 * log10( 1 + sqr( epsilon * sqr( J0 ))));
SJ0 = 1 / epsilon ;
SJ1 = J1 * (recurs18( SJ0 ));
SJ2 = J2 * (recurs18( SJ1 ));
SJ3 = J3 * (recurs18( SJ2 ));
S4 = 2 * SJ3 ;
S5 = exp( log( J4 / S4 + sqrt( ( J4 / S4 ) * ( J4 / S4 ) + 1 )) / (double)( M ) );
for( i = 0; i <= M/2 ; i++ ) {
double theta;
theta = (double)i * PIE / (double)(M);
cset( P5[i], (S5 * cos( theta) ), (S5 * sin( theta) ) );
recurs20( P5[i], A4, P4[i]);
recurs20( P4[i], A3, P3[i]);
recurs20( P3[i], A2, P2[i]);
recurs20( P2[i], A1, P1[i]);
recurs20( P1[i], A0, P0[i]);
if( fabs( P0[i].imag / P0[i].real ) < 0.0000001 ) {
printf( "P %f ; \n", P0[i].real);
}
else printf( "PC %f %f ; \n", P0[i].real, P0[i].imag );
}
/* Calcualte the scaling constant */
cset( num, 1.0, 0.0);
for( i = 0; i <= M/2 ; i++ ) {
if( fabs( P0[i].imag / P0[i].real ) < 0.0000001 ) { cmul( num, P0[i], num); }
else {
cmul( num, P0[i], num);
cmpl( P0[i]);
cmul( num, P0[i], num);
}
}
den = 1.0 ;
for( i = 1; i <= M/2 ; i++ ) {
den *= sqr( Y0[i] );
}
K_coeff = fabs(num.real / den);
printf( "K %f ; \n", K_coeff );
}
#undef scanf
/* A lousy fucking patch for micro-soft C. It doesn't
* seem to read floating point numbers properly.
*/
my_scanf( str, reg)
char *str;
union {
int x;
double y;
char t;
} *reg;
{
double atof();
if( !strcmp( str, "%f") ) {
char temp[40];
scanf( "%s", temp);
reg->y = atof( temp);
}
else scanf( str, reg);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.