巴特沃斯滤波器c语言x

1、巴特沃斯滤波器C语言 日期：1.模拟滤波器的设计1 .1巴特沃斯滤波器的次数根据给定的参数设计模拟滤波器 ,然后进行变数变换,求取数字滤波器的方法,称为滤波器的间接设计.做为数字滤波器的设计根底的模拟滤波器,称Z为原型滤波器.这里,我们首先介绍的是最简单最根底的原型滤波器,巴特沃斯低通滤波器.由于HR滤波器不具 有线性相位特性,因此不必考虑相位特性,直接考虑其振幅特性.+焙严在这里,N是滤波器的次数,Qc是截止频率.从上式的振幅特性可以看出,这个是单 调递减的函数,其振幅特性是不存在纹波的.设计的时候,一般需要先计算跟所需要设计参数 相符合的次数N.首先,就需要先由阻带频率,计算出阻带衰减4二

2、-201.歸凤丿业|将巴特沃斯低通滤波器的振幅特性,直接带入上式 ,那么有儿二一20 lo贏= 101.弘最后,可以解得次数N为N_llogZ10 寻-1) 2应密1.2巴特沃斯滤波器的传递函数巴特沃斯低通滤波器的传递函数,可由其振幅特性的分母多项式求得.其分母多项式当然,这里的N只能为正数,因此,假设结果为小数,那么舍弃小数,向上取整.2.V根据S解开,可以得到极点.这里,为了方便处理,我们分为两种情况去解这个方程.当N 为偶数的时候,1+这里,使用了欧拉公式N为奇数的时候,加+ = °斤二0丄22N 1O同样的,当1-=0同样的,这里也使用了欧拉公式.归纳以上,极点的解为5 k

3、= 02卫N1M + 1 兀n :even number2N J(k exp j一odd number .k = O.L25.,2N-1 N 上式所求得的极点,是在S平面内,在半径为Qc的圆上等间距的点,其数量为2N个.为了使得其 HR滤波器稳定,那么,只能选取极点在S平面左半平面的点.选定了稳定的极点之后,其模拟滤波 器的传递函数就可由下式求得.1.3巴特沃斯滤波器的实现C语言首先,是次数的计算.次数的计算,我们可以由下式求得.4?lgg止口其对应的C语言程序为:C pp v phh c opy1. N = C e il(0 .5*(1 o g 1 0( p 0 w ( 10, S topb

4、and_attenuatbnA0 ) 1 ) /2.bg1 0 ( S t 0 pba n d/C 0 to田)；然后是极点的选择,这里由于涉及到复数的操作,我们就声明一个复数结构体就可以了.最重要的是,极点的计算含有自然指数函数,这点对于计算机来讲,不是太方便,所以,我们将其替换为三角函数,门七+ (1/2)门.cos兀+ jQsm° Nc Nck 冗 口 kjvQ cos + /Q sm , c N c N 'W 化 N: even number 5N: odd number * = 0丄2 = 0?L2.这样的话,实部与虚部就还可以分开来计算.其代码实现为Epp v i

5、 ew p 1 a h c o p y1. ty p e d ef st juct2. 3. do u bl e R e a L P art；4. doub 1 e Im a g _ Part ;5. C0 M PLEX ;6.7.6. COMPLE X po 1 esN；9.7. fcr ( k = 0<= (2 疥)-1 ): k +)8. 9. if (C o t of f * co s ( ( k+dk) * (p述)<0)10. 11. pol escoun t R ea 1 _ p a r t= -C ot o3cos(k+ dk)*(p i /N)12. po 1 e

6、 s c oun t . In a g_P a r t = C otof f *s i n «k+dk)*(p i /N );16co u n 廿+ ；17. if ( cou nt二二 N) bre a k;18. 19. 计算出稳定的极点z后,就可以进行传递函数的计算了.传递的函数的计算,就像 下式一样h 左 x n 二-这里,为了得到模拟滤波器的系数 ,需要将分母乘开.很显然,这里的极点不一定是整数,或者来说,这里的乘开需要做复数运算.其复数的乘法代码如下,fcpp vi e w pl a in c op C om pl ex _M u 1 t pie (COMP

7、 LEX a ,C 0 M P LEX b,2.d o u b 1 e 駅 es_R e al, double* R es_ I m ag)3.4.5.*(R <? s _R e a D =( a.R e al pa r R eaLpart) - (a I mag_P a rt*(bin ag_P ar t);6.*(Resm ag)= ( a. In a g_P art ) *(b . R eaLpart) + (a.Re alpa r t)* bin ag_Par$;7. retur n (int ) 1 ；8. 右了乘法代码Z后,我们现在简单的情况下,看看其如何计算其滤波器系数.我

8、们做如下假 设N = 2, a = g +皿4 =入 +Jb,这个时候,其传递函数为ResOReaLP a rt Res 01 I m ag_ P artF比® 飞-购冷-矚+加将其乘开,其大致的关系就像下列图所示一样.Rtf Inaw lai61010_ fli « 2 a + b« *- 10计算的关系一目了然,这样的话,实现就简单多了.高阶的情况下也一样 ,重复这种计算就可以To其代码为fcp p v i ew phhc o pyp oles 0 R eal_pa r t;p oles£) I m ag_P art ；3. Rest! R e a

9、1 _ p a r 1;4 ResUlIn ag_P a r t= 0;5.6 for (count _1 = 0 ； c ount_l< N -1:count_1 +)7.&f 0 r fc0 unt二 0 ;c ount = co unt._1+ 2 :count+)9.10.if(0 = = count)11.12.Com plex_ M uhp 1e ( Re s co Un t 1 polescount_1 + U13.&(R es_ S avecount R eal_pa r t ),14.& (Re s_S av ecou ntl.In ag_P a

10、 rt);15.16.e1 S e if (count_l + 2)二二c ount)17.18.4-R es_S ave fcoun £R ealpart+ = Re s c ount- 11R ea Lpa19.v,Res_S a v e fco untI m a g_Part+二 Rescount llln ag _Pa rt;20.21.slse22.23.Com piex _M u 1 tple (Res fcount , p 0 les Poun 匸 1 + 124& (Res_S ave fco unt R ea 1 _ par t ),25.& (R

11、es_S av ecoun t 1I m ag_ P art) ；26.1Res _S ave fcoun t R e al_ part4二 R e s c ount-1 1 R ea1 _part;27.Res_Sav efcou n t In ag_Pa rt + 二Res count-1.In ag_Par t :28. 29. 30. * ( b+N)二 *(a+N)；到此,我们就可以得到一个模拟滤波器巴特沃斯低通滤波器了.92 .双1次z变换2.1双1次z变换的原理我们为了将模拟滤波器转换为数字滤波器的,可以用的方法很多. 这里着重说说双1次Z变换.我们希望通过双 1次Z变换,建立一

12、个S平面到Z平面的映射关系,将模拟滤波器转换为数字滤波器.和Z前的例子一样,我们假设有如下模拟滤波器的传递函数.H(沪宀 (s + a)将其做拉普拉斯逆变换,可得到其时间域内的连续微分方程式,響+砒2斶)其屮,X (0表示输入,y(t)表示输出.然后我们需要将其离散化,假设其采样周期是方程去近似的替代微分方程,可以得到下面结果当 b&) -一 1)十彳+ 歹 5 -1)卜彳 k(n)+xs -1)JL然后使用z变换,再将其化简.可得到如下结果ri + z1从而,我们可以得到了 S平面到Z平面的映射关系,即由于所有的高阶系统都可以视为一阶系统的并联,所以,这个映射关系在高阶系统中,也是成

13、立的.然后,将关系式z ejc: and s = 3 + jQ带入上式,可得2 1严s =:T 1 +严2 1 一 costo - /sin cdT 1 +cosco + j sin to2 l-(cos2c?j + sin2 co) - j2sin(yT2 + 2cos9到这里,我们可以就可以得到 Q与3的对应关系了.这里的Q与3的对应关系很重要.我们最终的目的设计的是数字滤波器,所以,设计时候给的参数必定是数字滤波器的指标.而我们通过间接设计设计HR滤波器时候,首先是要设计模拟滤波器,再通过变换,得到数字滤波器.那么,我们首先需要做的,就是将数 字滤波器的指标,转换为模拟滤波器的指标,基于

14、这个指标去设计模拟滤波器.另外, 这里的采样时间T的取值很随意,为了方便计算,一般取1 s就可以.2 .2双1次z变换的实现C语言我们设计好的巴特沃斯低通滤波器的传递函数如下所示.)=口aQs + aYs 斗+.43 + 0片我们将其进行双1次 z变换,我们可以得到如下式子f-i可以看岀,我们还是需要将式子乘开,进行合并同类项,这个跟之前说的算法相差不大.其代码'帕.IPpp view phhco p y1f o r ( C o un t=0 £ ou n tK=N c o u n t + )I2-3.for Coun t _ Z = 0 £ ount_Z <

15、= N ; C o u n 匸Z+)4.5.ResCount_Z=0 ；6.ReS _ S a v e C 0 un t_Z = 0 ；I 7.8.R es_S ave0 = 1 :9.f or(C ount_ 1 二 0 : C ou n 匸 1 < N - ount£ ount_1 +)10.11.fo r (C ount_2 =0 ; C o u n t _2 < = C o u n t_ 1+ 1 ;C ount_2+)12.13.i f C ou n t_2 = =0 ) R esCoun t _2 + 二 R e S _ S avent_2 14.i15.16

16、.2.el s e i f (C o u n t 2 二二(C o u n t _1 + 1)& & (C ount 1 =0)Re s C ount_ 2 += -R es_S av eCount_2 -1;els e R e s IC o u n 匸 2+=Res_Sa ve C ount_2 -R es_S ave lC ou nt _ 2 -1 ；f o r C o u nt_Z = 0 ； C ou n t _Z二 N £o u nt_ Z + +)Re £

17、s_S ave Co unt_ Z =Reslo u n 匸Z ；ResC ount_Z J = 0；f or (C ount _1 = (N-C ou n t ) ； C ount _ 1< N f oun t_l+ )ifor C oun t _2 二 0;C ount2<=Cou n t 1 + 1 ； Coun t 2+ + )if C o u nt_2 =0 ) R e sfou nt_2 + = R e s _S ave C ount_2;e 1 se if ( ( C ou n t 2 二二(C o u nt 1 + 1 )&&Q ount 1 !=0

18、)R e s C ount_ 2+二 Res _ S ave Count-2-11;Res;C o u n 匸2 + 二 R es_S ave C o un匸2 + R es _ Sa ve C ount_ 2-1;6.rf or C ount_Z = 0；C o unt_Z V 二 N :Count _Z+)Re s _ Sa ve C 0 unt_Z =R e s C ount_Z ;R esC ou n 匸 Z =o ；0.41.t) *fcr C ount_Z 二 0 ount_Z<=N £ o unt_Z +)42.43.(az

19、+C ount_ Z) += pow (2, NCount )*(*G s +C ounR es_ S a ve 1 oun匸Z ；(bz+C 0 un t _ Z ) + = (*(b s + C ou nt) *R e s_S ave £ oun telseZ44.45.到此,我们就已经实现了一个数字滤波器.3. I E滤波器的间接设计代码C语言fcp p v i ew phhc o p y1.# h c lide < st d i o.h>2.tthcli d e <m a th h>3.# i n elide V ma lbc.h>4.ttinc

20、l u de <stri ngh>5.6.7.ttdefinePi (double)3. 1 41592 6 )&9.10.s tr u c tD ESIGN._S P E CFEATD N11. 12d 0 uble C ot 0 ff;13. d 0 u b le S t opb a nd ；14. do u b le S topb a nd _at t e n u atbn;15. ;1617. t ypede f str uct18. 20. double In ag _P art;21. JCOMPLEX;5. in tCe£L( d

21、oub 1 e input)26. 27. i f (h put !二(int)input) r eiur n ( h t ) i nput) +1 ；28. else return ( ( i nt) i nput)；I29. 30.31.32.htC 0 m plex_M u 1 t ph(C 0 MPLEXa,C0MP LEX b33.d0 u ble * R es_R GJ al, d 0 u b le *Res_Inag)34.35.36.* (R es_R eaD=QR eal_.part)* ( b R ea 1 _ p ar t) -( a. I m ag_Part)*( b

22、. In ag_P art);37. *(R es_lm ag)二 6.Im ag_P art) * bR e ap a rt) + (a.R e a 1 _ p a rt) *(b Im ag_ P art)；38. r e turn ( i n t ) 1 ;39. 1.42.i n tBu t t0id ( d O ubleC o toff ,43.dou b 1e S topb and,44.d oubleS t O pband_att en u atb n )45.46.in tN :47.48.printf(HW c=%1fBad/S ec巾,Cot of f

23、);49.pri ntf(V s = % Ifiad/4 en ,S t 0 pband);50.p r i ntf("As = % 1f d B,St op band_a t t e n u a t bn)；51.p rin tf( f n )；52.53. N = C e i 1(0.5* ( loglO ( pow (10, S t opb a n d_ at t enuatb n /10) -1)54. log1 0( S t o pband/C ot 0 田 )；55.56.55. ret u m ( ht )N ;56. 59.60.61.htBut t er (iit

24、N, do ubl e Co toff,62.63.do u b 1 e* b)9.doubleu n t 二 0 ,count _1 = 0 ;COMPLEX pobsN;ht coC 0 MPLEX ResN +l,ResS ave IN71. if ( W % 2)=二 0 ) dk 二 0 5；72. else d k 二 0 ;74. for(k 二 0女<=(2 粗)一1 ) :k+)75. 76. if C otof* c os ( (k+dk ) * ( pi/N) <0 )double 也78.o 1 es fcoun t .

25、R e al_ pa r t = -C o t o ff * c o s «k+dk)*8.89.90.述；PartF -Ccount±+ ;i f fc o u n t =二 N) bre a k;prhtf(n );f o r ( c ou nt = 0 : c o u n tN;r i ntf( H (% 1 f ) +( % 疋?91. pr in t f(,z- c ount+-+)(k+dk) * ©述);n,z,-po 1 es cou n tl R e apartp oles count.In

26、 ag_P art);5.96.97.R es01R eaL Part=po 1 es 01R e a Lpa r t;01.102.Res0 Im a g_Pa r IF poleslOllm aR e S 1 . R ea 1 _ part二 1 ;R es1 .1 m a g_Par tg _Pa r t；for (cou n t _1 = 0; cou n匸 1< NT ;co u n 匸 1+)for fcou nt = 0 ;count<:z coun t_l + 2 ;co u n廿+)103.104.if © =二

27、 count)12.113.else(i if (fco u n匸 1+ 2)二=co u n t )Re s_ S ave tou ntlR eaLpar t + 二 ResB 0105.C om p lex_Mu kpl e (R ej S fcount , p O Jes co u nt_1 + 1 ,106.& ( R es_S a vecount Real_part),107.&(Res_S avec0 u n t lbag_Part);108./p r htf( /ZR es_Sa ve :( %1 f) +& If D VT,R

28、e s_ S a0 .R eal_pa r tR es_ S ave 0 丄 In ag_ Part);115.unt- 11R eal_part ;114.Re s_Save fcountJ.Inag_Part + = Re sfcount-1 Im ag_P art;+ 1,116.e 1 se117.118.C om pex_M uitip le ( R e s co unt, polesco unt_ 1119.& ( R es_S av e c o u nt.R e a 1 _p a r t ),120.&(R es_Sav e fco u n tI mag _121

29、.25.R e sfcou n t 1 R eaLpart,R es fcoun t -1. In a g _P art);unt-1126.a rt;127.128.fcount.R49.150.152.%ri n tf(tl. ReapartR es_SR es S aR eaLpart; R es_S a ve :( % 10 +ave co u nt In ag_Pa r t);0 u nt).R eal_par t

30、 + 二 Res coR es_S ave fco u nt. I m ag_ Part + 二 R esfcount-ills _S ave :隔 10 + &ma g_ PSavee aLpar t ,R es_ S ave fcountJ.In a/ / pr i ntf(z,Theref o r ( c o u n t = 0 : c o u ntV 二S fcount.R eapart= Ri ntf* G+N)二 *G+Ng _P art):N ； c ounV-+)e s S avecount. R eaLpa r t;e s fcoun t J.In ag _ Part

31、F R es_S avN -count)二pr i n t f(,zbs = D；0 unt.In ag _P art;R es fc o untlR e aLpart;d i s p 1 ay-for (count二 0: c o unt二 N ;coun t +)printf,* (b +coun t)；54.pr i n tf" ；155.156.p ri n tf( a s 二了);157.fcr (poun t=0:coun t < = N ; c oun廿+)158.159.pri ntf(總 1 f ,*(a +c

32、ou n t );160.161.pr intf("'巾；162.163.p rhtf( n ；164.165.r eiur n(int ) 1;69.i ntB ilh e ar( ht N,170.d o u ble * a s, double*bs,171.d ouble轮 z , d oubh*bz)172.173.i n t C ou n t = 0,COun 匸 1 = 0,Count _2=0 , C o u n 1174.d 0 u b JeR es N + 1;175.doub 1 e Res_Sav e N + 1 ;177.fc

33、r ( COuntZ 二0 f ountZ <= N ; C ount_ Z+)178.179.*(az+C ount_Z)二 0 ;180.*Gz+Co un t _ Z)二 0 ；181.183.184.fcr C ount二0 C ount<=N C0 un廿+)185.186.f 0r ( C ou n t_ Z = 0；Coun t_Z<=N ； Count_Z + +)187.188.ResC o u n t _Z二 0 ;189.R es_Sa v e C ount_Z = 0 ;190.191.R e s_S ave 0= 1;193.f 0r C ount_

34、l - 0； C ou n匸 1 < N C o u nt ； C o u nt194.195.for C oun t _ 2=0；Cou n t _ 2<= C ou n t_l + l C1 + +)oun t 2 + + )Z 二 0;192.196.197.if ( C o unt_2 二二 0)198.199.R es Co unt 2 + 二 R e s200.Zferintf ( Res%d_% 1CountS ave C ou n 匸 2 ;2 ,R es C o u201.202.215.nt_2 )；203.1 =0)204.205.206._2 ,R es C

35、 o10.211.e L C o untelsei f (C oun t 2C ou n t _1 + 1) && C oun 匸R es C o+二-R es/ / p rhtf ( ,ZRun匸2 );Re2 -1else2+二 RS av e to u nt_2 -1;n", C o unte s_S ave C oun t _2 - R e s_S a212.2 ,R esC oun匸2);t f ( n Res % d %n,z , C ount213.214.216./pri n tfCRe217.0 unt_Z =0 £

36、;t_Z<= N ;C oun t _ Z + + )218.219.e s S ave 1 ounZ 二 R esQ ount_220.221.rhtfCI% d% F,Count_Z,Res_Sa ve C oun t222.223./ht f ( n )；224.225.227.for C ount_=0J-C ount);C o u nt_l < N £ ou n 匸 1+)229.for C o u n t_2 =0 ; C oun t _2二 C oun t _ 1 + 1 ； C o u231.nt 2+)if C ount_2 二二 0)227.228.2

37、34.s C oun t 2)235.236.237.untl = 0)44.245.Co55.ZD;256.257.258.259.260.261.262.263.264.t 2 += R/hC ount_2 ,R es S3)；r i n tf( H R e s 1% d % 1 foun t_2 ,R eseount)*i f ( (C oun t_ 2 = = (Count_l +1)&&e s Co unt_ 2 +二 R

38、 es S aveC o u nt_2 -1;/ /prht f("Res風 d%unt_ 2 );R es2 -1C ount_ 2 +Res Save ount_ 2 +R e s_ S aver i n t f ( z, Res隅 d% IfC o u n t_2 ,R esIC ounRe/ prin t f (,ZR e s :,r ( C o u nt_ Z 二 0 C oun Z <= N £ ounZ+ +)S _S ave C o u npZ二 RestRest ou n t Z二 0 ;/ferht f (隔 d %1 f ,z ,C ount_Z

39、,Res _ S ave C 0 un tfo/ P rntfC/ p rh 社("Res : 0;r C ount_ Z0 ； C ount_Z< = N C 0 un t _Z + +)* ( a z+ C ount_ Z) + = p ow(2, N-CRes _S a v e C ount Z ;nt)265.Z；z + C ou n t_Z ) + 二(*( b s+ C o unt) *R esS ave C ount266.Zt>rhtf(z,% If",*©z+Coun；268./冷 r i ntf ( W)267.269.270.27

40、1.272.273.274.fcr ( C o un t_Z 二 N;Co u n匸 Z >=0o un 匸 Z)275.276.z +C ou n t_Z )=(*(b z+C o u nt_ Z)/毬 Gz+ 0 );277.*(az +C 0 u nt_ Z )二(* ( az +C ount_ Z ) / ( * ( a z+0);278.279.280.281./ 一 ./a i s p Jay-/ /282.P rntf( ,?bz =");283.f o r ( C ount_Z 二 0 ;C o u n t_Z <=N f o u nt_ Z + + )284.285.printf(羯 If" , *( bz + C ou n t_Z)；286.1287.Pr htf( n ")；288.pr htf ("az = D;289.for ( C ount_ Z 二 0; C o u n t _Z < = N ; Co unZ + +)290.291.Printf( f",* (az+ C o u nt_ Z );292.293.294.prntf(" W):prin tf(-n "