带通滤波器

this is the sixth passage of HEXO
SEVO 实验室

带通滤波器

巴特沃斯滤波器

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少，趋向负无穷大。
一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。

设计步骤

如设计一个数字低通滤波器，其技术指标为：
通带临界频率fp ，通带内衰减小于rp；
阻带临界频率fs，阻带内衰减大于αs；采样频率为FS
1、将指标变为角频率 wp=fp2pi;ws= fs2pi;
2、将数字滤波器的频率指标{Wk}由wk=（2/T）tan（Wk/2）转换为模拟滤波器的频率指标{wk}，由于是用双线性不变法设计，故先采取预畸变。

3、将高通指标转换为低通指标，进而设计高通的s域模型
4、归一化处理

由以上三式计算出N，查表可得模拟低通滤波器的阶数，从而由下式确定模拟高通滤波器的参数。

仿真程序的设计与调试

数字域指标变换成模拟域指标
其程序为：

1
2
3
4
5
fp = 400 fs= 300;
Rp = 1; Rs = 20;
wp =fp*2*pi;
ws =fs*2*pi;
FS=1000;T=1/FS；

程序执行结果为：wp=2.5133e+003 ws=1.8850e+003 与实际计算结果相符。

数字域频率进行预畸变

其程序为:

1
2
3
4
wp2=2*tan（Wp/2）/T;
ws2=2*tan（Ws/2）/T;
经过预畸变，可以发现频率变为： wp2= 6.1554e+003
ws2= 2.7528e+003

模拟滤波器的设计

其程序为:

1
2
3
4
5
6
7
8
9
10
%设计模拟滤波器
［N，Wn］ = buttord（wp2，ws2，Rp，Rs，‘s’）
［z，p，k］=buttap（N）; %创建Buttord低通滤波器原型
［Bap，Aap］=zp2tf（z，p，k）; %由零极点转换为传递函数的形式
figure（1） freqs（Bap，Aap）; %模拟低通滤波器的频率响应
TItle（‘模拟滤波器（低通原型）的频率响应’）
［Bbs，Abs］=lp2hp（Bap，Aap，Wn）; %模拟低通变高通
figure（2）
freqs（Bbs，Abs）;
TItle（‘模拟滤波器的频率响应’）

程序执行后可以发现其频率响应为： N=4，其波形如下图

模拟滤波器的频率响应
由上图分析可得：其符合高通的一般特征，与预期的效果一样。 而在此条件下，Butterworth滤波器低通原型的波形如下图。

模拟滤波器（低通原型）的频率响应
在设计的过程中，涉及一个频率变换的问题，即将模拟低通原型变为高通，其函数及用法如下：
［b，a］=lp2hp（Bap，Aap，Wn）;
功能：把模拟滤波器原型转换成截至频率为 Wn 的高通滤波器。 其中，Bap，Aap分别为低通传递函数的分子向量和分母向量；
b，a分别为高通传递函数的分子向量和分母向量。

模拟滤波器变成数字滤波器

其程序为：

1
［Bbz，Abz］=bilinear（Bbs，Abs，FS）; %用双线性变换法设计数字滤波器 freqz（Bbz，Abz，512，FS）;

程序运行的结果为：如下图

数字滤波器的频率响应
由于使用的是双线性不变法设计的，其相位为非线性。此处主要是基于要获得严格的频率响应，以及较准确地控制截止频率的位置，故画出了详细的幅频响应。（如下图）

详细的幅频响应
分析该图可知其在0.6（即300Hz）处的衰减为40dB，而在0.8（即400Hz）处的衰减极小，应小于1dB。由此可见，此设计符合要求设计的参数。
而在调试的过程中发现：通带衰减越小，滤波器的性能越好 阻带衰减越大，滤波器的性能越好 其曲线也越陡峭，选择性越好，当然所用的滤波器阶数也越高。
当阻带衰减变为40dB（之前为20dB），通带不变时，其波形如下图。对比上图可知，其在阻带临界频率处衰减变为了40dB，曲线变陡峭了。

详细的幅频响应（阻带衰减为40dB）
当通带变为5dB时，阻带不变时，其波形如下图。对比图3-3可知，其在通带处的衰减变为了5dB，曲线平滑了一些。

详细的幅频响应（通带衰减为5dB）

理论计算数字滤波器的仿真
1
2
3
4
5
6
7
8
9
10
11
wp=0.8*pi;
ws=0.6*pi;
OmegaP=2*1000*tan（wp/2）;
OmegaS=2*1000*tan（ws/2）;
lamdas=OmegaP/OmegaS;
N=0.5*log10（（10.^（20/10）-1）/（10.^（1/10）-1））/log10（lamdas）;
%笔算的结果为N=3.6947；故取N=4 %
此处为计算高通的传递函数 Wn= 4.8890e+003 az=［0 0 0 0 1］;
bz=［1 2.613 3.414，2.613，1］; ［Bbs，Abs］=lp2hp（az，bz，Wn） %用双线性不变法处理
［Bbz，Abz］=bilinear（Bbs，Abs，1000）; 其运行结果为：N=3.6947；图形如下图
![](/images/巴特沃斯滤波器12.png)
理论计算的滤波器的幅频响应

综上所述，本滤波器以四阶即实现了预期的设计目标：采样频率为1000Hz，通带临界频率fp =400Hz，通带内衰减小于1dB（αp=1）；阻带临界频率fs=300Hz，阻带内衰减大于20dB（αs=25），其在通带内的性能更好。

MATLAB应用

滤波器设计目标：设计一个1Hz截止频率的2阶低通巴特沃斯数字滤波器，并转化成C语言函数。（国标里提的要求）
滤波器指标：指标：截止频率Fc = 1Hz，阶数N=2，低通巴特沃斯滤波器，采样频率Fs = 15Hz。

一、Matlab计算滤波器系数

Matlab计算巴特沃斯低通滤波器系数过程如下：
①根据给定的通带截止频率、通带截止增益、阻带截止频率、阻带截止增益，利用buttord函数计算巴特沃斯滤波器所需的最小阶数和截止频率。
②根据上述计算得到的阶数，利用buttap函数计算出巴特沃斯滤波器原型。
③利用lp2lp函数，将原型滤波器转换成目标截止频率的滤波器。
④利用脉冲响应不变法（impinvar函数）或是双线性变换法（bilinear函数）将模拟滤波器转换为数字滤波器。数字滤波器形式为z域有理函数，分子分母系数即为滤波器系数。
我这里选用的是脉冲响应不变法，因为计算得到的滤波器比较简单，运算速度比较快。
（从左到右：滤波器原型、模拟滤波器、数字滤波器）
设计过程matlab源码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Fs = 15;        %采样频率
Nn = 12800;
N = 2;          %阶数
Wc = 1*2*pi;    %截止频率
[z,p,k] = buttap(N);  %计算巴特沃斯滤波器原型
[Bap,Aap] = zp2tf(z,p,k);  %转换成多项式模式
[b,a] = lp2lp(Bap,Aap,Wc);  %根据截止频率计算模拟巴特沃斯滤波器系数
[bz,az] = impinvar(b,a,Fs);  %用脉冲响应不变法离散化
figure(1)
[H,W] = freqz(bz,az,Nn,Fs);  %绘制频率特性曲线
subplot(2,1,1)
plot(W,20*log10(abs(H)));
grid on;
subplot(2,1,2)
plot(W,180/pi*unwrap(angle(H)));
grid on;

二、Matlab计算验证
先在Matlab中验证滤波函数。先编写带噪声的输入函数，然后经过滤波器函数后，观察滤波效果。其中滤波器函数写法为：

Filter函数为Matlab自带函数，其算法为：

其中，a即为z域传递函数的分母系数，b为分子系数。例如本应用中：

则算法为

1
az(1)*y(k) = bz(1)*x(k) + bz(2)*x(k-1) – az(2)*y(k-1) – az(3)*y(k-2)

Matlab中得到的结果如下（信号频率0.1Hz，噪声频率6Hz）：

三、C语言函数编写与验证
将上述算法翻译成C语言，写入单片机中。利用信号源输出各种波形，单片机AD采样进去之后，对采样点进行滤波处理，将原始数据和滤波后的数据发送到上位机进行绘图，得到图像对比如下：

C函数源码如下：

1
2
3
4
5
6
7
const float bz[2] = {0,0.128580115806658};  //分子
const float az[3] = {1,-1.42252474659021,0.553007125840971};  //分母
float Data_Output[DATA_LENTH];  //输出数据
float* but_filter(unsigned int len, float* x)  //len为输入数据数组长度，x为输入数据数组指针
{unsigned int i = 2;static float init[2] = {0,0};  //初值，一开始设为0
  if(len<2)                   //如果长度小于2，直接返回return Data_Output;Data_Output[0] = init[0];   //赋初值Data_Output[1] = init[1];for(i = 2;i < len;i++){Data_Output[i] = bz[0]*x[i] + bz[1]*x[i-1] - az[1]*Data_Output[i-1] - az[2]*Data_Output[i-2];/*算法为a1*y(k) = b1*x(k) + b2*x(k-1) - a(2)*y(k-1) - a(3)*y(k-2)*//*由于a1 = 1，故不做除法*/}init[0] = Data_Output[len-2];  //考虑到会被连续调用，此次的终值作为下次的初值init[1] = Data_Output[len-1];return Data_Output; 
}