软件滤波方法浅谈 (2)

   使用 FIR 滤波有很多优点，比如，总是可以通过先构造滤波器的幅频特性，然后逆变换（比如利用Inverse Fourier Transform) 得到离散的冲击响应序列。不过问题是，如果FIR 冲击响应不收敛，或收敛时间很长，这就需要大量的运算。即便是收敛快的，简单的 sinc 函数，为了保证效果，至少需要 16 点以上的序列。为减少运算量，迭代就会被考虑，也就是说，滤波器的离散公式为：

        y(k+1) = a0 * y(k) + a1 * y(k-1) + … + an*y(k-n)

            + b0*x(k) + b1*x(k-1) + … + bm*x(k-m)                                            (1)

  对照 FIR 的离散公式:

         y(k+1) = b0*x(k) + b1*x(k-1) + … + bm*x(k-m) 

       where the impulse response series H = [b0, b1, ..., bm]  

 

由于 y 的不断迭代，公式 (1) 中的每个 y 包含了以前所有输入x 的信息，所以这种滤波器被称作无限冲击响应 (IIR) 滤波器  (Infinite Impulse Response Filter).

 

IIR的设计需要一定的理论知识，比如Laplace and Z transform. 好在一些强大的工具比如 matlab可以把我们弱智化, 只要我们画出幅频特性，计算机就会自动给出公式。不过，了解一些理论知识，总是没有坏处的。

 

先看看最简单，最常用的滤波器 ---- RC 滤波器, 不知道的RC 滤波器的工程师赶快买个电阻撞死算了。RC filter 的传递函数是：

          H(s) = Wc / (s + Wc)                                                                           (2)

         Where Wc = 1/RC is the desired bandwidth.

 

由于冲击信号的拉氏变换为1，所以一个线性系统的传递函数就是系统的冲击响应。如果要把公式 (2) 离散化，那么我们先看看熟悉的积分器离散公式，假设 Ts 是采样时间：

           y(k) = y (k-1) + (x(k) +x(k-1)) * Ts /2                                                  (3)                                               

把公式 (3) 变换到 z 域:

        Y = Y * z + (X + X *z) * Ts/2

        Y = Ts/2 * (1+z)/(1-z) * X

则积分器的z域传递函数是:

        H(z) = Ts/2 * (1+z)/(1-z)                                                                      (4)

       

 而积分器的 s域传递函数是H(s) = 1/ s, 所以可以得出

         s = 2/Ts * (1-z)/(1+z)                                                                         (5)

 

把公式 (5) 带入公式 (2)，可得:

          H(z) = Wc*Ts/(2+Wc*Ts)*(1 + z) / (z – (2-Wc*Ts)/(2+Wc*Ts))            (6)

 

这样，就可从(6) 中得到 RC 滤波器的迭代公式：

         y(k) = (2-Wc*Ts)/(2+Wc*Ts) * y(k-1)

              + Wc*Ts/(2+Wc*Ts)* (x(k) + x(k-1))                                              (7)

由于 Wc, Ts 都是常数，公式 （7）的运算量并不大。

由上可以看出，只要有传递函数， IIR 的迭代公式推导并不难。 许多滤波器比如 chebyshov，butterworth，相关滤波器书上都有传递函数，通过 matlab 还是自己动手推导，都不是难题，就不啰嗦了。