巴特沃兹滤波器 C语言

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