#include <math.h>
 #include <stdio.h>
 #define M 40
 #define pi 3.141592653589793
   float bessel(float x,float eps)
   { int n;
     float f,fac,fx;

     f=1.0;    n=1;
     fx=x/2;   fac=1.0;
     fx=fx*fx;
     while(fx>=eps)
       { f=f+fx;
         n=n+1;
         fac*=n;
         fx=pow(x/2,n)/fac;
         fx=fx*fx;
         }
     return f;
     }

   //////////////////////////////////////////////////////////////
   // w為3dB截止頻率;要求w<1,1對應采樣頻率的一般;濾波器階數n;  //
   // b返回濾波器n+1個抽頭系數;                                //
   // type整型變量標志位用于確定所選的窗函數                   //
   //   type: 1:  巴特利特(Bartlett)窗                         //
   //         2:  漢寧(Hanning)窗                              //
   //         3:  哈明(Hamming)窗                              //
   //         4:  布萊克曼(Blackman)窗                         //
   //         5:  凱塞(Kaiser)窗(默認beta=4,誤差精度eps=1e-20) //
   //         其他:矩形窗                                      //
   // 函數返回0為正常返回;返回-1為異常返回:原因w>=1 or w<=0    //
   //////////////////////////////////////////////////////////////
 int fir1(float w,int n,float b[],int type)
 { int i;
   float alpha,m;

   if(w>=1||w<=0)
    { printf("\n W must be in the range of (0,pi]!\n");
      return(-1);
      }
   alpha=n/2.0;
   for(i=0; i<=(n+1)/2; i++)
    { //m=pi*(i-alpha);
      //if(m==0) b[i]=w;
      //else     b[i]=sin(w*m)/m;
      m=pi*(i-alpha+1e-20);
      b[i]=sin(w*m)/m;
      }
   alpha=pi/alpha;
   switch(type)
    { case 1 : for(i=0; i<=(n+1)/2; i++)
                  if(i<=n/2) b[i]*=i*alpha/pi;
                  else       b[i]*=2-i*alpha/pi;
               break;
      case 2 : for(i=0; i<=(n+1)/2; i++)
                  b[i]*=0.5*(1-cos(i*alpha));
               break;
      case 3 : for(i=0; i<=(n+1)/2; i++)
                  b[i]*=0.54-0.46*cos(i*alpha);
               break;
      case 4 : for(i=0; i<=(n+1)/2; i++)
                  b[i]*=0.42-0.50*cos(i*alpha)
                            +0.08*cos(2*i*alpha);
               break;
      case 5 : double ni,bes,be;
               ni=2.0/n;  bes=bessel(4,1e-20);
               for(i=0; i<=(n+1)/2; i++)
                 { be=1-i*ni;
                   be=sqrt(1-be*be);
                   b[i]*=bessel(4*be,1e-20)/bes;
                   }
               break;
      default: {   };
      }
   for(i=(n+1)/2+1; i<=n; i++)
      b[i]=b[n-i];
   return(0);
   }

   void main( )
   { float b[M+1],w;
     int rtn,i;

     w=0.25;
     rtn=fir1(w,M,b,3);
     printf("\n\n");
     if(rtn==0)
      for(i=0; i<=M; i++)
       printf("%.4f\t",b[i]);
       
     scanf("%f",&w);
     }