手把手教系列之梳状滤波器设计实现
[导读]：前面一篇文章关于IIR/移动平均滤波器设计的文章。本文来聊一聊陷波滤波器，该滤波器在混入谐波干扰时非常有用，算法简单，实现代价低。本文来一探其在机理、应用场景。
注：尽量在每篇文章写写摘要，方便阅读。信息时代，大家时间都很宝贵，如此亦可节约粉丝们的宝贵时间。
前文所说学习的倡导2W1H原则，思来想来并不全面，本文决定从What Why Where When How几个维度展开。我称之为4W1H学习法，借鉴管理学领域中的5W1H，起源于1932年，美国政治学家拉斯维尔提出“5W分析法”，后延伸出5W1H法。有兴趣的可以找来阅读，题外话技术人员读一些方法论管理学方面的书籍于做人做事个人认为是非常有益的。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传
梳状滤波器之What？
在信号处理中，梳状滤波器是通过向其自身添加信号的延迟而实现的，从而造成增强或削弱干扰的滤波器。 梳状滤波器的频率响应由一系列规则间隔的凹口组成，从而呈现出梳状外观。其大体拓扑形式如下：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传
梳状滤波器有着大量不同形式的传递函数，其作用是对周期性信号增强或削弱周期性信号，本文主要介绍其中一种形式的Z传递函数
其中：
其信号流图如下：
梳状滤波器英文称为comb(梳子) filter，这个名字真是无与伦比的绝！为何？谈到滤波器一定会重点关注其对幅频响应曲线，梳状滤波器，正是描述其幅频响应的。而幅频响应从本质上讲是描述系统各频率能量的放大或者衰减。本文中谈到的滤波器就是一个系统，对其输入能量按频率不同进行放大或者衰减，从而起到过滤作用。
梳状滤波器之Why？
前面说到梳状滤波器其幅频响应样子和梳子长的很像，为啥长的像，来一探究竟：
其频率响应为：
现以采样率20000Hz，10阶，阻带带宽50Hz为例。其幅频响应曲线如下：
相频响应曲线为：
从幅频响应曲线可看出，其形状真是如梳子形状，当阶数越大，其齿数越多。这种形式的梳状滤波器对于2KHz，4KHz，6KHz，8KHz，10KHz信号是衰减作用，也即阻止该频率通过，阻带宽度为50Hz。前面谈到了我们可以对某些频率信号加强，而其他的信号衰减吸收。这里引申出其互补滤波器，其Z传递函数刚好有一点形式上的对称性：
其中b为：
此互补滤波器其幅频响应曲线为：
这两个滤波器从其幅频响应曲线刚好是互补对称的。至此，从原理上已基本明了为什么称其为梳状滤波器。
梳状滤波器就其本质也是一种IIR滤波器，因为滤波器的输出反馈参与了滤波运算。
梳状滤波器之Where？
费这么大劲研究梳状滤波器，须的知道在什么地方可以去使用它，事实上梳状滤波器有着大量的应用场合：
- 级联积分梳状（CIC）滤波器，通常用于插值和抽取操作期间的抗混叠，这些操作会更改离散时间系统的采样率。
- PAL和NTSC电视解码器的芯片（也可能是软件滤波）实现的2D和3D梳状滤波器用以降低网点伪影干扰。
- 音频信号处理，包括延迟，镶边和数字波导合成中大量应用。 例如，如果将延迟设置为几毫秒，则可以使用梳状滤波器对圆柱空腔或振动弦中的声驻波的影响进行建模。
- 在天文学中，天体梳滤波器有望将现有光谱仪的精度提高近一百倍。
- 在声学方面，梳齿滤波会以某些不需要的方式出现。 例如，当两个扬声器在距收听者不同距离处播放同一信号时，会对信号产生梳状滤波效果。 在任何封闭的空间中，听众会听到直接声音和反射声音的混合声音。 由于反射的声音路径较长，因此构成了直接声音的延迟版本，并创建了一个梳状滤波器，使两者在听众处合并。
- 仪器仪表也常用梳状滤波器消除谐波干扰，或者选频放大。
- …
梳状滤波器之How？
本文依然借助于matlab的fdatool进行设计示例：
将其重要设置标注如上，其他的与IIR一文类似，就不罗嗦举例了。将重要参数输入，点击设计就轻松设计出相应的滤波器参数了。这里以1000Hz采样率，40Hz带宽，20阶为例，设计出滤波器参数如下：
% Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22.
% Generated on: 05-Apr-2020 13:40:29
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure : Direct-Form II
% Numerator Length : 21
% Denominator Length : 21
% Stable : Yes
% Linear Phase : No
Numerator:
0.38970091175151578
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-0.38970091175151578
Denominator:
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0.22059817649696845
C语言实现非常简单，由于梳状滤波器本质上是IIR滤波器，所以也可以直接利用前文ARM的库函数实现部署。由前面Z传递函数，很容易推导出其差分方程如下：
KaTeX parse error: No such environment: align* at position 8:
\begin{̲a̲l̲i̲g̲n̲*̲}̲
y(n)&=bx(n)-bx…
其互补滤波器差分方程为：
KaTeX parse error: No such environment: align* at position 8:
\begin{̲a̲l̲i̲g̲n̲*̲}̲
y(n)&=bx(n)-bx…
依据上面公式非常容易设计出C代码，这里将浮点数实现共享，如需定点实现也非常容易，代码如下：
#include <stdio.h>
#include <math.h>
#include <string.h>
/*长度应为阶数+1*/
#define CMF_RANK 20
typedef double E_SAMPLE;
typedef enum _E_CMF_TYPE{
CMF_TYPE_STOP_BANDS,
CMF_TYPE_PASS_BANDS
}E_CMF_TYPE;
/*定义移动平均寄存器历史状态*/
typedef struct _t_CMF
{
E_SAMPLE x[CMF_RANK];
E_SAMPLE y[CMF_RANK];
double b;
double r;
E_CMF_TYPE type;
int index;
}t_CMF;
void comb_filter_init(t_CMF * pCmf,double rho,E_CMF_TYPE type)
{
memset(pCmf,0,sizeof(t_CMF));
pCmf->index = CMF_RANK-1;
pCmf->r = rho;
pCmf->type = type;
if(type==CMF_TYPE_STOP_BANDS)
pCmf->b = (1+rho)/2;
else
pCmf->b = (1-rho)/2;
}
E_SAMPLE comb_filter(t_CMF * pCmf,E_SAMPLE xn)
{
E_SAMPLE yn=0;
int n_N;
int i=0;
n_N = pCmf->index;
if(pCmf->type == CMF_TYPE_STOP_BANDS)
{
/*y[n] = bx[n]-bx[n-N]+ry[n-N]*/
yn = pCmf->b*(xn-pCmf->x[n_N])+pCmf->r*pCmf->y[n_N];
}
else
{
/*y[n] = bx[n]+bx[n-N]+ry[n-N]*/
yn = pCmf->b*(xn+pCmf->x[n_N])+pCmf->r*pCmf->y[n_N];
}
/*存储yn为下次迭代准备*/
pCmf->y[n_N] = yn;
pCmf->x[n_N] = xn;
if(pCmf->index==0)
pCmf->index = CMF_RANK-1;
else
pCmf->index--;
return yn;
}
#define SAMPLE_RATE 1000.0f
#define SAMPLE_SIZE 512
#define PI 3.415926f
int main()
{
E_SAMPLE rawSin[SAMPLE_SIZE];
E_SAMPLE outSin[SAMPLE_SIZE];
t_CMF cmf;
FILE *pFile=fopen("./simulationSin.csv","wt+");
if(pFile==NULL)
{
printf("simulationSin.csv opened failed");
return -1;
}
for(int i=0;i<SAMPLE_SIZE;i++)
{
rawSin[i] = 10*sin(2*PI*25*i/SAMPLE_RATE);//+rand()%10;
rawSin[i] += 10*sin(2*PI*50*i/SAMPLE_RATE);
}
/*初始化*/
comb_filter_init(&cmf,-0.22059817649696845,CMF_TYPE_STOP_BANDS);
/*滤波*/
for(int i=0;i<SAMPLE_SIZE;i++)
{
outSin[i]=comb_filter(&cmf,rawSin[i]);
}
for(int i=0;i<SAMPLE_SIZE;i++)
{
fprintf(pFile,"%f,",rawSin[i]);
}
fprintf(pFile,"\n");
for(int i=0;i<SAMPLE_SIZE;i++)
{
fprintf(pFile,"%f,",outSin[i]);
}
/*初始化*/
comb_filter_init(&cmf,-0.22059817649696845,CMF_TYPE_PASS_BANDS);
/*滤波*/
for(int i=0;i<SAMPLE_SIZE;i++)
{
outSin[i]=comb_filter(&cmf,rawSin[i]);
}
/*存储数据*/
fprintf(pFile,"\n");
for(int i=0;i<SAMPLE_SIZE;i++)
{
fprintf(pFile,"%f,",outSin[i]);
}
fclose(pFile);
return 0;
}
同样利用excel,生成波形如下：
可见,经过梳状滤波器过滤后，50Hz谐波已经被过滤掉，25Hz 保留下来，而经过其互补滤波器后，25Hz被过滤，其50Hz被保留。如看时域波形不直观，可利用Python代码进行FFT展开分析：
梳妆滤波后FFT谱线图：
互补梳状滤波器过滤后FFT谱线图：
总结：
- 梳妆滤波器本质上是一种IIR滤波器
- 梳妆滤波器在滤除谐波上，利用极小代价就可以比较好的滤除谐波干扰
- 其互补滤波器在时域时会失真，使用时需要考虑
- 如需要考虑计算效率，可以考虑定点优化，但精度会下降。
- 梳妆滤波器在2D/3D图像处理广为应用，如有兴趣可深入研究
版权声明：所有文章版权归嵌入式客栈所有，如商业使用，须嵌入式客栈授权。欢迎关注微信公众号，内容更丰富。
推荐阅读