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 struct{double real,image;}COMPLEX;struct rptr{double *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 title[80],tmp[20];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,h[i]);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;i<=nz;i++)printf("\na[-]=_f,b[-]=_f",i,ptr->a[i],i,ptr->b[i]);printf("\n\nPress 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;i<=nz;i++){hwdb1.real+=ptr->a[i]*cos(i*jw);hwdb2.real+=ptr->b[i]*cos(i*jw);hwdb1.image+=ptr->a[i]*sin(i*jw);hwdb2.image+=ptr->b[i]*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;hwdb[k]=10*log10(amp1/amp2);if(hwdb[k]<-200) hwdb[k]=-200;if(k_==9) printf("*");}strcpy(title,"Transfer Property");if(type==1) strcat(title,"(BW Type) N=");if(type==2) strcat(title,"(CB Type) N=");itoa(ns,tmp,10);strcat(title,tmp);strcpy(tmp,"PI(rad)");draw_image(hwdb,N,title,"The Attenuation(dB)","0",tmp,0);free(ptr->b);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((fp>fr)||(*ap>*ar)||(fs<2*fr)){do{sound(1000);delay(200);nosound();printf("Invalid input! Please Reinput:");scanf("%lf,%lf,%lf,%lf,%lf",&fp,ap,&fr,ar,&fs);}while((fp>fr)||(*ap>*ar)||(fs<2*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)||(fs<2*fp)){do{sound(1000);delay(200);nosound();printf("Invalid input! Please Reinput:");scanf("%lf,%lf,%lf,%lf,%lf",&fp,ap,&fr,ar,&fs);}while((fp*ar)||(fs<2*fr));}*wp=fabs(1/tan(PI*fp/fs));*ws=fabs(1/tan(PI*fr/fs));*f1=1.0;*(f1+1)=1.0;*f2=-1.0;*(f2+1)=1.0;*nf=1;}/****************************************************************************/
void bandpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf){double fp1,fp2,fr1,fr2,fs,wp1,wp2,wr1,wr2,cw0,pwp1,pwp2,pws1,pws2;printf("\nPlease input the Fp1,Fp2,Ap,Fr1,Fr2,,Ar,Fs value");printf("\nFp1,Fp2,Ap:Passband frequency(Hz) and MAXAttenuation(dB)");printf("\nFr1,Fr2,Ar:Stopband frequency(Hz) and MINAttenuation(dB)");printf("\nFs is the sample frequency(Hz) bandpass filter");printf("\nInput parameters Fp1,Fp2,Ap,Fr1,Fr2,Ar,Fs:");scanf("%lf,%lf,%lf,%lf,%lf,%lf,%lf",&fp1,&fp2,ap,&fr1,&fr2,ar,&fs);if((fp1>fp2)||(fp1fr2)||(*ap>*ar)||(fs<2*fr2)){do{sound(1000);delay(200);nosound();printf("Invalid input! Please Reinput:");scanf("%lf,%lf,%lf,%lf,%lf,%lf,%lf",&fp1,&fp2,ap,&fr1,&fr2,ar,&fs);}while((fp1>fp2)||(fp1fr2)||(*ap>*ar)||(fs<2*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 pp[20];double xs[8][8]={{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;i<*n;i++){y=(2.0*i+1.0)*PI/(2.0*(*n));(b+i)->real=-wc*FNSH1(x)*sin(y);(b+i)->image=-wc*FNCH1(x)*cos(y);}}else{c=(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;i<*n;i++){p=PI*(0.5+(2.0*i+1.0)/2.0/(*n));(b+i)->real=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;i<*n;i++){for(k=0;k cs1=(b1+k)->real-(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;k<=i+1;k++){(b1+k)->real=(b2+k)->real;(b1+k)->image=(b2+k)->image;(b2+k)->real=0.0;(b2+k)->image=0.0;}}}for(i=0;i<=*n;i++)h[i]=(b1+i)->real;for(i=0;i<=*n;i++)printf("\nz[-]=_f",i,h[i]);printf("\nz[0]=_f,\nz[1]=_f,\nz[2]=_f,\nz[3]=_f,\nz[4]=_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++)pp[i]=xs[*n-1][i];for(i=0;i<=*n-1;i++)h[i]=pp[i]/pow(wc,i);free(b);/* free(b1);free(b2);*/return h;}/********************************************************************************
/********************************************************************************bsf-有理分式变换子程序**********************************************************************************/struct rptr * bsf(double *c,int ni,double *f1,double *f2,int nf,struct rptr *ptr,int *no){int i,j,ii,nd,ne,ng;double *d=NULL,*e=NULL,*g=NULL;ptr->a=pnpe(f2,nf,ni,ptr->a,no); /*calcute the numberator coefficients */ptr->b=(double *)calloc(*no+2,sizeof(double));if(ptr->b==NULL){printf("\nNot enough memory to allocate!");exit(0);}for(i=0;i<=ni;i++){ /* calculate the denominator coefficients */d=pnpe(f1,nf,i,d,&nd);e=pnpe(f2,nf,ni-i,e,&ne);g=ypmp(d,nd,e,ne,g,&ng);for(j=0;j<=ng;j++)ptr->b[j]+=c[i]*g[j];free(d);free(e);free(g);}return ptr;}/*********************************************************************************//*pnpe-*//**********************************************************************************/double *pnpe(double *a,int m,int n,double *b,int *mn){int i,j,k,nk;double *c;*mn=m*n;c=(double *)calloc(*mn+3,sizeof(double));b=(double *)calloc(*mn+3,sizeof(double));if(b==NULL){printf("\nNot enough memory to allocate!");exit(0);}if(n==0) {*b=1.00;free(c);return b;}else{for(i=0;i<=m;i++) b[i]=a[i];if(n==1) {free(c);return b;}else {nk=m;for(i=1;i for(j=0;j<=m;j++)for(k=0;k<=nk;k++) c[k+j]+=a[j]*b[k];nk+=m;for(k=0;k<=nk;k++) {b[k]=c[k];c[k]=0.0;}}}}free(c);return b;}/******************************************************************************
本文来源:https://www.2haoxitong.net/k/doc/dcec8263ddccda38376bafd6.html
文档为doc格式