ellf.c

巴特沃斯、切比雪夫I和椭圆滤波器设计的源程序！非常难得！ 简洁的C语言编写

/* ellf.c
 * 
 * Read ellf.doc before attempting to compile this program.
 */


#include <stdio.h>

/* size of arrays: */
#define ARRSIZ 50


/* System configurations */
#include "mconf.h"


extern double PI, PIO2, MACHEP, MAXNUM;

static double aa[ARRSIZ];
static double pp[ARRSIZ];
static double y[ARRSIZ];
static double zs[ARRSIZ];
cmplx z[ARRSIZ];
static double wr = 0.0;
static double cbp = 0.0;
static double wc = 0.0;
static double rn = 8.0;
static double c = 0.0;
static double cgam = 0.0;
static double scale = 0.0;
double fs = 1.0e4;
static double dbr = 0.5;
static double dbd = -40.0;
static double f1 = 1.5e3;
static double f2 = 2.0e3;
static double f3 = 2.4e3;
double dbfac = 0.0;
static double a = 0.0;
static double b = 0.0;
static double q = 0.0;
static double r = 0.0;
static double u = 0.0;
static double k = 0.0;
static double m = 0.0;
static double Kk = 0.0;
static double Kk1 = 0.0;
static double Kpk = 0.0;
static double Kpk1 = 0.0;
static double eps = 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 m1 = 0.0;
static double m1p = 0.0;
static double cang = 0.0;
static double sang = 0.0;
static double bw = 0.0;
static double ang = 0.0;
double fnyq = 0.0;
static double ai = 0.0;
static double pn = 0.0;
static double an = 0.0;
static double gam = 0.0;
static double cng = 0.0;
double gain = 0.0;
static int lr = 0;
static int nt = 0;
static int i = 0;
static int j = 0;
static int jt = 0;
static int nc = 0;
static int ii = 0;
static int ir = 0;
int zord = 0;
static int icnt = 0;
static int mh = 0;
static int jj = 0;
static int jh = 0;
static int jl = 0;
static int n = 8;
static int np = 0;
static int nz = 0;
static int type = 1;
static int kind = 1;

static char wkind[] =
{"Filter kind:\n1 Butterworth\n2 Chebyshev\n3 Elliptic\n"};

static char salut[] =
{"Filter shape:\n1 low pass\n2 band pass\n3 high pass\n4 band stop\n"};

#ifdef ANSIPROT
extern double exp ( double );
extern double log ( double );
extern double cos ( double );
extern double sin ( double );
extern double sqrt ( double );
extern double fabs ( double );
extern double asin ( double );
extern double atan ( double );
extern double atan2 ( double, double );
extern double pow ( double, double );
extern double cabs ( cmplx *z );
extern void cadd ( cmplx *a, cmplx *b, cmplx *c );
extern void cdiv ( cmplx *a, cmplx *b, cmplx *c );
extern void cmov ( void *a, void *b );
extern void cmul ( cmplx *a, cmplx *b, cmplx *c );
extern void cneg ( cmplx *a );
extern void csqrt ( cmplx *z, cmplx *w );
extern void csub ( cmplx *a, cmplx *b, cmplx *c );
extern double ellie ( double phi, double m );
extern double ellik ( double phi, double m );
extern double ellpe ( double x );
extern int ellpj ( double, double, double *, double *, double *, double * );
extern double ellpk ( double x );
int getnum ( char *line, double *val );
double cay ( double q );
int lampln ( void );
int spln ( void );
int xfun ( void );
int zplna ( void );
int zplnb ( void );
int zplnc ( void );
int quadf ( double, double, int );
double response ( double, double );
#else
double exp(), log(), cos(), sin(), sqrt();
double ellpk(), ellik(), asin(), atan(), atan2(), pow();
double cay(), cabs();
double response();
int lampln(), spln(), xfun(), zplna(), zplnb(), zplnc(), quadf();
#define fabs(x) ( (x) < 0 ? -(x) : (x) )
#endif

int main()
{
char str[80];

dbfac = 10.0/log(10.0);

top:

printf( "%s ? ", wkind );	/* ask for filter kind */
gets( str );
sscanf( str, "%d", &kind );
printf( "%d\n", kind );
if( (kind <= 0) || (kind > 3) )
	exit(0);

printf( "%s ? ", salut );	/* ask for filter type */
gets( str );
sscanf( str, "%d", &type );
printf( "%d\n", type );
if( (type <= 0) || (type > 4) )
	exit(0);

getnum( "Order of filter", &rn ); /* see below for getnum() */
n = rn;
if( n <= 0 )
	{
specerr:
	printf( "? Specification error\n" );
	goto top;
	}
rn = n;	/* ensure it is an integer */
if( kind > 1 ) /* not Butterworth */
	{
	getnum( "Passband ripple, db", &dbr );
	if( dbr <= 0.0 )
		goto specerr;
	if( kind == 2 )
		{
/* For Chebyshev filter, ripples go from 1.0 to 1/sqrt(1+eps^2) */
		phi = exp( 0.5*dbr/dbfac );

		if( (n & 1) == 0 )
			scale = phi;
		else
			scale = 1.0;
		}
	else
		{ /* elliptic */
	eps = exp( dbr/dbfac );
	scale = 1.0;
	if( (n & 1) == 0 )
		scale = sqrt( eps );
	eps = sqrt( eps - 1.0 );
		}
	}

getnum( "Sampling frequency", &fs );
if( fs <= 0.0 )
	goto specerr;

fnyq = 0.5 * fs;

getnum( "Passband edge", &f2 );
if( (f2 <= 0.0) || (f2 >= fnyq) )
	goto specerr;

if( (type & 1) == 0 )
	{
	getnum( "Other passband edge", &f1 );
	if( (f1 <= 0.0) || (f1 >= fnyq) )
		goto specerr;
	}
else
	{
	f1 = 0.0;
	}

if( f2 < f1 )
	{
	a = f2;
	f2 = f1;
	f1 = a;
	}
if( type == 3 )	/* high pass */
	{
	bw = f2;
	a = fnyq;
	}
else
	{
	bw = f2 - f1;
	a = f2;
	}
/* Frequency correspondence for bilinear transformation
 *
 *  Wanalog = tan( 2 pi Fdigital T / 2 )
 *
 * where T = 1/fs
 */
ang = bw * PI / fs;
cang = cos( ang );
c = sin(ang) / cang; /* Wanalog */
if( kind != 3 )
	{
	wc = c;
/*printf( "cos( 1/2 (Whigh-Wlow) T ) = %.5e, wc = %.5e\n", cang, wc );*/
	}


if( kind == 3 )
	{ /* elliptic */
	cgam = cos( (a+f1) * PI / fs ) / cang;
	getnum( "Stop band edge or -(db down)", &dbd );
	if( dbd > 0.0 )
		f3 = dbd;
	else
		{ /* calculate band edge from db down */
		a = exp( -dbd/dbfac );
		m1 = eps/sqrt( a - 1.0 );
		m1 *= m1;
		m1p = 1.0 - m1;
		Kk1 = ellpk( m1p );
		Kpk1 = ellpk( m1 );
		q = exp( -PI * Kpk1 / (rn * Kk1) );
		k = cay(q);
		if( type >= 3 )
			wr = k;
		else
			wr = 1.0/k;
		if( type & 1 )
			{
			f3 = atan( c * wr ) * fs / PI;
			}
		else
			{
			a = c * wr;
			a *= a;
			b = a * (1.0 - cgam * cgam) + a * a;
			b = (cgam + sqrt(b))/(1.0 + a);
			f3 = (PI/2.0 - asin(b)) * fs / (2.0*PI);
			}
		}
switch( type )
	{
	case 1:
		if( f3 <= f2 )
			goto specerr;
		break;

	case 2:
		if( (f3 > f2) || (f3 < f1) )
			break;
		goto specerr;

	case 3:
		if( f3 >= f2 )
			goto specerr;
		break;

	case 4:
		if( (f3 <= f1) || (f3 >= f2) )
			goto specerr;
		break;
	}
ang = f3 * PI / fs;
cang = cos(ang);
sang = sin(ang);
	
if( type & 1 )
	{
	wr = sang/(cang*c);
	}
else
	{
	q = cang * cang  -  sang * sang;
	sang = 2.0 * cang * sang;
	cang = q;
	wr = (cgam - cang)/(sang * c);
	}

if( type >= 3 )
	wr = 1.0/wr;
if( wr < 0.0 )
	wr = -wr;
y[0] = 1.0;
y[1] = wr;
cbp = wr;

if( type >= 3 )
	y[1] = 1.0/y[1];

if( type & 1 )
	{
	for( i=1; i<=2; i++ )
		{
		aa[i] = atan( c * y[i-1] ) * fs / PI ;
		}
	printf( "pass band %.9E\n", aa[1] );
	printf( "stop band %.9E\n", aa[2] );
	}
else
	{
	for( i=1; i<=2; i++ )
		{
		a = c * y[i-1];
		b = atan(a);
		q = sqrt( 1.0 + a * a  -  cgam * cgam );
#ifdef ANSIC
		q = atan2( q, cgam );
#else
		q = atan2( cgam, q );
#endif
		aa[i] = (q + b) * fnyq / PI;
		pp[i] = (q - b) * fnyq / PI;
		}
	printf( "pass band %.9E %.9E\n", pp[1], aa[1] );
	printf( "stop band %.9E %.9E\n", pp[2], aa[2] );
	}
lampln();	/* find locations in lambda plane */
if( (2*n+2) > ARRSIZ )
	goto toosml;
	}

/* 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 )
 *
 *
 * Band edges
 *            cos( Wcenter T) - cos( Wdigital T )
 *  Wanalog = -----------------------------------
 *                        sin( Wdigital T )
 */

if( kind == 2 )
	{ /* Chebyshev */
	a = PI * (a+f1) / fs ;
	cgam = cos(a) / cang;
	a = 2.0 * PI * f2 / fs;
	cbp = (cgam - cos(a))/sin(a);
	}
if( kind == 1 )
	{ /* Butterworth */
	a = PI * (a+f1) / fs ;
	cgam = cos(a) / cang;
	a = 2.0 * PI * f2 / fs;
	cbp = (cgam - cos(a))/sin(a);
	scale = 1.0;
	}

spln();		/* find s plane poles and zeros */

if( ((type & 1) == 0) && ((4*n+2) > ARRSIZ) )
	goto toosml;

zplna();	/* convert s plane to z plane */
zplnb();
zplnc();
xfun(); /* tabulate transfer function */
goto top;

toosml:
printf( "Cannot continue, storage arrays too small\n" );
goto top;
}


int lampln()
{

wc = 1.0;
k = wc/wr;
m = k * k;
Kk = ellpk( 1.0 - m );
Kpk = ellpk( m );
q = exp( -PI * rn * Kpk / Kk );	/* the nome of k1 */
m1 = cay(q); /* see below */
/* Note m1 = eps / sqrt( A*A - 1.0 ) */
a = eps/m1;
a =  a * a + 1;
a = 10.0 * log(a) / log(10.0);
printf( "dbdown %.9E\n", a );
a = 180.0 * asin( k ) / PI;
b = 1.0/(1.0 + eps*eps);
b = sqrt( 1.0 - b );
printf( "theta %.9E, rho %.9E\n", a, b );
m1 *= m1;
m1p = 1.0 - m1;
Kk1 = ellpk( m1p );
Kpk1 = ellpk( m1 );
r = Kpk1 * Kk / (Kk1 * Kpk);
printf( "consistency check: n= %.14E\n", r );
/*   -1
 * sn   j/eps\m  =  j ellik( atan(1/eps), m )
 */
b = 1.0/eps;
phi = atan( b );
u = ellik( phi, m1p );
printf( "phi %.7e m %.7e u %.7e\n", phi, m1p, u );
/* consistency check on inverse sn */
ellpj( u, m1p, &sn, &cn, &dn, &phi );
a = sn/cn;
printf( "consistency check: sn/cn = %.9E = %.9E = 1/eps\n", a, b );
u = u * Kk / (rn * Kk1);	/* or, u = u * Kpk / Kpk1 */
return 0;
}




/* calculate s plane poles and zeros, normalized to wc = 1 */
int spln()
{
for( i=0; i<ARRSIZ; i++ )
	zs[i] = 0.0;
np = (n+1)/2;
nz = 0;
if( kind == 1 )
	{
/* Butterworth poles equally spaced around the unit circle
 */
	if( n & 1 )
		m = 0.0;
	else
		m = PI / (2.0*n);
	for( i=0; i<np; i++ )
		{	/* poles */
		lr = i + i;
		zs[lr] = -cos(m);
		zs[lr+1] = sin(m);
		m += PI / n;
		}	
	/* high pass or band reject
	 */
	if( type >= 3 )
		{
		/* map s => 1/s
		 */
		for( j=0; j<np; j++ )
			{
			ir = j + j;
			ii = ir + 1;
			b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
			zs[ir] = zs[ir] / b;
			zs[ii] = zs[ii] / b;
			}
		/* The zeros at infinity map to the origin.
		 */
		nz = np;
		if( type == 4 )
			{
			nz += n/2;
			}
		for( j=0; j<nz; j++ )
			{
			ir = ii + 1;
			ii = ir + 1;
			zs[ir] = 0.0;
			zs[ii] = 0.0;
			}
		}
	}
if( kind == 2 )
	{
	/* For Chebyshev, find radii of two Butterworth circles
	 * See Gold & Rader, page 60
	 */
	rho = (phi - 1.0)*(phi+1);  /* rho = eps^2 = {sqrt(1+eps^2)}^2 - 1 */
	eps = sqrt(rho);
	/* sqrt( 1 + 1/eps^2 ) + 1/eps  = {sqrt(1 + eps^2)  +  1} / eps
	 */
	phi = (phi + 1.0) / eps;
	phi = pow( phi, 1.0/rn );  /* raise to the 1/n power */
	b = 0.5 * (phi + 1.0/phi); /* y coordinates are on this circle */
	a = 0.5 * (phi - 1.0/phi); /* x coordinates are on this circle */
	if( n & 1 )
		m = 0.0;
	else
		m = PI / (2.0*n);
	for( i=0; i<np; i++ )
		{	/* poles */
		lr = i + i;
		zs[lr] = -a * cos(m);
		zs[lr+1] = b * sin(m);
		m += PI / n;
		}	
	/* high pass or band reject
	 */
	if( type >= 3 )
		{
		/* map s => 1/s
		 */
		for( j=0; j<np; j++ )
			{
			ir = j + j;
			ii = ir + 1;
			b = zs[ir]*zs[ir] + zs[ii]*zs[ii];
			zs[ir] = zs[ir] / b;
			zs[ii] = zs[ii] / b;
			}
		/* The zeros at infinity map to the origin.
		 */
		nz = np;
		if( type == 4 )
			{
			nz += n/2;
			}
		for( j=0; j<nz; j++ )
			{
			ir = ii + 1;
			ii = ir + 1;
			zs[ir] = 0.0;
			zs[ii] = 0.0;
			}
		}
	}
if( kind == 3 )
	{
	nz = n/2;
	ellpj( u, 1.0-m, &sn1, &cn1, &dn1, &phi1 );
	for( i=0; i<ARRSIZ; i++ )