滤波器的使用及算例

“前一篇文章我们讲解了离散傅立叶变换的公式、推导及应用方法，本文我们将基于离散傅立叶变换来进行滤波器的讲解，并举例说明频域滤波和时域滤波的异同”

01

—

频域滤波：双边谱计算

图1第一幅图是横坐标时间，纵坐标加速度的时域离散信号（工程中采集到的最原始的信号）。

该信号由直流分量为6（频率为0），频率40Hz幅值为8的余弦信号，频率100Hz幅值为6的余弦信号，频率400Hz幅值为4的余弦信号组成。

该信号的采样频率为1024Hz，采样点数为512个点，那么诸如时间分辨率等其他参数均可如图1中的计算结果。

图1

图1第二幅图是第一幅图离散傅立叶变换除以N后的结果，从图中可以看出：0频幅值是6，40Hz幅值是4（和其以Fs/2为中心对称的984Hz幅值是4）；100Hz幅值是3（和其以Fs/2为中心对称的924Hz幅值是3）；400Hz幅值是2（和其以Fs/2为中心对称的624Hz幅值是2）；

该图未做后续处理，以中心频率Fs/2为轴成对称，称为双边谱（实际信号中不含984Hz，924Hz，624Hz频率成分）。

所谓滤波器，即设计一个方法，滤掉不要的频率成分，保留感兴趣的频率成分。图2是设计一个滤波器，滤掉图1时域信号中400Hz的频率成分。

图2第二幅图是4阶300Hz低通巴特沃斯滤波器的幅频传递特性：

1）低通滤波，300Hz以前传递率约为1，即0dB；

2）300Hz传递率为0.707，即衰减-3dB，300Hz以后大幅衰减。

需要注意的是图2中的频谱图均为N个点。

图2

将图2中第一幅图（复数）乘以图2中第二幅图（复数）后，得到图3中第一幅图（复数），即滤波后的双边谱（图3第一幅图），然后对其进行离散反傅立叶变换，接着乘以N，得到的结果只取实部，即得到滤波后的时域信号（图3第二幅图）。

图3

从图3第二幅图可以看出，滤掉400Hz频率成分后，信号显得更加光滑。至此信号滤波完成。

总结下来的步骤为：

1）时域信号傅立叶变换得到双边谱（复数）；

2）双边谱乘以滤波器频响函数（复数，全频段即到Fs）；

3）复数相乘的结果进行反傅立叶变换，取实部。

*注意，如果步骤1）除以N，那么在步骤3）处应该乘以N。

02

—

频域滤波：单边谱计算

图1第二幅图双边谱我们在信号分析中不常见。常见到的是只显示到最高分析频率为Fs/2的单边谱（符合采样定理：采样频率Fs要大于最高分析频率的2倍以上），如图4第二幅图。

图4

得到单边谱的步骤如下：

1）时域信号傅立叶变换后除以N；

2）第0个数，不变；

3）从第1个数到第N/2-1个数均乘以2；

4）只显示第0个数到第N/2-1个数（共N/2个数）。

按照以上步骤处理后的单边谱结果才能准确反映各频率的幅值及相位成分（乘以2并不改变相位）。

单边谱是N/2个点，相应滤波器的频率点也是N/2个点，如图5

图5

将图5中第一幅图乘以图5中第二幅图后，并且将N/2到N-1个点都赋值为0后，得到图6中第一幅图（因为频谱是多少个点，反傅立叶变换后仍然是多少个点，为了保证时域信号仍然是N个点，所有对频谱信号的后半部分采取赋值为0 的方法）。然后对图6第一幅图进行离散反傅立叶变换，然后乘以N，得到的结果取实部，即得到滤波后的信号，如图6第二幅图。

图6

总结下来的步骤为：

1）时域信号傅立叶变换得到双边谱（复数）；

2）将双边谱转换为单边谱（复数）;

2.1) 第0个数保持不变；

2.2) 第1个～N/2-1个数乘以2；且只取0～N/2-1个数（即前N/2个数）

3）单边谱乘以滤波器频响函数（复数，全频段到Fs/2）；

4）复数相乘的结果从N/2～N-1个数赋值为0；

5）对这N个复数进行反傅立叶变换，取实部。

*注意，如果步骤1）除以N，那么在步骤5）处应该乘以N。

03

—

滤波器的构造

滤波器的构造，实际上是设计一个频谱的传递函数（不管是频域滤波还是时域滤波都是该传递函数），下面的公式是一个4阶滤波器的传递函数，自变量是频率f（或角频率w）。

对于3阶滤波器：a，b分别有4个值（a0,a1,a2,a3；b0,b1,b2,b3);

对于6阶滤波器：a，b分别有7个值（a0,a1, ... , a6；b0,b1, ... ,b6);

图7是分别是3阶，4阶，6阶300Hz低通，Butterworth滤波器。3阶，4阶，6阶的区别在于：阶次越高，衰减区域越陡；同时分子b, 分母a的个数也不一样。

图7

04

—

时域滤波器

前面构造的滤波器参数b, a 对于时域滤波很重要，我们仍然拿4阶Butterworth滤波器举例：

这个过程是z变换的转换过程。对于离散信号，z的-1次方表示前移一个数，z的-n次方表示前移n个数。即：

对于时域滤波器来说，x信号是滤波前，y信号是滤波后，则滤波过程即如图8，图9：

从图8可以看出，前4个数的计算没有用全 b,a

图8

从图9可以看出，从第5个数开始往后，所有的计算均用全了 b,a，此后的计算在编程中可以用循环语句实现：

图9

最终的滤波效果见图10，结合图8，图9的计算过程及图10的计算结果，可以看出：时域滤波是对x离散信号的实时滤波，也就是x来一个数据点，滤波计算就得到同一时刻，一个y的数据点。

图10

因为图8中前n个数没有用全b,a,所以前n个数的滤波效果不是很好。这里的n即是滤波器的阶次。如图11，6阶滤波器高频滤波效果较4阶要好，但是它会导致b,a的个数比4阶的多，因而前6个数滤波效果不是很好（4阶滤波器是前4个数滤波效果不是很好）。

工程中其实并不太在意前几个离散数据，因此，由于时域滤波具有计算快，实时性等特点，工程应用中多用时域滤波。

图11

05

—

总结

频域滤波，双边谱计算：

过程：傅立叶变换后不需要作任何处理，直接乘以滤波器频响传递函数（复数，双边谱），然后反傅立叶变换，取实部，得到结果（图12）。

优点：直接关注滤波结果，方法简单。

缺点：双边谱不能直观体现滤波前信号频谱特征，频域滤波是得到一段信号后再进行处理，所以有较大延时。

图12

频域滤波，单边谱计算：

过程：傅立叶变换后需要0频保持，双边谱转换成单边谱，乘以滤波器频响传递函数（复数，单边谱），后半部分补0，然后反傅立叶变换，取实部，得到结果（图13）。

优点：单边谱能直观体现信号频谱特征（即直流分量，各频率幅值及相位）。

缺点：双边谱转为单边谱，然后后半部分补0，较烧脑，频域滤波是得到一段信号后再进行处理，所以有较大延时。

图13

时域滤波：

过程：构造滤波器，计算分子分母的b,a, 然后通过一系列的乘法和加法，最终得到结果（图14）。

优点：计算机对乘和加计算迅速，能实现没有延时的实时滤波。

缺点：前n个数没有实现完全滤波（见图8的推导）。

图14

补充：针对前n个数没有实现完全滤波的缺点，有些方法采用倒序再滤波一次（图15），然后和正序滤波结果取算数平均。个人认为，这样就失去了时域滤波的实时性，因为工程中我们也不太关注前n个点，所以该方法也不常用。

图15

以上是个人在滤波器应用中的一些算法总结，不管是频域滤波还是时域滤波，关键是构造滤波器的频响函数及计算参数。

同样的滤波器：频域滤波采用滤波器的频响传函H；时域滤波采用构成滤波器传函的参数 b, a。

滤波器常用的一些场景是：高通滤波滤掉直流分量和趋势项；低通滤波滤掉信号高频毛刺；带通或带阻滤波提取或删掉特别的频率段信号等……

希望本文对您能有所帮助。如有疑问，也可将您的问题发消息提问。谢谢！