巴特沃斯濾波器c語言

1、1. 模擬濾波器的設(shè)計(jì)      1.1巴特沃斯濾波器的次數(shù)        根據(jù)給定的參數(shù)設(shè)計(jì)模擬濾波器，然后進(jìn)行變數(shù)變換，求取數(shù)字濾波器的方法，稱為濾波器的間接設(shè)計(jì)。做為數(shù)字濾波器的設(shè)計(jì)基礎(chǔ)的模擬濾波器，稱之為原型濾波器。這里，我們首先介紹的是最簡單最基礎(chǔ)的原型濾波器，巴特沃斯低通濾波器。由于IIR濾波器不具有線性相位特性，因此不必考慮相位特性，直接考慮其振幅特性。       在這里，N是濾波器的次數(shù)，c是截止頻率。從上式的振幅特性可以看出，這個是單調(diào)遞減的函數(shù)，其振幅特性是不存在

2、紋波的。設(shè)計(jì)的時候，一般需要先計(jì)算跟所需要設(shè)計(jì)參數(shù)相符合的次數(shù)N。首先，就需要先由阻帶頻率，計(jì)算出阻帶衰減將巴特沃斯低通濾波器的振幅特性，直接帶入上式，則有最后，可以解得次數(shù)N為當(dāng)然，這里的N只能為正數(shù)，因此，若結(jié)果為小數(shù)，則舍棄小數(shù)，向上取整。      1.2巴特沃斯濾波器的傳遞函數(shù)         巴特沃斯低通濾波器的傳遞函數(shù)，可由其振幅特性的分母多項(xiàng)式求得。其分母多項(xiàng)式根據(jù)S解開，可以得到極點(diǎn)。這里，為了方便處理，我們分為兩種情況去解這個方程。當(dāng)N為偶數(shù)的時候，這里，使用了歐拉公式。同樣的，當(dāng)N為奇數(shù)的時候

3、，同樣的，這里也使用了歐拉公式。歸納以上，極點(diǎn)的解為上式所求得的極點(diǎn)，是在s平面內(nèi)，在半徑為c的圓上等間距的點(diǎn)，其數(shù)量為2N個。為了使得其IIR濾波器穩(wěn)定，那么，只能選取極點(diǎn)在S平面左半平面的點(diǎn)。選定了穩(wěn)定的極點(diǎn)之后，其模擬濾波器的傳遞函數(shù)就可由下式求得。       1.3巴特沃斯濾波器的實(shí)現(xiàn)（C語言）          首先，是次數(shù)的計(jì)算。次數(shù)的計(jì)算，我們可以由下式求得。         其對應(yīng)的C語言程序?yàn)閏pp view plaincop

4、y1. N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /   2.             log10 (Stopband/Cotoff) );           然后是極點(diǎn)的選擇

5、，這里由于涉及到復(fù)數(shù)的操作，我們就聲明一個復(fù)數(shù)結(jié)構(gòu)體就可以了。最重要的是，極點(diǎn)的計(jì)算含有自然指數(shù)函數(shù)，這點(diǎn)對于計(jì)算機(jī)來講，不是太方便，所以，我們將其替換為三角函數(shù)，這樣的話，實(shí)部與虛部就還可以分開來計(jì)算。其代碼實(shí)現(xiàn)為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.           計(jì)算出穩(wěn)定的極點(diǎn)之后，就可以進(jìn)行傳遞函數(shù)的計(jì)算了。傳遞的函數(shù)的計(jì)算，就像下式一樣這里，為了得到模擬濾波

9、器的系數(shù)，需要將分母乘開。很顯然，這里的極點(diǎn)不一定是整數(shù)，或者來說，這里的乘開需要做復(fù)數(shù)運(yùn)算。其復(fù)數(shù)的乘法代碼如下，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.   有了乘法代碼之后，我們現(xiàn)在簡單的情況下，看看其如何計(jì)算其濾波器系數(shù)。我們做如下假設(shè)這個時候，其傳遞函數(shù)為將其乘開，其大致的關(guān)系就像下圖所示一樣。計(jì)算的關(guān)系一目了然，這樣的話，實(shí)現(xiàn)就簡單多了。高階的情況下也一樣，重復(fù)這種計(jì)算就可以了。其代碼為

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變換的原理        我們?yōu)榱藢⒛M濾波器轉(zhuǎn)換為數(shù)字濾波器的，可以用的方法很多。這里著重說說雙1次

21、z變換。我們希望通過雙1次z變換，建立一個s平面到z平面的映射關(guān)系，將模擬濾波器轉(zhuǎn)換為數(shù)字濾波器。        和之前的例子一樣，我們假設(shè)有如下模擬濾波器的傳遞函數(shù)。將其做拉普拉斯逆變換，可得到其時間域內(nèi)的連續(xù)微分方程式，其中，x(t)表示輸入，y(t)表示輸出。然后我們需要將其離散化，假設(shè)其采樣周期是T，用差分方程去近似的替代微分方程，可以得到下面結(jié)果然后使用z變換，再將其化簡?？傻玫饺缦陆Y(jié)果從而，我們可以得到了s平面到z平面的映射關(guān)系，即由于所有的高階系統(tǒng)都可以視為一階系統(tǒng)的并聯(lián)，所以，這個映射關(guān)系在高階系統(tǒng)中，也是成立的。然后，將關(guān)系式帶入上式，

22、可得到這里，我們可以就可以得到與的對應(yīng)關(guān)系了。         這里的與的對應(yīng)關(guān)系很重要。我們最終的目的設(shè)計(jì)的是數(shù)字濾波器，所以，設(shè)計(jì)時候給的參數(shù)必定是數(shù)字濾波器的指標(biāo)。而我們通過間接設(shè)計(jì)設(shè)計(jì)IIR濾波器時候，首先是要設(shè)計(jì)模擬濾波器，再通過變換，得到數(shù)字濾波器。那么，我們首先需要做的，就是將數(shù)字濾波器的指標(biāo)，轉(zhuǎn)換為模擬濾波器的指標(biāo)，基于這個指標(biāo)去設(shè)計(jì)模擬濾波器。另外，這里的采樣時間T的取值很隨意，為了方便計(jì)算，一般取1s就可以。       2.2雙1次z變換的實(shí)現(xiàn)（C語言）    &

23、#160;    我們設(shè)計(jì)好的巴特沃斯低通濾波器的傳遞函數(shù)如下所示。     我們將其進(jìn)行雙1次z變換，我們可以得到如下式子可以看出，我們還是需要將式子乘開，進(jìn)行合并同類項(xiàng)，這個跟之前說的算法相差不大。其代碼為。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、;    到此，我們就已經(jīng)實(shí)現(xiàn)了一個數(shù)字濾波器。3.IIR濾波器的間接設(shè)計(jì)代碼（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.