win.c

FIR滤波器系数计算 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 <beta  <10) \n");
    scanf("%lf",&beta);

}
pi=4.0*atan(1.0);
 if((n%2)==0)
 	{
        n2=n/2-1;
        mid=1;
 	}
 else
 	{
       n2=n/2;
	 mid=0;
 	}
 delay=n/2.0;	
 wc1=2.0*pi*fln;
 if(band>=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<sum*(1.0e-8))
	break;
}
return (sum);
}