firlhbsp.cpp

数个关于滤波器的产生的C程序

 #include <math.h>
 #include <stdio.h>
 #define M 38
 #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) //
   //         其他:矩形窗                                      //
   //  flag整型变量标志位确定滤波器的类型                      //
   //    flag: 0:低通滤波器 1:高通滤波器                       //
   // 函数返回 0为正常返回;                                    //
   //     返回-1为异常返回:原因w>1                             //
   //     返回-2为异常返回:flag!=0 or 1                        //
   //////////////////////////////////////////////////////////////
 int fir1(float w,int n,float b[],int type,int flag)
 { int i;
   float alpha,m;

   if(w>=1||w<=0)
    { printf("\n W must be in the range of (0,1)!\n");
      return(-1);
      }
   alpha=n/2.0;
   if(flag==0)                                     //  LP FILTER
    for(i=0; i<=(n+1)/2; i++)
     { //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;
       }
   else if(flag==1)                                //  HP FILTER
         for(i=0; i<=(n+1)/2; i++)
           { //if(m==0) b[i]=1-w;
             //else b[i]=(sin(m)-sin(w*m))/m;
             m=pi*(i-alpha+1e-20);
             b[i]=(sin(m)-sin(w*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;
      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);
   }
   ///////////////////////////////////////////////////////////////
   // w1,w2为通带的下上3dB截频;n为滤波器的阶数                  //
   // n返回滤波器的n+1个抽头系数                                //
   // type整型变量标志位用于确定所选的窗函数                    //
   //   type: 1:   巴特利特(Bartlett)窗                         //
   //         2:   汉宁(Hanning)窗                              //
   //         3:   哈明(Hamming)窗                              //
   //         4:   布莱克曼(Blackman)窗                         //
   //         5:   凯塞(Kaiser)窗(默认beta=4,误差精度eps=1e-20) //
   //         其他:矩形窗                                       //
   //  flag标志位用于确定滤波器的类型                           //
   //    flag: 0:带通滤波器 1:带阻滤波器                        //
   // 函数返回0为正常返回－1为异常返回异常原因                  //
   //    0<w1<w2<1 不满足;1对应于采样频率的一般                 //
   ///////////////////////////////////////////////////////////////
 int fir1(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++)
     { //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++)
             { //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-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.12;
     rtn=fir1(w,M,b,3,0);
     printf("\n\n");
     if(rtn==0)
      for(i=0; i<=M; i++)
       { printf("%.4f\t",b[i]);
         if((i+1)%6==0) printf("\n");
         }
       
     scanf("%f",&w);
     }