巴特沃斯濾波器c語言

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是截止頻率。從上式的振幅特性可以看出，這個(gè)是單調(diào)遞減的函數(shù)，其振幅特性是不存在紋波的。設(shè)計(jì)的時(shí)候，一般需要先計(jì)算跟所需要設(shè)計(jì)參數(shù)相符合的次數(shù)N。首先，就需要先由阻帶頻率，計(jì)算出阻帶衰減將巴特沃斯低通濾波器的振幅特性，直接帶入上

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

3、器的傳遞函數(shù)就可由下式求得。2/131.3 巴特沃斯濾波器的實(shí)現(xiàn)（C 語言）首先，是次數(shù)的計(jì)算。次數(shù)的計(jì)算，我們可以由下式求得。其對(duì)應(yīng)的C 語言程序?yàn)閏ppview plaincopyN = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /log10 (Stopband/Cotoff) );然后是極點(diǎn)的選擇，這里由于涉及到復(fù)數(shù)的操作，我們就聲明一個(gè)復(fù)數(shù)結(jié)構(gòu)體就可以了。最重要的是，極點(diǎn)的計(jì)算含有自然指數(shù)函數(shù)，這點(diǎn)對(duì)于計(jì)算機(jī)來講，不是太方便，所以，我們將其替換為三角函數(shù)，這樣的話，實(shí)部與虛部就還可以分開來計(jì)算。其代碼實(shí)現(xiàn)為cpp

4、view plaincopy1.typedefstructdouble Real_part;double Imag_Part; COMPLEX;6.7.COMPLEX polesN;for (k = 0;k = (2*N)-1) ; k+)if (Cotoff*cos(k+dk)*(pi/N) 0)polescount.Real_part = -Cotoff*cos(k+dk)*(pi/N);15.polescount.Imag_Part= -Cotoff*sin(k+dk)*(pi/N);count+;17.if(count = N)break ;3/13計(jì)算出穩(wěn)定的極點(diǎn)之后，就可以進(jìn)行傳遞

5、函數(shù)的計(jì)算了。傳遞的函數(shù)的計(jì)算，就像下式一樣這里，為了得到模擬濾波器的系數(shù)，需要將分母乘開。很顯然，這里的極點(diǎn)不一定是整數(shù)，或者來說，這里的乘開需要做復(fù)數(shù)運(yùn)算。其復(fù)數(shù)的乘法代碼如下，cppview plaincopyint Complex_Multiple(COMPLEX a,COMPLEX b,2.double*Res_Real,double*Res_Imag)3.*(Res_Real) = (a.Real_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);*(Res_Imag)= (a.Imag_Part)*(b.Real_part) +

6、(a.Real_part)*(b.Imag_Part);return ( int )1;有了乘法代碼之后，我們現(xiàn)在簡單的情況下，看看其如何計(jì)算其濾波器系數(shù)。我們做如下假設(shè)這個(gè)時(shí)候，其傳遞函數(shù)為將其乘開，其大致的關(guān)系就像下圖所示一樣。計(jì)算的關(guān)系一目了然，這樣的話，實(shí)現(xiàn)就簡單多了。高階的情況下也一樣，重復(fù)這種計(jì)算就可以了。其代碼為cppview plaincopyRes0.Real_part = poles0.Real_part;Res0.Imag_Part= poles0.Imag_Part;Res1.Real_part = 1;4/13Res1.Imag_Part= 0;for (count_

7、1 = 0;count_1 N-1;count_1+)for (count = 0;count = count_1 + 2;count+)if (0 = count)Complex_Multiple(Rescount, polescount_1+1,13.&(Res_Savecount.Real_part),14.&(Res_Savecount.Imag_Part);16.elseif (count_1 + 2) = count)Res_Savecount.Real_part += Rescount - 1.Real_part;Res_Savecount.Imag_Part += Rescou

8、nt - 1.Imag_Part;elseComplex_Multiple(Rescount, polescount_1+1,24.&(Res_Savecount.Real_part),25.&(Res_Savecount.Imag_Part);26.1Res_Savecount.Real_part += Rescount - 1.Real_part;Res_Savecount.Imag_Part += Rescount - 1.Imag_Part;*(b+N) = *(a+N);到此，我們就可以得到一個(gè)模擬濾波器巴特沃斯低通濾波器了。2.雙 1 次 z 變換2.1 雙 1 次 z 變換的原理

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

10、，這個(gè)映射關(guān)系在高階系統(tǒng)中，也是成立的。然后，將關(guān)系式帶入上式，可得到這里，我們可以就可以得到 與 的對(duì)應(yīng)關(guān)系了。這里的 與 的對(duì)應(yīng)關(guān)系很重要。我們最終的目的設(shè)計(jì)的是數(shù)字濾波器，所以，設(shè)計(jì)時(shí)候給的參數(shù)必定是數(shù)字濾波器的指標(biāo)。而我們通過間接設(shè)計(jì)設(shè)計(jì)IIR 濾波器時(shí)候，首先是要設(shè)計(jì)模擬濾波器，再通過變換，得到數(shù)字濾波器。那么，我們首先需要做的，就是將數(shù)字濾波器的指標(biāo)，轉(zhuǎn)換為模擬濾波器的指標(biāo)，基于這個(gè)指標(biāo)去設(shè)計(jì)模擬濾波器。另外，這里的采樣時(shí)間T 的取值很隨意，為了方便計(jì)算，一般取1s 就可以。2.2 雙 1 次 z 變換的實(shí)現(xiàn)（ C 語言）我們?cè)O(shè)計(jì)好的巴特沃斯低通濾波器的傳遞函數(shù)如下所示。我們將其進(jìn)

11、行雙1 次 z 變換，我們可以得到如下式子6/13可以看出，我們還是需要將式子乘開，進(jìn)行合并同類項(xiàng)，這個(gè)跟之前說的算法相差不大。其代碼為。cppview plaincopyfor (Count = 0;Count=N;Count+)3.for (Count_Z = 0;Count_Z = N;Count_Z+)5.ResCount_Z = 0;6.Res_SaveCount_Z = 0;8.Res_Save 0 = 1;9.for (Count_1 = 0; Count_1 N-Count;Count_1+)11.for(Count_2 = 0; Count_2 = Count_1+1;Cou

12、nt_2+)12.13.if (Count_2 = 0)ResCount_2 += Res_SaveCount_2;14.elseif(Count_2 = (Count_1+1)&(Count_1 != 0)15.ResCount_2 += -Res_SaveCount_2 - 1;16.elseResCount_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.23.f

13、or (Count_1 = (N-Count); Count_1 N;Count_1+)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(

14、Count_Z = 0;Count_Z= N;Count_Z+)34.35.Res_SaveCount_Z = ResCount_Z ;36.ResCount_Z = 0;37.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;到此，我們就已經(jīng)實(shí)現(xiàn)了一個(gè)數(shù)字濾波器。7/133.IIR 濾波器的間接設(shè)計(jì)代碼（C 語言）cpp

15、view plaincopy#include #include #include #include 7.#definepi(double)3.1415926)8.9.struct DESIGN_SPECIFICATIONdouble Cotoff;double Stopband;double Stopband_attenuation;16.17.typedefstructdouble Real_part;double Imag_Part; COMPLEX;22.23.24.int Ceil( double input)27.if (input != (int )input)return(int

16、)input) +1;else return ( int )input);30.31.32.intComplex_Multiple(COMPLEX a,COMPLEX b33.,double *Res_Real, double *Res_Imag)34.*(Res_Real) = (a.Real_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);*(Res_Imag)= (a.Imag_Part)*(b.Real_part) + (a.Real_part)*(b.Imag_Part);return ( int )1;40.41.42.int B

17、uttord(double Cotoff,43.doubleStopband,44.doubleStopband_attenuation)int N;48.printf(Wc =%lfrad/sec n,Cotoff);49.printf(Ws =%lfrad/sec n,Stopband);50.printf(As =%lfdB n,Stopband_attenuation);51.printf(-n);8/1352.N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /54.log10 (Stopband/Cotoff

18、) );55.56.return ( int )N;59.60.61.int Butter(int N, double Cotoff,62.double*a,63.double*b)64.65.doubledk = 0;66.intk = 0;67.intcount = 0,count_1 = 0;COMPLEX polesN;COMPLEX ResN+1,Res_SaveN+1;if (N%2) = 0) dk = 0.5;else dk = 0;73.for (k = 0;k = (2*N)-1) ; k+)76.if (Cotoff*cos(k+dk)*(pi/N) 0)78.poles

19、count.Real_part = -Cotoff*cos(k+dk)*(pi/N);polescount.Imag_Part= -Cotoff*sin(k+dk)*(pi/N);count+;81.if(count = N)break ;85.printf(Pk =n);for (count = 0;count N ;count+)88.printf(%lf) + (%lf i) n,-polescount.Real_part89.,-polescount.Imag_Part);91.printf(-n);92.Res0.Real_part = poles0.Real_part;Res0.I

20、mag_Part= poles0.Imag_Part;Res1.Real_part = 1;Res1.Imag_Part= 0;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);9/13108./print

21、f( Res_Save : (%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.&(Res_S

22、avecount.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_part;126.Res_Save

23、count.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.Rescount.Imag_Part= Res_Savecoun

24、t.Imag_Part;138.139.*(a + N - count) = Rescount.Real_part;140.141.142./printf(There! n );46.*(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.printf(as = );157.for (count = 0;count = N ;co

25、unt+)158.159.printf(%lf , *(a+count);160.10/13161.printf( n);162.163.printf(-n);164.165.return( int) 1;169.int Bilinear(int N,170.double*as,double*bs,171.double*az,double*bz)173.int Count = 0,Count_1 = 0,Count_2 = 0,Count_Z = 0;174.doubleResN+1;175.doubleRes_SaveN+1;176.177.for(Count_Z = 0;Count_Z =

26、 N;Count_Z+)178.179.*(az+Count_Z) = 0;180.*(bz+Count_Z) = 0;84.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;190.191.Res_Save 0 = 1;192.193.for(Count_1 = 0; Count_1 N-Count;Count_1+)194.195.for (Count_2 = 0; Cou

27、nt_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)&(Count_1 != 0)204.205.ResCount_2 += -Res_SaveCount_2 - 1;206./printf( Res%d %lf n , Count_2 ,ResCount_2);207.208.209.else21

28、0.211.ResCount_2 += Res_SaveCount_2 - Res_SaveCount_2 - 1;212./printf( Res%d %lf n , Count_2 ,ResCount_2);16./printf( Res : );217.for (Count_Z = 0;Count_Z= N;Count_Z+)11/13218.219.Res_SaveCount_Z = ResCount_Z ;220.ResCount_Z = 0;221./printf( %d %lf ,Count_Z, Res_SaveCount_Z);222.223./pr

29、intf( n );27.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.ResCount_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