#include<xxgc.h>
/*****IIR Digital Filter Desgin Process***/
int orderBtw(int,double,double,double,double,double,...);
int orderBtw(bandType,db1,db2,fs,f1,f2,f3,f4)
int bandType;
double db1,db2,fs,f1,f2,f3,f4;
{
double pi,omga1,omga2,omga3,omga4,p1,p2,p;
int N;
if(bandType<1
||bandType>4
||db1<0
||db2<db1
||f1<0
||f2<f1
||fs<f2)
{printf("input argument error in orderBtw()\n");
getch();
exit(1);}
pi=M_PI;/*4*atan(1.0);*/
/*********for low_pass and high_pass***********/
omga1=tan(pi*f1/fs);
omga2=tan(pi*f2/fs);
if(bandType<3)p=omga2/omga1;
/*********for band_pass and band_stop***********/
else{
if(f3<f2
||f4<f3
||fs<f4)
{printf("Input argument error in orderBtw()\n");
getch();exit(1);}
omga3=tan(pi*f3/fs);
omga4=tan(pi*f4/fs);
if(bandType==3)
{
p1=(omga1*omga1-omga3*omga2)/(omga1*(omga3-omga2));
p2=(omga4*omga4-omga3*omga2)/(omga4*(omga3-omga2));
p1=-p1;
p=(p1<p2)?p1:p2;
}
else if(bandType==4)
{
p1=-omga2*(omga4-omga1)/(omga2*omga2-omga4*omga1);
p2=-omga3*(omga4-omga1)/(omga3*omga3-omga4*omga1);
p2=-p2;
p=(p1<p2)?p1:p2;}
}
p1=pow(10,db1/10)-1;
p2=pow(10,db2/10)-1;
p=0.5*log10(p2/p1)/log10(p);
N=ceil(p);
return N;
}
int getBtw(double b[][2],int order)
/*notice:please declare double d[N],N=filter_order */
{
int k,M,L;
double pi;
pi=4*atan(1.0);
M=order/2;L=M;
if(order!=M*2){b[M][0]=0;b[M][1]=1;L++;}
for(k=0;k<M;k++){
b[k][0]=1;b[k][1]=(-2)*cos((2*k+order+1)*pi/(2*order));}
return L;
}
double getc23(double [2][3],
int,int,double,double,double,double,...);
double getc23(c,bandType,N,db1,fs,f1,f2,f3,f4)
int bandType,N;
double c[2][3],db1,fs,f1,f2,f3,f4;
{
double fc,pi,*b,namda;
pi=4*atan(1.0);
if(bandType==LOWPASS)
{c[0][0]=1;c[0][1]=-c[0][0];c[0][2]=0.0;
c[1][0]=tan(pi*f1/fs);c[1][1]=c[1][0];c[1][2]=0.0;
fc=f1;}
else if(bandType==HIGHPASS)
{c[0][0]=tan(pi*f2/fs);c[0][1]=c[0][0];c[0][2]=0.0;
c[1][0]=1;c[1][1]=-1;c[1][2]=0.0;
fc=f2;}
else if(bandType==BANDPASS)
{double U,L;
U=tan(pi*f3/fs);L=tan(pi*f2/fs);
c[0][0]=1+U*L;c[0][1]=2*(U*L-1);c[0][2]=1+U*L;
c[1][0]=U-L;c[1][1]=0;c[1][2]=L-U;
fc=f2;}
else if(bandType==BANDSTOP)
{double U,L;
U=tan(pi*f4/fs);L=tan(pi*f1/fs);
c[1][0]=1+U*L;c[1][1]=2*(U*L-1);c[1][2]=1+U*L;
c[0][0]=U-L;c[0][1]=0;c[0][2]=L-U;
fc=f1;}
namda=pow(10,db1/10)-1;
namda=pow(1/namda,0.5/N);
c[1][0]=namda*c[1][0];
c[1][1]=namda*c[1][1];
c[1][2]=namda*c[1][2];
return fc;
}
void BtwAFtoDF(double H[][2][5],int L,
double b[][2],double c[2][3]
)
/*jiang s=(a+b/z+c/zz)/(d+e/z+f/zz) dairu
h(s)=1/(m*ss+n*s+1) de xinshi zhong
dedao:
h(1/z)=[dd+2de/z+(ee+2df)/zz+2ef/zzz+ff/zzzz]
/{m[aa+2ab/z+(bb+2ac)/zz+2bc/zzz+cc/zzzz]
+[dd+2de/z+(ee+2df)/zz+2ef/zzz+ff/zzzz]
+n[ad+(ae+bd)/z+(cd+af+be)/zz+(bf+ce)/zzz+cf/zzzz]}
=[]/
[maa+dd+nad
+(2mab+2de+nae+nbd)/z
+(mbb+2mac+2df+ee+n(cd+af+be))/zz
+(2mbc+2ef+nbf+nce)/zzz
+(mcc+ff+ncf)/zzzz]
****/
{
int i;
for(i=0;i<L;i++)
{H[i][0][0]=c[1][0]*c[1][0];
H[i][0][1]=2*c[1][0]*c[1][1];
H[i][0][2]=c[1][1]*c[1][1]+2*c[1][0]*c[1][2];
H[i][0][3]=2*c[1][1]*c[1][2];
H[i][0][4]=c[1][2]*c[1][2];
H[i][1][0]=b[i][0]*c[0][0]*c[0][0]+c[1][0]*c[1][0]+b[i][1]*c[0][0]*c[1][0];
H[i][1][1]=2*b[i][0]*c[0][0]*c[0][1]+2*c[1][0]*c[1][1]+b[i][1]*(c[0][0]*c[1][1]+c[0][1]*c[1][0]);
H[i][1][2]=b[i][0]*(c[0][1]*c[0][1]+2*c[0][0]*c[0][2])+2*c[1][0]*c[1][2]
+c[1][1]*c[1][1]+b[i][1]*(c[0][2]*c[1][0]+c[0][0]*c[1][2]+c[0][1]*c[1][1]);
H[i][1][3]=2*b[i][0]*c[0][1]*c[0][2]+2*c[1][1]*c[1][2]
+b[i][1]*(c[0][1]*c[1][2]+c[0][2]*c[1][1]);
H[i][1][4]=b[i][0]*c[0][2]*c[0][2]+c[1][2]*c[1][2]+b[i][1]*c[0][2]*c[1][2];
}}
double pzValue(double *a,int N,double x)
{
int i;
double psumI,psumR;
psumR=a[0];psumI=0.0;
for(i=1;i<N;i++){psumR=psumR+a[i]*cos(x*i);
psumI=psumI-a[i]*sin(x*i);}
return(sqrt(psumR*psumR+psumI*psumI));
}
double Btw20lgHz(double f,double *fs,double H[][2][5],int *L)
{
double wT,sum;
int i;
wT=8*f*atan(1.0)/fs[0];
for(sum=0,i=0;i<L[0];i++)
{sum=sum+20*log10(pzValue(&(H[i][0][0]),5,wT)/pzValue(&(H[i][1][0]),5,wT));}
return sum;}
/**********************(1)*************************************/
/***To do here,according to book p74(20-1),
draw Analog Buttworth Filter's figure ***/
double Btw20lgHs(double f,double *C,int *N)
{/***program your function here****/
return -10*log10(1+pow(*C*f**C*f,*N));
}
/****************************************************************/
void IIR(int,double,double,double,double,double,...);
void IIR(bandType,db1,db2,fs,f1,f2,f3,f4)
int bandType;
double db1,db2,fs,f1,f2,f3,f4;
{int order,L,i,NK;
double fc,b[10][2],c[2][3],H[10][2][5],d[3];
order=orderBtw(bandType,db1,db2,fs,f1,f2,f3,f4);
L=getBtw(b,order);
/****fc=****/
getc23(c,bandType,order,db1,fs,f1,f2,f3,f4);
BtwAFtoDF(H,L,b,c);
/***************************(2)***********************************/
/*to do here,please use printf() display Digital Filter's coeaffence */
/*****************************************************************/
{int gd=VGA,gm=1;double C;
initgraph(&gd,&gm,"");
window2("lowpass",-1.,5.,1000.,-60.,"hz","db",BLUE,BLUE);
plotxy2(GREEN,Btw20lgHz,&fs,H,&L);
/***************************(3)************************************/
/***please modify here draw Analong Buttworth Filter **/
/*****C=?? *******/
C=pow(pow(10,0.1*db1)-1,0.5/order);
plotxy2(RED,Btw20lgHs,&C,&order);
/*********************************************************/
setcolor(LIGHTRED);
if(bandType==LOWPASS)
{rectangle2(0.0,0.0,f1,-db1);rectangle2(0.0,0.0,f2,-db2);}
else if(bandType==HIGHPASS)
{rectangle2(0.0,.0,f2,-db1);rectangle2(0.0,0.0,f1,-db2);}
else if(bandType==BANDSTOP)
{rectangle2(f1,0.0,f4,-db1);rectangle2(f2,0.0,f3,-db2);}
else if(bandType==BANDPASS)
{rectangle2(f2,0.0,f3,-db1);rectangle2(f1,0.0,f4,-db2);}
getch();
closegraph();
}
}
main()
{IIR(1,3.0,18.0,2000.,200.,300.,400.,500.);
IIR(2,3.0,18.0,2000.,200.,300.,400.,500.);
IIR(3,3.0,18.0,2000.,200.,300.,400.,500.);
IIR(4,3.0,18.0,2000.,200.,300.,400.,500.);
getch();
/************(4)***************************/
/*****to do here, Realize IIR filter ****/
/******************************************/
}
评论0