#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);
}
