数字滤波器设计之⼀：巴特沃斯（Butterworth ）滤波器 我们常说的经典滤波器是根据傅⾥叶分析和变换设计出来的，只允许⼀定频率范围内的信号成分正常通过，⽽阻⽌另⼀部分频率成分通过。按照最佳逼近特性或者滤波通带特性分类，主要为巴特沃斯滤波器（Butterworth）、切⽐雪夫滤波器（Chebyshev）、贝塞尔滤波器（Bessel）和椭圆滤波器（Elliptic）四种。每种MATLAB都有相应的函数，⽤起来也⽐较⽅便，但是却缺少C/C++的程序，于是⾃⼰仔细研究了每种滤波器的特性和原理，并且部分滤波器实现了C语⾔的代码化，接下来的时间会对这些滤波器的原理和C语⾔的实现进⾏介绍。 该系列均以低通滤波器为原型来介绍，其他类型的滤波器可以由低通滤波器通过频率变换转换得到，这⾥不过多介绍。低通滤波器的主要性能指标有以下⼏个：通带截⽌频率fp、阻带截⽌频率fs、通带衰减 ( Ap)、阻带衰减 ( As) 以及归⼀化频率时需要⽤到的-3dB的转折频率fc。
1. Butterworth 滤波器原理
Butterworth滤波器因其在通带内的幅值特性具有最⼤平坦的特性⽽闻名，是四种经典滤波器中最简单的，巴特沃斯滤波器只需要两个
参数表征，滤波器的阶数N和-3dB处的截⽌频率。其幅度平⽅函数为：
N是滤波器的阶数，从幅度平⽅函数可以看出，N阶滤波器有2N个极点，⽽且这2N个极点均布在⼀个圆上，圆的半径为，称之为
Butterworth圆，Butterworth滤波器系统是⼀个线性系统，要使其稳定，其极点必须位于S平⾯的左半平⾯，所以取左半平⾯内的N个极点作为滤波器的极点，滤波器就是稳定的了，求出极点之后，计算模拟滤波器的系数as、bs，然后通过双线性变换（不懂得⾃⾏查书）由模拟域到数字域，求出系数az和bz 。最后通过差分⽅程就可以计算滤波结果了。 2. C 语⾔实现
A．求阶数
玻璃垫片
公式为：
代码：
B．求极点
公式：
代码：
1N = ceil(0.5*( log10 (( pow (10, Stopband_attenuation/10) - 1)/2 ( pow (10, Passband_attenuation/10) - 1)) / log10 (Stopband/Passband) ));
1for(k = 0;k <= ((2*N)-1) ; k++)
2 {
山人全息码
3 if(Cutoff*cos((2*k+N-1)*(pi/(2*N))) < 0.0)
4 {
5 poles[count].x = -Cutoff*cos((2*k+N-1)*(pi/(2*N)));
6 poles[count].y = -Cutoff*sin((2*k+N-1)*(pi/(2*N)));
7 count++;
8 if (count == N) break;
9 }
10 }
C．计算模拟滤波器系数
公式：
代码（这部分需要⾃⼰推导，多算⼏步就到规律了）：
1 Res[0].x = poles[0].x;
2 Res[0].y = poles[0].y;
3
4 Res[1].x = 1;
5 Res[1].y= 0;
6
7 for(count_1 = 0;count_1 < N-1;count_1++)//N个极点相乘次数
8 {
9 for(count = 0;count <= count_1 + 2;count++)
10 {
11 if(0 == count)
12 {
13 Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] );
14 }
15
16 else if((count_1 + 2) == count)
17 {
18 Res_Save[count].x += Res[count - 1].x;
19 Res_Save[count].y += Res[count - 1].y;
20 }
21 else
22 {
23 Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] );
24 Res_Save[count].x += Res[count - 1].x;
25 Res_Save[count].y += Res[count - 1].y;
26 }
27 }
28
29 for(count = 0;count <= N;count++)//Res[i]=a(i),i越⼤次数越⾼
30 {
31 Res[count].x = Res_Save[count].x;
32 Res[count].y = Res_Save[count].y;
33
34 *(a + N - count) = Res[count].x;
35 }
36 }
37
38 *(b+N) = *(a+N);
D．双线性变换
⽤下式进⾏替换 中的变量，得到 ，然后类似上⾯的计算⽅法计算Z域的系数az和bz，其中T为采样周期，但是因为在计算中会被约去，所以简化计算，这⾥取1。
公式：
代码：
1int Count = 0,Count_1 = 0,Count_2 = 0,Count_Z = 0;
2 double *Res, *Res_Save;
3 Res = new double[N+1]();
4 Res_Save = new double[N+1]();
5 memset(Res, 0, sizeof(double)*(N+1));
6 memset(Res_Save, 0, sizeof(double)*(N+1));
7
8 for(Count_Z = 0;Count_Z <= N;Count_Z++)
9 {
10 *(az+Count_Z) = 0;
11 *(bz+Count_Z) = 0;
12 }
13
14
15 for(Count = 0;Count<=N;Count++)
16 {
17 for(Count_Z = 0;Count_Z <= N;Count_Z++)
18 {
19 Res[Count_Z] = 0;
20 Res_Save[Count_Z] = 0;
21 }
22 Res_Save [0] = 1;
23
24 for(Count_1 = 0; Count_1 < N-Count;Count_1++)//计算（1-Z^-1）^N-Count的系数,
25 { //Res_Save[]=Z^-1多项式的系数，从常数项开始
26 for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
27 {
28 if(Count_2 == 0)
29 {规划沙盘
30 Res[Count_2] += Res_Save[Count_2];
31 }
32
33 else if((Count_2 == (Count_1+1))&&(Count_1 != 0))
34 {
35 Res[Count_2] += -Res_Save[Count_2 - 1];
36 }
37
38 else
39 {
40 Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];
41 }
42 }
43
44 for(Count_Z = 0;Count_Z<= N;Count_Z++)
45 {
46 Res_Save[Count_Z] = Res[Count_Z] ;
47 Res[Count_Z] = 0;
48 }
49 }
50
51 for(Count_1 = (N-Count); Count_1 < N;Count_1++)//计算(1-Z^-1)^N-Count*（1+Z^-1）^Count的系数,
52 { //Res_Save[]=Z^-1多项式的系数，从常数项开始
53 for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
54 {
55 if(Count_2 == 0)
55 if(Count_2 == 0)
56 {
57 Res[Count_2] += Res_Save[Count_2];
58 }
59
60 else if((Count_2 == (Count_1+1))&&(Count_1 != 0))
61 {
62 Res[Count_2] += Res_Save[Count_2 - 1];
63 }
64
立式导热油加热器65 else
66 {
67 Res[Count_2] += Res_Save[Count_2] + Res_Save[Count_2 - 1];
68 }
69 }
70
71 for(Count_Z = 0;Count_Z<= N;Count_Z++)
72 {
73 Res_Save[Count_Z] = Res[Count_Z] ;
74 Res[Count_Z] = 0;
75 }
76 }
77
78 for(Count_Z = 0;Count_Z<= N;Count_Z++)
79 {
80 *(az+Count_Z) += pow(2,N-Count) * (*(as+Count)) * Res_Save[Count_Z];
81 *(bz+Count_Z) += (*(bs+Count)) * Res_Save[Count_Z];
82 }
83
84 }//最外层for循环
85
86 for(Count_Z = N;Count_Z >= 0;Count_Z--)
87 {
88 *(bz+Count_Z) = (*(bz+Count_Z))/(*(az+0));
89 *(az+Count_Z) = (*(az+Count_Z))/(*(az+0));
90 }
E．由由差分⽅程计算滤波结果
采⽤滤波器直接II型结构，可以减少⼀半的中间缓存内存，具体代码如下：
1double FiltButter(double *pdAz, //滤波器参数表1
2 double *pdBz, //滤波器参数表2
3 int nABLen, //参数序列的长度
4 double dDataIn,//输⼊数据
5 double *pdBuf) //数据缓冲区
6{
7 int i;
8 int nALen;
9 int nBLen;
10 int nBufLen;
11 double dOut;
12
13 if( nABLen<1 )return 0.0;
14 //根据参数,⾃动求取序列有效长度
15 nALen = nABLen;
16 for( i=nABLen-1; i; --i )
17 {
18 if( *(pdAz+i) != 0.0 )//从最后⼀个系数判断是否为0
滤菌器19 {
20 nALen = i+1;
21 break;
22 }
23 }
24 //printf("%lf ", nALen);
25 if( i==0 ) nALen = 0;
26
27 nBLen = nABLen;
28 for( i=nABLen-1; i; --i )
29 {
30 if( *(pdBz+i) != 0.0 )
31 {
32 nBLen = i+1;
33
34 break;
35 }
36 }
37 //printf("%lf ", nBLen);
38 if( i==0 ) nBLen = 0;
39 //计算缓冲区有效长度
40 nBufLen = nALen;
41 if( nALen < nBLen)
42 nBufLen = nBLen;
43
44 //滤波: 与系数a卷乘
破真空阀
45 dOut = ( *pdAz ) * dDataIn; // a(0) * x(i)
46
47 for( i=1; i<nALen; i++) // a(i) * w(n-i),i=1toN
48 {
49 dOut -= *(pdAz+i) * *(pdBuf + (nBufLen - 1) - i);
50 }
51
52 //卷乘结果保存为缓冲序列的最后⼀个
53 *(pdBuf + nBufLen - 1) = dOut;
54
55 //滤波: 与系数b卷乘
56 dOut = 0.0;
57 for( i=0; i<nBLen; i++) // b(i) * w(n-i)
58 {
59 dOut += *(pdBz+i) * *(pdBuf + (nBufLen - 1) - i);
60 }
61
62 //丢弃缓冲序列中最早的⼀个数, 最后⼀个数清零
63 for( i=0; i<nBufLen-1; i++)
64 {
65 *(pdBuf + i) = *(pdBuf + i + 1);
66 }
67 *(pdBuf + nBufLen - 1) = 0;
68
69 //返回输出值
70 return dOut;
71}
VC6.0开发环境滤波结果如下：
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议