IIR滤波器设计程序.docx
《IIR滤波器设计程序.docx》由会员分享，可在线阅读，更多相关《IIR滤波器设计程序.docx（12页珍藏版）》请在冰豆网上搜索。
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;kcs1=（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;
}
/********************************************************************************
/**************************************************************
升级会员 