5.IIR数字滤波器

---

1.简介

在做音频eq的时候，iir滤波器是必不可少的，所以这节主要学习一下iir滤波器，当然，不会偏于理论，而是着重怎么用。

iir滤波器主要是用来处理信号的频率的，这是第一个最基本的认识。

摆一下基本的公式：

系统函数：

差分方程：

再来说一下该滤波器的结构：
a.直接I型
b.直接II型(典范型)
c.级联型
d.并联型

bcd三种类型都是从直接I型转化而来的，这是第二个最基本的认识。

iir滤波器对于叠加了噪声的信号，也就是信号和噪声的频谱相互混叠是无能为力的，这就是其用在eq上一个基本的原因，这是第三个最基本的认识。

具体再使用一个例子来表述如何使用该公式，由于使用matlab仿真条件下，反馈符号是负号，故需要将代入一个负号。
假设是二阶滤波器，并且分子，分母如下：


样本点如下：

则计算一下输出值：

使用matlab编写上述程序，如下所示：

clc;
clear all;
format long;
b =  [0.6,0.8,1.0];
a = [1,0.5,1.5];
x = [1,2.5,1.5,0.5];
y=filter(b,a,x)

有一点要注意，在做dsp开发的时候，a,b的排序并不一定是按照，有可能是，这一点需要查看具体的API即可。

2.关于iir滤波器系数

可以使用matlab来计算出iir滤波器的系数，当然，也有其他的工具可以达成这样的功能，例如：一个直接计算的在线工具

3.二阶滤波器例子

3.1 简单的一段2阶滤波器

再来写一点例子程序，例如c代码:

void test_iir_raw()
{
    int i,j,k;
    int len = 4;
    float in[4]={0.0,0.1,0.2,0.3};
    float geq_para[5]={0.1,0.3,0.5,0,0};
    float a[3],b[3];
    float geq_state[2]={0.0};
    float geq_state_in[2]={0.0};
    float ym;
    float g=1.0;
    float tmpdata[4];

    a[0] = 1;
    a[1] = geq_para[1];
    a[2] = geq_para[2];
    b[0] = geq_para[0];
    b[1] = 0;
    b[2] = -b[0];

    for (k = 0; k < len; k++)
    {
        ym = geq_para[0] * (in[k] - geq_state_in[1]) - geq_para[2] * geq_state[1] - geq_para[1] * geq_state[0];
        geq_state[1] = geq_state[0];
        geq_state[0] = ym;
        geq_state_in[1] = geq_state_in[0];
        geq_state_in[0] = in[k];
        tmpdata[k] += ym * g;
    }
    for(i = 0 ; i < 4 ; i ++){
        printf("tmpdata:%f\n",tmpdata[i]);
    }
}

再例如matlab代码：

clc;
clear all;
format long;
b =  [0.1,  0,   -0.1];
a = [1,   0.3,   0.5];
xk=[0.0,0.1,0.2,0.3];
y=filter(b,a,xk)

可以很简单的分析到：

geq_state可以代表y(n)的临时变量.geq_state_in代表的是x(n)的临时变量.

我再给出结果:

tmpdata:0.000000
tmpdata:0.010000
tmpdata:0.017000
tmpdata:0.009900

最后，在说明一下本例子中的一个临时结果:

geq_state[0]=0.009900,geq_state[1]=0.017000
geq_state_in[0]=0.300000,geq_state_in[1]=0.200000

3.2 简单的二段2阶滤波器

这是一个并联的二阶滤波器，相当于一段音频做两次二阶滤波器效果，对应的是在每一个点上加上上一次做的结果值，我列出c代码：

void test_iir_raw()
{
    int i,j,k;
    int len = 4;
    float in[4]={0.0,0.1,0.2,0.3};
    float geq_para[GEQ_BANDS][5]={{0.1,0.3,0.5,0,0},{0.1,0.3,0.5,0,0}};
    float a[3],b[3];
    float geq_state[GEQ_BANDS][2]={0.0};
    float geq_state_in[GEQ_BANDS][2]={0.0};
    float ym;
    float g=1.0;
    float tmpdata[4];

    for(j = 0 ; j < GEQ_BANDS ; j ++){
        for (k = 0; k < len; k++)
        {
            ym = geq_para[j][0] * (in[k] - geq_state_in[j][1]) - geq_para[j][2] * geq_state[j][1] - geq_para[j][1] * geq_state[j][0];
            geq_state[j][1] = geq_state[j][0];
            geq_state[j][0] = ym;
            geq_state_in[j][1] = geq_state_in[j][0];
            geq_state_in[j][0] = in[k];
            tmpdata[k] += ym * g;
        }
    }

    for(i = 0 ; i < 4 ; i ++){
        printf("tmpdata:%f\n",tmpdata[i]);
    }
}

对应的matlab结果就是上次给出的值中每一项乘以2即可.至此，这个部分结束.

4.定点

首先给出biquad滤波器的差分函数形式：

注意:
Matlab里的计算就是按照上面的式子计算的，但是STM32F4DSP库里的系数a1，a2是取反的。

再贴出来biquad滤波器的实现代码：

        for (k = 0; k < len; k++)
        {
            ym = geq_para[j][0] * (in[k] - geq_state_in[j][1]) - geq_para[j][2] * geq_state[j][1] - geq_para[j][1] * geq_state[j][0];
            geq_state[j][1] = geq_state[j][0];
            geq_state[j][0] = ym;
            geq_state_in[j][1] = geq_state_in[j][0];
            geq_state_in[j][0] = in[k];
            tmpdata[k] += ym * g;
        }

为了让定点化来的顺利，可以在第一步的时候只定点化ym，也就是说：

            geq_state[j][1] = geq_state[j][0];
            geq_state[j][0] = ym;
            geq_state_in[j][1] = geq_state_in[j][0];
            geq_state_in[j][0] = in[k];
            tmpdata[k] += ym * g;

这部分不变，而先将：

ym = geq_para[j][0] * (in[k] - geq_state_in[j][1]) - 
geq_para[j][2] * geq_state[j][1] - geq_para[j][1] * geq_state[j][0];

各个部分进行定点：

long long int tmp1 = geq_para[j][0] * (1 << 26);
long long int tmp2 = geq_para[j][2] * (1 << 26);
long long int tmp3 = geq_para[j][1] * (1 << 26);
tmp1 = (tmp1 * tmp0) >> 26;
tmp2 = (tmp2 * tmpp1) >> 26;
tmp3 = (tmp3 * tmpp2) >> 26;
ym = (tmp1 - tmp2 - tmp3)*1.0 / (1 << 26);

之后再将整体浮点数替换，从而真正实现定点化(64位)。为了方便移植到硬件上，还需要进一步定点化。

可以按照如下方法做进一步处理：

long long int tmp0 = (in[i][k])*(1 << 28) - geq_state_in_int[i][j][1];
long long int tmp1 = geq_para[j][0] * (1 << 26);
long long int tmp2 = geq_para[j][2] * (1 << 26);
long long int tmp3 = geq_para[j][1] * (1 << 26);

long long int tmpp1 = geq_state_int[i][j][1] ;
long long int tmpp2 = geq_state_int[i][j][0] ;

tmp1 = (tmp1 * tmp0) >> 26;
tmp2 = (tmp2 * tmpp1) >> 26;
tmp3 = (tmp3 * tmpp2) >> 26;

ym_int = tmp1 - tmp2 - tmp3;
geq_state_int[i][j][1] = geq_state_int[i][j][0];
geq_state_int[i][j][0] = ym_int;
geq_state_in_int[i][j][1] = geq_state_in_int[i][j][0];
geq_state_in_int[i][j][0] = in[i][k] * (1 << 28);
tmpdata_int[k] += ym_int * g;