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