巴特沃斯滤波器c语言.docx
《巴特沃斯滤波器c语言.docx》由会员分享，可在线阅读，更多相关《巴特沃斯滤波器c语言.docx（29页珍藏版）》请在冰点文库上搜索。
巴特沃斯滤波器c语言
巴特沃斯滤波器ｃ语言
————————————————————————————————作者：
————————————————————————————————日期:
ﻩ
1.模拟滤波器的设计
１.1巴特沃斯滤波器的次数
根据给定的参数设计模拟滤波器,然后进行变数变换，求取数字滤波器的方法，称为滤波器的间接设计。
做为数字滤波器的设计基础的模拟滤波器,称之为原型滤波器。
这里，我们首先介绍的是最简单最基础的原型滤波器，巴特沃斯低通滤波器。
由于IIR滤波器不具有线性相位特性,因此不必考虑相位特性，直接考虑其振幅特性。
在这里，Ｎ是滤波器的次数,Ωc是截止频率。
从上式的振幅特性可以看出，这个是单调递减的函数,其振幅特性是不存在纹波的。
设计的时候，一般需要先计算跟所需要设计参数相符合的次数Ｎ。
首先，就需要先由阻带频率，计算出阻带衰减
将巴特沃斯低通滤波器的振幅特性，直接带入上式,则有
最后,可以解得次数N为
当然，这里的N只能为正数,因此，若结果为小数，则舍弃小数，向上取整。
1.2巴特沃斯滤波器的传递函数
巴特沃斯低通滤波器的传递函数，可由其振幅特性的分母多项式求得。
其分母多项式
根据Ｓ解开，可以得到极点。
这里，为了方便处理，我们分为两种情况去解这个方程。
当Ｎ为偶数的时候,
这里，使用了欧拉公式
。
同样的，当Ｎ为奇数的时候，
同样的,这里也使用了欧拉公式。
归纳以上,极点的解为
上式所求得的极点,是在s平面内,在半径为Ωc的圆上等间距的点，其数量为2N个。
为了使得其IIＲ滤波器稳定,那么，只能选取极点在S平面左半平面的点。
选定了稳定的极点之后，其模拟滤波器的传递函数就可由下式求得。
1.3巴特沃斯滤波器的实现（C语言）
首先,是次数的计算。
次数的计算，我们可以由下式求得。
其对应的C语言程序为
[ｃpp］ ｖiewplainｃopy
1.Ｎ ＝ Ceil（０.5*（ ｌoｇ１０ （ ｐｏw （10, Stopband_attenuation/1０） － １） /
2. log１０ （Ｓtｏpbaｎd/Cｏtoff） ））;
然后是极点的选择，这里由于涉及到复数的操作，我们就声明一个复数结构体就可以了。
最重要的是，极点的计算含有自然指数函数,这点对于计算机来讲,不是太方便，所以，我们将其替换为三角函数,
这样的话，实部与虚部就还可以分开来计算。
其代码实现为
[cpp] vｉew ｐｌａinｃoｐy
1.tyｐeｄef sｔruct
2.{
3. doｕblｅ Rｅａl_ｐart；
4. doubｌe Imaｇ＿Ｐaｒｔ;
5.} COＭPLEX;
6.
7.
8.COMPLEＸ ｐoｌes[N]；
9.
10.for（ｋ = 0;k <= （（２*N）-1） ； k＋+）
11.｛
12. if（Ｃoｔofｆ＊coｓ（（k+dk）＊（pi/N）） < ０）
13. {
14. ｐoｌｅs[ｃouｎt］.Reaｌ_ｐaｒt = -Cotｏff*ｃos（（k＋dk）*（pｉ/N））;
15. poｌeｓ［ｃounｔ]．Imａg_Pａｒt＝ -Cotofｆ*sｉn（（k+dk）*（pｉ/N））;
16. coｕｎt++；
17. if （couｎt == N） breａk;
18. ｝
19.}
计算出稳定的极点之后，就可以进行传递函数的计算了。
传递的函数的计算，就像下式一样
这里，为了得到模拟滤波器的系数,需要将分母乘开。
很显然，这里的极点不一定是整数，或者来说,这里的乘开需要做复数运算。
其复数的乘法代码如下，
[cpp］ viｅwplａiｎｃoｐy
1.int Complｅx＿Ｍuｌｔiple（CＯMPＬEX ａ,CＯＭＰLＥX b,
2. ｄoｕbｌｅ *Res_Rｅal，dｏuｂle ＊Ｒes_Ｉmag）
3.
4.{
5. ＊（Ｒeｓ_Rｅａl） ＝ （a.Rｅal＿paｒt）*（b．Real_part） - （ａ.Iｍag_Pａrt）*（b.Imag_Ｐａｒt）;
6. *（Reｓ_Iｍag）= （ａ.Imａg_Part）*（b．Real_ｐａｒt） ＋ （ａ.Rｅal＿paｒt）＊（b.Imag_Part）;
7. returｎ （int）１；
8.}
有了乘法代码之后,我们现在简单的情况下，看看其如何计算其滤波器系数。
我们做如下假设
这个时候，其传递函数为
将其乘开，其大致的关系就像下图所示一样。
计算的关系一目了然,这样的话,实现就简单多了。
高阶的情况下也一样,重复这种计算就可以了。
其代码为
[cpｐ] vｉewplaincｏpy
1. Ｒes[０].Real_pａrt ＝ ｐoles［０].Ｒeal_paｒt;
2. Ｒｅｓ[0].Ｉmag_Ｐart= ｐoles[0］．Ｉｍag_Part；
3. Res[1]．Rｅaｌ＿pａｒt = 1;
4. Res[1].Imag_Pａｒt＝ 0;
5.
6.for（count＿1 ＝ ０；ｃount_1 ＜ N-1;coｕnt＿1++）
7.｛
8. fｏr（cｏunt = 0;ｃount <＝ coｕnｔ_1 ＋ 2;count++）
9. {
10. if（0 =＝ count）
11. {
12. Complex＿Ｍultipｌｅ（Ｒｅｓ［coｕnｔ], poles［count_１＋1],
13. &（Ｒes_Ｓave［count］．Real_paｒｔ）,
14. &（Rｅs_Savｅ[couｎt].Imag_Pａrt））;
15. }
16. ｅlｓe if（（count_1 + 2） == ｃount）
17. {
18. Res_Save[couｎt].Real＿ｐart ＋＝ Ｒeｓ[ｃount - 1].Reａl_paｒt;
19. Reｓ_Saｖe[cｏunt］.Iｍaｇ_Part += Rｅs[count － 1].Imaｇ_Ｐaｒt;
20. ｝
21. else
22. {
23. Compleｘ_Muｌtiple（Rｅs[counｔ]， pｏles[couｎt_1+1]，
24. &（Res＿Save[cｏｕnt]．Reaｌ＿parｔ）,
25. &（Res_Ｓavｅ［counｔ].Ｉmag＿Ｐart））；
26.1 Res＿Savｅ[counｔ］.Rｅal＿part ＋= Reｓ[ｃount - １].Ｒeaｌ_part;
27. Res_Savｅ[couｎｔ]．Imag_Pａｒｔ += Ｒes［count - 1］．Imag_Parｔ;
28. ｝
29. }
30. *（b+N） = *（a+Ｎ）;
到此,我们就可以得到一个模拟滤波器巴特沃斯低通滤波器了。
ﻫ
２.双１次ｚ变换
２.1双1次z变换的原理
我们为了将模拟滤波器转换为数字滤波器的,可以用的方法很多。
这里着重说说双1次ｚ变换。
我们希望通过双1次ｚ变换，建立一个ｓ平面到ｚ平面的映射关系,将模拟滤波器转换为数字滤波器。
和之前的例子一样，我们假设有如下模拟滤波器的传递函数。
将其做拉普拉斯逆变换，可得到其时间域内的连续微分方程式，
其中,ｘ（t）表示输入,y（t）表示输出。
然后我们需要将其离散化，假设其采样周期是T，用差分方程去近似的替代微分方程，可以得到下面结果
然后使用z变换,再将其化简。
可得到如下结果
从而，我们可以得到了s平面到z平面的映射关系，即
由于所有的高阶系统都可以视为一阶系统的并联，所以,这个映射关系在高阶系统中，也是成立的。
然后，将关系式
带入上式，可得
到这里,我们可以就可以得到Ω与ω的对应关系了。
这里的Ω与ω的对应关系很重要。
我们最终的目的设计的是数字滤波器,所以，设计时候给的参数必定是数字滤波器的指标。
而我们通过间接设计设计IIR滤波器时候，首先是要设计模拟滤波器，再通过变换，得到数字滤波器。
那么，我们首先需要做的，就是将数字滤波器的指标，转换为模拟滤波器的指标，基于这个指标去设计模拟滤波器。
另外，这里的采样时间T的取值很随意,为了方便计算，一般取１ｓ就可以。
２.2双1次z变换的实现（C语言）
我们设计好的巴特沃斯低通滤波器的传递函数如下所示。
我们将其进行双１次z变换,我们可以得到如下式子
可以看出,我们还是需要将式子乘开,进行合并同类项,这个跟之前说的算法相差不大。
其代码为。
[cpp] viewplaincoｐy
1.ｆｏr（Ｃｏuｎｔ ＝ 0;Couｎt<=N;Coｕnｔ++）
2. {
3. for（Counｔ＿Z ＝ ０;Count_Ｚ ＜= N;Ｃoｕｎt_Z++）
4. ｛
5. Rｅs[Ｃoｕｎｔ_Z] = 0；
6. Rｅｓ_Ｓaｖe[Cｏunｔ_Z］ = 0；
7. ｝
8. Res_Save ［0] ＝ １;
9. ｆor（Ｃount＿１ = ０; Ｃouｎt_1 < N-Count;Count_１++）
10. {
11. foｒ（Count_2 = ０; Cｏｕnｔ_2 <＝ Ｃｏuｎt_１＋1;Count_2++）
12. {
13. iｆ（Couｎt_2 =＝ ０） Res[Ｃounｔ_2］ ＋= Ｒeｓ_Ｓave[Coｕnt_2］;
14. elｓｅ iｆ（（Coｕｎt＿2 == （Ｃｏuｎｔ_1+1））&＆（Count＿１ !
= 0））
15. Rｅs［Ｃount_２] += -Res_Savｅ[Ｃｏｕnｔ_2 - 1］;
16. elsｅ Rｅs[Coｕnt_2] += Res_Sａve［Count_2] - Rｅs_Save[Couｎｔ＿2 - １]；
17. ｆｏr（Cｏｕnt_Z ＝ 0；Couｎｔ_Z<= N;Coｕnt_Ｚ＋+）
18. {
19. Rｅs_Save［Ｃｏunt_Ｚ] ＝ Res[Coｕｎt_Z］ ；
20. Ｒes［Count_Z] = 0；
21. }
22. }
23. ｆor（Count＿1 = （N-Ｃouｎt）； Counｔ＿１ < N;Cｏuｎt_1++）
24. {
25. ｆor（Counｔ_2 = 0; Coｕnt＿2 <= Couｎｔ_1+1；Counｔ_2＋+）
26. {
27. if（Coｕnt_2 == ０） Reｓ[Coｕnt_2］ +＝ Ｒeｓ_Save［Count_２];
28. eｌse if（（Couｎt_２ == （Ｃoｕnt＿1+1））&&（Count＿１ !
＝ 0））
29. Rｅｓ[Ｃount_２] += Res＿Save［Ｃoｕnｔ_2 - 1];
30. else
31. Ｒｅs［Cｏuｎt_2] += Res_Savｅ[Cｏunt_2] + Reｓ_Ｓaｖe[Ｃount＿2 - 1];
32. }
33. ｆor（Count_Z ＝ 0；Cｏunt_Z＜= Ｎ;Count＿Z++）
34. {
35. Rｅｓ＿Saｖe[Cｏunt＿Ｚ] = Rｅs［Ｃount_Z］ ;
36. Res［Couｎt_Ｚ] ＝ ０；
37. }
38. ｝
39. for（Count_Z = 0;Count_Z<＝ N;Cｏunt＿Ｚ++）
40. {
41. *（az+Count_Ｚ） += pow（２,Ｎ－Coｕnｔ） * （*（aｓ+Counｔ）） *
42. Res＿Ｓａve[Count_Z］；
43. *（bz+Cｏunｔ_Ｚ） +＝ （*（bｓ+Ｃouｎt）） * Rｅs_Save[Counｔ＿Ｚ］;
44. ｝
45. }
到此,我们就已经实现了一个数字滤波器。
3.ＩIR滤波器的间接设计代码（C语言）
[cpｐ] vｉew plaincｏｐｙ
1.＃inｃlude ＜stｄｉo.h>
2.#incluｄe
3.#ｉｎclude ＜mａlloｃ.ｈ>
4.#inclｕde 5.
6.
7.#define ｐi （（double）3.１41592６）
8.
9.
10.ｓtrｕct DＥSＩGN_SＰＥCIFICAＴIOＮ
11.{
12. ｄｏuble Cotｏff;
13. dｏuｂle Sｔopbａnd；
14. doｕｂle Stopbａnd＿atｔeｎｕation;
15.｝;
16.
17.ｔypedeｆ ｓtｒuct
18.{
19. dｏuｂlｅ Real_part；
20. double Imag＿Pａｒｔ;
21.} COMPLEX;
22.
23.
24.
25.iｎt Ceil（doubｌｅ input）
26.{
27. iｆ（inｐut ！
= （ｉｎt）iｎpuｔ） ｒeturｎ （（inｔ）ｉnput） ＋1；
28. else return （（ｉnｔ）ｉnput）；
29.}
30.
31.
32.int Cｏmplex_Mｕｌｔiple（CＯMPLＥX a,COMＰLEX b
33. ,ｄoｕble ＊Res_Rｅal，dｏuｂle *Res_Imaｇ）
34.
35.{
36. ＊（Ｒes_Real） ＝ （a.Ｒeal＿part）＊（ｂ.Reaｌ_ｐaｒt） - （a.Ｉmag＿Part）*（ｂ．Imag_Pａrt）;
37. *（Ｒes_Imag）= （a．Ｉｍag_Part）＊（b.Rｅal_pａrt） ＋ （a.Reａｌ＿pａrt）＊（b.Ｉmag＿Part）；
38. rｅturn （iｎｔ）１;
39.}
40.
41.
42.iｎt Buｔｔｏrd（dｏuble Ｃｏtｏfｆ,
43. douｂlｅ Stopbａｎd,
44. ｄouble Ｓtｏpband_atｔenｕatioｎ）
45.{
46. iｎt N；
47.
48. printf（＂Wc = ％lｆ [rad/ｓec] \n" ，Coｔoｆｆ）;
49. prｉｎｔf（"Ws = %lf ［rad/sｅc] \ｎ" ，Stｏpband）;
50. pｒｉｎtf（"Ａs = %lｆ [ｄB] \n" ,Stｏｐband_aｔtｅnｕaｔion）；
51. ｐrinｔf（"-－---－-－----－-----－----－－－-－－----------－-----－-------－--\n＂ ）;
52.
53. N = Ceｉl（0.5*（ log10 （ pow （10, Ｓｔopbａｎd＿atｔenuatioｎ/10） - 1） ／
54. log１0 （Ｓｔｏpband/Cotｏff） ））;
55.
56.
57. retｕrn （int）N;
58.}
59.
60.
61.int Butｔｅｒ（int Ｎ， doｕblｅ Cotｏｆｆ,
62. double *a,
63. doｕｂｌｅ ＊b）
64.{
65. doublｅ dk ＝ ０;
66. ｉｎt k ＝ 0；
67. int coｕnｔ = ０,count＿1 = 0;
68. COMPLEX poles[N];
69. ＣOＭPLEX Ｒes[Ｎ+1],Res＿Save[N＋1］;
70.
71. if（（N％2） ＝= 0） dk = ０．5;
72. elｓe ｄk = 0;
73.
74. fｏr（k = 0;k <＝ （（2*N）-１） ; ｋ++）
75. ｛
76. if（Cotofｆ*ｃos（（k+dk）*（pi/N）） < ０）
77. ｛
78. pｏｌes[counｔ].Rｅal＿paｒｔ = -Cｏｔｏff＊cｏs（（k+dk）*（ｐi/N））;
79. ｐoleｓ[cｏｕnt].Imag＿Part= -Cｏｔoff*sin（（k+dk）＊（pi/N））;
80. count++;
81. ｉf （cｏuｎt ＝= N） breａk;
82. }
83. }
84.
85. printf（"Pk = \ｎ＂ ）;
86. fｏr（ｃouｎt ＝ 0；ｃｏuｎt < N ;ｃount++）
87. {
88. pｒｉntf（＂（%lｆ） + （%lf i） ＼n" ,-poｌes［couｎt].Ｒｅal_part
89. ,-ｐoles[count].Imag＿Part）;
90. ｝
91. pｒinｔf（"－---－----－-－---－----------－--－---－-------－------------－-\n" ）；
92.
93. Res[0].Real_ｐaｒｔ ＝ poｌes[0].Reａl_paｒt;
94. Res［0]．Imaｇ_Paｒt= poles[0].Imaｇ_Ｐaｒt；
95.
96. Reｓ[1]．Reaｌ＿part = １;
97. Res［1］.Iｍａg＿Parｔ= ０；
98.
99. ｆｏr（couｎｔ_1 ＝ 0;ｃｏｕnt_1 ＜ N-1;coｕｎt_1++）
100. {
101. for（couｎt ＝ 0;count <= couｎt_1 + 2;coｕnt++）
102. {
103. if（0 ＝= count）
104. {
105. Coｍｐlex_Mｕltiplｅ（Reｓ[counｔ], ｐｏleｓ［coｕnt＿1+1］，
106. ＆（Ｒｅs_Sａve［count］.Ｒeal_part）,
107. &（Reｓ_Save[ｃｏｕｎｔ].Imag_Part））;
108. //pｒintf（ "Res_Sａve :
（%ｌf） ＋ （%lf i） \n" ，Rｅs＿Sａｖｅ［０].Real_paｒt,Res_Ｓave[0].Imag_Ｐaｒt）;
109. ｝
110.
111. eｌｓe if（（coｕn
升级会员 