IIR滤波器设计程序.docx
《IIR滤波器设计程序.docx》由会员分享,可在线阅读,更多相关《IIR滤波器设计程序.docx(12页珍藏版)》请在冰豆网上搜索。
![IIR滤波器设计程序.docx](https://file1.bdocx.com/fileroot1/2023-1/21/e882a856-66c1-4a62-a26a-2827e19889ef/e882a856-66c1-4a62-a26a-2827e19889ef1.gif)
IIR滤波器设计程序
IIRD_DF.c-双线性变换法设计IIR数字滤波器
InfiniteImpulseResponseDigitalFilterDesignByDoubleConverting
====================================================
包括1.Lowpass 2.Highpass 3.bandpass三部分,如果有特殊需要
把其他部分删除掉即可。
***********************************************************/
/*#include*/
#include
#include
#include
#include
#include
#include
#include
/*COMPLEXSTRUCTURE*/
typedefstruct{
doublereal,image;
}COMPLEX;
structrptr{
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(doubleap,doubleas,doublewp,doublews,int*n,double*h,int*type);
structrptr*bsf(double*c,intni,double*f1,double*f2,intnf,structrptr*ptr,int*no);
double*pnpe(double*a,intm,intn,double*b,int*mn);
double*ypmp(double*a,intm,double*b,intn,double*c,int*mn);
voidlowpass_input(double*wp,double*ws,double*ap,double*ar,double*f1,double*f2,int*nf);
voidhighpass_input(double*wp,double*ws,double*ap,double*ar,double*f1,double*f2,int*nf);
voidbandpass_input(double*wp,double*ws,double*ap,double*ar,double*f1,double*f2,int*nf);
voiddraw_image(double*x,intm,char*title1,char*title2,char*xdis1,char*xdis2,intdis_type);
/***********************************************************************************************/
voidmain(void)
{
double*f1,*f2,*hwdb;
double*h=NULL;
intN,nf,ns,nz,i,j,k,ftype,type;
COMPLEXhwdb1,hwdb2;
doublewp,ws,ap,as,jw,amp1,amp2;
chartitle[80],tmp[20];
structrptr*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("\nNotenoughmemorytoallocate!
");
exit(0);
}
printf("\n1.Lowpass2.Highpass3.bandpass");
printf("\nPleaseselectthefiltertype:
");
scanf("%d",&ftype);
switch(ftype){
case1:
lowpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);
break;
case2:
highpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);
break;
case3:
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("\nTheanalogfilterdenominatorcoefficientsofHa(s):
");
for(i=0;i<=ns;i++)
printf("\nb[%2d]=%16f",i,h[i]);
ptr=bsf(h,ns,f1,f2,nf,ptr,&nz);
printf("\nThedigitalfiltercoefficientsofH(z):
");
printf("\n(aisnumeratorcoefficient.bisdenominatorcoefficient.)");
for(i=0;i<=nz;i++)
printf("\na[%2d]=%10f,b[%2d]=%16f",i,ptr->a[i],i,ptr->b[i]);
printf("\n\nPressanykeytocalculatethefilterresponseofH(z)...");
getch();
printf("\nWaittingforcalculating...");
/*Calculatethemagnitude-frequencyresponse*/
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%10==9)printf("*");
}
strcpy(title,"TransferProperty");
if(type==1)strcat(title,"(BWType)N=");
if(type==2)strcat(title,"(CBType)N=");
itoa(ns,tmp,10);
strcat(title,tmp);
strcpy(tmp,"PI(rad)");
draw_image(hwdb,N,title,"TheAttenuation(dB)","0",tmp,0);
free(ptr->b);
free(ptr->a);
free(h);
free(hwdb);
free(f2);
free(f1);
}
/***************************************************************/
voidlowpass_input(double*wp,double*ws,double*ap,double*ar,double*f1,double*f2,int*nf)
{
doublefp,fr,fs;
printf("\nPleaseinputtheFp,Ap,Fr,Ar,Fsvalue");
printf("\nFp,Ap:
Passbandfrequency(Hz)andMAXAttenuation(dB)");
printf("\nFr,Ar:
Stopbandfrequency(Hz)andMINAttenuation(dB)");
printf("\nFsisthesamplefrequency(Hz)Lowpassfilter");
printf("\nInputparametersFp,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("Invalidinput!
PleaseReinput:
");
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;
}
/**********************************************************************************************/
voidhighpass_input(double*wp,double*ws,double*ap,double*ar,double*f1,double*f2,int*nf)
{
doublefp,fr,fs;
printf("\nPleaseinputtheFp,Ap,Fr,Ar,Fsvalue");
printf("\nFp,Ap:
Passbandfrequency(Hz)andMAXAttenuation(dB)");
printf("\nFr,Ar:
Stopbandfrequency(Hz)andMINAttenuation(dB)");
printf("\nFsisthesamplefrequency(Hz)highpassfilter");
printf("\nInputparametersFp,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("Invalidinput!
PleaseReinput:
");
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;
}
/****************************************************************************/
voidbandpass_input(double*wp,double*ws,double*ap,double*ar,double*f1,double*f2,int*nf)
{
doublefp1,fp2,fr1,fr2,fs,wp1,wp2,wr1,wr2,cw0,pwp1,pwp2,pws1,pws2;
printf("\nPleaseinputtheFp1,Fp2,Ap,Fr1,Fr2,,Ar,Fsvalue");
printf("\nFp1,Fp2,Ap:
Passbandfrequency(Hz)andMAXAttenuation(dB)");
printf("\nFr1,Fr2,Ar:
Stopbandfrequency(Hz)andMINAttenuation(dB)");
printf("\nFsisthesamplefrequency(Hz)bandpassfilter");
printf("\nInputparametersFp1,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("Invalidinput!
PleaseReinput:
");
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(doubleap,doubleas,doublewp,doublews,int*n,double*h,int*type)
{
inti,j,k;
doublea,c,e,p,q,x,y,wc,cs1,cs2;
COMPLEX*b;
doublepp[20];
doublexs[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("\nTYPE1.Butterworth2.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("\nNotenoughmemorytoallocate!
");
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("\nNotenoughmemorytoallocate!
");
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("\nTheorderofprototypefilteris:
%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("\nNotenoughmemorytoallocate!
");
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[%2d]=%16f",i,h[i]);
printf("\nz[0]=%16f,\nz[1]=%16f,\nz[2]=%16f,\nz[3]=%16f,\nz[4]=%16f",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);*/
returnh;
}
/********************************************************************************
/**************************************************************