|
|
实现FIR滤波器这么简单，为什么各种数字信号处理教材堆砌那么多公式、引入那么多复杂概念，让人望而却步？
关注者
1,206
被浏览
434,800
77 个回答
Mr.看海
好好学习，天天向上
今天开始，本专栏将开通一个关于滤波的新专题。
数字滤波器中最基础的莫过于FIR和IIR这两个类型，而且这两类方法也有着系统的理论基础。相关课程中虽然对其有着非常详尽且专业的讲解，但是对于部分同学可能不那么容易理解。
本篇将从FIR入手，形象地讲解有限冲激响应数字滤波的理论知识，相信跟着阅读完，你一定能有所收获。
一、任务描述
如下图，我们仿真生成了两段信号。其中一条是未加入噪声的 纯净信号 ，另外一条是加入白噪声后的 含噪声信号 。其中纯净信号使用了一段正弦信号。我们的目的就是将 含噪声信号 中的噪声部分滤除，得到一条尽量还原 纯净信号 的滤波结果。
二、一种直观的解决方案
从上图容易看出，含噪声信号是在纯净信号上下波动的，所以很容易想到，如果取某几个含噪声信号的平均值，作为滤波结果，则可以抵消掉噪声在上下随机波动中的干扰效果。假如我们每3个点做一次平均，实现过程见下图（为了便于观察，使用前0.25秒数据演示）：
能够看出，经过该方法处理后的数据相对于纯净信号的波动幅值得到了一定程度的抑制，也就是可以取得一定程度的滤波效果。
这是一种最简单的FIR滤波器，也可以叫做滑动平均滤波方法。
三、另外一种理解角度——引入权重系数概念
上边讲到的滑动平均，是将3个相邻信号求取平均值，这也可以理解为这3个值分别乘以权重系数1/3再求和，那么上述滤波过程可以用下边这张动图来演示，应该是比较直观的了：
看到这里大家想到什么了吗。是的——是 卷积 。
四、重温卷积的概念
这里大家可能存在两个疑问：为什么上图能联系到卷积？卷积怎样影响到滤波算法的构造？
1.先解答第一个疑惑，为什么上图能联系到卷积。
在之前的文章里比较详细地讲解了卷积的概念： Mr.看海：这篇文章能让你明白卷积 ，大家在读完本篇文章后如果感兴趣可以再跳转过去看，现在我们继续讲。
在那篇文章里，我引用了一个形象的例子：
比如说老板命令张三干活，张三却到楼下打台球去了，后来被老板发现，他非常气愤，扇了张三一巴掌（注意，这就是输入信号，脉冲），于是张三的脸上会渐渐地（贱贱地）鼓起来一个包，张三的脸就是一个系统，而鼓起来的包就是张三的脸对巴掌的响应，好，这样就和信号系统建立起来意义对应的联系。下面还需要一些假设来保证论证的严谨：假定张三的脸是线性时不变系统，也就是说，无论什么时候老板打张三一巴掌，打在张三脸的同一位置（这似乎要求张三的脸足够光滑，如果张三长了很多青春痘，甚至整个脸皮处处连续处处不可导，那难度太大了，我就无话可说了哈哈），张三的脸上总是会在相同的时间间隔内鼓起来一个相同高度的包来，并且假定以鼓起来的包的大小作为系统输出。好了，那么，下面可以进入核心内容——卷积了！
如果张三每天都到地下去打台球，那么老板每天都要扇张三一巴掌，不过当老板打张三一巴掌后，一天后就消肿了，所以时间长了，张三甚至就适应这种生活了……如果有一天，老板忍无可忍，以0.5秒的间隔开始不间断的扇张三的过程，这样问题就来了，第一次扇张三鼓起来的包还没消肿，第二个巴掌就来了，张三脸上的包就可能鼓起来两倍高，老板不断扇张三，脉冲不断作用在张三脸上，效果不断叠加了，这样这些效果就可以求和了，结果就是张三脸上的包的高度随时间变化的一个函数了（注意理解）；如果老板再狠一点，频率越来越高，以至于都辨别不清时间间隔了，那么，求和就变成积分了。可以这样理解，在这个过程中的某一固定的时刻，张三的脸上的包的鼓起程度和什么有关呢？和之前每次打张三都有关！但是各次的贡献是不一样的，越早打的巴掌，贡献越小，所以这就是说，某一时刻的输出是之前很多次输入乘以各自的衰减系数之后的叠加而形成某一点的输出，然后再把不同时刻的输出点放在一起，形成一个函数，这就是卷积，卷积之后的函数就是张三脸上的包的大小随时间变化的函数。本来张三的包几分钟就可以消肿，可是如果连续打，几个小时也消不了肿了，这难道不是一种平滑过程么？反映到剑桥大学的公式上，f(a)就是第a个巴掌，g(x-a)就是第a个巴掌在x时刻的作用程度，乘起来再叠加就ok了，大家说是不是这个道理呢？我想这个例子已经非常形象了，你对卷积有了更加具体深刻的了解了吗？
转自GSDzone论坛
来源：人人网
看一下下图，红线的方波代表“打了一巴掌”，蓝线是肿起来的鼓包的高度。假设鼓包一天的时间（24小时）可以完全回复，也就是高度会变为0。
需要注意，“鼓包的高度”不是随着时间的推移无限趋近于0，而是在1天后就变成了0，也就是对于巴掌这个“冲激”的响应是“有限的”，这就是有限冲激响应。在这里我们假设“鼓包的高度”随时间的推移用g(t)表示，很明显g(0)=1，而g(1)=0。
如果进一步增加扇巴掌的频次，如果时间间隔小于24小时，那么第二次的鼓包就会叠加上上一次没有消除的鼓包高度，对应鼓包高度会像下图这样：
假如我们只求整时刻（即0.1的整数倍）的鼓包情况，为了更便于理解，干脆将“鼓包的高度”这个冲激响应也离散化，这并不影响我们对于整时刻的求解：
对于每一个时刻，其鼓包高度的计算是前几次为消除鼓包高度的叠加，对于0.3天这个时刻，此时鼓包的高度就是：
0.3天这一刻刚刚打了一巴掌，刚刚出炉的鼓包高度很明显是1，用公式表示就是g(0)
0.2天这一刻打的巴掌，距离0.3天已经过去0.1天，鼓包高度还残余g(0.1)
0.1天这一刻打的巴掌，距离0.3天已经过去0.2天，鼓包高度还残余g(0.2)
0天这一刻打的巴掌，距离0.3天已经过去0.3天，鼓包高度还残余g(0.3)
所以第0.3天这一刻的鼓包高度就是 g(0)+g(0.1)+g(0.2)+g(0.3)
现在再贴近现实一些，老板也很难控制每次巴掌力度的一致，如果他用了两倍的力气，那鼓包的高度很容易理解就是2了（因为是线性时不变系统），这里我们引入一个变量f(t)，表示t时刻老板的力度。那么对于0.3天这个时刻，此时鼓包的高度就变成：
0.3天这一刻刚刚打了一巴掌，此时力气值是f(0.3)，用鼓包高度就是f(0.3)g(0)
0.2天这一刻打的巴掌，此时力气值是f(0.2)，距离0.3天已经过去0.1天，鼓包高度还残余f(0.2)g(0.1)
0.1天这一刻打的巴掌，此时力气值是f(0.1)，距离0.3天已经过去0.2天，鼓包高度还残余f(0.1)g(0.2)
0天这一刻打的巴掌，此时力气值是f(0)，距离0.3天已经过去0.3天，鼓包高度还残余f(0)g(0.3)
所以第0.3天这一刻的鼓包高度就是 f(0.3)g(0)+f(0.2)g(0.1)+f(0.1)g(0.2)+f(0)g(0.3)
如果老板2天之内使用了如下的“连环掌”：
那张三脸上的鼓包将会是怎样呢。
请开始表演：
看完上图有没有觉得有点眼熟，是的！和图3简直一毛一样！
所以综上可以得到一个直观的结论： 当一组数像滑窗一样滑过另外一组数时，将对应的数据相乘并求和得到一组新的数，这个过程必然和卷积有着莫大的关系。
在后边可能会讲到的卷积神经网络中，只是将一维数据扩展到了二维平面，但这种滑窗求和的基本思想是一脉相承的，乃至于该卷积思想可以推广到多维空间当中。就像下图：
不过这里我们不做过多扩展了。
我们只需要知道： FIR滤波就是在时域上卷积的过程。
2.再解答第二个疑惑，上述卷积怎样影响到滤波算法的构造。
卷积有一个重要性质： 时域的卷积等于频域相乘。
这个性质可以说和滤波是息息相关的，因为 滤波的目的就是想要得到某特定的频率段。
而频域相乘恰恰能够实现这个筛选的目的。
回到文章最开始的例子，想要对含噪声信号滤波，那我们不妨先看看他的频谱：
可以看到我们想要保留的频率段范围基本在3Hz以下，高于3Hz的都可以抛弃，很自然地，我们可以想到一个办法，即构造如下一组频域信号：
此时将图11中的频谱和图12相乘，高于3Hz的数值就全部变为了0，此时频谱就成了：
再将此频域数值进行傅里叶逆变换，得到时域数值：
可以看到我们近乎完美的实现了滤波。一句话表述上述滤波过程：设计一个频域滤波器（将想要保留的频率段赋值为1，其他频率段赋值为0），将其与含噪声信号的频谱 在频域上相乘 ，再将乘积做傅里叶逆变换，即可实现滤波，这种滤波器叫 频域滤波器。
虽然得到了较好的滤波结果，但是事情到这里还没完，不要忘记这篇文章是讲FIR滤波器的~
注意看上边那段话中的加黑字体： 在频域上相乘 。回想起来了吗？
时域的卷积等于频域相乘！
所以，上述在频域空间的一大堆操作，可以简单地转化成在频域上的卷积操作，具体来说，就是 将含噪声信号与低通滤波器的傅里叶逆变换值进行卷积，这个过程就是FIR滤波。
我用一张图表示他们之间的关系：
到这里你应该已经理解了， FIR滤波器的本质是设计一组系数，这组系数实际就是滤波器的IFFT（冲激响应，还记得张三的鼓包吗）离散化以后的结果。
至此，FIR到底是怎样一个概念应该明白了吧。
五、MATLAB中滤波器的设计方法
MATLAB为滤波设置了种工具，比如图形化的设计工具filterDesginer、根据差分方程直接设计滤波器的filter函数、根据滤波目的进行设计的lowpass函数、highpass函数等等。
不过我最终选用了 designfilt 函数进行FIR、IIR滤波器设计实现。因为该函数兼具了方法的 全面性 和 统一性 ，而且相对于图形界面，纯代码的形式在很多场景下也更便于调用。
5.1 使用designfilt函数进行滤波
这里我们举一个例子来说明designfilt函数的用法。
准备滤波的对象数据如下，我们对该数据进行低通滤波：
1.导入数据。
我们将该数据命名为data，在程序中使用如下代码进行导入：
%% 1.导入数据
load data.mat此步骤也可以是导入其他格式的数据文件，比如Excel、csv、txt、一些音频格式等等；也可以不导入数据而是使用仿真数据，那么此处就对应改成数据仿真的计算公式。
2.参数设置
下面我们设计一个FIR低通滤波器，所以type设置成了“lowpassfir”，此外根据滤波器类型和设计方法的不同，可以选择：'lowpassfir' 、 'lowpassiir' 、'highpassfir'、 'highpassiir'、 'bandpassfir' 、'bandpassiir'、'bandstopfir'、'bandstopiir' 。
在这个例子里边我们需要指定四个参数，分别是滤波器阶数、低通截止频率、数据采样频率、设计方法。
type = 'lowpassfir';
order = 50; %滤波器阶数
CutoffFrequency = 75; %低通截止频率
Fs = 500; %数据采样频率
DesignMethod = 'window' %设计方法3.实现滤波并画图
下边这行代码是将上述参数导入designfilt函数，得到设计完成的滤波器d。
d =designfilt(type,'FilterOrder',order,...'CutoffFrequency',Fst,'SampleRate',Fs，'DesignMethod',DesignMethod);然后调用filter函数，实现滤波，filter入口参数为滤波器d和待滤波数据data。
xf = filter(d,data); %xf为滤波后数据此时绘制一下滤波前数据和滤波后数据的对比图：
可以看到滤波结果是有相位延迟的，对于相位延迟，在 之前的文章 中也介绍过。由于FIR滤波的相位延迟是线性的，可以通过平移进行延迟补偿(不用担心，IIR滤波器也可以相位补偿)：
4. 关于designfilt函数的其他规范集
这里要提到一个“规范集”的概念，可以看一下下图：
上图是lowpassfir的规范集，对于其他七种方法，都有这样一张规范集列表。所谓的“规范集”简单来说就是designfilt函数支持不同的滤波器设计方法（这里的设计方法不是指FIR或者IIR，而是更加细分的“window、butter”等等），而在不同的设计方法中，需要输入的各种参数都会有所不同。而上表就是列出在不同的设计方法下需要设置哪些参数。
我们在上边程序中用到的规范集是下边红框框出来的这一行：
如果在你的研究需求中，要使用其他规范集，就需要参考到这张表了。
5.2 代码的傻瓜化封装——一行代码实现FIR/IIR滤波
通过上述方法结合官方教程文档，我们可以实现滤波器的设计，但是对于新手来说多少还是有一些门槛的，而且对于很多应用场景，我们并不需要考虑那么多种规范集形式。
有没有一种更 简便易用 ，兼顾 专业性 ，最好还可以 多功能分析 和 画图 的实现方法呢？
现在有了——笔者按照专栏以往的的代码风格，对FIR、IIR算法进行了二次封装，实现了“一行代码”完成滤波的效果。我将封装后的函数命名为 funFirIirFliter 。
比如还是对于上边的例子，只需要像下边这样设置参数并调用funFirIirFliter函数：
resp = 'lowpassfir'; %选择滤波器类型
option.order = 50; %设置滤波器阶数
option.CutoffFrequency = 75; %设置滤波器截止频率
option.Fs = 500; %设置信号采样频率
compFlag = 'on'; %compFlag：是否进行延迟补偿，'on'为开启补偿，'off'为关闭补偿。
%% 3.调用funFirIirFliter实现滤波
dataOut = funFirIirFliter(data,resp,option,compFlag);就可以得到如下结果：
图1是设计得到的滤波器响应图；图2是滤波效果图，其中包括滤波前信号及频谱、滤波后信号及频谱、滤波前后信号对比图（可以用于比较相位滞后情况）。
也就是说只要向funFirIirFliter函数导入待分析的数据、选择滤波器类型（如'firlowpass'等）、进行option相关参数设置（对于不同的滤波器类型设置大概3-4个参数），还可以选择是否进行相位补偿（通过设置compFlag参数为'on'或'off'）。然后调用一行“funFirIirFliter(data,resp,option,compFlag);”，就可以实现上述所有功能了。
对于funFirIirFliter函数更详细的注释说明如下：
function [dataOut,d] = funFirIirFliter(data,resp,option,compFlag)
% 进行FIR或者IIR滤波，可以实现滤波器设计、滤波、画图
% 如需对下述代码进行二次开发，可参考此官方帮助文档：https://ww2.mathworks.cn/help/signal/ref/designfilt.html?s_tid=doc_ta
% 输入：
% data：待滤波数据
% resp：滤波器响应类型和滤波类型，具体包含：
% -'lowpassfir' FIR低通滤波器，采用'window'设计方法对应的规范集，需要设置的参数包括：
% - option.order 滤波器阶数
% - option.CutoffFrequency 滤波器截止频率
% - option.Fs 采样频率
% -'lowpassiir' IIR低通滤波器，采用最小阶设计对应的规范集
% - option.Fs 采样频率
% - option.PassbandFrequency 通带频率
% - option.StopbandFrequency 阻带频率
% - option.DesignMethod 设计方法，选项： 'butter'、'cheby1'、'cheby2'、'ellip'
% -'highpassfir' FIR高通滤波器，采用'window'设计方法对应的规范集，需要设置的参数包括：
% - option.order 滤波器阶数
% - option.CutoffFrequency 滤波器截止频率
% - option.Fs 采样频率
% -'highpassiir' IIR高通滤波器，采用最小阶设计对应的规范集
% - option.Fs 采样频率
% - option.PassbandFrequency 通带频率
% - option.StopbandFrequency 阻带频率
% - option.DesignMethod 设计方法，选项： 'butter'、'cheby1'、'cheby2'、'ellip'
% -'bandpassfir' FIR带通滤波器，采用'window'设计方法对应的规范集，需要设置的参数包括：
% - option.order 滤波器阶数
% - option.CutoffFrequency1,... %通带起始频率
% - option.CutoffFrequency2,... %通带结束频率
% - option.Fs 采样频率
% -'bandpassiir' IIR带通滤波器，采用最小阶设计对应的规范集
% - option.Fs 采样频率
% - option.StopbandFrequency1, ... % 阻带起始频率
% - option.PassbandFrequency1, ... % 通带起始频率
% - option.PassbandFrequency2, ... % 阻带截止频率
% - option.StopbandFrequency2, ... % 通带截止频率
% - option.DesignMethod 设计方法，选项： 'butter'、'cheby1'、'cheby2'、'ellip'
% -'bandstopfir' FIR带阻滤波器，采用'window'设计方法对应的规范集，需要设置的参数包括：
% - option.order 滤波器阶数
% - option.CutoffFrequency1,... %通带起始频率
% - option.CutoffFrequency2,... %通带结束频率
% - option.Fs 采样频率
% -'bandstopiir' IIR带阻滤波器
% - option.order 滤波器阶数
% - option.CutoffFrequency1,... %通带起始频率
% - option.CutoffFrequency2,... %通带结束频率
% - option.Fs 采样频率