巴特沃斯滤波器c语言

1、1. 模拟滤波器的设计      1.1巴特沃斯滤波器的次数        根据给定的参数设计模拟滤波器，然后进行变数变换，求取数字滤波器的方法，称为滤波器的间接设计。做为数字滤波器的设计基础的模拟滤波器，称之为原型滤波器。这里，我们首先介绍的是最简单最基础的原型滤波器，巴特沃斯低通滤波器。由于IIR滤波器不具有线性相位特性，因此不必考虑相位特性，直接考虑其振幅特性。       在这里，N是滤波器的次数，c是截止频率。从上式的振幅特性可以看出，这个是单调递减的函数，其振幅特性是不存在

2、纹波的。设计的时候，一般需要先计算跟所需要设计参数相符合的次数N。首先，就需要先由阻带频率，计算出阻带衰减将巴特沃斯低通滤波器的振幅特性，直接带入上式，则有最后，可以解得次数N为当然，这里的N只能为正数，因此，若结果为小数，则舍弃小数，向上取整。      1.2巴特沃斯滤波器的传递函数         巴特沃斯低通滤波器的传递函数，可由其振幅特性的分母多项式求得。其分母多项式根据S解开，可以得到极点。这里，为了方便处理，我们分为两种情况去解这个方程。当N为偶数的时候，这里，使用了欧拉公式。同样的，当N为奇数的时候

3、，同样的，这里也使用了欧拉公式。归纳以上，极点的解为上式所求得的极点，是在s平面内，在半径为c的圆上等间距的点，其数量为2N个。为了使得其IIR滤波器稳定，那么，只能选取极点在S平面左半平面的点。选定了稳定的极点之后，其模拟滤波器的传递函数就可由下式求得。       1.3巴特沃斯滤波器的实现（C语言）          首先，是次数的计算。次数的计算，我们可以由下式求得。         其对应的C语言程序为cpp view plaincop

4、y1. N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /   2.             log10 (Stopband/Cotoff) );           然后是极点的选择

5、，这里由于涉及到复数的操作，我们就声明一个复数结构体就可以了。最重要的是，极点的计算含有自然指数函数，这点对于计算机来讲，不是太方便，所以，我们将其替换为三角函数，这样的话，实部与虚部就还可以分开来计算。其代码实现为cpp view plaincopy1. typedef struct   2.   3.     double Real_part;  4.     double Imag_Part;

6、0; 5.  COMPLEX;  6.   7.   8. COMPLEX polesN;  9.   10. for(k = 0;k <= (2*N)-1)  k+)  11.   12.     if(Cotoff*cos(k+dk)*(pi/N) < 0)  13. &#

7、160;     14.         polescount.Real_part = -Cotoff*cos(k+dk)*(pi/N);  15.     polescount.Imag_Part= -Cotoff*sin(k+dk)*(pi/N);        16.   

8、      count+;  17.         if (count = N) break;  18.       19.           计算出稳定的极点之后，就可以进行传递函数的计算了。传递的函数的计算，就像下式一样这里，为了得到模拟滤波

9、器的系数，需要将分母乘开。很显然，这里的极点不一定是整数，或者来说，这里的乘开需要做复数运算。其复数的乘法代码如下，cpp view plaincopy1. int Complex_Multiple(COMPLEX a,COMPLEX b,  2.                  double *Res_Real,double *Res_Imag

10、)  3.       4.   5.        *(Res_Real) =  (a.Real_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);  6.        *(Res_Imag)=  (a.Imag_P

11、art)*(b.Real_part) + (a.Real_part)*(b.Imag_Part);      7.      return (int)1;   8.   有了乘法代码之后，我们现在简单的情况下，看看其如何计算其滤波器系数。我们做如下假设这个时候，其传递函数为将其乘开，其大致的关系就像下图所示一样。计算的关系一目了然，这样的话，实现就简单多了。高阶的情况下也一样，重复这种计算就可以了。其代码为

12、cpp view plaincopy1.  Res0.Real_part = poles0.Real_part;   2.  Res0.Imag_Part= poles0.Imag_Part;  3.  Res1.Real_part = 1;   4.  Res1.Imag_Part= 0;  5.   6. for(count_1 = 0;cou

13、nt_1 < N-1;count_1+)  7.   8.   for(count = 0;count <= count_1 + 2;count+)  9.     10.       if(0 = count)  11.    12.  

14、0;            Complex_Multiple(Rescount, polescount_1+1,  13.                    &(Res_Savecount.Real_part),  14. &#

15、160;                  &(Res_Savecount.Imag_Part);  15.         16.       else if(count_1 + 2) = count

16、)  17.         18.             Res_Savecount.Real_part  += Rescount - 1.Real_part;  19.    Res_Savecount.Imag_Part += Rescount

17、0;- 1.Imag_Part;  20.                21.   else   22.     23.               Complex_Multiple(Re

18、scount, polescount_1+1,  24.                    &(Res_Savecount.Real_part),  25.                 

19、;   &(Res_Savecount.Imag_Part);                 26. 1   Res_Savecount.Real_part  += Rescount - 1.Real_part;  27.    Res_Savecount.I

20、mag_Part += Rescount - 1.Imag_Part;  28.    29.     30.  *(b+N) = *(a+N);  到此，我们就可以得到一个模拟滤波器巴特沃斯低通滤波器了。2.双1次z变换      2.1双1次z变换的原理        我们为了将模拟滤波器转换为数字滤波器的，可以用的方法很多。这里着重说说双1次

21、z变换。我们希望通过双1次z变换，建立一个s平面到z平面的映射关系，将模拟滤波器转换为数字滤波器。        和之前的例子一样，我们假设有如下模拟滤波器的传递函数。将其做拉普拉斯逆变换，可得到其时间域内的连续微分方程式，其中，x(t)表示输入，y(t)表示输出。然后我们需要将其离散化，假设其采样周期是T，用差分方程去近似的替代微分方程，可以得到下面结果然后使用z变换，再将其化简。可得到如下结果从而，我们可以得到了s平面到z平面的映射关系，即由于所有的高阶系统都可以视为一阶系统的并联，所以，这个映射关系在高阶系统中，也是成立的。然后，将关系式带入上式，

22、可得到这里，我们可以就可以得到与的对应关系了。         这里的与的对应关系很重要。我们最终的目的设计的是数字滤波器，所以，设计时候给的参数必定是数字滤波器的指标。而我们通过间接设计设计IIR滤波器时候，首先是要设计模拟滤波器，再通过变换，得到数字滤波器。那么，我们首先需要做的，就是将数字滤波器的指标，转换为模拟滤波器的指标，基于这个指标去设计模拟滤波器。另外，这里的采样时间T的取值很随意，为了方便计算，一般取1s就可以。       2.2双1次z变换的实现（C语言）    &

23、#160;    我们设计好的巴特沃斯低通滤波器的传递函数如下所示。     我们将其进行双1次z变换，我们可以得到如下式子可以看出，我们还是需要将式子乘开，进行合并同类项，这个跟之前说的算法相差不大。其代码为。cpp view plaincopy1. for(Count = 0;Count<=N;Count+)  2.              3. 

24、0;         for(Count_Z = 0;Count_Z <= N;Count_Z+)  4.               5.              

25、60;   ResCount_Z = 0;  6.              Res_SaveCount_Z = 0;    7.               8.    

26、;             Res_Save 0 = 1;  9.           for(Count_1 = 0; Count_1 < N-Count;Count_1+)  10.     

27、0;         11.             for(Count_2 = 0; Count_2 <= Count_1+1;Count_2+)  12.              

28、;     13.                     if(Count_2 = 0)  ResCount_2 += Res_SaveCount_2;      14.     

29、0;   else if(Count_2 = (Count_1+1)&&(Count_1 != 0)    15.                            ResCount_2&#

30、160;+= -Res_SaveCount_2 - 1;   16.         else  ResCount_2 += Res_SaveCount_2 - Res_SaveCount_2 - 1;  17.            for(Cou

31、nt_Z = 0;Count_Z<= N;Count_Z+)  18.                 19.                      Res_SaveCo

32、unt_Z  =  ResCount_Z   20.                    ResCount_Z  = 0;  21.              

33、;             22.               23.         for(Count_1 = (N-Count); Count_1 < N;Count_1+) &

34、#160;24.               25.                         for(Count_2 = 0; Count_2 <= Cou

35、nt_1+1;Count_2+)  26.                   27.                      if(Count_2 = 0) 

36、;ResCount_2 += Res_SaveCount_2;     28.                  else if(Count_2 = (Count_1+1)&&(Count_1 != 0)    29.    

37、;                    ResCount_2 += Res_SaveCount_2 - 1;  30.                  else

38、60;   31.                        ResCount_2 += Res_SaveCount_2 + Res_SaveCount_2 - 1;     32.    

39、60;          33.                   for(Count_Z = 0;Count_Z<= N;Count_Z+)  34.          

40、           35.                        Res_SaveCount_Z  =  ResCount_Z   36.    &#

41、160;               ResCount_Z  = 0;  37.                     38.       

42、;        39.             for(Count_Z = 0;Count_Z<= N;Count_Z+)  40.           41.       

43、0;             *(az+Count_Z) += pow(2,N-Count) * (*(as+Count) *  42.        Res_SaveCount_Z;  43.          

44、60;      *(bz+Count_Z) +=  (*(bs+Count) * Res_SaveCount_Z;               44.              45.   

45、;    到此，我们就已经实现了一个数字滤波器。3.IIR滤波器的间接设计代码（C语言）cpp view plaincopy1. #include <stdio.h>  2. #include <math.h>  3. #include <malloc.h>  4. #include <string.h>  5.   6.   7. #de

46、fine     pi     (double)3.1415926)  8.   9.   10. struct DESIGN_SPECIFICATION  11.   12.     double Cotoff;     13.     doubl

47、e Stopband;  14.     double Stopband_attenuation;  15. ;  16.   17. typedef struct   18.   19.     double Real_part;  20.     double Imag_Pa

48、rt;  21.  COMPLEX;  22.   23.   24.   25. int Ceil(double input)  26.   27.      if(input != (int)input) return (int)input) +1;  28.    

49、0; else return (int)input);   29.   30.   31.   32. int Complex_Multiple(COMPLEX a,COMPLEX b  33.                    

50、60;                 ,double *Res_Real,double *Res_Imag)  34.       35.   36.        *(Res_Real) =  (a.Rea

51、l_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);  37.        *(Res_Imag)=  (a.Imag_Part)*(b.Real_part) + (a.Real_part)*(b.Imag_Part);      38.      return (int)1

52、;   39.   40.   41.   42. int Buttord(double Cotoff,  43.                  double Stopband,  44.       &#

53、160;          double Stopband_attenuation)  45.   46.    int N;  47.   48.    printf("Wc =  %lf  rad/sec n" ,Cotoff); &#

54、160;49.    printf("Ws =  %lf  rad/sec n" ,Stopband);  50.    printf("As  =  %lf  dB n" ,Stopband_attenuation);  51.    printf("-n"

55、 );  52.        53.    N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /   54.              

56、;       log10 (Stopband/Cotoff) );  55.      56.      57.    return (int)N;  58.   59.   60.   61. int Butter(int N, dou

57、ble Cotoff,  62.                double *a,  63.                double *b)  64.   65.  

58、0;  double dk = 0;  66.     int k = 0;  67.     int count = 0,count_1 = 0;  68.     COMPLEX polesN;  69.     C

59、OMPLEX ResN+1,Res_SaveN+1;  70.   71.     if(N%2) = 0) dk = 0.5;  72.     else dk = 0;  73.   74.     for(k = 0;k <= (2*N)

60、-1)  k+)  75.       76.          if(Cotoff*cos(k+dk)*(pi/N) < 0)  77.            78.        

61、        polescount.Real_part = -Cotoff*cos(k+dk)*(pi/N);  79.           polescount.Imag_Part= -Cotoff*sin(k+dk)*(pi/N);        80.   &#

62、160;           count+;  81.             if (count = N) break;  82.            83.  &#

63、160;     84.   85.      printf("Pk =   n" );     86.      for(count = 0;count < N count+)  87.    

64、0;   88.            printf("(%lf) + (%lf i) n" ,-polescount.Real_part  89.                   &#

65、160;                   ,-polescount.Imag_Part);  90.        91.      printf("-n" );  92.    

66、60;   93.      Res0.Real_part = poles0.Real_part;   94.      Res0.Imag_Part= poles0.Imag_Part;  95.   96.      Res1.Real_part = 1;   97

67、.      Res1.Imag_Part= 0;  98.   99.     for(count_1 = 0;count_1 < N-1;count_1+)  100.       101.          for(count 

68、= 0;count <= count_1 + 2;count+)  102.            103.               if(0 = count)  104.     

69、0;        105.                     Complex_Multiple(Rescount, polescount_1+1,  106.           &

70、#160;                        &(Res_Savecount.Real_part),  107.                  &

71、#160;                 &(Res_Savecount.Imag_Part);  108.                /printf( "Res_Save : (%lf) +&#

72、160;(%lf i) n" ,Res_Save0.Real_part,Res_Save0.Imag_Part);  109.                 110.   111.               else

73、 if(count_1 + 2) = count)  112.                 113.                      Res_Sa

74、vecount.Real_part  += Rescount - 1.Real_part;  114.                 Res_Savecount.Imag_Part += Rescount - 1.Imag_Part;    115.    

75、                    116.             else   117.             

76、0; 118.                      Complex_Multiple(Rescount, polescount_1+1,  119.                 &

77、#160;                  &(Res_Savecount.Real_part),  120.                        &

78、#160;           &(Res_Savecount.Imag_Part);  121.   122.                        /printf( "Res

79、60;         : (%lf) + (%lf i) n" ,Rescount - 1.Real_part,Rescount - 1.Imag_Part);  123.                 /print

80、f( "Res_Save : (%lf) + (%lf i) n" ,Res_Savecount.Real_part,Res_Savecount.Imag_Part);  124.                   125.       

81、          Res_Savecount.Real_part  += Rescount - 1.Real_part;  126.                 Res_Savecount.Imag_Part += Rescount 

82、- 1.Imag_Part;  127.               128.                 /printf( "Res_Save : (%lf) + (%lf i) n&

83、quot; ,Res_Savecount.Real_part,Res_Savecount.Imag_Part);  129.                   130.               131.    &#

84、160;        /printf("There n" );  132.            133.   134.          for(count = 0;count <= N;c

85、ount+)  135.            136.                Rescount.Real_part = Res_Savecount.Real_part;   137.       &#

86、160;          Rescount.Imag_Part= Res_Savecount.Imag_Part;  138.                    139.         

87、60;   *(a + N - count) = Rescount.Real_part;  140.            141.            142.         

88、0;/printf("There! n" );  143.            144.           145.   146.      *(b+N) = *(a+N);  147.   14

89、8.      /-display-/  149.      printf("bs =  " );     150.      for(count = 0;count <= N count+)  151.    &

90、#160;   152.            printf("%lf ", *(b+count);  153.        154.      printf("  n" );  155.  

91、60;156.      printf("as =  " );     157.      for(count = 0;count <= N count+)  158.        159.     &

92、#160;      printf("%lf ", *(a+count);  160.        161.      printf("  n" );  162.   163.      printf("-n

93、" );  164.   165.      return (int) 1;  166.   167.   168.   169. int Bilinear(int N,   170.              

94、;    double *as,double *bs,  171.                  double *az,double *bz)  172.   173.       int Count =&

95、#160;0,Count_1 = 0,Count_2 = 0,Count_Z = 0;  174.       double ResN+1;  175.     double Res_SaveN+1;   176.   177.         fo

96、r(Count_Z = 0;Count_Z <= N;Count_Z+)  178.       179.                  *(az+Count_Z)  = 0;  180.      

97、       *(bz+Count_Z)  = 0;  181.       182.   183.       184.     for(Count = 0;Count<=N;Count+)  185.     

98、         186.           for(Count_Z = 0;Count_Z <= N;Count_Z+)  187.               188.   &#

99、160;              ResCount_Z = 0;  189.              Res_SaveCount_Z = 0;    190.      

100、60;        191.              Res_Save 0 = 1;  192.       193.           for(Count_1 =

101、 0; Count_1 < N-Count;Count_1+)  194.               195.             for(Count_2 = 0; Count_2 <= Count_1+1;Co

102、unt_2+)  196.                   197.                      if(Count_2 = 0)  

103、  198.                    199.                        ResCount_2 += Re

104、s_SaveCount_2;  200.                      /printf( "Res%d %lf  n" , Count_2 ,ResCount_2);  201.      

105、0;               202.   203.                  else if(Count_2 = (Count_1+1)&&(Count_1 != 0) 

106、   204.                    205.                        ResCount_2 +=&#

107、160;-Res_SaveCount_2 - 1;     206.                               /printf( "Res%d %lf  n&qu

108、ot; , Count_2 ,ResCount_2);  207.                     208.   209.                 

109、60;else    210.                    211.                        ResCoun

110、t_2 += Res_SaveCount_2 - Res_SaveCount_2 - 1;  212.                     /printf( "Res%d %lf  n" , Count_2 ,ResCount

111、_2);  213.                                   214.            

112、0;  215.   216.                     /printf( "Res : ");  217.               

113、;for(Count_Z = 0;Count_Z<= N;Count_Z+)  218.                 219.                      

114、Res_SaveCount_Z  =  ResCount_Z   220.                    ResCount_Z  = 0;  221.            

115、;     /printf( "%d  %lf  " ,Count_Z, Res_SaveCount_Z);       222.                 223.     

116、0;         /printf(" n" );  224.               225.               226.   22

117、7.         for(Count_1 = (N-Count); Count_1 < N;Count_1+)  228.               229.            

118、60;        for(Count_2 = 0; Count_2 <= Count_1+1;Count_2+)  230.                   231.        &#

119、160;             if(Count_2 = 0)    232.                    233.        

120、60;               ResCount_2 += Res_SaveCount_2;  234.                      /printf( "Res%

121、d %lf  n" , Count_2 ,ResCount_2);  235.                      236.   237.            &

122、#160;     else if(Count_2 = (Count_1+1)&&(Count_1 != 0)    238.                    239.       

123、60;                ResCount_2 += Res_SaveCount_2 - 1;  240.                      

124、0;        /printf( "Res%d %lf  n" , Count_2 ,ResCount_2);  241.                      242.  

125、0;243.                  else    244.                    245.       

126、60;                ResCount_2 += Res_SaveCount_2 + Res_SaveCount_2 - 1;  246.                  

127、    /printf( "Res%d %lf  n" , Count_2 ,ResCount_2);  247.                             

128、;      248.               249.   250.                    /   printf( "

129、Res : ");  251.               for(Count_Z = 0;Count_Z<= N;Count_Z+)  252.                 253. &#

130、160;                    Res_SaveCount_Z  =  ResCount_Z   254.                  &

131、#160; ResCount_Z  = 0;  255.                 /printf( "%d  %lf  " ,Count_Z, Res_SaveCount_Z);       256.  

132、;               257.                /printf(" n" );  258.           &#

133、160;   259.   260.   261.              /printf( "Res : ");  262.         for(Count_Z = 0;Count_Z<= N;Count_

134、Z+)  263.           264.                     *(az+Count_Z) +=  pow(2,N-Count)  *  (*(as+Count) *

135、0;Res_SaveCount_Z;  265.              *(bz+Count_Z) +=  (*(bs+Count) * Res_SaveCount_Z;         266.          

136、            /printf( "  %lf  " ,*(bz+Count_Z);           267.              268.

137、         /printf(" n" );  269.   270.       271.   272.   273.         274.      for(Count_Z =&

138、#160;N;Count_Z >= 0;Count_Z-)  275.        276.           *(bz+Count_Z) =  (*(bz+Count_Z)/(*(az+0);  277.           *(az+Count_Z) =  (*(az+Count_Z)/(*(az+0);  278.        279.