www.pudn.com > FIR.rar > win.c

 
 
#include "math.h" 
#include "stdio.h" 
 
void firwin(int n,int band,double fln,double fhn,double wn,double h[]); 
static double window(int type,int n,int i,double beta); 
static double kaiser(int i,int n,double beta); 
static double bessel0(double x ); 
 
 
main() 
{ 
int i,j,n,n2,band,wn; 
double fl,fh,fs,freq; 
double h[100]; 
char fname[40]; 
FILE *fp; 
 
printf("选择滤波器类型\n"); 
printf(" 1 --低通滤波器; 2 -- 高通滤波器\n"); 
printf(" 3 --带通滤波器; 4 -- 带阻滤波器\n"); 
do{ 
scanf("%d",&band); 
if(band>4) 
     { 
	printf("输入错误请重新输入\n"); 
	band=0; 
     } 
} 
while(band==0); 
 
printf("输入阶数 n\n"); 
scanf("%d",&n); 
fh=0.0; 
switch(band) 
{ 
case 1: printf("输入通带上边界频率 fl\n");scanf("%lf",&fl);break; 
case 2: printf("输入通带下边界频率 fl\n");scanf("%lf",&fl);break; 
case 3: printf("输入通带下边界频率 fl\n"); 
            scanf("%lf",&fl); 
            printf("输入上边界频率 fh\n"); 
            scanf("%lf",&fh);break; 
case 4: printf("输入通带下边界频率 fl\n"); 
           scanf("%lf",&fl); 
           printf("输入上边界频率 fh\n"); 
           scanf("%lf",&fh);break; 
default: ; 
} 
 
printf("输入取样频率 fs\n"); 
scanf("%lf",&fs); 
 
printf("选择窗函数类型 \n"); 
printf(" 1 -- retangular;  2 -- tapered rectangular\n"); 
printf(" 3 --triangular;    4 -- Hanning\n"); 
printf(" 5 -- Hamming ;   6 -- Blackman\n"); 
printf(" 7 -- Kaiser\n"); 
scanf("%d",&wn); 
 
printf("输入文件名 保存滤波系数\n"); 
scanf("%s",fname); 
if((fp=fopen(fname,"w"))==NULL) 
{ 
printf("不能打开文件\n"); 
exit(1); 
} 
fl=fl/fs; 
fh=fh/fs; 
firwin(n, band, fl, fh,  wn,  h ); 
printf("FIR 滤波系数表 \n\n\n"); 
 
 
n2=n/2; 
for(i=0;i<=n2;i++) 
{ 
j=n-i; 
printf(" h(%2d) = %12.8lf ; h(%2d)=%12.8lf \n",i,h[i],j,h[i]); 
fprintf(fp,"h(%2d) = %12.8lf ; h(%2d)=%12.8lf \n",i,h[i],j,h[i]); 
} 
fclose(fp); 
scanf("%d",&band); 
} 
 
void firwin(int n,int band,double fln,double fhn,double wn,double h[]) 
{ 
int i,n2,mid; 
double s,pi,wc1,wc2,beta,delay; 
 
 
beta=0; 
if(wn==7) 
{ 
printf("input beta parameter of Kaiser window ( 3 =3) 
 	wc2=2.0*pi*fhn; 
switch(band) 
{ 
case 1:  
	   for(i=0;i<=n2;i++) 
              { 
              s=i-delay; 
	        h[i]=(sin(wc1*s)/(pi*s))*window(wn,n+1,i,beta); 
		 h[n-i]=h[i];	 
              } 
            if(mid==1) 
			h[n/2]=wc1/pi; 
			 
          break; 
case 2: 
	   for(i=0;i<=n2;i++) 
	   	{ 
	   	s=i-delay; 
		h[i]=(sin(pi*s)-sin(wc1*s))/(pi*s); 
		h[i]=h[i]*window(wn,n+1,i,beta); 
		h[n-i]=h[i]; 
	   	} 
	   if(mid==1) 
	   	  h[n/2]=1.0-wc1/pi; 
	   break; 
case 3:  
	   for(i=0;i<=n2;i++) 
	   	{ 
	   	s=i-delay; 
		h[i]=(sin(wc2*s)-sin(wc1*s))/(pi*s); 
		h[i]=h[i]*window(wn,n+1,i,beta); 
		h[n-i]=h[i]; 
	   	} 
	   if(mid==1) 
	   	   h[n/2]=(wc2-wc1)/pi; 
	   break; 
case 4:  
	    for(i=0;i<=n2;i++) 
	    	{ 
	    	s=i-delay; 
		h[i]=(sin(wc1*s)+sin(pi*s)-sin(wc2*s))/(pi*s); 
		h[i]=h[i]*window(wn,n+1,i,beta); 
		h[n-i]=h[i]; 
	    	} 
		if(mid==1) 
			h[n/2]=(wc1+pi-wc2)/pi; 
		break; 
default: ;		 
} 
} 
 
static double window(int type,int n,int i,double beta) 
{ 
int k; 
double pi,w; 
 
pi=4.0*atan(1.0); 
w=1.0; 
switch(type) 
{ 
case 1: 
	w=1.0; 
	break; 
case 2: 
	 k=(n-2)/10; 
	 if(i<=k) 
	 	w=0.5*(1.0-cos(i*pi/(k+1))); 
        if(i>n-k-2) 
		w=0.5*(1.0-cos((n-i-1)*pi/(k+1)));	 
	break; 
case 3: 
	  w=1.0-fabs(1.0-2*i/(n-1.0)); 
	  break; 
case 4: 
	 w=0.5*(1.0-cos(2*i*pi/(n-1))); 
	 break; 
case 5: 
	  w=0.54-0.46*cos(2*i*pi/(n-1)); 
	  break; 
case 6: 
	w=0.42-0.5*cos(2*i*pi/(n-1))+0.08*cos(4*i*pi/(n-1)); 
	break; 
case 7: 
	 w=kaiser(i,n,beta); 
	  break; 
default: ;	   
} 
return (w); 
} 
 
  
static double kaiser(int i,int n,double beta) 
{ 
double a,w,a2,b1,b2,beta1; 
 
b1=bessel0(beta); 
a=2.0*i/(double)(n-1)-1.0; 
a2=a*a; 
beta1=beta*sqrt(1.0-a2); 
b2=bessel0(beta1); 
w=b2/b1; 
return (w); 
} 
 
static double bessel0(double x ) 
{ 
int i; 
double d,y,d2,sum; 
 
y=x/2.0; 
d=1.0; 
sum=1.0; 
for(i=1;i<=25;i++) 
{ 
d=d*y/i; 
d2=d*d; 
sum=sum+d2; 
if(d2