手把手教系列之梳状滤波器设计实现

时间 2020-05-23

标签手把手系列滤波器设计实现繁體版

原文原文链接

[导读]：前面一篇文章关于IIR/移动平均滤波器设计的文章。本文来聊一聊陷波滤波器，该滤波器在混入谐波干扰时很是有用，算法简单，实现代价低。本文来一探其在机理、应用场景。算法

注：尽可能在每篇文章写写摘要，方便阅读。信息时代，你们时间都很宝贵，如此亦可节约粉丝们的宝贵时间。微信

前文所说学习的倡导2W1H原则，思来想来并不全面，本文决定从What Why Where When How几个维度展开。我称之为4W1H学习法，借鉴管理学领域中的5W1H，起源于1932年，美国政治学家拉斯维尔提出“5W分析法”，后延伸出5W1H法。有兴趣的能够找来阅读，题外话技术人员读一些方法论管理学方面的书籍于作人作事我的认为是很是有益的。
函数

梳状滤波器之What？

在信号处理中，梳状滤波器是经过向其自身添加信号的延迟而实现的，从而形成加强或削弱干扰的滤波器。梳状滤波器的频率响应由一系列规则间隔的凹口组成，从而呈现出梳状外观。其大致拓扑形式以下：
学习

梳状滤波器有着大量不一样形式的传递函数，其做用是对周期性信号加强或削弱周期性信号，本文主要介绍其中一种形式的Z传递函数\(H(Z)=b\frac{1-Z^{-N}}{1-{\rho}N^{-N}}\)优化

其中：\(b=\frac{1+\rho}{2}\)spa

其信号流图以下：
设计

梳状滤波器英文称为comb(梳子) filter，这个名字真是无与伦比的绝！为什么？谈到滤波器必定会重点关注其对幅频响应曲线，梳状滤波器，正是描述其幅频响应的。而幅频响应从本质上讲是描述系统各频率能量的放大或者衰减。本文中谈到的滤波器就是一个系统，对其输入能量按频率不一样进行放大或者衰减，从而起到过滤做用。excel

梳状滤波器之Why？

前面说到梳状滤波器其幅频响应样子和梳子长的很像，为啥长的像，来一探究竟：code

其频率响应为：orm

\[H(e^{j\omega})=b\frac{1-e^{j\omega{N}}}{1-\rho{e^{j\omega{N}}}} \]

现以采样率20000Hz，10阶，阻带带宽50Hz为例。其幅频响应曲线以下：

相频响应曲线为：

![](E:\blog\embInn\DSP\comb Filter\pic\phase12.png)

从幅频响应曲线可看出，其形状真是如梳子形状，当阶数越大，其齿数越多。这种形式的梳状滤波器对于2KHz，4KHz，6KHz，8KHz，10KHz信号是衰减做用，也即阻止该频率经过，阻带宽度为50Hz。前面谈到了咱们能够对某些频率信号增强，而其余的信号衰减吸取。这里引伸出其互补滤波器，其Z传递函数恰好有一点形式上的对称性：

\[H(Z)=b\frac{1+Z^{-N}}{1-{\rho}N^{-N}} \]

其中b为：

\[b=\frac{1-\rho}{2} \]

此互补滤波器其幅频响应曲线为：

这两个滤波器从其幅频响应曲线恰好是互补对称的。至此，从原理上已基本明了为何称其为梳状滤波器。

梳状滤波器就其本质也是一种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传递函数，很容易推导出其差分方程以下：

\[\begin{align*} y(n)&=bx(n)-bx(n-N)+\rho{y(n-N)}\\&=\frac{1+\rho}{2}[x(n)-x(n-N)]+\rho{y(n-N)} \end{align*} \]

其互补滤波器差分方程为：

\[\begin{align*} y(n)&=bx(n)-bx(n-N)+\rho{y(n-N)}\\&=\frac{1-\rho}{2}[x(n)+x(n-N)]+\rho{y(n-N)} \end{align*} \]

依据上面公式很是容易设计出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图像处理广为应用，若有兴趣可深刻研究

文章出自微信公众号：嵌入式客栈，更多内容，请关注本人公众号,严禁商业使用，违法必究