低通滤波器的设计与DSP实现

前言

低通滤波器是很常见的工具。比如，在电机控制DSP的AD采样之前，会有大约600Hz的硬件滤波器，用于滤除高频噪声，一般采用两级运放的有源滤波器。再比如，在EQEP的输入信号端口，会有阻容硬件滤波器，用于消除长距信号线可能存在的噪声。在电机控制器内部一般倒是不太会用到滤波器，因为加入低通滤波器有可能造成相位延时，导致控制系统的性能改变。不过，在一些新能源应用和在线监测场景下，还是会用到这个工具，毕竟虽然理论性比较强，但是它实现起来很简单，而且效果明显。

设计目标

这篇文章采用例子的形式，来实现一个数字滤波器的设计以及在线应用。设计的过程利用matlab实现，最终推导得到结果后，可以在DSP中实现。参考郑君里的《信号与系统》，现在的数字滤波器都是基于模拟滤波器得到的，而模拟滤波器中最常用的就是巴特沃斯滤波器和切比雪夫I型滤波器。这里的应用场景下，我选择巴特沃斯低通滤波器。

在采样频率2000Hz的情况下，设计一个截止频率为5Hz的2阶数字低通滤波器。这里截止频率即：

The cutoff frequency is the frequency at which the magnitude response of the filter is 1 / sqrt(2) = 0.7071

实现过程

首先，设计一个模拟滤波器，在matlab中可以使用butter这个函数，这个函数应用很多，这里我采用三个输入参数的形式。

[z,p,k] = butter(2,5*2*pi,'s');

其中，第1个参数规定了滤波器的阶数；第2个参数在模拟滤波器设计过程中，是截止频率的弧度形式，即2π倍的频率；第3个参数设为’s’，对应模拟滤波器设计。同时，输出也有很多，这里使用传递函数分子分母的形式，得到如下结果。因为我们关注工程实现，这里不去追问这个结果是如何得到，但是《信号与系统》的书里有的。大概意思是根据截止频率和巴特沃斯滤波器的基本结构进行求解。

下面，要将s函数转化为z函数，一般有两种方法：冲激响应不变法和双线性变换法。由于第1种方法会有高频混叠，实际上最常用的是双线性变换法。利用matlab中的bilinear函数，可以直接得到，免去了代入计算的麻烦，得到如下结果。

然后，只要把上式变成差分方程就可以在DSP上用了，这部分[1]给了手推公式，这里针对我们的应用，整理一下。

最后，只要实现一下上面的公式，就可以实现在线的低通滤波了。为了方便，这里在m文件中模拟一个在线的过程，构造了一个带噪声的2+0.4cos(2π50t)+0.3cos(2π100t)信号，结果如下。同时还给了离散函数的幅频特性，没错，我不喜欢dB形式的图，我喜欢线性，哈哈哈。

阶数一般化

在上节中，为了将z函数变化为差分方程，我们推导了2阶的形式，那么3阶呢？高阶呢？从下面可以看到其中的规律，我相信matlab里应该有对应的一般化函数，但是这里我决定自己完善一下上面的例子，于是得到了更为一般化的结果。

根据扩展，我们可以得到如下所示的滤波结果。可以看到，随着阶数增加，滤波延时和过冲会增加。但是，并不是意味着阶数可以一直增加，在本例中，8阶情况下，系统会出现不收敛的情况，大概是离散精度之类的问题，这个属于系统系的研究范畴了。

小结

对了，补充一点，这里实现的滤波器，应当是属于IIR，即无限冲击响应滤波器范畴。至于FIR，我并没有看懂[doge]。

参考与代码

[1] https://www.cnblogs.com/xklzw/p/13439873.html

clearvars; clc; close all;

% 设计指标
Fs = 2000; order = 3; Wn = 5*2*pi;

% 构造目标信号
Ls = 2000;
t = (0:1/Fs:1/Fs*(Ls-1))';
signal = 2+0.4*cos(2*pi*50*t)+0.3*cos(2*pi*100*t)+rand(Ls,1);

% designs a lowpass analog Butterworth filter 
[num,den] = butter(order,Wn,'s');
% 模拟传递函数
H = tf(num,den);

% 双线性变换
[numd,dend] = bilinear(num,den,Fs);
% 数字传递函数
Hd = tf(numd,dend,1/Fs);

% 滤波器测试 order阶
b = numd; a = dend;
yd = zeros(1,order); ud = zeros(1,order+1);
for k=1:1:Ls                     
    ud = [signal(k),ud(1:order)];             
    yd = [1/a(1)*sum(b.*ud) - 1/a(1)*sum(a(2:order+1).*yd) , yd(1:order-1)];               
    filt_signal(k) = yd(1);
end

figure(1);
plot(t,signal);hold on; grid on;xlabel('Time (s)'); ylabel('Magnitude'); title('Butterworth Low Pass Filter')
plot(t,filt_signal,'linewidth',2);hold on; grid on;

% figure(2);
% P=bodeoptions;
% % P.FreqUnits='rad/s';
% P.FreqUnits='Hz';
% P.FreqScale='linear';
% % P.MagScale='log';
% P.MagUnits='abs'; 
% bode(Hd,P);grid on;