数字滤波器设计之二：贝塞尔（Bessel）滤波器

1.前言

经典滤波器按照最佳逼近特性可分为巴特沃斯（Butterworth）、切比雪夫（Chebyshev）、贝塞尔（Bessel）和椭圆（Elliptic）滤波器四种。滤波器按照频带又可分为低通、高通、带通和带阻四种，其中高通、带通和带阻滤波器均可由低通滤波器由频率变换得到，本文仅考虑低通滤波器的分析与设计。N阶巴特沃斯滤波器由于在0频率处的2N-1阶导数都等于零，所以称为最大平坦近似，在通带和阻带幅值响应均单调递减；切比雪夫滤波器可分为切比雪夫Ⅰ型滤波器和切比雪夫Ⅱ型滤波器，Ⅰ型滤波器幅值响应具有通带等波纹变化而阻带单调变化，Ⅱ型滤波器则相反，在通带单调变化而阻带等波纹变化。椭圆滤波器的幅值响应在通带和阻带都是等波纹变化。贝塞尔滤波器具有最大平坦的群时延。上述前三种滤波器对于阶跃响应均具有很大的过冲，只有贝塞尔滤波器对阶跃响应有很小的过冲。当被测信号有很多阶跃信号时，应该采用贝塞尔滤波器进行滤波。

2.Bessel滤波器的原理

2.1Bessel滤波器的幅值响应

Bessel滤波器的传递函数如下：

所以Bessel模拟低通滤波器的幅值特性函数为：

其中为Bessel多项式，可以通过递推公式得到，如下所示：

...

将Bessel多项式写成一般n次多项式的形式为

由公式可以得到的各个系数，如表1所示

表1

1	1
2	3	3
3	15	15	6
4	105	105	45	10
5	945	945	420	105	15
6	10395	10395	4725	1260	210	21
7	135135	135135	62370	17325	3150	378	28
8	2027025	2027025	945945	270270	51975	6930	630	36

Bessel滤波器不像Butterworth、Chebyshev和Elliptic滤波器的幅值响应，没有简单的方法和公式确定 $B_{n}(s)=0$ 的根，即Bessel幅值响应的极点。我们可以通过计算机数值计算的方法得到Bessel幅值响应的极点，如表2所示。

表2

1	-1.0000000
2	-1.5000000	±j0.8660254
3	-2.3221854	-1.8389073	±j1.7543810
4	-2.8962106	±j0.8672341	-2.1037894	±j2.6574180
5	-3.6467386	-3.3519564	±j1.7426614	-2.3246743	±j3.5710229
6	-4.2483594	±j0.8675097	-3.7357084	±j2.6262723	-2.5159322	±j4.4926730
7	-4.9717869	-4.7582905	±j1.7392861	-4.0701392	±j3.5171740	-2.6856769	±j5.4206941
8	-5.5878860	±j0.8676144	-2.8389840	±j6.3539113	-4.3682892	±j4.4144425	-5.2048408	±j2.6161751

2.2 Bessel滤波器的群时延

Bessel滤波器是具有线性相位响应的线性滤波器，在整个通频带具有恒定的群时延。下图1为2到10阶的Bessel滤波器的群时延。

3. Bessel滤波器的C语言实现

一个N阶数字滤波器输入和输出的关系可以用N次差分方程表示，如下所示，

用传递函数表示如下

当系数全部为零时表示有限冲激响应滤波器，系数不全为零时表示无限冲激响应滤波器。若求得系数和，便可由输入信号递推得到。流程图如下图所示：

Bessel滤波器程序设计流程图对应步骤如下：

1.由滤波器的性能指标确定滤波器的阶数N

Bessel滤波器不像巴特沃斯和切比雪夫滤波器，有公式来计算滤波器的阶数，Bessel滤波器的阶数只能通过给定的性能指标通过已有的设计曲线查找得到，主要通过时延D和通带衰减来确定阶数N，本文所设计的程序是直接人为给定滤波器的阶数。

2.由阶数N查表确定模拟滤波器的系数和

Bessel滤波器的极点没有简单的计算公式，只能通过数值方法计算得到，所以滤波器的系数通过计算极点和零点得到，所以本文采用查表的方法直接根据滤波器的阶数确定滤波器的系数和。则滤波器的传递函数为：

3.s域频率变换

原型低通滤波器可以通过s域的频率变换得到模拟低通、高通、带通和带阻滤波器，本文主要考虑低通到低通的变换，只需做如下变换：

替换结果为

4.双线性变换求数字滤波器系数和

双线性变换将模拟滤波器的传递函数转换到数字滤波器的传递函数，相对冲击响应不变法，可以避免频率混叠，双线性变化法将s域的虚轴映射到z域的单位圆上，将s域中左半平面的极点映射到单位圆内，是一种保角映射。具体变换方法为用z的分式替换s：

$s=\frac{2}{T}\cdot \frac{1-z^{-1}}{1+z^{-1}}$

变换结果整理如下：

C语言编程实现步骤如下：

初始化、和中间缓存数组为0；
计算展开式的各项系数；
合并同类项的系数；
最高次项系数归一化。

5.由差分方程递推计算滤波结果

一个无限冲激响应系统的差分方程所示，实现结构有直接型、并联型、级联型，其中最简单直观的是直接型，直接型分为直接Ⅰ型和直接Ⅱ型，用流图表示方法如下：

从图3和图4可以看出直接Ⅱ型的滤波器结构相比直接Ⅰ型减少一半的延时单元，程序实现上所需的存储空间也减少一半。基于C语言的编程步骤如下：

$y\left ( i \right )=\sum_{k=0}^{N}b_{k}w_{k}$

4.滤波结果

一个一个输入待滤波数据；
初始化所有中间缓存状态变量为0；
由数组计算当前状态变量：
利用状态变量和数组计算滤波后的输出：
刷新中间状态变量：

取一段含有噪声的信号，对其进行Bessel滤波，设置滤波器阶数为4，最大平坦相位频率500Hz，滤波结果在MATLAB中图形化如下图5所示，可以看出滤波后的信号基本没有高频噪声而且曲线光滑，只是滤波后信号相对滤波前信号有一定的时延，这是Bessel滤波器的固有属性。设计的滤波器对不对，本文通过与MATLAB中自带的bessel滤波器函数对信号的滤波结果进行对比，可以看出MATLAB的滤波结果和本文所设计的Bessel滤波器的滤波结果完全吻合，如图6所示，说明本文所设计的滤波器是合理的。