IIR数字滤波器设计C程序.doc

IIR数字滤波器设计C程序 - 包括图形显示转帖/*IIRD_DF.c-双线性变换法设计IIR数字滤波器Infinite Impulse Response Digital Filter Design By Double Converting=包括1.Lowpass 2.Highpass 3.bandpass三部分，如果有特殊需要把其他部分删除掉即可。*/*#include */#include #include #include #include #include #include #include /*COMPLEX STRUCTURE*/typedef structdouble real,image;COMPLEX;struct rptrdouble *a;double *b;#define PI (double)(4.0*atan(1.0)#define FNSSH(x) log(x+sqrt(x*x+1)#define FNCCH(x) log(x+sqrt(x*x-1)#define FNSH1(x) (exp(x)-exp(-x)/2#define FNCH1(x) (exp(x)+exp(-x)/2/*/double *bcg(double ap,double as,double wp,double ws,int *n,double *h,int *type);struct rptr *bsf(double *c,int ni,double *f1,double *f2,int nf,struct rptr *ptr,int *no);double *pnpe(double *a,int m,int n,double *b,int *mn);double *ypmp(double *a,int m,double *b,int n,double *c,int *mn);void lowpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf);void highpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf);void bandpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf);void draw_image(double *x,int m,char *title1,char *title2,char *xdis1,char *xdis2,int dis_type);/*/void main(void)double *f1,*f2,*hwdb;double *h=NULL;int N,nf,ns,nz,i,j,k,ftype,type;COMPLEX hwdb1,hwdb2;double wp,ws,ap,as,jw,amp1,amp2;char title80,tmp20;struct rptr *ptr=NULL;/* printf(%lf,PI);getch();*/N=500;f1=(double *)calloc(4,sizeof(double);f2=(double *)calloc(4,sizeof(double);hwdb=(double *)calloc(N+2,sizeof(double);if (hwdb=NULL)printf(n Not enough memory to allocate!);exit(0);printf(n1.Lowpass 2.Highpass 3.bandpass);printf(nPlease select the filter type:);scanf(%d,&ftype);switch(ftype)case 1:lowpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);break;case 2:highpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);break;case 3:bandpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);break;default:lowpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);h=bcg(ap,as,wp,ws,&ns,h,&type);printf(nThe analog filter denominator coefficients of Ha(s):);for(i=0;i=ns;i+)printf(nb-=f,i,hi);ptr=bsf(h,ns,f1,f2,nf,ptr,&nz);printf(nThe digital filter coefficients of H(z):);printf(n(a is numerator coefficient. b is denominator coefficient.);for(i=0;iai,i,ptr-bi);printf(nnPress any key to calculate the filter response of H(z).);getch();printf(nWaitting for calculating.);/* Calculate the magnitude-frequency response*/for(k=0;k=N;k+)jw=k*PI/N;hwdb1.real=0;hwdb1.image=0;hwdb2.real=0;hwdb2.image=0;for(i=0;iai*cos(i*jw);hwdb2.real+=ptr-bi*cos(i*jw);hwdb1.image+=ptr-ai*sin(i*jw);hwdb2.image+=ptr-bi*sin(i*jw);amp1=(pow(hwdb1.real,2)+pow(hwdb1.image,2);amp2=(pow(hwdb2.real,2)+pow(hwdb2.image,2);if(amp1=0) amp1=1e-90;if(amp2=0) amp2=1e-90;hwdbk=10*log10(amp1/amp2);if(hwdbkb);free(ptr-a);free(h);free(hwdb);free(f2);free(f1);/*/void lowpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf)double fp,fr,fs;printf(nPlease input the Fp,Ap,Fr,Ar,Fs value);printf(nFp,Ap:Passband frequency(Hz) and MAXAttenuation(dB);printf(nFr,Ar:Stopband frequency(Hz) and MINAttenuation(dB);printf(nFs is the sample frequency(Hz) Lowpass filter);printf(nInput parameters Fp,Ap,Fr,Ar,Fs:);scanf(%lf,%lf,%lf,%lf,%lf,&fp,ap,&fr,ar,&fs);if(fpfr)|(*ap*ar)|(fsfr)|(*ap*ar)|(fs2*fr);*wp=tan(PI*fp/fs);*ws=tan(PI*fr/fs);*f1=-1.0;*(f1+1)=1.0;*f2=1.0;*(f2+1)=1.0;*nf=1;/*/void highpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf)double fp,fr,fs;printf(nPlease input the Fp,Ap,Fr,Ar,Fs value);printf(nFp,Ap:Passband frequency(Hz) and MAXAttenuation(dB);printf(nFr,Ar:Stopband frequency(Hz) and MINAttenuation(dB);printf(nFs is the sample frequency(Hz) highpass filter);printf(nInput parameters Fp,Ap,Fr,Ar,Fs:);scanf(%lf,%lf,%lf,%lf,%lf,&fp,ap,&fr,ar,&fs);if(fp*ar)|(fs2*fp)dosound(1000);delay(200);nosound();printf(Invalid input! Please Reinput:);scanf(%lf,%lf,%lf,%lf,%lf,&fp,ap,&fr,ar,&fs);while(fp*ar)|(fsfp2)|(fp1fr2)|(*ap*ar)|(fsfp2)|(fp1fr2)|(*ap*ar)|(fs2*fr2);wp1=2*PI*fp1/fs;wr1=2*PI*fr1/fs;wp2=2*PI*fp2/fs;wr2=2*PI*fr2/fs;cw0=sin(wp1+wp2)/(sin(wp1)+sin(wp2);pwp1=fabs(cw0-cos(wp1)/sin(wp1);pws1=fabs(cw0-cos(wr1)/sin(wr1);pwp2=fabs(cw0-cos(wp2)/sin(wp2);pws2=fabs(cw0-cos(wr2)/sin(wr2);if(fabs(pws1-pwp1) *wp=pws1;*ws=pws1;else*wp=pwp2;*ws=pws2;*f1=1.0;*(f1+1)=-2.0*cw0;*(f1+2)=1.0;*f2=-1.0;*(f2+1)=0.0*cw0;*(f2+2)=1.0;*nf=2;/*bcg-Chebyshev 和 Butterworth 型模拟原型传输函数生成子程序即程序得到系统函数 H(s).输出格式为：*/double *bcg(double ap,double as,double wp,double ws,int *n,double *h,int *type)int i,j,k;double a,c,e,p,q,x,y,wc,cs1,cs2;COMPLEX *b;double pp20;double xs88=1.0,1.0,1.41421356,1.0,2.0,2.0,1.0,2.61312593,3.41421356,2,61312593,1.0,3.23606789,5.23606789,5.23606789,3.23606789,1.0,3.86370331,7.46410162,9.14162017,7.46410162,3.86370331,1.0,4.49395921,10.09783468,14.59179389,14.59579389,10.09783468,4.49395921,1.0,5.12583090,13.13707118,21.84615097,25.68835593,21.84615097,13.13707118,5.12583090;printf(nTYPE 1.Butterworth 2.Chebyshev?(1/2):);scanf(%d,type);if(*type=2)c=sqrt(pow(10,as/10.0)-1.0);e=sqrt(pow(10,ap/10.0)-1.0);*n=(int)(fabs(FNCCH(c/e)/FNCCH(ws/wp)+0.99999);b=(COMPLEX*)calloc(*n+2,sizeof(COMPLEX);if(b=NULL)printf(nNot enough memory to allocate!);exit(0);wc=wp;a=pow(wc,(*n)/(e*pow(2.0,(*n)-1);q=1/e;x=FNSSH(q)/(*n);for(i=0;ireal=-wc*FNSH1(x)*sin(y);(b+i)-image=-wc*FNCH1(x)*cos(y);elsec=(pow(10.0,(0.1*ap)-1.0)/(pow(10.0,(0.1*as)-1.0);*n=(int)(fabs(log(c)/log(wp/ws)/2.0)+0.99999);b=(COMPLEX *)calloc(*n+2,sizeof(COMPLEX);if(b=NULL)printf(nNot enough memory to allocate!);exit(0);wc=wp/pow(pow(10.0,0.1*ap)-1.0,1.0/(2.0*(*n);a=pow(wc,(double)(*n);for(i=0;ireal=wc*cos(p);(b+i)-image=wc*sin(p);printf(nThe order of prototype filter is:%d,*n);/* b1=(COMPLEX*)calloc(*n+2,sizeof(COMPLEX);b2=(COMPLEX*)calloc(*n+2,sizeof(COMPLEX);h=(double *)calloc(*n+2),sizeof(double);if(h=NULL)printf(nNot enough memory to allocate!);exit(0);b1-real=-(b-real);b1-image=-(b-image);(b1+1)-real=1.0;(b1+1)-image=0.0;if(*n!=1)for(i=1;ireal-(b1+k+1)-real*(b+i)-real;cs2=(b1+k)-image-(b1+k+1)-real*(b+i)-image;(b2+k+1)-real=cs1+(b1+k+1)-image*(b+i)-image;(b2+k+1)-image=cs2-(b1+k+1)-image*(b+i)-real;b2-real=-(b1-real*(b+i)-real-b1-image*(b+i)-image);b2-image=-(b1-real*(b+i)-image+b1-image*(b+i)-real);(b2+i+1)-real=(b1+i)-real);(b2+i+1)-image=(b1+i)-image);for(k=0;kreal=(b2+k)-real;(b1+k)-image=(b2+k)-image;(b2+k)-real=0.0;(b2+k)-image=0.0;for(i=0;ireal;for(i=0;i=*n;i+)printf(nz-=f,i,hi);printf(nz0=f,nz1=f,nz2=f,nz3=f,nz4=f,1.0,2.6131/wc,3.4142/pow(wc,2),2.6131/pow(wc,3),1/a);*/for(i=0;i=*n-1;i+)ppi=xs*n-1i;for(i=0;ia=pnpe(f2,nf,ni,ptr-a,no); /*calcute the numberator coefficients */pt