巴特沃斯滤波器c语言

1、模拟滤波器的设计1.1 巴特沃斯滤波器的次数根据给定的参数设计模拟滤波器，然后进行变数变换，求取数字滤波器的方法，称为滤波器的间接设计。做为数字滤波器的设计基础的模拟滤波器，称之为原型滤波器。这里，我们首先介绍的是最简单最基础的原型滤波器，巴特沃斯低通滤波器。由于 IIR 滤波器不具有线性相位特性，因此不必考虑相位特性，直接考虑其振幅特性。在这里， N 是滤波器的次数，c是截止频率。从上式的振幅特性可以看出，这个是单调递减的函数，其振幅特性是不存在纹波的。设计的时候，一般需要先计算跟所需要设计参数相符合的次数N。首先，就需要先由阻带频率，计算出阻带衰减将巴特沃斯低通滤波器的振幅特性，直接带入上

2、式，则有最后，可以解得次数N 为当然，这里的N 只能为正数，因此，若结果为小数，则舍弃小数，向上取整。1.2 巴特沃斯滤波器的传递函数巴特沃斯低通滤波器的传递函数，可由其振幅特性的分母多项式求得。其分母多项式1/13根据 S 解开，可以得到极点。这里，为了方便处理，我们分为两种情况去解这个方程。当N 为偶数的时候，这里，使用了欧拉公式。同样的，当N 为奇数的时候，同样的，这里也使用了欧拉公式。归纳以上，极点的解为上式所求得的极点，是在s 平面内，在半径为c的圆上等间距的点，其数量为2N 个。为了使得其IIR 滤波器稳定，那么，只能选取极点在S 平面左半平面的点。选定了稳定的极点之后，其模拟滤波

3、器的传递函数就可由下式求得。2/131.3 巴特沃斯滤波器的实现（C 语言）首先，是次数的计算。次数的计算，我们可以由下式求得。其对应的C 语言程序为cppview plaincopyN = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /log10 (Stopband/Cotoff) );然后是极点的选择，这里由于涉及到复数的操作，我们就声明一个复数结构体就可以了。最重要的是，极点的计算含有自然指数函数，这点对于计算机来讲，不是太方便，所以，我们将其替换为三角函数，这样的话，实部与虚部就还可以分开来计算。其代码实现为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计算出稳定的极点之后，就可以进行传递

5、函数的计算了。传递的函数的计算，就像下式一样这里，为了得到模拟滤波器的系数，需要将分母乘开。很显然，这里的极点不一定是整数，或者来说，这里的乘开需要做复数运算。其复数的乘法代码如下，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;有了乘法代码之后，我们现在简单的情况下，看看其如何计算其滤波器系数。我们做如下假设这个时候，其传递函数为将其乘开，其大致的关系就像下图所示一样。计算的关系一目了然，这样的话，实现就简单多了。高阶的情况下也一样，重复这种计算就可以了。其代码为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);到此，我们就可以得到一个模拟滤波器巴特沃斯低通滤波器了。2.双 1 次 z 变换2.1 双 1 次 z 变换的原理

9、我们为了将模拟滤波器转换为数字滤波器的，可以用的方法很多。这里着重说说双1 次 z 变换。我们希望通过双1 次 z 变换，建立一个s 平面到 z 平面的映射关系，将模拟滤波器转换为数字滤波器。和之前的例子一样，我们假设有如下模拟滤波器的传递函数。将其做拉普拉斯逆变换，可得到其时间域内的连续微分方程式，其中， x(t) 表示输入， y(t) 表示输出。然后我们需要将其离散化，假设其采样周期是T ，用差分方程去近似的替代微分方程，可以得到下面结果然后使用z 变换，再将其化简。可得到如下结果5/13从而，我们可以得到了s 平面到 z 平面的映射关系，即由于所有的高阶系统都可以视为一阶系统的并联，所以

10、，这个映射关系在高阶系统中，也是成立的。然后，将关系式带入上式，可得到这里，我们可以就可以得到 与 的对应关系了。这里的 与 的对应关系很重要。我们最终的目的设计的是数字滤波器，所以，设计时候给的参数必定是数字滤波器的指标。而我们通过间接设计设计IIR 滤波器时候，首先是要设计模拟滤波器，再通过变换，得到数字滤波器。那么，我们首先需要做的，就是将数字滤波器的指标，转换为模拟滤波器的指标，基于这个指标去设计模拟滤波器。另外，这里的采样时间T 的取值很随意，为了方便计算，一般取1s 就可以。2.2 双 1 次 z 变换的实现（ C 语言）我们设计好的巴特沃斯低通滤波器的传递函数如下所示。我们将其进

11、行双1 次 z 变换，我们可以得到如下式子6/13可以看出，我们还是需要将式子乘开，进行合并同类项，这个跟之前说的算法相差不大。其代码为。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;到此，我们就已经实现了一个数字滤波器。7/133.IIR 滤波器的间接设计代码（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