手把手教系列之IIR滤波器设计

时间  2020-05-20
标签 手把手 系列 iir 滤波器 设计 繁體版
原文   原文链接

[导读]：在嵌入式系统中常常须要采集模拟信号，采集模拟信号的信号链中不免引入干扰，那么如何滤除干扰呢？今天就来个一步一步描述如何设计部署一个IIR滤波器到你的系统。算法

何为IIR滤波器？

无限冲激响应（IIR： Infinite Impulse Response）是一种适用于许多线性时不变系统的属性，这些系统的特征是具备一个冲激响应h（t），该冲激响应h（t）不会在特定点上彻底变为零，而是无限期地持续。 这与有限冲激响应（FIR： Finite Impulse Response）系统造成对比，在有限冲激响应（FIR）系统中，对于某个有限T，在时间t> T时，冲激响应确实刚好变为零。 线性时不变系统的常见示例是大多数电子和数字滤波器。 具备此属性的系统称为IIR系统或IIR滤波器。编程

这是常见的教科书式数学严谨定义，不少人看到这一下就蒙了，能说人话吗？数组

线性时不变系统理论俗称LTI系统理论，源自应用数学，直接在核磁共振频谱学、地震学、电路、信号处理和控制理论等技术领域运用。它研究的是线性、非时变系统对任意输入信号的响应。虽然这些系统的轨迹一般会随时间变化（例如声学波形）来测量和跟踪，可是应用到图像处理和场论时，LTI系统在空间维度上也有轨迹。所以，这些系统也被称为线性非时变平移，在最通常的范围理论给出此理论。在离散（即采样）系统中对应的术语是线性非时变平移系统。由电阻、电容、电感组成的电路是LTI系统的一个很好的例子。好比一个运放系统在必定频带范围内知足信号的时域叠加，输入一个100Hz和200Hz正弦信号，输出频率是这两种信号的线性叠加。微信

用数学对LTI系统描述：函数

线性：输入\(x_1(t)\),产生响应\(y_1(t)\),而输入\(x_2(t)\),产生响应\(y_2(t)\),那么缩放和加和输入\(a_1x_1(t)+a_2x_2(t)\),产生缩放、加和的响应\(a_1y_1(t)+a_2y_2(t)\),其中\(a_1\)和\(a_2\)是标量，对于任意的有：输入\(\sum_{k=0}^{N}a_kx_k(t)\),产生响应\(\sum_{k=0}^{N}a_ky_k(t)\)测试

时不变性：指若是将系统的输入信号延迟\(\tau\)秒，那么获得的输出响应也相应延时\(\tau\)秒。用数学描述，也即若是输入\(x_1(t)\),产生响应\(y_1(t)\),而输入\(x_1(t+\tau)\)，产生响应\(y_1(t+\tau)\)。这么描述仍是不易懂，来个图，有图有真相：优化

假定一个信号放大电路对100Hz正弦信号放大2倍：ui

image-20200326205428645

则输出为：spa

image-20200326205811819

而对200Hz的正弦信号，其放大倍数为1.7倍。（作过运放电路设计的朋友应该有经验，在其同频带其放大倍数每每并不平坦，也即幅频响应在频带内不平坦，这是比较常见的）。命令行

也即输入为：

image-20200326210319201

响应为:

image-20200326210332571

那么若是输入100Hz和200Hz的时域叠加信号，则其输入为：

image-20200326210550908

则其响应为：

image-20200326210603870

由这些图可看出，输入信号的形状保持不变，输出为对应输入的线性时域叠加。

对于时不变，就不用图描述了，在一个真实电路中，若是输入延迟必定时间，则响应对应延迟相同时间输出。

上面这么多文字只是为了描述在什么场合可使用IIR滤波器对信号进行数字滤波。总结而言，就是在线性时不变系统中适用。换言之，在大多数电路系统中咱们均可以尝试采用IIR滤波器进行数字滤波。

那么究竟什么是IIR滤波器呢？从数字信号处理的书籍中咱们能看到这样的Z变换信号流图：

image-20200326221206401

image-20200326214242009表示延迟一拍,在数字系统中表示对于输入信号而言,即为上一次采样值,对于输出而言,即为上一次的输出值。

在时域中对于上述流图，用时域描述即为：

image-20200326221258617

若是熟悉Z变换，则Z变换传递函数为：

image-20200326221348262

上述数字滤波器，若是从编程的角度来看，x(n-1),表示上一次的信号，多是来自ADC的上次采样，而y(n-1)则为上一次滤波器的输出值，对应就比较好理解x(n-N)就表示前第n次输入样本信号，而y(n-M)则为前第M次滤波器的输出。

说了这么多，只是为了更好的理解概念，只有概念理解正确，才能用争取。 概念理解这对工程师而言，很是之重要。

如何设计呢？

打开fdatool

MATLAB提供了很是容易使用的FDATool帮助咱们设计数字滤波器，真正精彩的地方开始了，让咱们拭目以待究竟如何一步一步设计并实施一个IIR滤波器。

首先打开MATLAB,在命令行中敲fdatool

image-20200326233121867

弹出窗体就是fdatool了，以下：

在设计具体,有几个相关概念须要澄清：

Fs：采样率，单位为Hz,真实部署在系统中，请务必确保样本是按照恒定采样率进行采样，不然将得不到想要的效果。

Fpass: 通频带，单位为Hz,即系统中指望经过的最高频率。

Fstop: 截至频率，即幅频响应的-3dB处的频率，这个如不理解，请自行查阅相关书籍。

分贝dB: 这是一个无单位反应输出与输入倍数的一个术语。电学中分贝与放大倍数的转换关系为：

  • A(V)(dB)=20lg(Vo/Vi)；电压增益,Vo 为输出电压，Vi为输入电压
  • A(I)(dB)=20lg(Io/Ii)；电流增益,Io 为输出电流，Ii为输入电流
  • A(p)(dB)=10lg(Po/Pi)；功率增益，Po 为输出功率，Pi为输入功率

滤波器类型：这里有Butterworth(巴特沃斯)、Chebyshev Type I,Chebyshev Type II、(切比雪夫)、Elipic 等可选。

image-20200327002802811

  • 巴特沃斯 Butterworth，也被称做最大平坦滤波器。巴特沃斯滤波器的特色是通频带内的频率响应曲线最大限度平坦，没有纹波。
  • 切比雪夫 Chebyshev，是在通带或阻带上频率响应幅度等波纹波动的滤波器。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。
  • 椭圆 Elliptic，椭圆滤波器是在通带和阻带等波纹的一种滤波器。
  • ...这里就不一一介绍了，有兴趣能够去查信号处理书籍。

就其特色，这里对其中几种略做介绍：

  • 巴特沃斯具备最平坦的通带。
  • 椭圆滤波器衰减最快，可是通带、阻带都有波纹。
  • 切比雪夫滤波器衰减比巴特沃斯快，但比椭圆滤波器慢，波纹区域可选择。

假设咱们须要设计一个IIR滤波器，采样率为32000Hz, 有用信号频率在2000Hz内，设计IIR滤波器对信号进行数字滤波。这里为节省算力，咱们指定滤波器的阶数，也即传递函数中N/M中的最大值，通常而言N大于M。

这里指定阶数为8阶，类型指定为巴特沃斯型IIR滤波器，输入阶数8阶，采样率32000Hz,截至频率设置为2000Hz而后点击Design Filter以下图所示：

其相频响应曲线以下：

楚辞以外，咱们还能够将幅频与相频曲线放在一个频率坐标上去看设计结果：

导出滤波器参数，这里咱们选择，

选择ASCII码

而后就获得了一个文件，保存2kHz_LPF.fcf，文件名随你喜欢。

文件内容以下：

% Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22.
% Generated on: 27-Mar-2020 21:27:06

% Coefficient Format: Decimal

% Discrete-Time IIR Filter (real)                            
% -------------------------------                            
% Filter Structure    : Direct-Form II, Second-Order Sections
% Number of Sections  : 4                                    
% Stable              : Yes                                  
% Linear Phase        : No                                   

                                                            
SOS Matrix:                                                  
1  2  1  1  -1.7193929141691948  0.8610574795347461          
1  2  1  1  -1.5237898734101736  0.64933827386370635         
1  2  1  1  -1.4017399331200424  0.51723237044751591         
1  2  1  1  -1.3435020629061745  0.45419615396638446         
                                                             
Scale Values:                                                
0.035416141341387819                                         
0.031387100113383172                                         
0.028873109331868367                                         
0.027673522765052503

至此设计工做就结束了，立刻进入滤波器的部署测试阶段。

部署测试滤波器

到这里，没有经验的朋友可能会说，这么一堆参数我该咋用呢？

须要本身去写前面描述的计算公式吗？固然你也能够这么作，这里就不写了,ARM的CMSIS库已经帮你们设计好了种类繁多的数字信号处理函数实现了，并且通过了测试，这里直接拿来用便可。有兴趣本身写也不难，只要理解Z传递函数概念内涵，很是容易实现。这里咱们采用32位浮点实现函数：arm_biquad_cascade_df1_f32。

该函数位于：CMSIS\DSP\Source\FilteringFunctions\arm_biquad_cascade_df1_init_f32.c中

以及CMSIS\DSP\Source\FilteringFunctions\arm_biquad_cascade_df1_f32.c

咱们来看一看这个函数：

arm_biquad_cascade_df1_init_f32.c：

/*
*做用      :初始化滤波器
*S        :指向浮点SOS级联结构的实例。
*numStages:滤波器中二阶SOS的数量
*pCoeffs  :滤波器参数指针,参数按下列顺序存储
*          {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...}
*pState   :历史状态缓冲区指针
*/
void arm_biquad_cascade_df1_init_f32(
        arm_biquad_casd_df1_inst_f32 * S,
        uint8_t numStages,
  const float32_t * pCoeffs,
        float32_t * pState)
{
  /* Assign filter stages */
  S->numStages = numStages;

  /* Assign coefficient pointer */
  S->pCoeffs = pCoeffs;

  /* Clear state buffer and size is always 4 * numStages */
  memset(pState, 0, (4U * (uint32_t) numStages) * sizeof(float32_t));

  /* Assign state pointer */
  S->pState = pState;
}

arm_math.h 定义了须用到的结构体，对于本例相关的结构体为arm_biquad_casd_df1_inst_f32

typedef struct
{
  unsigned int numStages; /*2阶节的个数，应为2*numStages.        */
  float *pState;          /*状态系数数组指针，数组长度为4*numStages*/
  float *pCoeffs;         /*系数数组指针， 数组的长度为5*numStages.*/
} arm_biquad_casd_df1_inst_f32;

滤波器具体滤波函数为arm_biquad_cascade_df1_f32

/**
 *  *S       ：指向浮点Biquad级联结构的实例.
 *  *pSrc    ：指向输入数据块。
 *  *pDst    ：指向输出数据块。
 *  blockSize：每次调用要处理的样本数。
 *  返回值    ：无.
 */
void arm_biquad_cascade_df1_f32(
  const arm_biquad_casd_df1_inst_f32 * S,
  float * pSrc,
  float * pDst,
  unsigned int blockSize)
{
  float *pIn = pSrc;                         /*源指针     */
  float *pOut = pDst;                        /*目的指针    */
  float *pState = S->pState;                 /*状态指针    */
  float *pCoeffs = S->pCoeffs;               /*参数指针    */
  float acc;                                 /*累加器      */
  float b0, b1, b2, a1, a2;                  /*滤波器参数   */
  float Xn1, Xn2, Yn1, Yn2;                  /*滤波器状态变量*/
  float Xn;                                  /*临时输入     */
  unsigned int sample, stage = S->numStages; /*循环计数     */

  do
  {
    /* Reading the coefficients */
    b0 = *pCoeffs++;
    b1 = *pCoeffs++;
    b2 = *pCoeffs++;
    a1 = *pCoeffs++;
    a2 = *pCoeffs++;

    Xn1 = pState[0];
    Xn2 = pState[1];
    Yn1 = pState[2];
    Yn2 = pState[3];

    sample = blockSize >> 2u;

    while(sample > 0u)
    {
      /* 读第一个输入 */
      Xn = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn2 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);

      /* Store the result in the accumulator in the destination buffer. */
      *pOut++ = Yn2;

      /* 每次计算输出后，状态都应更新. */
      /* 状态应更新为:  */
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */

      /* Read the second input */
      Xn2 = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn1 = (b0 * Xn2) + (b1 * Xn) + (b2 * Xn1) + (a1 * Yn2) + (a2 * Yn1);

      /* 将结果存储在目标缓冲区的累加器中. */
      *pOut++ = Yn1;

      /* 每次计算输出后，状态都应更新. */
      /* 状态应更新为:  */
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */

      /*读第三个输入 */
      Xn1 = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn2 = (b0 * Xn1) + (b1 * Xn2) + (b2 * Xn) + (a1 * Yn1) + (a2 * Yn2);

      /* 将结果存储在目标缓冲区的累加器中. */
      *pOut++ = Yn2;

      /* 每次计算输出后，状态都应更新. */
      /* 状态应更新为: */
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      /* 读第四个输入 */
      Xn = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn1 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn2) + (a2 * Yn1);

      /* 将结果存储在目标缓冲区的累加器中. */
      *pOut++ = Yn1;

      /* 每次计算输出后，状态都应更新. */
      /* 状态应更新为:  */
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      Xn2 = Xn1;
      Xn1 = Xn;

      /* 递减循环计数器 */
      sample--;
    }

    /* 若是blockSize不是4的倍数，
    *请在此处计算任何剩余的输出样本。
    *不使用循环展开. */
    sample = blockSize & 0x3u;

    while(sample > 0u)
    {
      /* 读取输入 */
      Xn = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      acc = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);

      /* 将结果存储在目标缓冲区的累加器中. */
      *pOut++ = acc;

      /* 每次计算输出后，状态都应更新。 */
      /* 状态应更新为:    */
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      Xn2 = Xn1;
      Xn1 = Xn;
      Yn2 = Yn1;
      Yn1 = acc;

      /* d递减循环计数器 */
      sample--;
    }

    /*  将更新后的状态变量存储回pState数组中 */
    *pState++ = Xn1;
    *pState++ = Xn2;
    *pState++ = Yn1;
    *pState++ = Yn2;

    /*第一阶段从输入缓冲区到输出缓冲区.     */
    /*随后的numStages在输出缓冲区中就地发生*/
    pIn = pDst;

    /* 重置输出指针 */
    pOut = pDst;

    /* 递减循环计数器 */
    stage--;

  } while(stage > 0u);
}

开始测试：

#include <stdio.h>
#include <math.h>
/*
SOS Matrix:                                                  
1  2  1  1  -1.7193929141691948  0.8610574795347461          
1  2  1  1  -1.5237898734101736  0.64933827386370635         
1  2  1  1  -1.4017399331200424  0.51723237044751591         
1  2  1  1  -1.3435020629061745  0.45419615396638446         
                                                             
Scale Values:                                                
0.035416141341387819                                         
0.031387100113383172                                         
0.028873109331868367                                         
0.027673522765052503  
作以下转换：
1.缩放
[1  2  1] * 0.035416141341387819
[1  2  1] * 0.031387100113383172
[1  2  1] * 0.028873109331868367
[1  2  1] * 0.027673522765052503
获得：
[0.035416141341387819  2*0.035416141341387819  0.035416141341387819]
[0.031387100113383172  2*0.031387100113383172  0.031387100113383172] 
[0.028873109331868367  2*0.028873109331868367  0.028873109331868367] 
[0.027673522765052503  2*0.027673522765052503  0.027673522765052503]
2.舍掉第四列参数
3.将后两列分别乘以-1，即：
0.035416141341387819  2*0.035416141341387819  0.035416141341387819  -1.7193929141691948  0.8610574795347461          
0.031387100113383172  2*0.031387100113383172  0.031387100113383172  -1.5237898734101736  0.64933827386370635         
0.028873109331868367  2*0.028873109331868367  0.028873109331868367  -1.4017399331200424  0.51723237044751591         
0.027673522765052503  2*0.027673522765052503  0.027673522765052503  -1.3435020629061745  0.45419615396638446 
这样就获得了滤波器系数组了
*/
#define IIR_SECTION 4                /*见前面设计输出为4个SOS块*/
static float iir_state[4*IIR_SECTION];/*历史状态缓冲区         */
const float iir_coeffs[5*IIR_SECTION]={
    0.035416141341387819,2*0.035416141341387819,0.035416141341387819,1.7193929141691948,-0.8610574795347461,    0.031387100113383172,2*0.031387100113383172,0.031387100113383172,1.5237898734101736,-0.64933827386370635,    0.028873109331868367,2*0.028873109331868367,0.028873109331868367,1.4017399331200424,-0.51723237044751591,    0.027673522765052503,2*0.027673522765052503,0.027673522765052503,1.3435020629061745,-0.45419615396638446
};
static arm_biquad_casd_df1_inst_f32 S;
/*假定采样512个点*/
#define BUF_SIZE 512
#define PI 3.1415926
#define SAMPLE_RATE  32000 /*32000Hz*/
int main()
{
    float raw[BUF_SIZE];
    float raw_4k[BUF_SIZE];
    float raw_out[BUF_SIZE];

    float raw_noise[BUF_SIZE];
    float raw_noise_out[BUF_SIZE];

    arm_biquad_casd_df1_inst_f32 S;
    FILE *pFile=fopen("./simulation.csv","wt+");
    if(pFile==NULL)
    {
        printf("file opened failed");
        return -1;
    }

    for(int i=0;i<BUF_SIZE;i++)
    {
    /*模拟9000Hz输入信号         */
        raw[i] = 0.5*1024.0/3*sin(2*PI*800*i/32000.0f)+rand()%50;
        raw_4k[i] = 0.5*1024.0/3*sin(2*PI*4000*i/32000.0f);
    /*模拟9000Hz +11000Hz叠加输入*/
        raw_noise[i] = raw[i] + raw_4k[i];
    }

    arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]);
    arm_biquad_cascade_df1_f32(&S, raw, raw_out, BUF_SIZE);

    for(int i=0;i<BUF_SIZE;i++)
    {
       fprintf(pFile,"%f,",raw[i]);
    }

    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_4k[i]);
    }
    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_out[i]);
    }

    /*从新初始化*/
    arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]);
    arm_biquad_cascade_df1_f32(&S, raw_noise, raw_noise_out, BUF_SIZE);

    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_noise[i]);
    }

    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_noise_out[i]);
    }

    fclose(pFile);

    return 0;
}

利用csv文件，将模拟数据存储，直接用excel打开，将行数据生成曲线图以下：

image-20200517125200485

  • 第一幅图，为800Hz信号混入随机噪声的波形
  • 第二幅图，为4000Hz信号，对假定系统为无用干扰信号
  • 第三幅图， 为800Hz 混入随机噪声过滤后，已经很好的还原有用信号频率
  • 第四幅图， 为800Hz信号混入随机噪声，同时叠加4000Hz干扰的波形，对系统而言，从时域中，明显可见，有用信号已经彻底扭曲
  • 第五幅图，为800Hz信号混入随机噪声，同时叠加4000Hz干扰的输入，通过该低通滤波器后的波形，与第三幅图基本同样，已经很是好的滤除了干扰信号。

总结：

  • IIR滤波器在线性时不变系统中能够很好的解决工程中通常噪声问题

  • 若是须要设计带通、高通滤波器其步骤基本相似，只是滤波器的参数以及SOS块个数可能不同而已

  • 须要提醒的时，IIR的相频响应不线性，若是系统对相频响应有严格要求，就须要采用其余的数字滤波器拓扑形式了

  • 实际应用中，若是阶数不高时，如今算力强劲的单片机或者DSP以及能够直接使用浮点处理。

  • 若是对处理速度有严格的实时要求，须要在极短期进行滤波处理，能够考虑下降阶数，或采用定点IIR滤波算法实现。也或者将文中函数进行汇编级优化。

文章出自微信公众号：嵌入式客栈，关注公众号，获取更新更多内容

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程