//-------------------------------------------------------
// Digital Filter Coefficient Generator using
// the Bilinear Coefficient Formula.
// Peter J. Pupalaikis
//-------------------------------------------------------
#include <math.h>
#include <stdio.h>
#define PI 3.141592

/* Prewarp function */

double prewarp(double fc, double fs){

    double Fc;
    
    Fc = (fs/PI)*tan(PI*fc/fs);

    return Fc;
}    

// simple factorial function
static double f(double i)
{
if (i==0)
return 1;
else return i*f(i-1);
}
#define max(a,b) (((a)>(b))?(a):(b))
#define min(a,b) ((a<b)?a:b)
// computes a digital filter coefficient from the analog
// coefficients using the bilinear coefficient formula
double s_to_z
(
double *a, // s domain coefficients
int N, // index of last coefficient (order)
int n, // z domain coeffient desired
double fs // sample rate
)
{
double z_coef = 0.0;
for (int i=0; i< N; i++)
{
double acc = 0.0;
for (int k=max(n-N+i,0); k<=min(i,n); ++k)
{
acc += ( ( f(i)*f(N-i) ) /
( f(k)*f(i-k)*f(n-k)*f(N-i-n+k) ) ) *
( (k & 1) ? -1.0 : 1.0 );
} z_coef += (
a[i] *
pow(2.0*fs,i) *
acc);
}
return z_coef;
}

/* ----------------Use the BLT----------------------*/
#define N 4
int main(){


double *den=new double[N];
double *denw=new double[N];
double fs =1.0,wc,fc, wcpw,fcpw;
double A;

den[0] = 1.0; den[1] = 2.0;den[2] = 2.0;den[3] = 1.0;//Butterworth lp
wc = 1;
fc = wc/(2.0*PI);

//Prewarp analog filter critical frequencies
fcpw = prewarp(fc,fs);
wcpw = 2.0*PI*fcpw;
for(int j =0; j<N;j++){denw[j] = wcpw*den[j]; }
printf("Prewarped f cutoff: %f\n",fcpw);
printf("Prewarped w cutoff: %f\n",wcpw);
for(int n =0; n<N;n++){printf("Denominator coeficient %i: %f\n",n,denw[n]); }
//Estimate Digital FIlter Coefs
for(int k=0; k<N;k++){
 A = s_to_z(denw,N,k,fs);
 printf("%f\n",A);   
}


return 0;
}

Prewarped f cutoff: 0.173894
Prewarped w cutoff: 1.092605
Denominator coeficient 0: 1.092605
Denominator coeficient 1: 2.185210
Denominator coeficient 2: 2.185210
Denominator coeficient 3: 1.092605
22.944705
-4.370420
-10.926050
13.111260