在探究完滤波器原理之后,又面临一个问题就是数字滤波器如何实现的问题,因为在实际应用过程中,如果不接触硬件的话,低通滤波器一般都是通过编程实现的,具体代码应该怎么编写,在应用过程中又应该注意什么问题呢?这里值得探究一下。
首先我们回顾一下低通滤波器的原理:
从自动控制原理的角度,一阶低通滤波器实际上是一个一阶惯性环节,低通滤波器的传递函数表达形式如下所示,即为一阶惯性环节的传递函数表达式,截止频率为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;