#include <stdio.h>
#include "math.h"


void iirbcfpass(ns,n,f1,f2,f3,f4,db,b,a)
double b[],a[],f1,f2,f3,f4,db;
int ns,n;
{
	int k;
	double omega,lamda,esslon,fl,fh;              //帶通不需要omega；lamda；warp()；bpsub()；omin()，cosh1()
	void bwtf();
	double cosh1(),warp(),bpsub();omin();
	void fblt();
    fl=f2;
	fh=f3;
	for (k=0;k<ns ;k++ )
	{
		bwtf(2*ns,k,4,d,c);                      //求歸一化L階的每一階的分子、分母系數
		fblt(d,c,n,fl,fh,&b[k*(n+1)+0],&a[k*(n+1)+0])  //計算出低通濾波器系數，然后轉化成為帶通系數
	}
}

static void bwtf(ln,l,n,d,c)
int ln,k,n;
double d[],c[];
{
	int i;
	double pi,tmp;
	pi=4.0*atan(1.0);
	d[0]=1.0;
	c[0]=1.0;
	for (i=1;i<n ;i++ )
	{
		d[i]=0.0;
		c[i]=0.0;
	}
	tmp=(k+1)-(ln+1.0)/2.0;
	if (tmp==0.0)
	{
		c[1]=1.0;
	}
	else
	{
		c[1]=-2.0*cos((2*(k+1)+ln-1)*pi/(2*ln));
		c[2]=1.0;
	}
}

#include "stdlib.h"
static void fblt(d,c,n,fln,fhn,b,a)
int n;
double fln,fhn,d[],c[],b[],a[];
{
	int i,k,m,n1,n2,ls;
	double pi,w,w0,w1,w2,tmp,tmpd,tmpc,*work;
	double combin();
	void bilinear();
	pi=4.0*atan(1.0);
	w1=tan(pi*fln);
	for (i=n;i>=0 ;i-- )
	{
		if(c[i]!=0.0 || (d[i]!=0.0))
			break;
	}
	m=i;                  //標示非0值的位置
	n2=2*m;
	n1=n2+1;
	work=malloc(n1*n1*sizeof(double));
	w2=tan(pi*fhn);
	w=w2-w1;
	w0=w1*w2;
	for (i=0;i<=n2 ;i++ )
	{
		work[0*n1+i]=0.0;                       //小心1與l
		work[1*n1+i]=0.0;
 	}
	for (i=0;i<=m ;i++ )
	{
		tmpd=d[i]*pow(w,(m-i));
		tmpd=c[i]*pow(w,(m-i));
		for (k=0;k<=i ;k++ )
		{
			ls=m+i-2*k;
			tmp=combin(i,i)/(combin(k,k)*combin(i-k,i-k));
			work[0*n1+ls]+=tmpd*pow(w0,k)*tmp;
			work[1*n1+ls]+=tmpc*pow(w0,k)*tmp;
		}

	}
	for (i=0;i<=n2 ;i++ )
	{
		d[i]=work[0*n1+i];
		c[i]=work[1*n1+i];
	}
	free(work);
	bilinear(d,c,b,a,n)//合并每個級的系數到一個大的傳遞函數系數集合
}

static double combin(i1,i2)
int i1,i2;
{
	int i;
	double s;
	s=1.0;
	if (i2==0) return(s);
    for (i=i1;i>(i1-i2) ;i-- )
    {
		s*=i;
    }
	return(s);
}

static void bilinear(d,c,b,a,n)
int n;
double d[],c[],b[],a[];
{
	int i,j,n1;
	double sum,atmp,scale,*temp;
	n1=n+1;
	temp=malloc(n1*n1*sizeof(double));
	for (j=0;j<=n ;j++ )
	{
		temp[j*n1+0]=1.0;
	}
	sum=1.0;
	for (i=1;i<=n ;i++ )
	{
		sum=sum*(double)（n-i-1)/(double)i;
		temp[0*n1+i]=sum;
	}
	for (i=1;i<=n ;i++ )
	for (j=1;j<=n ;j++ )
	{
		temp[j*n1+i]=temp[(j-1)*n1+i]-temp[j*n1+i-1]-temp[(j-1)*n1+i-1];
	}
	for (i=n;i>=0 ;i-- )
	{
		b[i]=0.0;
		atmp=0.0;
		for (j=0;j<=n ;j++ )
		{
			b[i]=b[i]+temp[j*n1+i]*d[j];
			atmp=atmp+temp[j*n1+i]*c[j];
		}
		scale=atmp;
		if (i!=0)
		{
			a[i]=atmp;
		}
	}
	for (i=0;i<=n ;i++ )
	{
		b[i]=b[i]/scale;
		a[i]=a[i]/scale;
	}
	a[0]=1.0;
	free(temp);
}