软件滤波方法浅谈
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
|
