步态数据去噪滤波器之巴特沃斯滤波器

对信号进行去渲染、水平倾斜校正和低通滤波（巴特沃斯滤波器，4阶；截止频率10Hz）

Butterworth filter

定义：巴特沃斯滤波器是电子滤波器的一种。巴特沃斯滤波器的特点是通频带的频率响应曲线最平滑。这种滤波器最先由英国工程师斯替芬·巴特沃斯（StephenButterworth）在1930年发表在英国《无线电工程》期刊的一篇论文中提出的。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。

巴特沃斯滤波器特性

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波特图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少，趋向负无穷大。

一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频。单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。

巴特沃斯滤波器原理

巴特沃斯型滤波器在现代设计方法设计的滤波器中，是最为有名的滤波器，由于它设计简单，性能方面又没有明显的缺点，又因它对构成滤波器的元件Q值较低，因而易于制作且达到设计性能，因而得到了广泛应用。其中，巴特沃斯滤波器的特点是通频带的频率响应曲线最平滑。

滤波器的截止频率的变换是通过先求出待设计滤波器的截止频率与基准滤波器的截止频率的比值M，再用这个M去除滤波器中的所有元件值来实现的，其计算公式如下：M=待设计滤波器的截止频率/基准滤波器的截止频率。

滤波器的特征阻抗的变换是通过先求出待设计滤波器的特征阻抗与基准滤波器的特征阻抗的比值K，再用这个K去乘基准滤波器中的所有电感元件值和用这个K去除基准滤波器中的所有电容元件值来实现的。

巴特沃斯滤波器与其它滤波器比较

下图是巴特沃斯滤波器（左上）和同阶第一类切比雪夫滤波器（右上）、第二类切比雪夫滤波器（左下）、椭圆函数滤波器（右下）的频率响应图。

巴特沃斯滤波器优点

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波得图上，从某一边界角频率开始，振幅随着角频率的增加而逐渐减少，趋向负无穷大。

一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝，三阶巴特沃斯滤波器的衰减率为每倍频18分贝，如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数、振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。

巴特沃斯滤波器是滤波器的一种设计分类，类同于切比雪夫滤波器，它有高通，低通，带通，高通，带阻等多种滤波器。它在通频带内外都有平稳的幅频特性，但有较长的过渡带，在过渡带上很容易造成失真，我在调用MATLAB里的巴特沃斯滤波器做仿真时，信号总会在第一个周期略微有些失真，但往后的幅频特性就非常的好。

巴特沃斯滤波器主要参数介绍

（1）［N，wc］=buttord（wp，ws，RP，As，’s’）

　　该格式用于计算巴特沃斯模拟滤波器的阶数N和3db截止频率wc。Wp、ws和wc是实际模拟角频率（rad\s）。Rp和As为通带最大衰减和最小衰减。

（2）［Z，P，k］=buttap（N）

　　该格式用于计算N阶巴特沃斯归一化模拟低通原型滤波器系统函数的零、极点和增益因子，返回长度为N的列向量Z和P，分别给出N个零点和极点的位置，K表示滤波器增益。

　　（3）Y=filter（b，a，x）

　　式中b表示系统传递函数的分子多项式的系数矩阵；a表示系统传递函数的分母多项式的系数矩阵；x表示输入序列；filter表示输出序列。IIR函数实现的直接形式。

（4）［b，a］=butter（N，wc，‘ftype’）

　　计算N阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a。

　　说明：调用参数N和wc分别为巴特沃斯数字滤波器的阶数和3dB截止频率的归一化值，一般是调用buttord格式（1）计算N和wc。系数b、a是按照z-1的升幂排列。

　　（5）［B，A］=butter（N，Ωc，‘ftype’，‘s’）

　　计算巴特沃斯模拟滤波器系统函数的分子、分母多项式系数向量。

　　说明：调用参数N和Ωc分别为巴特沃斯模拟滤波器的阶数和3dB截止频率（实际角频率），可调用buttord（2）格式计算N和Ωc。系数B、A按s的正降幂排列。

　　tfype为滤波器的类型：

　　◇ftype=high时，高通；Ωc只有1个值。

　　◇ftype=stop时，带阻；Ωc=［Ωcl，Ωcu］，分别为带阻滤波器的通带3dB下截止频率和上截止频率。

　　◇ftype缺省时：若Ωc只有1个值，则默认为低通；若Ωc有2个值，则默认为带通；其通带频率区间Ωcl《Ω《Ωcu。

　　（6）［H，w］=freqz（b，a，N）

　　b和a分别为离散系统的系统函数分子、分母多项式的系数向量，返回量H则包含了离散系统频响在0～pi范围内N个频率等分点的值（其中N为正整数），w则包含了范围内N个频率等分点。调用默认的N时，其值是512。可以先调用freqz（）函数计算系统的频率响应，然后利用abs（）和angle（）函数及plot（）函数，绘制出系统的频响曲线。

　　（7）lp2lp函数

　　［bt，at］=lp2lp（b，a，w0）

　　该函数用于实现由低通模拟原型滤波器至低通滤波器的频率变换，可以用传递函数和状态空间进行转换，但无论哪种形式，其输入必须是模拟滤波器原型。

　　（8）［bz，az］=impinvar（b，a，fs）

　　把具有［b，a］模拟滤波器传递函数模型转换为采样频率为fs的数字滤波器的传递函数模型［bz，az］，如果在函数中没有确定频率fs时，函数默认为1Hz.

　　低通巴特沃斯数字滤波器的实例

　　一、Matlab计算滤波器系数

　　Matlab计算巴特沃斯低通滤波器系数过程如下：

　　①根据给定的通带截止频率、通带截止增益、阻带截止频率、阻带截止增益，利用buttord函数计算巴特沃斯滤波器所需的最小阶数和截止频率。（由于要求中已给出截止频率，故这步省去）

　　（图是网上扒拉的，指标与本设计不符）

　　②根据上述计算得到的阶数，利用buttap函数计算出巴特沃斯滤波器原型。

　　③利用lp2lp函数，将原型滤波器转换成目标截止频率的滤波器。

　　④利用脉冲响应不变法（impinvar函数）或是双线性变换法（bilinear函数）将模拟滤波器转换为数字滤波器。数字滤波器形式为z域有理函数，分子分母系数即为滤波器系数。

　　我这里选用的是脉冲响应不变法，因为计算得到的滤波器比较简单，运算速度比较快。

　　（从左到右：滤波器原型、模拟滤波器、数字滤波器）

　　设计过程matlab源码如下：

　　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）））;

　　gridon;

　　subplot（2，1，2）

　　plot（W，180/pi*unwrap（angle（H）））;

　　gridon;

　　二、Matlab计算验证

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

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

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

　　则算法为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函数源码如下：

　　constfloatbz［2］={0，0.128580115806658};//分子

　　constfloataz［3］={1，-1.42252474659021，0.553007125840971};//分母

　　floatData_Output［DATA_LENTH］;//输出数据

　　float*but_filter（unsignedintlen，float*x）//len为输入数据数组长度，x为输入数据数组指针

　　{

　　unsignedinti=2;

　　staticfloatinit［2］={0，0};//初值，一开始设为0

　　if（len《2）//如果长度小于2，直接返回

　　returnData_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］;

　　returnData_Output;

　　}