主旨
在生理信号处理中,经常有滤除低频信号的需求,例如在分析心率相关问题时排除掉呼吸等因素造成的基线漂移等低频干扰。在这个场景的工程实践中,移动平均滤波器由于设计实现简单,经常被用到。但是移动平均滤波器的频谱分析缺乏类似于FIR或IIR相关定量资料,在工作中时常凭经验确定阶数,严谨性不足。本文试图采用理论分析+数值计算方法给出移动平均滤波器频谱的一般规律,用于指导日常工作
代码位置
本文中的数值计算使用Matlab实现,代码和输出表格放在如下位置
https://github.com/wangyaobupt/MovingAverageFilter
移动平均滤波器介绍
考虑时间上无限长度的输入信号 x(n),n=−∞,..,0,1,...∞
移动平均滤波器的输出为
y(n)=12m+1∑k=n−mn+mx(n−k)
其中m为正整数,2m+1称为滤波器阶数
工程上不可能有无限长度信号,对于输入信号 x(n),n=0,1,..N ,当 k<0 时取 x(k)=0 。此滤波器引入m个点的延迟,即x(0)对应的输出点为y(m). y(0)/y(1)/…/y(m-1)是滤波器非稳态输出,应该抛弃。
数学推导频谱
考虑输入信号为冲击函数
h(n)={10if n=0otherwise.
移动平均滤波器的输出为
y(n)=12m+1∑k=n−mn+mh(n−k)
按照傅里叶变换,滤波器频域响应函数为
H(ejΩ)=12m+1∑k=n−mn+me−jkΩ=12m+1{(ejmΩ+e−jmΩ)+(ej(m−1)Ω+e−j(m−1)Ω)+...+(ej(Ω+e−jΩ)+1}=12m+1{1+2∑k=1mcos(kΩ)}
其中
Ω=ωfs=2πffs
上述两个式子综合起来,移动平均滤波器的频率响应可以写为如下形式
H(ej∗2πf)=12m+1{1+2∑k=1mcos(k∗2πf)}
数值计算频谱特征
先考察几个典型阶数的滤波器幅度-频率响应曲线,考虑采样率为每秒125个点
本节提到的所有图片都可以在我的GitHub MovingAverageFilter工程的picture目录下载到清晰版本
- m=1, 3阶移动平均滤波器
- m=10, 21阶移动平均滤波器
- m=61, 123阶移动平均滤波器
从上述图形可以看出,随着滤波器阶数的上升,通带截止频率(以3dB功率点考虑)持续降低,对于3阶移动平均滤波器,通带截止频率在20Hz附近,而对于123阶移动平均滤波器,通带截止频率在0.45Hz附近。同时,阶数上升对于阻带的衰减能力也持续上升。
那么,阶数与通带截止频率的定量关系如何呢,我们继续通过数值计算寻找答案。下图中横坐标为移动平均滤波器阶数,纵坐标为截止频率,通过这张图,可以清楚的看出随着滤波器阶数上升,截止频率以指数速度降低。
为了工程上便于使用,我将上述曲线保存成表格,设计移动平均滤波器时根据截止频率要求查表即可。
表格位置
https://github.com/wangyaobupt/MovingAverageFilter/blob/master/Code/MovingAverageFilter_PassBandUpperLimit_VS_order.csv
结论
移动平均滤波器可以实现低通滤波器的效果,其滤波器截止频率和衰减特征虽然不易给出直观的数学表达式,但是可以通过数值计算给出符合工程应用的表格,用于指导日常工作。