#include <math.h>
 #include <stdio.h>
 #define M 67
 #define pi 3.141592653589793
   /////////////////////////////////////////////////////
   // w1,w2為通帶的下上3dB截頻;n為濾波器的階數        //
   // n返回濾波器的n+1個抽頭系數                      //
   // type整型變量標志位用于確定所選的窗函數          //
   //   type: 1:   巴特利特(Bartlett)窗               //
   //         2:   漢寧(Hanning)窗                    //
   //         3:   哈明(Hamming)窗                    //
   //         4:   布萊克曼(Blackman)窗               //
   //         其他:矩形窗                             //
   //  flag標志位用于確定濾波器的類型                 //
   //    flag: 0:帶通濾波器 1:帶阻濾波器              //
   // 函數返回0為正常返回－1為異常返回異常原因        //
   //    0<w1<w2<1 不滿足;1對應于采樣頻率的一般       //
   /////////////////////////////////////////////////////
 int fir2(float w1,float w2,int n,float b[],int type,int flag)
 { int i;
   float alpha,m;

   if(w1>=w2||w2>=1||w1<=0)
    { printf("\n  w1 and w2 must be as follow:");
      printf("\n     0<=w1<w2<1\n");
      return(-1);
      }
   alpha=n/2.0;
   if(flag==0)                  // BP FILTER
     for(i=0; i<=(n+1)/2; i++)
      { //m=pi*(i-alpha);
        //if(m==0) b[i]=w2-w1;
        //else b[i]=(sin(w2*m)-sin(w1*m))/m;
        m=pi*(i-alpha+1e-20);
        b[i]=(sin(w2*m)-sin(w1*m))/m;
        }
   else if(flag==1)            // BS FILTER
         { if(n%2 !=0)
            { printf("\n N must be even!\n"); n+=1; }
           for(i=0; i<=(n+1)/2; i++)
            { //m=pi*(i-alpha);
              //if(m==0) b[i]=w1+1-w2;
              //else b[i]=(sin(w1*m)+sin(m)-sin(w2*m))/m;
              m=pi*(i-alpha+1e-20);
              b[i]=(sin(w1*m)+sin(m)-sin(w2*m))/m;
              }
           }
        else
          { printf("\n ERROR!Flag must be 0 or 1!EXIT!\n");
            return(-2);
            }
   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;
       default: {   };
       }
   for(i=(n+1)/2+1; i<=n; i++)
     b[i]=b[n-i];
   return(0);
   }


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

     w1=0.25/pi;
     w2=0.5/pi;
     rtn=fir2(w1,w2,M,b,3,0);
     if(rtn==0)
       for(i=0; i<=M; i++)
         printf("%.4f\t",b[i]);
         
     scanf("%f",&w1);
     }