C Code for IIR Filter Design - Practical Analog and Digital Filter Design - page 255

/*====================================================
Bilinear_Transform() - use bilinear transform to
convert transfer function from s-domain to z-domain
Prototype: int Bilinear_Transform(Filt_Params *FP,
double fsamp);
Return: error value
Arguments: FP - ptr to struct holding filter params
====================================================*/
int Bilinear_Transform(Filt_Params *FP)
{ int i,j,start, /* loop counters and index */
num_quads; /* number of quad factors */
double f2,f4, /* 2 * fsamp, and 4 * fsamp^2 */
N0,N1,N2, /* numerator temp variables */
D0,D1,D2; /* denominator temp variables */
if( (!FP->acoefs) || (!FP->bcoefs) )
{ return ERR_NULL;}
/* determine some constants */
f2 = 2 * FP->fsamp;
f4 = f2 * f2;
num_quads = (FP->order + 1)/2;
/* handle first-order factor if present */
start = 0;
if(FP->order % 2)
{ N0 = FP->acoefs[2] + FP->acoefs[1] * f2;
N1 = FP->acoefs[2] - FP->acoefs[1] * f2;
D0 = FP->bcoefs[2] + FP->bcoefs[1] * f2;
D1 = FP->bcoefs[2] - FP->bcoefs[1] * f2;
FP->acoefs[0] = 1.0;
FP->acoefs[1] = N1 / N0;
FP->acoefs[2] = 0.0;
FP->bcoefs[0] = 1.0;
FP->bcoefs[1] = D1 / D0;
FP->bcoefs[2] = 0.0;
FP->gain *= (N0 / D0);
start = 1;
}
/* Handle quadratic factors. */
for(i = start; i < num_quads ;i++)
{ j = 3 * i;
N0 = FP->acoefs[j]*f4 + FP->acoefs[j+1]*f2 + FP->acoefs[j+2];
N1 = 2 * (FP->acoefs[j+2] - FP->acoefs[j]*f4);
N2 = FP->acoefs[j]*f4 - FP->acoefs[j+1]*f2 + FP->acoefs[j+2];
D0 = FP->bcoefs[j]*f4 + FP->bcoefs[j+1]*f2 + FP->bcoefs[j+2];
D1 = 2 * (FP->bcoefs[j+2] - FP->bcoefs[j]*f4);
D2 = FP->bcoefs[j]*f4 - FP->bcoefs[j+1]*f2 + FP->bcoefs[j+2];
FP->acoefs[j] = 1.0;
FP->acoefs[j+1] = N1 / N0;
FP->acoefs[j+2] = N2 / N0;
FP->bcoefs[j] = 1.0;
FP->bcoefs[j+1] = D1 / D0;
FP->bcoefs[j+2] = D2 / D0;
FP->gain *= (N0 / D0);
}
return ERR_NONE;
}
Listing G.2 Bilinear_Transform function.