Digital Signal Processing Reference
In-Depth Information
/*====================================================
Calc_Butter_Coefs() - calcs normal Butterworth coefs
Prototype: int Calc_Butter_Coefs(Filt_Params *FP);
Return: error value
Arguments: FP - ptr to struct holding filter params
====================================================*/
int Calc_Butter_Coefs(Filt_Params *FP)
{ int m,a,b; /* Loop counter and indices*/
double R,epsilon, /* Intermediate values */
theta, /* Angle location of poles */
sigma,omega; /* Real/imag pos of poles */
/* Check for NULL ptrs and zero order. */
if( !FP->acoefs)
{ return ERR_NULL;}
if( !FP->bcoefs)
{ return ERR_NULL;}
if(FP->order <= 0)
{ return ERR_VALUE;}
/* Make calculations of necessary constants. */
epsilon = sqrt( pow(10.0,-0.1*FP->apass1) - 1.0 );
R = pow(epsilon,-1.0/FP->order);
/* Initialize gain to 1.0. Start indices at 0 */
FP->gain = 1.0;
a = 0; b = 0;
/* Handle odd order if necessary. */
if(FP->order % 2)
{ FP->acoefs[a++] = 0.0;
FP->acoefs[a++] = 0.0;
FP->acoefs[a++] = R;
FP->bcoefs[b++] = 0.0;
FP->bcoefs[b++] = 1.0;
FP->bcoefs[b++] = R;
}
/* Handle all quadratic terms. */
for(m = 0;m < FP->order/2;m++)
{ /* Calc angle first, then real and imag pos. */
theta = PI*(2*m + FP->order +1) / (2 * FP->order);
sigma = R * cos(theta);
omega = R * sin(theta);
/* Set the quadratic coefs. */
FP->acoefs[a++] = 0.0;
FP->acoefs[a++] = 0.0;
FP->acoefs[a++] = sigma*sigma+omega*omega;
FP->bcoefs[b++] = 1.0;
FP->bcoefs[b++] = -2 * sigma;
FP->bcoefs[b++] = sigma*sigma+omega*omega;
}
return ERR_NONE;
}
Listing D.4
Calc_Butter_Coefs
function.
The
Ellip_Integral
function is shown in Listing D.5 and uses the
arithmetic-geometric mean method of determining the complete elliptic integral,
as defined in (2.62). It takes as an argument the modulus
k
and returns the value of
the complete elliptic integral.
MAX_TERMS
and
ERR_SMALL
have been defined
