温故而知新：IIR滤波器设计的方法，幅频计算和参数理解

LTI系统、$\delta$、差分方程和FIR、IIR的由来

“学而时习之、不亦乐乎… 温故而知新、可以为师矣。”很多经典的知识理论，每次复习回顾都会有新的收获和体会。

LTI线性时不变系统

顾名思义，LTI系统线性时不变两个特性，这令系统设计比较自适应和匹配滤波器等复杂方法简单多了。那么最这个系统最准确的（特性）抽象定义就是，当单位脉冲响应函数$\delta (n)$作为系统的输入，输出信号$y(n)$就是这个系统的响应，即$h(n)$。这里需要强调的是，$\delta (n)$虽然只是一个单位冲击，但$y(n)$有可能是一串响应，像“阿拉伯的故事，伯的故事，的故事”，有点余音绕梁的感觉。 从数学的角度，可以看成系统对冲击输入，有不同路径的延时加权，然后累加后输出一个结果，有些路径0延时，有些路径延时比较长。数字信号系统通常用时钟（节拍）来衡量这个延时，有的一拍，有的很多拍，下图是借用z变换延时算子来表示的LTI系统。

一般来说，这个延时系统中，延时久的抽头加权值越小，犹如实际中回音会越来越小一样。相反的话系统就不收敛了，导致系统最后会爆炸，犹如核子响应，这是无法控制的，暂且不去深究了。

用差分方程来理解FIR的意义

  \space  对LTI最直观的表征就是差分方程，这很好理解，软件也很好实现。上图的例子可以用$$y(n)=ax(n)+bx(n-1)+cx(n-2)$$来表述。很显然这是一个FIR滤波器，因为冲激响应经过滤波器后，最多能影响两个延时，再多的延时冲击对系统来说已经熟视无睹了，这就是所谓的有限冲激响应滤波器FIR的由来。当然FIR滤波器还有很多特性，不是本篇要展开的。

IIR 无限冲激带来的诱惑和危机

无限冲激就是说$\delta (n)$在系统里一直存在，就像你感染了某些病毒，刚开始的时候表现的感冒症状很严重（瞬态冲击），但即使症状痊愈了，他们在你身体里还存在，只能选择和他们并存，有朝一日卷土重来未可知。IIR系统也确实如此，它对系统的改变会很迅猛，同时因为在系统中一直有残留响应，也为系统的不稳定性埋下伏笔。举例来说，对前面的差分增加一个递归抽头，$$y(n)=ax(n)+bx(n-1)+cx(n-2)+dy(n-1)$$
很显然，一旦系统受到$\delta (n)$冲激，那么只要$d$不为0，这种冲激随着时间推移，会一直影响系统输出$y(n)$，$d$这个权系数很大，或者与其他权系数的关系设计不合理，系统一工作，就会出现非常大的响应结果（不收敛、发散），对于数字设计来说，大部分情况是灾难，数字信号处理要寻找线性时不变的稳定系统。IIR虽然有风险，但它在运算效率和系统响应方面的优势又令人向往，所以IIR一直在系统中扮演重要的角色。

IIR系统设计方法

所谓的滤波器设计，实质上就是寻找一组系数，与系统输入（iir还要加上输出）做卷积（乘累加），常见的设计方法是脉冲响应不变方法和双线性变换方法。但两种方法的基础都是基于模拟滤波器的laplace变换而来，所以模拟原型滤波器也是设计数字滤波器的基础。参考文档里罗列了巴特沃兹到贝塞尔等滤波器的介绍博客，他们为IIR设计提供了原型，但如果需要灵活改变参数，实时设计的话，还需要对滤波器的数学公式和滤波器的指标参数有更深刻的理解。引用数字信号处理 原萍 著 中国铁道出版社文中的一段话，从原型模拟滤波器到数字滤波器的设计，就是寻找s域到z域的映射：
（1）H（z）的频率响应要能模仿H（s）的频率响应，即s平面的虚轴应映射到z平面的单位圆上。（2）H（s）的因果稳定性映射成H（z）后保持不变，即s平面的左半平面应映射到z平面的单位圆以内。
我觉得作者这本书对于IIR设计部分理解的很深入，对其中的butterworth滤波器部分摘抄整理一下，加深理解。butterworth低通滤波器的经典幅频表达：$$|H(j\Omega ){|}^2=\frac{1}{1+(\frac{\Omega }{{\Omega }_c}{)}^{2N}},N=1,2,3,....$$
或者$$|H(j\Omega )|=\frac{1}{\sqrt{1+(\frac{\Omega }{{\Omega }_c}{)}^{2N}}},N=1,2,3,....$$
当$\Omega ={\Omega }_c$的时候，幅度为$\Omega =0$时的$1/\sqrt{2}$，即所谓的$-3db$截止频率。如何根据这个幅频表达推算转换函数是一个有挑战的事情，也是模拟滤波器设计的核心。长话短说，令$s/\Omega$代替$s$，取butterworth圆上复平面的极点（另一侧为共轭）作为传递函数极点，得到$$H_a(s)=\frac{{\Omega }_c^N}{{\prod }_{i=0}^{N-1}(s-s_i)}$$
N=3的时候，得出三阶滤波器原型：
$$H_a(s)=\frac{1}{(\frac{s}{{\Omega }_c}{)}^3+2(\frac{s}{{\Omega }_c})+(\frac{s}{{\Omega }_c})+1}$$
而对于获取滤波器阶数的办法，利用了一组不等式，即通带内衰减不能小于${\alpha }_p$，阻带内的纹波不能大于${\alpha }_s$，在通带和阻带上用（不）等式来表示：
α p ⩾ − 20 l o g ∣ H ( j Ω ) ∣ ⩾ − 10 l o g 1 1 + ( Ω p Ω c ) 2 N αp⩾−20log|H(jΩ)|⩾−10log11+(ΩpΩc)2N

αp​​⩾−20log∣H(jΩ)∣⩾−10log1+(Ωc​Ωp​​)2N1​​
$$\begin{array}{rl}{\alpha }_s& ⩽-20log|H(j\Omega )|\\ & ⩽-10log\frac{1}{1+(\frac{{\Omega }_s}{{\Omega }_c}{)}^{2N}}\end{array}$$
两个不等式做消减，令${\Omega }_c$消去，最后得出
$$(\frac{{\Omega }_p}{{\Omega }_s}{)}^{2N}⩽\frac{10^{-0.1{\alpha }_p-1}-1}{10^{-0.1{\alpha }_s-1}-1}$$
推导出$$N⩾\frac{log[(10^{-0.1{\alpha }_p-1}-1)/(10^{-0.1{\alpha }_s-1}-1)]}{2log\frac{{\Omega }_p}{{\Omega }_s}}$$这个结果我推导不出来（有些疑问），暂且按照原文的意思，主要为了理解如何求参，有了它和通带阻带要求，就可以估算需要几阶滤波了，然后根据之前的公式计算：
$${\Omega }_c=\frac{{\Omega }_p}{(10^{-0.1{\alpha }_p-1}-1{)}^{1/2N}}$$或者$${\Omega }_c=\frac{{\Omega }_s}{(10^{-0.1{\alpha }_s-1}-1{)}^{1/2N}}$$
而通过查表可以很容易得到归一化$p=s/\Omega$巴特沃兹原型滤波器。

那么只要求出截止频率${\Omega }_c$，即可反向求出你所需要的滤波器原型了。在The Butterworth Low-Pass Filter一文中讲到了matlab计算butterworth原型滤波器的方法，该文也绘制了不同阶数的幅频响应。

条条大路通罗马，就看哪条最适合你了。

参考文档

Signal Processing Fundamentals, implementations and applications
数字信号处理 原萍 著 中国铁道出版社
巴特沃斯滤波器详解
切比雪夫滤波器
椭圆滤波器
贝塞尔滤波器
IIR-冲激响应不变法各个步骤详解
Understanding Butterworth Filter Poles and Zeros
Butterworth Filter Construction along with its Applications
A BUTTERWORTH-FILTER COOKBOOK
The Butterworth Low-Pass Filter
hebyshev Filters
Matlab-style IIR filter design