高阶IIR滤波的设计与实现

刚离开学校时发现在学校学的很多东西在工作中都用不上。随着时间的推移，当自己需要用这些东西时却发现当初所学的知识都已经还给了老师。 对我来说，数字滤波器的设计就是例子。

前段时间花了大概有一两个月的时间来重温当初在学校学习过的数字信号处理相关内容，从欧拉公式到传递函数等，一点一滴慢慢的捡起来。以下就是我这一两个月的工作成果^_^

由模拟滤波器设计数字IIR滤波器： 脉冲响应不变法、双线性变换法

我主要关注双线性变换法。

在选择模拟滤波器原型时，有如下几种选择： Butterworth/Chebyshev I/Chebyshev II/Elliptic/...

N阶IIR LPF：

首先，我们要找到模拟滤波器原型的幅度平方函数。这个在wiki上很容易找到。不要问我怎么来的，我也不知道。

以Chebyshev I为例：

* 1
* |H(jw)|^2 = --------------------------
* sqrt(1 + e^2*Tn^2(w/wc))

由于有s = jw

* 1
* H(s)*H(-s) = -----------------------
* 1 + e^2*Tn^2(-j*s/wc)

由上面的公式即可得到模拟滤波器的极点，别问我怎么解出来的，wiki上都有

* pk = (+-sinh(A)*sin(B) + jcosh(A)*cos(B))*wc
*
* = Rp + j*Ip
*
* 1 1
* where A = ---*arsinh(---),
* n e
*
* PI 2*k-1
* B = ---*-------, k = 1, 2, ... n/2
* 2 n

当然，上式中只显示的n/2个极点均处于实轴之上，别一半极点位于实轴之下，且二者成共轭。为了在代码中写注释方便，共轭我就写成了pk~

到这里，所有与模拟滤波器类型相关的部分都处理完了。接下来的处理对各种模拟滤波器原型都适用：

代码里面实现数字滤波器的关键是要把高阶IIR滤波器转换成1阶或2阶滤波器，这里可以使用共轭复数的加法和乘法可以消去虚部的特性让n阶IIR滤波器转换成n/2个2阶结。

继续算啦

* n/2 1
* H(s) = MUL -----------------------
* k=1 (1 - s/pk)*(1 - s/pk~)
*
* n/2 1
* = MUL ------------------------------------------
* k=1 1 - ((pk+pk~)/(pk*pk~))*s + s^2/(pk*pk~)
*
* n/2 1
* = MUL ---------------------------------------------
* k=1 1 - ((2*Rp)/(Rp^2+Ip^2))*s + s^2/(Rp^2+Ip^2)
*
* n/2 1
* = MUL -------------------
* k=1 1 + A1*s + A2*s^2

到这里就已经很简单了，接下来就是双线性变将模拟滤波器转换成数字滤波器：

* n/2 1 + 2*z^-1 + z^-2
* H(z) = MUL --------------------------------------------
* k=1 (1+A1+A2) + (2-2*A2)*z^-1 + (1-A1+A2)*z^-2

接下来用代码实现就easy啦。

LPF OK了，HPF用频率转换方法也很easy。BPF有些麻烦，需要用4阶结。

======================

接下来是一些概念：

自然频率： Hz

规一化的数字频率： 自然频率/采样率

数字频率： 2*PI*规一化的数字频率

模拟频率： tan(PI*规一化的数字频率)

给大家猜猜： 上面的wc是什么频率

=============================

最后是一点点小体会：

在做公式推导时，多用变量代替某一阶段的结果，这样写代码时清晰可见。而且当发现公式推导出错时，不会造成整个代码重写。当然，在做性能优化时又另当别论啦。