数字滤波器的实现——低通滤波器再探究

在探究完滤波器原理之后，又面临一个问题就是数字滤波器如何实现的问题，因为在实际应用过程中，如果不接触硬件的话，低通滤波器一般都是通过编程实现的，具体代码应该怎么编写，在应用过程中又应该注意什么问题呢？这里值得探究一下。

首先我们回顾一下低通滤波器的原理：

从自动控制原理的角度，一阶低通滤波器实际上是一个一阶惯性环节，低通滤波器的传递函数表达形式如下所示，即为一阶惯性环节的传递函数表达式，截止频率为1000rad/s。

 那如果放在代码里面如何实现呢？

第一步：求取滤波器的差分方程，需要按照下述流程。S域传递函数——Z域离散函数——差分方程。

假设当前S域传递函数为：

通过后向差分方法，对其离散化，还有前向差分和双线性变换法，其中前向差分法不稳定，双线性变换精度更高，这里先探究一种；后向差分法中s算子的表达式为：

式中T为离散化过程中的单个步长时间，通常在单片机中是采样周期或为中断周期。可以得到离散域的传递函数：

​​​​​​由于在信号传递过程中：

即：

将其带入整合得到的输入输出关系式中，可以得到差分方程表达式：

若传递函数为文中开始的形式：

式中T为离散化过程中的单个步长时间。再次执行一下上面的推导过程；

经过推导后可以得出差分方程表达式为：

第二步：编程实现

 有了差分方程，下一步就是对其进行编程实现，插入一个function模块，代码段为：

function y = fcn(u)
%#codegen
persistent yk yk_1 Tsc wc
if isempty(yk)
    yk = 0;
    yk_1 = 0;
    Tsc =0.0001;
    wc =100;
end
 sum = (1+wc*Tsc);
    alpha = 1 / sum
    beta = wc*Tsc / sum
    yk = alpha*yk_1 + beta*u;
    yk_1 = yk;
y = yk;

如果按照文中开头给的传递函数，基波为幅值为100，角频率为10rad/s的正弦波，噪声是幅值为10，角频率为1000rad/s的正弦波。低通滤波器的截止频率100rad/s，采样时间为0.001s。对其进行仿真实验，下图第一行为原始波形，第二行为连续域传递函数滤波后效果，第三行为数字滤波器滤波后效果，可以看到，通过上述函数实现了对高频噪声的滤除，实现了指定截止频率的低通滤波器，代码具备实际应用价值。

附件：绘制截止频率为 100 rad/s的低通滤波器在连续域与离散域的bode图。

% transfor function
sys = tf([100], [1 100]);

% discretize
ts = 0.0001; % 采样周期
dsys = c2d(sys, ts, 'i'); % 转化为差分方程

% extract
[num, den] = tfdata(dsys,'v'); % 提取差分方程系数
sys =sys
dsys = dsys
opts = bodeoptions;
opts.FreqUnits = 'rad/s';
opts.XLim = [0.01, 10000];
opts.Grid = 'on';
 
bode(sys,dsys,opts);

 第三步：通用性低通滤波器的代码撰写

实际上在公司的代码里面，是不会直接对截止频率进行给定的，当然我认为直接给截止频率更方便一些，为了适应工作环境，下面附上按照滤波常数 Tf 和采样时间 Tsc 来设计低通滤波器的函数。

function y = fcn(u)
%#codegen
persistent yk yk_1 Tsc Tf
if isempty(yk)
    yk = 0;
    yk_1 = 0;
    Tsc =0.0001;
    Tf =0.01;
end
    alpha =  Tsc / (Tsc+Tf)
    beta = Tf / (Tsc+Tf)
    yk = (u * Tsc + yk_1*Tf )/(Tsc+Tf);
    yk_1 = yk;
y = yk;

 截止频率 wc  和 滤波常数 Tf 之间的比例关系如下所示：

 按照滤波时间常数 Tf 设计的滤波器结果如下图所示，输入和噪声与上文一致，从图中可以看出已实现了低通滤波器的功能。

任意低通滤波器传递函数转换为代码的一种方法：

% transfor function
sys = tf([100], [1 100]);
 
% discretize
ts = 0.0001; % 采样周期
dsys = c2d(sys, ts, 'i') % 转化为差分方程

转化为这个结果

基于这个结果推导一下公式：

 输入代码：例如SMO的反电动势滤波，输入为Zalpha，输出为Ealpha。由于代码是顺序执行，等号右侧的Ealpha实际上是上一次进入中断算出的Ealpha值。

Ealpha = 0.99Ealpha + 0.01Zalpha; 

为了便于显示当然也可以写成：这样在最后把Ealpha赋值给另一个变量，更加直观易读，但是也会多出一个变量定义。

Ealpha = 0.99Ealpha_old + 0.01Zalpha; 
.....//其他代码
Ealpha_old = Ealpha;