原创 二次IIR滤波器

 在声音信号处理中，经常需要设计滤波器，二次IIR滤波器由于其参数设计简单，运算量少，运用的最为广泛。这里介绍几种常用的滤波器的参数计算。

Peak  Equalizer
EQ----大家使用的音乐播放器中就有EQ的设置，如千千静听，cooledit等软件上面都有。

Peak   Equalizer的设计参数为中心频率f，最大增益db，形状参数Q。下面这个Python小程序通过这三个参数计算滤波器的参数a和b，取样频率为44100Hz。 Python是一种开源语言，这里我就不介绍了，我也只是略知一二。转化成C语言还是比较简单的。与C不同的地方我会加上注释。

Toggle Line Number
01    def peak(f, db, Q):
02    A = 10.0**(db/40.0)                        //A = 10.0^(db/40.0)
03    w0= 2*math.pi*f/44100.0             //44100=FS采样频率，你也可以取48000.
04    alpha = math.sin(w0)/2/Q
05    b0 = 1+alpha*A
06    b1 = -2*math.cos(w0)
07    b2 = 1-alpha*A
08    a0 = 1+alpha/A
09    a1 = -2*cos(w0)
10    a2 = 1-alpha/A
11    b0 /= a0
12    b1 /= a0
13    b2 /= a0
14    a1 /= a0
15    a2 /= a0
16    a0 /= a0
17    return [b0,b1,b2],[a0,a1,a2]
  






高、低通滤波器
二次高、低通滤波器的滤波倾斜度为-12dB/oct ，只有一个参数：滤波频率f。之后你会遇到要求设计-24db/0ct,-48db/oct斜率的滤波器，以后会慢慢介绍。

Toggle Line Number
01    def lpf(f):
02    a = [1.0, math.sqrt(2), 1.0]                //可以理解为三个元素的数组a
03    b = [0.0, 0.0, 1.0]
04    ad = [0.0] * 3                                      //开三个元素的数组
05    bd = [0.0] * 3
06    fs = 44100.0
07    f = 2*math.pi*f
08    b[2] *= f*f
09    a[1] *= f
10    a[2] *= f*f
11    ad[0] = a[2] - a[1]*2*fs + 4*fs*fs*a[0]
12    ad[1] = a[2]*2 - 8*fs*fs*a[0]
13    ad[2] = a[2] + 2*fs*a[1] + 4*fs*fs*a[0]
14    bd[0] = b[2] - b[1]*2*fs + 4*fs*fs*b[0]
15    bd[1] = b[2]*2 - 8*fs*fs*b[0]
16    bd[2] = b[2] + b[1]*2*fs + 4*fs*fs*b[0]
17    t = ad[2]
18    return [bd[2]/t, bd[1]/t, bd[0]/t], [ad[2]/t, ad[1]/t, ad[0]/t]

Toggle Line Number
01    def hpf(f):
02    a = [1.0, math.sqrt(2), 1.0]
03    b = [0.0, 0.0, 1.0]
04    ad = [0.0] * 3
05    bd = [0.0] * 3
06    fs = 44100.0
07    f = 2*math.pi*f
08    b[0],b[2] = b[2],b[0]               //t=b[0];b[0]=b[2];b[2]=t;
09    a[1] *= f
10    a[2] *= f*f
11    ad[0] = a[2] - a[1]*2*fs + 4*fs*fs*a[0]
12    ad[1] = a[2]*2 - 8*fs*fs*a[0]
13    ad[2] = a[2] + 2*fs*a[1] + 4*fs*fs*a[0]
14    bd[0] = b[2] - b[1]*2*fs + 4*fs*fs*b[0]
15    bd[1] = b[2]*2 - 8*fs*fs*b[0]
16    bd[2] = b[2] + b[1]*2*fs + 4*fs*fs*b[0]
17    t = ad[2]
18    return [bd[2]/t, bd[1]/t, bd[0]/t], [ad[2]/t, ad[1]/t, ad[0]/t]









shelving滤波器   
Shelving滤波器提升或者抑制高频或者低频的增益，有三个参数： 中心频率f，增益db，类型type。type为0时，对低频进行增益处理，为1时对高频进行增益处理。

Toggle Line Number
01def shelving(f, db, type):
02    if type==1:
03        f = 22050 - f
04    Wb = 2.0*f/44100.0*math.pi
05    K = math.tan(Wb/2)
06    g = 10.0 ** (db/20.0)
07    V = g ** (1.0/2.0) - 1
08    a0 = 1 + 2*K + K*K
09    a1 = 2*K*K - 2
10    a2 = 1 - 2*K + K*K
11    b0 = a0
12    b1 = a1
13    b2 = a2
14    b0 = b0 + 2*V*K*(K+1)
15    b1 = b1 + 2*V*K*(2*K)
16    b2 = b2 + 2*V*K*(K-1)
17    b0 = b0 + V*V*K*K*1
18    b1 = b1 + V*V*K*K*2
19    b2 = b2 + V*V*K*K*1
20    b0 = b0 / a0
21    b1 = b1 / a0
22    b2 = b2 / a0
23    a1 = a1 / a0
24    a2 = a2 / a0
25    a0 = a0 / a0
26    if type==1:
27        b1 = -b1
28        a1 = -a1
29    return [b0,b1,b2],[a0,a1,a2]

  

二次IIR滤波器的传递函数如下，
               b0 + b1*z^-1 + b2*z^-2
H(z) = ------------------------------------- 
                a0 + a1*z^-1 + a2*z^-2

它看上去有6个参数，其实独立的只有5个，把上式中每个系数都除以a0，可以得到

             (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
H(z) = --------------------------------------------------------        注意：有些传递函数分母是1-a1*z^-1+a2*z^-2,所以大家要根据实际情况进行处理。
             1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

         在上一篇文章 二次IIR滤波器公式介绍中，最后结论都是5个参数。传递函数可以直接转换为如下的直接1型计算公式：

y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]        DSP处理应用的公式。

         这里的n表示当前时刻(取样点)，x是IIR滤波器的输入，y是输出，这个计算公式表明：n时刻的输出是n，n-1，n-2时刻的输入和n-1，n-2时刻的输出的线性组合。

         实现这个公式的C语言程序不会比公式本身长多少，所以在写具体的计算程序之前，我们还是来看看如何设计二次IIR均衡器的系数b0，b1，b2，a0，a1，a2。

         二次IIR均衡器有中心频率f0，增益dBgain和Q(Q值决定频率响应的宽窄)。

数字滤波器还有一个额外的参数：取样频率Fs。
给出上面的4个参数，就可以按照下面的公式求出b0，b1，b2，a0，a1，a2了。
A  = sqrt( 10^(dBgain/20) )
w0 = 2*pi*f0/Fs
alpha = sin(w0)/(2*Q)                                      
b0 =   1 + alpha*A
b1 =  -2*cos(w0)
b2 =   1 - alpha*A
a0 =   1 + alpha/A
a1 =  -2*cos(w0)
a2 =   1 - alpha/A        这个公式可以根据频率响应的特征列出方程组求得，由于二次IIR的频率响应的通用公式很容易通过其传递函数求得，因此只需要确定频率响应几个特征值，就可以确定传递函数的系数了。

       下 面我们来看看如何用C语言实现二次IIR滤波器。上文给出直接1型计算公式就是计算二次IIR滤波器全部算法了。不过在DSP具体实现中，为了能够实时地 处理声音信号，经常使用乒乓缓存和块处理等技术。无论是输入还是输出流，都有一个乒乓缓存与之对应。当系统正在写入或者读出乒缓存的时候，我们的程序就处 理乓缓存中的数据，当系统完成写入或者读出之后，就立刻转到乓缓存继续，而我们则开始对乒缓存处理。通过这样的交替，使得整个输入输出数据流可以不间断地 工作。所以为了实时处理数据的，必须一次能够处理一个乒缓存或者乓缓存，简单地说就是一次能处理一块数据。下面先给出二次IIR滤波器块处理程序，然后再 做分析。
void iir2(short *buf, float *pars, int size, int inc, float *prev_buf){
    int i,j;
    float out;
    float *t;
    short *p=buf;
    t=prev_buf;
    for(i=0;i
        t[0]=(float)(*p);
        out=0;
        for(j=0;j<5;j++)
            out+=t[j]*pars[j];
        if(out>0x7fff) *p=0x7fff;
        else
        if(out<0x8000) *p=0x8000;
        else
        *p=(short)out;
        t[2]=t[1];
        t[1]=t[0];
        t[4]=t[3];
        t[3]=out;
        p+=inc;
    }
}      由 于要计算当前块的第一个输出的时候，必须使用前一个块的信息，所以我们使数组prev_buf来保存前一个块的信息。它的元素分别是：x0，x1，x2， y1，y2。x0表示当前的输入值，x1表示上一个输入值，x2则是上上个输入值，y1是上个输出值，y2是上上个输出值。
Pars是IIR滤波器的系数数组，由程序
for(j=0;j<5;j++)
    out+=t[j]*pars[j];
可知，系数数组的值是[(b0/a0), (b1/a0), (b2/a0), -(a1/a0), -(a2/a0)]。

       Buf 是滤波器的输入数据，它的长度为size（样本数）。一般的声音都是16bit的整型数据，所以它的类型是short *，而系数和中间的计算结果都用float保存。在实际的DSP程序中，为了获得最佳的计算速度，系数和中间结果也都使用short保存，这样的话就需要 做定点小数运算了。

        最后，由于声音数据有stereo（双声道）和mono（单声道）之分，所以增加一个inc参数，这样如果是stereo输入的话，就把inc设置为2，只对一个声道的数据进行处理。