IIR数字滤波器设计C程序文档格式.docx
《IIR数字滤波器设计C程序文档格式.docx》由会员分享,可在线阅读,更多相关《IIR数字滤波器设计C程序文档格式.docx(11页珍藏版)》请在冰点文库上搜索。
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);
}
\n1.Lowpass2.Highpass3.bandpass"
\nPleaseselectthefiltertype:
scanf("
%d"
&
ftype);
switch(ftype){
case1:
lowpass_input(&
wp,&
ws,&
ap,&
as,f1,f2,&
nf);
break;
case2:
highpass_input(&
case3:
bandpass_input(&
default:
h=bcg(ap,as,wp,ws,&
ns,h,&
type);
\nTheanalogfilterdenominatorcoefficientsofHa(s):
for(i=0;
i<
=ns;
i++)
\nb[-]=_f"
i,h[i]);
ptr=bsf(h,ns,f1,f2,nf,ptr,&
nz);
\nThedigitalfiltercoefficientsofH(z):
\n(aisnumeratorcoefficient.bisdenominatorcoefficient.)"
=nz;
\na[-]=_f,b[-]=_f"
i,ptr->
a[i],i,ptr->
b[i]);
\n\nPressanykeytocalculatethefilterresponseofH(z)..."
\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;
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,"
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);
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;
\nPleaseinputtheFp,Ap,Fr,Ar,Fsvalue"
\nFp,Ap:
Passbandfrequency(Hz)andMAXAttenuation(dB)"
\nFr,Ar:
Stopbandfrequency(Hz)andMINAttenuation(dB)"
\nFsisthesamplefrequency(Hz)Lowpassfilter"
\nInputparametersFp,Ap,Fr,Ar,Fs:
%lf,%lf,%lf,%lf,%lf"
fp,ap,&
fr,ar,&
fs);
if((fp>
fr)||(*ap>
*ar)||(fs<
2*fr)){
do{
sound(1000);
delay(200);
nosound();
Invalidinput!
PleaseReinput:
}while((fp>
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)
\nFsisthesamplefrequency(Hz)highpassfilter"
if((fp*ar)||(fs<
2*fp)){
}while((fp*ar)||(fs<
*wp=fabs(1/tan(PI*fp/fs));
*ws=fabs(1/tan(PI*fr/fs));
*f1=1.0;
*f2=-1.0;
/****************************************************************************/
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;
\nPleaseinputtheFp1,Fp2,Ap,Fr1,Fr2,,Ar,Fsvalue"
\nFp1,Fp2,Ap:
\nFr1,Fr2,Ar:
\nFsisthesamplefrequency(Hz)bandpassfilter"
\nInputparametersFp1,Fp2,Ap,Fr1,Fr2,Ar,Fs:
%lf,%lf,%lf,%lf,%lf,%lf,%lf"
fp1,&
fp2,ap,&
fr1,&
fr2,ar,&
if((fp1>
fp2)||(fp1fr2)||(*ap>
2*fr2)){
}while((fp1>
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)=-2.0*cw0;
*(f1+2)=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}};
\nTYPE1.Butterworth2.Chebyshev?
(1/2):
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){
\nNotenoughmemorytoallocate!
wc=wp;
a=pow(wc,(*n))/(e*pow(2.0,(*n)-1));
q=1/e;
x=FNSSH(q)/(*n);
*n;
y=(2.0*i+1.0)*PI/(2.0*(*n));
(b+i)->
real=-wc*FNSH1(x)*sin(y);
image=-wc*FNCH1(x)*cos(y);
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));
wc=wp/pow(pow(10.0,0.1*ap)-1.0,1.0/(2.0*(*n)));
a=pow(wc,(double)(*n));
p=PI*(0.5+(2.0*i+1.0)/2.0/(*n));
real=wc*cos(p);
image=wc*sin(p);
\nTheorderofprototypefilteris:
*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){
b1->
real=-(b->
real);
image=-(b->
image);
(b1+1)->
real=1.0;
image=0.0;
if(*n!
=1){
for(i=1;
kcs1=(b1+k)->
real-(b1+k+1)->
real*(b+i)->
real;
cs2=(b1+k)->
image-(b1+k+1)->
image;
(b2+k+1)->
real=cs1+(b1+k+1)->
image*(b+i)->
image=cs2-(b1+k+1)->
b2->
real=-(b1->
real-b1->
image=-(b1->
image+b1->
(b2+i+1)->
real=((b1+i)->
image=((b1+i)->
=i+1;
(b1+k)->
real=(b2+k)->
image=(b2+k)->
(b2+k)->
real=0.0;
=*n;
h[i]=(b1+i)->
\nz[-]=_f"
\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);
=*n-1;
pp[i]=xs[*n-1][i];
h[i]=pp[i]/pow(wc,i);
free(b);
/*free(b1);
free(b2);
returnh;
/********************************************************************************
bsf-有理分式变换子程序
*********************