📄 iirbcf.c
字号:
#include"math.h"
#include"stdio.h"
#include"stdlib.h"
//#include"gain.c"
void iirbcf(int ifilt,int band,int ns,int n,double f1,double f2,
double f3,double f4,double db,double b[],double a[])
{
int k;
double omega,lamda,epslon,fl,fh;
double d[5],c[5];
void chebyi(),chebyii(),bwtf();
double cosh1(),warp(),bpsub(),omin();
void fblt();
if((band==1)||(band==4)) fl=f1;
if((band==2)||(band==3)) fl=f2;
if(band<=3) fh=f3;
if(band==4) fh=f4;
if(ifilt<3)
{
switch(band)
{
case 1:
case 2:
{omega=warp(f2)/warp(f1);
break;
}
case 3:
{
omega=omin(bpsub(warp(f1),fh,fl),bpsub(warp(f4),fh,fl));
break;
}
case 4:
{
omega=omin(1.0/bpsub(warp(f2),fh,fl),1.0/bpsub(warp(f3),fh,fl));
}
}
lamda=pow(10.0,(db/20.0));
epslon=lamda/cosh(2*ns*cosh1(omega));
}
for(k=0;k<ns;k++)
{
switch(ifilt)
{
case 1:
{
chebyi(2*ns,k,4,epslon,d,c);
break;
}
case 2:
{
chebyii(2*ns,k,4,omega,lamda,d,c);
break;
}
case 3:
{
bwtf(2*ns,k,4,d,c);
break;
}
}
fblt(d,c,n,band,fl,fh,&b[k*(n+1)+0],&a[k*(n+1)+0]);
}
}
static double cosh1(double x)
{
double z;
z=log(x+sqrt(x*x-1.0));
return(z);
}
static double warp(double f)
{
double pi,z;
pi=4.0*atan(1.0);
z=tan(pi*f);
return(z);
}
static double bpsub(double om,double fh,double fl)
{
double z;
z=(om*om-warp(fh)*warp(fl))/((warp(fh)-warp(fl))*om);
return(z);
}
static double omin(double om1,double om2)
{
double z,z1,z2;
z1=fabs(om1);
z2=fabs(om2);
z=(z1<z2)?z1:z2;
return(z);
}
static void bwtf(int ln,int k,int n,double d[],double 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;
}
}
static void chebyi(int ln,int k,int n,double ep,double d[],double c[])
{
int i;
double pi,gam,omega,sigma;
pi=4.0*atan(1.0);
gam=pow(((1.0+sqrt(1.0+ep*ep))/ep),1.0/ln);
sigma=0.5*(1.0/gam-gam)*sin((2*(k+1)-1)*pi/(2*ln));
omega=0.5*(1.0/gam+gam)*cos((2*(k+1)-1)*pi/(2*ln));
for(i=0;i<=n;i++)
{
d[i]=0.0;
c[i]=0.0;
}
if(((ln%2)==1)&&((k+1)==(ln+1)/2))
{
d[0]=-sigma;
c[0]=d[0];
c[1]=1.0;
}
else
{
c[0]=sigma*sigma+omega*omega;
c[1]=-2.0*sigma;
c[2]=1.0;
d[0]=c[0];
if(((ln%2)==0)&&(k==0))
d[0]=d[0]/sqrt(1.0+ep*ep);
}
}
static void chebyii(int ln,int k,int n,double ws,double att,double d[],double c[])
{
int i;
double pi,gam,alpha,beta,sigma,omega,scln,scld;
pi=4.0*atan(1.0);
gam=pow((att+sqrt(att*att-1.0)),1.0/ln);
alpha=0.5*(1.0/gam-gam)*sin((2*(k+1)-1)*pi/(2*ln));
beta=0.5*(1.0/gam+gam)*cos((2*(k+1)-1)*pi/(2*ln));
sigma=ws*alpha/(alpha*alpha+beta*beta);
omega=-1.0*ws*beta/(alpha*alpha+beta*beta);
for(i=0;i<=n;i++)
{
d[i]=0.0;
c[i]=0.0;
}
if(((ln%2)==1)&&((k+1)==(ln+1)/2))
{
d[0]=-1.0*sigma;
c[0]=d[0];
c[1]=1.0;
}
else
{
scln=sigma*sigma+omega*omega;
scld=pow((ws/cos((2*(k+1)-1)*pi/(2*ln))),2);
d[0]=scln*scld;
d[2]=scln;
c[0]=d[0];
c[1]=-2.0*sigma*scld;
c[2]=scld;
}
}
static void fblt(double d[],double c[],int n,
int band,double fln,double fhn,double b[],double 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;
switch(band)
{
case 1:
case 2:
{
n2=m;
n1=n2+1;
if(band==2)
{
for(i=0;i<=m/2;i++)
{
tmp=d[i];
d[i]=d[m-i];
d[m-i]=tmp;
tmp=c[i];
c[i]=c[m-i];
c[m-i]=tmp;
}
}
for(i=0;i<=m;i++)
{
d[i]=d[i]/pow(w1,i);
c[i]=c[i]/pow(w1,i);
}
break;
}
case 3:
case 4:
{
n2=2*m;
n1=n2+1;
work=malloc(n1*n1*sizeof(double));
w2 = tan(pi*fhn);
w=w2-w1;
w0=w1*w2;
if(band==4)
{
for(i=0;i<=m/2;i++)
{
tmp=d[i];
d[i]=d[m-i];
d[m-i]=tmp;
tmp=c[i];
c[i]=c[m-i];
c[m-i]=tmp;
}
}
for(i=0;i<=n2;i++)
{
work[0*n1+i]=0.0;
work[1*n1+i]=0.0;
}
for(i=0;i<=m;i++)
{
tmpd=d[i]*pow(w,(m-i));
tmpc=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(int i1,int 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(double d[],double c[],double b[],double a[],int n)
{
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);
}
void gainc(double b[],double a[],int m,int n,
double x[],double y[],int len,int sign)
{
int i,k;
double ar,ai,br,bi,zr,zi,im,re,den,numr,numi,freq,temp;
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);
br=0.0;
bi=0.0;
for(i=m;i>0;i--)
{
re=br;
im=bi;
br=(re+b[i])*zr-im*zi;
bi=(re+b[i])*zi+im*zr;
}
ar=0.0;
ai=0.0;
for(i=n;i>0;i--)
{
re=ar;
im=ai;
ar=(re+a[i])*zr-im*zi;
ai=(re+a[i])*zi+im*zr;
}
br=br+b[0];
ar=ar+1.0;
numr=ar*br+ai*bi;
numi=ar*bi-ai*br;
den=ar*ar+ai*ai;
x[k]=numr/den;
y[k]=numi/den;
switch(sign)
{
case 1:
{
temp=sqrt(x[k]*x[k]+y[k]*y[k]);
y[k]=atan2(y[k],x[k]);
x[k]=temp;
break;
}
case 2:
{
temp=x[k]*x[k]+y[k]*y[k];
y[k]=atan2(y[k],x[k]);
x[k]=10.0*log10(temp);
}
}
}
}
/* ifilt:滤波器的类型(取值1，2，3 分别对应切比雪夫、逆切比雪夫和巴特沃兹滤波器);
band:滤波器的通带(取值1、2、3、4分别对应低通、高通、带通和带阻滤波器);
ns:滤波器的n阶节数;
fs:采样频率;
db:阻带衰减; fname:幅频响应的文件名
对于低通和高通滤波器 fc:通带边界频率 fr:阻带边界频率
对于带通和带阻滤波器 flc:通带下边界频率; fhc:通带上边界频率; fls:阻带下边界频率;
fhs:阻带上边界频率*/
main()
{
int i,k,n,ns,ifilt,band;
double a[50],b[50],x[300],y[300];
double f1,f2,f3,f4,fc,fr,fs,flc,fls,fhc,fhs,freq,db;
char fname[40];
FILE *fp;
printf("enter 1 for Chebyshev I,2 for Chebyshev II,3 for Butterworth\n ");
scanf("%d",&ifilt);
printf("enter 1 for lowpass,2 for highpass,3 for bandpass,4 for bandstop\n");
scanf("%d",&band);
n=(band<=2)?2:4;
printf("enter the number of filter section\n");
scanf("%d",&ns);
printf("enter sample frequence fs\n");
scanf("%f",&fs);
if(ifilt<=2)
{
switch (band)
{
case 1 :
case 2:
{
printf("enter passband edge frequency fc\n");
scanf("%f",&fc);
printf("enter stopband edge frequency fr\n");
scanf("%f",&fr);
if(band ==1)
{
f1=fc;
f2=fr;
}
else
{
f1=fr;
f2=fc;
}
f2=f4=0.0;
break;
}
case 3:
case 4:
{
printf("enter the lower passband edge frequency flc\n");
scanf("%f",&flc);
printf("enter the higher passband edge frequency fhc\n");
scanf("%f",&fhc);
printf("enter the lower stopband edge frequency fls\n");
scanf("%f",&fls);
printf("enter the higher stopband edge frequency fhs\n");
scanf("%f",&fhs);
if(band == 3)
{
f1=fls;
f2=flc;
f3=fhc;
f4=fhs;
}
else
{
f1=flc;
f2=fls;
f3=fhs;
f4=fhc;
}
}
}
printf("enter stopband attenuation (DB)\n");
scanf("%f",&db);
}
else
{
switch(band)
{
case 1:
case 2:
{
printf("enter passband edge frequency fc\n");
scanf("%f",&fc);
f1=f2=f3=f4=0.0;
if(band == 1)
{ f1=fc; }
else
{ f2=fc; }
break;
}
case 3:
case 4:
{
printf("enter the lower passband edge frequency flc\n");
scanf("%f",&flc);
printf("enter the higher passband edge frequency fhc\n");
scanf("%f",&fhc);
f1=f2=f3=f4=0.0;
if(band==3)
{
f2=flc;
f3=fhc;
}
else
{
f1=flc;
f4=fhc;
}
}
}
}
f1=f1/fs;
f2=f2/fs;
f3=f3/fs;
f4=f4/fs;
iirbcf(ifilt,band, ns,n,f1,f2,f3,f4,db,b,a);
for(k=0;k<ns;k++)
{
printf("\nsection %d\n\n",k+1);
for(i=0;i<=n;i++)
{
printf(" b[%d][%d]=%.10f",k,i,b[k*(n+1)+i]);
if(((i%2)==0)&&(i!=0))printf("\n");
}
printf("\n");
for(i=0;i<=n;i++)
{
printf("a[%d][%d]=%.10f",k,i,a[k*(n+1)+i]);
if(((i%2)==0)&&(i!=0))printf("\n");
}
printf("\n");
}
printf("\nenter file name of magnitude response\n");
scanf("%s",fname);
if((fp=fopen(fname,"w"))==NULL)
{
printf("cannot open file %s\n",fname);
exit(0);
}
gainc(b,a,n,ns,x,y,300,2);
for(i=0;i<300;i++)
{
freq=i*0.5/300.0;
fprintf(fp,"%f %f",freq,x[i]);
}
fclose(fp);
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -