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

87 阅读 0 评论 58 点赞

我是靠谱客的博主寒冷凉面，最近开发中收集的这篇文章主要介绍数字滤波器的实现——低通滤波器再探究，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

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

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

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

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

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

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

通过前向差分方法，对其离散化：（还有后向差分和双线性变换法，这里仅探究一种）

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

对其进行化简与整合可以得到输入和输出的关系式：

由于在信号传递过程中：

即：

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

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

式中Tf为一阶低通滤波器的截止频率。经过推导后可以得出差分方程表达式为：

第二步：编程实现

有了差分方程，下一步就是对其进行编程实现，插入一个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 之间的比例关系如下所示：