巴特沃斯濾波器c語言

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

2、特性，直接帶入上式，則有最后，可以解得次數(shù)n 為.當(dāng)然，這里的n 只能為正數(shù)，因此，若結(jié)果為小數(shù)，則舍棄小數(shù)，向上取整。1.2 巴特沃斯濾波器的傳遞函數(shù)巴特沃斯低通濾波器的傳遞函數(shù)，可由其振幅特性的分母多項式求得。其分母多項式根據(jù) s 解開，可以得到極點。這里，為了方便處理，我們分為兩種情況去解這個方程。當(dāng) n為偶數(shù)的時候，這里，使用了歐拉公式。同樣的，當(dāng)n 為奇數(shù)的時候，.同樣的，這里也使用了歐拉公式。歸納以上，極點的解為上式所求得的極點，是在s 平面內(nèi)，在半徑為c 的圓上等間距的點，其數(shù)量為2n 個。為了使得其iir 濾波器穩(wěn)定， 那么， 只能選取極點在s 平面左半平面的點。選定了穩(wěn)定的極

3、點之后，其模擬濾波器的傳遞函數(shù)就可由下式求得。 1.3 巴特沃斯濾波器的實現(xiàn)（c 語言）首先，是次數(shù)的計算。次數(shù)的計算，我們可以由下式求得。其對應(yīng)的 c 語言程序為cppview plaincopy1.n = ceil(0.5*( log10 ( pow (10, stopband_attenuation/10) - 1) / 2. log10 (stopband/cotoff) ); 然后是極點的選擇，這里由于涉與到復(fù)數(shù)的操作，我們就聲明一個復(fù)數(shù)結(jié)構(gòu)體就可以了。最重要的是，極點的計算含有自然指數(shù)函數(shù)，這點對于計算機來講，不是太方便，所以，我們將其替換為三角函數(shù)，.這樣的話，實部與虛部就還可以

4、分開來計算。其代碼實現(xiàn)為cppview plaincopy1.typedefstruct2. 3.double real_part; 4.double imag_part; 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. 14. polescount.real_part = -cotoff*cos(k+dk)*(pi/n); 15. polescount.imag_part= -cotoff*sin(k+dk)*(pi/n

5、); 16. count+; 17.if (count = n) break ; 18. 19. 計算出穩(wěn)定的極點之后，就可以進行傳遞函數(shù)的計算了。傳遞的函數(shù)的計算，就像下式一樣這里，為了得到模擬濾波器的系數(shù)，需要將分母乘開。很顯然，這里的極點不一定是整數(shù)，或者來說，這里的乘開需要做復(fù)數(shù)運算。其復(fù)數(shù)的乘法代碼如下，cppview complex_multiple(complex a,complex b, 2.double *res_real,double *res_imag) 3.4. 5. *(res_real) = (a.real_part)*(b.real

6、_part) - (a.imag_part)*(b.imag_part); 6. *(res_imag)= (a.imag_part)*(b.real_part) + (a.real_part)*(b.imag_part); 7.return ( int )1; 8. 有了乘法代碼之后，我們現(xiàn)在簡單的情況下，看看其如何計算其濾波器系數(shù)。我們做如下假設(shè)這個時候，其傳遞函數(shù)為將其乘開，其大致的關(guān)系就像下圖所示一樣。計算的關(guān)系一目了然，這樣的話，實現(xiàn)就簡單多了。高階的情況下也一樣，重復(fù)這種計算就可以了。其代碼為cppview plaincopy1. res0.real_part = poles0.r

7、eal_part; 2. res0.imag_part= poles0.imag_part; 3. res1.real_part = 1; .4. res1.imag_part= 0; 5.6.for (count_1 = 0;count_1 n-1;count_1+) 7. 8.for (count = 0;count = count_1 + 2;count+) 9. 10.if(0 = count) 11. 12. complex_multiple(rescount, polescount_1+1, 13. &(res_savecount.real_part), 14. &

8、(res_savecount.imag_part); 15. 16.elseif (count_1 + 2) = count) 17. 18. res_savecount.real_part += rescount - 1.real_part; 19. res_savecount.imag_part += rescount - 1.imag_part; 20. 21.else22. 23. complex_multiple(rescount, polescount_1+1, 24. &(res_savecount.real_part), 25. &(res_savecount.

9、imag_part); 26.1 res_savecount.real_part += rescount - 1.real_part; 27. res_savecount.imag_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次 z 變換。我們希望通過雙1 次 z 變換，建立一個s 平面到 z 平面的映射關(guān)系，將模擬濾波器轉(zhuǎn)換為數(shù)字

10、濾波器。和之前的例子一樣，我們假設(shè)有如下模擬濾波器的傳遞函數(shù)。將其做拉普拉斯逆變換，可得到其時間域內(nèi)的連續(xù)微分方程式，.其中， x(t)表示輸入， y(t) 表示輸出。然后我們需要將其離散化，假設(shè)其采樣周期是t，用差分方程去近似的替代微分方程，可以得到下面結(jié)果然后使用z 變換，再將其化簡。可得到如下結(jié)果從而，我們可以得到了s 平面到 z 平面的映射關(guān)系，即由于所有的高階系統(tǒng)都可以視為一階系統(tǒng)的并聯(lián)，所以， 這個映射關(guān)系在高階系統(tǒng)中，也是成立的。然后，將關(guān)系式帶入上式，可得到這里，我們可以就可以得到 與 的對應(yīng)關(guān)系了。.這里的 與 的對應(yīng)關(guān)系很重要。我們最終的目的設(shè)計的是數(shù)字濾波器，所以， 設(shè)計

11、時候給的參數(shù)必定是數(shù)字濾波器的指標(biāo)。而我們通過間接設(shè)計設(shè)計iir 濾波器時候， 首先是要設(shè)計模擬濾波器，再通過變換，得到數(shù)字濾波器。那么，我們首先需要做的，就是將數(shù)字濾波器的指標(biāo)，轉(zhuǎn)換為模擬濾波器的指標(biāo)，基于這個指標(biāo)去設(shè)計模擬濾波器。另外，這里的采樣時間t 的取值很隨意，為了方便計算，一般取1s 就可以。 2.2 雙 1 次 z 變換的實現(xiàn)（ c 語言）我們設(shè)計好的巴特沃斯低通濾波器的傳遞函數(shù)如下所示。我們將其進行雙1 次 z 變換，我們可以得到如下式子可以看出， 我們還是需要將式子乘開，進行合并同類項，這個跟之前說的算法相差不大。其代碼為。cppview plaincopy1.for (co

12、unt = 0;count=n;count+) 2. 3.for (count_z = 0;count_z = n;count_z+) 4. 5. rescount_z = 0; 6. res_savecount_z = 0; 7. 8. res_save 0 = 1; 9.for (count_1 = 0; count_1 n-count;count_1+) 10. 11.for (count_2 = 0; count_2 = count_1+1;count_2+) 12. 13.if (count_2 = 0) rescount_2 += res_savecount_2; 14.elsei

13、f (count_2 = (count_1+1)&(count_1 != 0) .15. rescount_2 += -res_savecount_2 - 1; 16.else rescount_2 += res_savecount_2 - res_savecount_2 - 1; 17.for (count_z = 0;count_z= n;count_z+) 18. 19. res_savecount_z = rescount_z ; 20. rescount_z = 0; 21. 22. 23.for (count_1 = (n-count); count_1 n;count_1

14、+) 24. 25.for (count_2 = 0; count_2 = count_1+1;count_2+) 26. 27.if (count_2 = 0) rescount_2 += res_savecount_2; 28.elseif (count_2 = (count_1+1)&(count_1 != 0) 29. rescount_2 += res_savecount_2 - 1; 30.else31. rescount_2 += res_savecount_2 + res_savecount_2 - 1; 32. 33.for (count_z = 0;count_z=

15、 n;count_z+) 34. 35. res_savecount_z = rescount_z ; 36. rescount_z = 0; 37. 38. 39.for (count_z = 0;count_z= n;count_z+) 40. 41. *(az+count_z) += pow(2,n-count) * (*(as+count) * 42. res_savecount_z; 43. *(bz+count_z) += (*(bs+count) * res_savecount_z; 44. 45. 到此，我們就已經(jīng)實現(xiàn)了一個數(shù)字濾波器。3.iir 濾波器的間接設(shè)計代碼（c 語言

16、）cppview plaincopy1.#include 2.#include .3.#include 4.#include 5.6.7.#define pi (double)3.1415926)8.9.10.struct design_specification 11. 12.double cotoff; 13.double stopband; 14.double stopband_attenuation; 15.; 16.17.typedefstruct18. 19.double real_part; 20.double imag_part; 21. complex;

17、5.int ceil(double input) 26. 27.if (input != (int )input) return (int )input) +1; 28.elsereturn (int )input); 29. 30.31.32.int complex_multiple(complex a,complex b 33. ,double *res_real,double *res_imag) 34.35. 36. *(res_real) = (a.real_part)*(b.real_part) - (a.imag_part)*(b.imag_part); 37. *(res_im

18、ag)= (a.imag_part)*(b.real_part) + (a.real_part)*(b.imag_part); 38.return ( int )1; 39. 40.41.42.int buttord(double cotoff, 43.double stopband, 44.double stopband_attenuation) .45. 46.int n; 47.48. printf(wc = %lf rad/sec n ,cotoff); 49. printf(ws = %lf rad/sec n ,stopband); 50. printf(as = %lf db n

19、 ,stopband_attenuation); 51. printf(-n ); 52.53. n = ceil(0.5*( log10 ( pow (10, stopband_attenuation/10) - 1) / 54. log10 (stopband/cotoff) ); 55.56.57.return ( int )n; 58. 59.60.61.int butter(int n, double cotoff, 62.double *a, 63.double *b) 64. 65.double dk = 0; 66.int k = 0; 67.int count = 0,cou

20、nt_1 = 0; 68. complex polesn; 69. complex 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)-1) ; k+) 75. 76.if (cotoff*cos(k+dk)*(pi/n) 0) 77. 78. polescount.real_part = -cotoff*cos(k+dk)*(pi/n); 79. polescount.imag_part= -cotoff*sin(k+dk)*(pi/n); 80. count

21、+; 81.if (count = n) break ; 82. 83. 84.85. printf(pk = n ); 86.for (count = 0;count n ;count+) 87. 88. printf(%lf) + (%lf i) n ,-polescount.real_part .89. ,-polescount.imag_part); 90. 91. printf(-n ); 92.93. res0.real_part = poles0.real_part; 94. res0.imag_part= poles0.imag_part; 95.96. res1.real_p

22、art = 1; 97. res1.imag_part= 0; 98.99.for (count_1 = 0;count_1 n-1;count_1+) 100. 101.for (count = 0;count = count_1 + 2;count+) 102. 103.if(0 = count) 104. 105. complex_multiple(rescount, polescount_1+1, 106. &(res_savecount.real_part), 107. &(res_savecount.imag_part); 108./printf( res_save

23、 : (%lf) + (%lf i) n ,res_save0.real_part,res_save0.imag_part);109. 110.111.elseif (count_1 + 2) = count) 112. 113. res_savecount.real_part += rescount - 1.real_part; 114. res_savecount.imag_part += rescount - 1.imag_part; 115. 116.else117. 118. complex_multiple(rescount, polescount_1+1, 119. &(

24、res_savecount.real_part), 120. &(res_savecount.imag_part); 121.122./printf( res : (%lf) + (%lf i) n ,rescount - 1.real_part,rescount - 1.imag_part);123./printf( res_save : (%lf) + (%lf i) n ,res_savecount.real_part,res_savecount.imag_part);124.125. res_savecount.real_part += rescount - 1.real_pa

25、rt; 126. res_savecount.imag_part += rescount - 1.imag_part; 127.128./printf( res_save : (%lf) + (%lf i) n ,res_savecount.real_part,res_savecount.imag_part);129.130. 131./printf(there n );132. 133.134.for (count = 0;count = n;count+) 135. 136. rescount.real_part = res_savecount.real_part; 137. rescou

26、nt.imag_part= res_savecount.imag_part; 138.139. *(a + n - count) = rescount.real_part; 140. 141.142./printf(there! n );143.144. 145.146. *(b+n) = *(a+n); 147.148./-display-/149. printf(bs = ); 150.for (count = 0;count = n ;count+) 151. 152. printf(%lf , *(b+count); 153. 154. printf( n ); 155.156. pr

27、intf(as = ); 157.for (count = 0;count = n ;count+) 158. 159. printf(%lf , *(a+count); 160. 161. printf( n ); 162.163. printf(-n );164.165.return ( int ) 1; 166. 167.168.169.int bilinear(int n, .170.double *as,double *bs, 171.double *az,double *bz) 172. 173.int count = 0,count_1 = 0,count_2 = 0,count

28、_z = 0; 174.double resn+1; 175.double res_saven+1; 176.177.for (count_z = 0;count_z = n;count_z+) 178. 179. *(az+count_z) = 0; 180. *(bz+count_z) = 0; 181. 182.183.184.for (count = 0;count=n;count+) 185. 186.for (count_z = 0;count_z = n;count_z+) 187. 188. rescount_z = 0; 189. res_savecount_z = 0; 1

29、90. 191. res_save 0 = 1; 192.193.for (count_1 = 0; count_1 n-count;count_1+) 194. 195.for (count_2 = 0; count_2 = count_1+1;count_2+) 196. 197.if (count_2 = 0) 198. 199. rescount_2 += res_savecount_2; 200./printf( res%d %lf n , count_2 ,rescount_2);201. 202.203.elseif (count_2 = (count_1+1)&(cou

30、nt_1 != 0) 204. 205. rescount_2 += -res_savecount_2 - 1; 206./printf( res%d %lf n , count_2 ,rescount_2);207. 208.209.else210. .211. rescount_2 += res_savecount_2 - res_savecount_2 - 1; 212./printf( res%d %lf n , count_2 ,rescount_2);213. 214. 215.216./printf( res : );217.for (count_z = 0;count_z= n

31、;count_z+) 218. 219. res_savecount_z = rescount_z ; 220. rescount_z = 0; 221./printf( %d %lf ,count_z, res_savecount_z); 222. 223./printf( n );224.225. 226.227.for (count_1 = (n-count); count_1 n;count_1+) 228. 229.for (count_2 = 0; count_2 = count_1+1;count_2+) 230. 231.if (count_2 = 0) 232. 233. r

32、escount_2 += res_savecount_2; 234./printf( res%d %lf n , count_2 ,rescount_2);235. 236.237.elseif (count_2 = (count_1+1)&(count_1 != 0) 238. 239. rescount_2 += res_savecount_2 - 1; 240./printf( res%d %lf n , count_2 ,rescount_2);241. 242.243.else244. 245. rescount_2 += res_savecount_2 + res_save

33、count_2 - 1; 246./printf( res%d %lf n , count_2 ,rescount_2);247. 248. .249.250./ printf( res : );251.for (count_z = 0;count_z= n;count_z+) 252. 253. res_savecount_z = rescount_z ; 254. rescount_z = 0; 255./printf( %d %lf ,count_z, res_savecount_z); 256. 257./printf( n );258. 259.260.261./printf( res : );262.for (count_z = 0;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.280.281./-display-/282. printf(bz = ); 283.for (