谈谈软件滤波

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