c实现巴特沃斯低通滤波器

#include <stdio.h> #include <math.h> #define pi 3.1415926 void bwtf(int ln,int k,int n,double d[],double c[])//n为每节的阶数，计算出级联型的模拟滤波器的每个因式的表达式 {int i; double tmp; d[0]=1.0;//d为每节的分子，C为每节的分母 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+ln+1)*pi/(2*ln));c[2]=1.0;} } //双线形变换法中，用模拟滤波器H（S）的系数表示数字滤波器H（Z）的系数即求解阶数为二的每节H（Z）表达式 void biliner(double d[],double c[],double b[],double a[]) { int c0,T=1; double r; c0=2/T;//c是变换常数，C=2/T（T=1），b[]表示分子，a[]表示分母 r=c[0]+c[1]*c0+c[2]*c0*c0; b[0]=(d[0]+d[1]*c0+d[2]*c0*c0)/r; b[1]=(2*d[0]-2*d[2]*c0*c0)/r; b[2]=(d[0]-d[1]*c0+d[2]*c0*c0)/r; a[0]=1; a[1]=(2*c[0]-2*c[2]*c0*c0)/r; a[2]=(c[0]-c[1]*c0+c[2]*c0*c0)/r; } void gainc(double b[],double a[],int n,int ns,double x[],double y[],int len,int sign) { int i,j,k,n1; double ar,ai,br,bi,zr,zi,im,re,den,numr,numi,freq,temp; double hr,hi,tr,ti; n1=n+1; for(k=0;k<len;k++) {freq=k*0.5/(len-1); zr=cos(-8.0*atan(1.0)*freq); zi=sin(-8.0*atan(1.0)*freq); x[k]=1.0; y[k]=0.0; for(j=0;j<ns;j++) {br=0.0; bi=0.0; for(i=n;i>0;i--) {re=br; im=bi; br=(re+b[j*n1+i])*zr-im*zi; bi=(re+b[j*n1+i])*zi+im+zr; } ar=0.0; ai=0.0; for(i=n;i>0;i--) {re=ar; im=ai; ar=(re+a[j*n1+i])*zr-im*zi; ai=(re+a[j*n1+i]*zi+im*zr); } br=br+b[j*n1+0]; ar=ar+1.0; numr=ar*br+ai*bi; numi=ar*bi-ai*br; den=ar*ar+ai*ai; hr=numr/den; hi=numi/den; tr=x[k]*hr-y[k]*hi; ti=x[k]*hi+y[k]*hr; x[k]=tr; y[k]=ti; } switch(sign) { case 1: { temp=sqrt(x[k]*x[k]+y[k]*y[k]); if(temp!=0.0) y[k]=atan2(y[k],x[k]); else y[k]=0.0; x[k]=temp; break; } case 2: {temp=x[k]*x[k]+y[k]*y[k]; if(temp!=0.0) {y[k]=atan2(y[k],x[k]);} else { temp=1.0e-40; y[k]=0.0; } x[k]=10.0*log10(temp); } } } } //假设所需设计的低通数字滤波器的要求是：通带内频率底于0.2*PIrad/sec时，容许的幅度误差为1DB， //即ωp=0.2*PI,αP=1DB，在0.3*PI到PI之间阻带的衰减大于15DB，即ωS=0.3*PI,αS=15DB main() {double ap,as,wp,ws,Wp,Ws,Wc,sp1,sp2,a[350],b[350],c[350],d[350],bz[350],az[350],x[300],y[300],N0,freq,m,p; int N,ln,k,n,ns,n1,i; printf("this is a lowpass filter of the butterworth\n"); printf("please enter the passband edge frequency wp and wsof the digital filter\n"); printf("wp="); wp=0.2; //scanf("%lf",&wp); printf("ws="); ws=0.3; // scanf("%lf",&ws); printf("please enter the tolerent error while the frequency<=wp\n"); printf("ap=");ap=1.0;//scanf("%lf",&ap); printf("please enter the tolerent error while the frequency>=ws\n"); printf("as=");as=15.0;//scanf("%lf",&as);// 以下为由数字滤波器的特性计算模拟滤波器的相对特性 Wp=2*tan(0.5*wp*pi);Ws=2*tan(0.5*ws*pi);//数字频率转模拟频率 sp1=sqrt((pow(10.0,0.1*ap)-1)/(pow(10.0,0.1*as)-1)); sp2=Ws/Wp; N0=-log10(sp1)/log10(sp2); if(N0>floor(N0)) N=(int)floor(N0)+1; else N=(int)floor(N0); printf("\nN=%d\n",N); m=pow(10.0,0.1*ap)-1; p=pow(m,(double)-1/(2*N)); Wc=p*Wp; ln=N;ns=ln/2;n=2; n1=n+1; for(k=1;k<=ns;k++) //ns表示节数 { bwtf(ln,k,n,d,c); biliner(d,c,b,a);//Hi(S)-Hi(Z) for(i=0;i<n1;i++) { bz[(k-1)*n1+i]=b[i]; az[(k-1)*n1+i]=a[i]; //printf("\n%lf,%lf\n",b[i],a[i]); } } gainc(bz,az,2,ns,x,y,300,2);//300为频率点的长度，2表示计算滤波器DB表示的幅频响应和相频响应 //FILE *fp; //fp=fopen("aram.dat","w"); //for(i=0;i<300;i++) //{freq=i*0.5/300.0;fprintf(fp,"%lf %lf\n",freq,x[i]);} //fclose(fp); for(i=0;i<300;i++) { freq=i*0.5/300.0; printf("%lf %lf \n",freq,x[i]); } getchar(); return 1; }