软件滤波方法浅谈 (1)

 

21icbb 上的程序匠人总结了10种软件滤波方法，简单实用，看得出是匠人自己对实际应用的归纳升华。不过如果能从理论角度分析一下， 更易于理解，因为 10 种方法大多数可以归结为一种做法：有限冲击响应滤波器 (FIR Filter) ，或者通俗的说是程序匠人所提到的“加权平均滤波器“. 

  

滤波器可以分为两种：IIR (Infinite Impulse Response) and FIR (Finite Impulse Response), 这里先谈谈 FIR. 首先介绍一下卷积 convolution, 原理很简单：如果一个线性系统（这里是滤波器）的冲击响应是H(t), 那么对于任何输入 x(t) , 可以通过convolution 计算出y(t). 如果在离散系统，就是由冲击响应序列：

   H = [h(0), h(1), h(2), …., h(N-1)]                                      (1)

以及输入序列：

       X = [x(k),   x(k-1) , …,   x(k-N-1) ]                                        (2)

 可以求出：

       Y(k) = Cov(H, X) = h(0)*x(k) + h(1)*x(k-1) + …. + h(N-1)*x(k-N-1)       (3)

 

公式 （3）即便用汇编 （TI DSP 甚至有专门的指令），在mcu, DSP, PC上的实现也都极其简单,  所以源代码就不提供了。现在软件滤波的设计问题变成了找到滤波的冲击响应序列的问题，或者说就是找到加权平均滤波器中的权重，如果从理论角度看 fir，那么要看看卷积一个极为重要的特性是：两个函数的时域的卷积，在频率域表现为对应的频率响应函数的乘积。即：

 y(t) = cov(x(t), h(t)) = inverseFrequency(X(f) * H(f))               (4)

  

这里的 X(f)*H(f) 就是滤波后的幅频特性曲线，因尔，可以构造一个合适的 H(f), 逆变换回 h(t), 就可以得到 fir.

假定H(f) 是标准的低通，如果用离散形式表示H(f), 就是：

 H = [1,  1,  ...,  1, 0,  0, ...,  0]                                            (5) 

我们可以推导出H(f) 逆变换的时域函数：

 h(t) = sin(t) / t                                                                    (6)

 

这个函数就是有名的 sinc 函数，又称窗口滤波器函数。

 

    再来看看极为常用的平均滤波公式：

        y(k) = ( x(k) + x(k-1) + … + x(k-N-1) ) / N                            (7)

 

   公式（7）可以导出平均滤波的冲击响应序列：

        H = [1,   1,   1, … , 1] / N                                                 (8)    

 

    用 matlab 画出 bode chart (幅频响应曲线)，那么可以清楚地看出是一个低通滤波器. 如果不使用对数坐标，进一步观察，可以发现幅频特性为：

       y(f) = sin(f) / f                                                                   (9)

 

 

 注意： 为了保证 fir 的增益为 1， fir 对输入各项加权后必须平均，即所有系数之和为 1

   K = (h(0) + h(1) + ... + h(N-1))

  

冲击响应序列为：

   H = [h(0), h(1), h(2), …., h(N-1)] / K     

 

 

看看一个例子：用 FIR 设计 RC 滤波器, 注意，这只是演示FIR过程，实际应用中 RC 滤波器不能真正的使用 FIR.

一个 RC 电路的冲击响应为：

 　　　　 h(t) = exp(-t / RC)                                                    (10)

 

公式 (10) 应该不会陌生，冲击响应离散公式为：

                 H(n) = exp( - n * Ts / RC)                                        (11)

 

假设　RC = 1 (R = 100K, C = 10 uF)，采样时间是 0.5 s, 则可以获得冲击响应序列：

           H = [h(0),  h(1),  h(2), …., h(23)] =

           [ 1.000000,   0.606531,  0.367879, 0.223130, 0.135335, 0.082085, 0.049787, ….. , 0.00001]

          

           K = h(0) + h(1) + ... + h(23) = 2.541478

 

最后我们得到了最终的 RC 滤波器的加权平均系数（冲击响应序列）：

           H = [0.393472， 0.238653， 0.144750， 0.087795，...., 0.000004]