一个设计低通巴特沃斯数字滤波器的实例

本人本科渣渣一个，前两天导师让我设计一个数字滤波器。由于本人基本没有数字信号处理基础，于是只能依靠百度和matlab，折腾了半天总算是摸索明白了。百度上有一些文章不靠谱，很容易误导别人，故在此发一篇博客。

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

滤波器指标：指标：截止频率Fc = 1Hz，阶数N=2，低通巴特沃斯滤波器，采样频率Fs = 15Hz。

一、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)));

grid on;

subplot(2,1,2)

plot(W,180/pi*unwrap(angle(H)));

grid on;

二、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函数源码如下：

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;

}