16、无限脉冲响应（IIR）与有限脉冲响应（FIR）滤波器的实现

无限脉冲响应（IIR）与有限脉冲响应（FIR）滤波器的实现

1. IIR滤波器的实现

1.1 并行实现示例

在一个并行实现示例中，有如下关系：
$w_3(n) = x(n) - 0.9w_3(n - 1)$
$y_3(n) = -0.21w_3(n)$
输出为$y(n) = y_1(n) + y_2(n) + y_3(n)$，其并行实现如图14.1所示。

1.2 转置直接形式I

转置定理可将直接形式I实现转换为转置IIR滤波器，如图14.12所示。在转置IIR滤波器中，可以注意到信号流、加法节点和分支节点的变化。当设置$a_i = 0$（$i = 1,2, \ldots, M - 1$）时，可得到FIR滤波器的转置结构。

直接形式和转置形式的IIR滤波器对量化误差都非常敏感，因为反馈部分的量化误差会在滤波器中反馈和累积，当滤波器阶数较高时，这种情况会变得更加严重。为了减少量化误差，一种更稳健的结构是将高阶IIR滤波器分解为二阶IIR滤波器（或双二阶）的级联。级联$K$个双二阶形成$2K$阶IIR滤波器的一般表达式为：
$H(z) = G\prod_{k = 1}^{K}\frac{b_{k0} + b_{k1}z^{-1} + b_{k2}z^{-2}}{1 + a_{k1}z^{-1} + a_{k2}z^{-2}}$
其中$G$是增益。例如，级联两个双二阶形成四阶IIR滤波器的情况如图14.12所示。级联单个一阶IIR滤波器与一系列滤波器转置也可称为流图反转，对单输入单输出（SISO）滤波器进行转置不会改变其传递函数，这可以由梅森增益公式或特勒根定理推导得出。转置直接形式I实现是一个二阶IIR数字滤波器，注意输入信号从右侧输入，输出在左侧。

1.3 转置直接形式II

图14.13展示了转置直接形式II结构。为了便于与原始结构进行比较，输入和输出信号保持“交换”，因此信号通常从右向左流动，而不是通常的从左向右。

下面通过几个例子来说明如何求$H(z)$：
- 例14.2 ：对于差分方程$y[n] = x[n] + x[n - 1]$，其脉冲响应为${h[n]} = {\ldots, 0, 1, 1, 0, \ldots}$，则$H(z) = \sum_{n = -\infty}^{\infty}h[n]z^{-n} = \sum_{n = 0}^{1}h[n]z^{-n} = 1 + z^{-1}$。
- 例14.3 ：对于递归差分方程$y[n] = a_0x[n] + a_1x[n - 1] - b_1y[n - 1]$，若$x[n] = z^n$，则$y[n] = H(z)z^n$，$y[n - 1] = H(z)z^{n - 1}$，代入差分方程可得：
$H(z)z^n = a_0z^n + a_1z^{n - 1} - b_1H(z)z^{n - 1}$
整理得$H(z) = \frac{a_0 + a_1z^{-1}}{1 + b_1z^{-1}}$（$z \neq -b_1$），当$z = -b_1$时，$H(z) = \infty$。
对于一般数字滤波器，其差分方程对应的$H(z)$为：
$H(z) = \frac{a_0 + a_1z^{-1} + a_2z^{-2} + \cdots + a_Nz^{-N}}{b_0 + b_1z^{-1} + b_2z^{-2} + \cdots + b_Mz^{-M}}$（$b_0 = 1$）

1.4 二阶数字滤波器的信号流图

对于二阶数字滤波器$H(z) = \frac{a_0 + a_1z^{-1} + a_2z^{-2}}{1 + b_1z^{-1} + b_2z^{-2}}$，其差分方程为$y[n] = a_0x[n] + a_1x[n - 1] + a_2x[n - 2] - b_1y[n - 1] - b_2y[n - 2]$，其信号流图如图14.14所示，这被称为“直接形式I”的双二阶部分。通过重新排列信号流图的两个部分，可以得到图14.15，它与图14.14具有相同的脉冲响应。进一步简化可得图14.16，这是“直接形式II”的双二阶部分，它具有最少的延迟框数量，被称为“规范”形式，其系统函数与“直接形式I”的信号流图相同，因此可以实现任何二阶双二阶系统函数。

1.5 用MATLAB实现陷波滤波器

假设要设计一个四阶“陷波”数字滤波器，以消除800Hz的不需要的正弦波，同时不严重影响信号的其余部分，采样率为$F_S = 10kHz$。可以使用MATLAB的“butter”函数，代码如下：

FS = 10000;
FL = 775;
FU = 825;
[a b] = butter(2, [FL FU]/(FS/2),'stop');
a = [0.98 -3.43 4.96 -3.43 0.98];
b = [1 -3.47 4.96 -3.39 0.96];
freqz(a, b);
freqz(a, b, 512, FS);
axis([0 FS/2 -50 5]);

最后两条MATLAB语句产生的频率响应（增益和相位）如图14.17所示。由于巴特沃斯带阻滤波器在两个截止频率$F_L = 775$和$F_U = 825$处的增益为 -3dB，因此陷波的“ -3dB频率带宽”为$25 + 25 = 50Hz$。

四阶数字滤波器的传递函数为：
$H(z) = \frac{0.98 - 3.43z^{-1} + 4.96z^{-2} - 3.43z^{-3} + 0.98z^{-4}}{1 - 3.47z^{-1} + 4.96z^{-2} - 3.39z^{-3} + 0.96z^{-4}}$
“直接形式II”实现的四阶陷波滤波器的信号流图如图14.18所示。但阶数大于二的“直接形式”IIR实现很少使用，因为系数值的舍入误差敏感性会很高，并且$z^{-1}$框中的“中间”信号范围会很大。在定点运算中会出现很大的困难，因此通常使用“级联双二阶部分”。可以使用以下代码将四阶传递函数$H(z)$转换为新形式：

[a b] = butter(2, [FL FU]/(FS/2),'stop');
[SOS G] = tf2sos(a,b)

MATLAB响应如下：

SOS =
    1 -1.753 1 1 -1.722 0.9776
    1 -1.753 1 1 -1.744 0.9785
G = 0.978

在MATLAB中，“SOS”代表“二阶部分”（即双二阶），函数“tf2SOS”将数组“a”和“b”中的系数转换为存储在数组“SOS”中的新系数集和常数$G$。数组SOS有两行，一行对应第一个双二阶部分，另一行对应第二个双二阶部分。在每一行中，前三个项指定非递归部分，后三个项定义递归部分。因此$H(z)$现在可以如图14.21所示实现。

1.6 无限脉冲响应滤波器的实现方法

设计数字IIR滤波器有两种方法：
- 模拟到数字滤波器设计 ：首先基于常用的模拟滤波器（如巴特沃斯、切比雪夫I和II以及椭圆滤波器）设计模拟低通滤波器，然后进行以下两种操作之一：
- 在模拟域中进行频率变换，然后应用模拟到数字（或s平面到z平面）映射。
- 应用模拟到数字映射，然后在数字域中进行频率变换。频率变换是使用一些映射函数将低通滤波器转换为具有不同截止频率的另一个低通、带通、高通或带阻滤波器的过程，该变换可以在模拟或数字域中进行。
- 双线性变换 ：一种更好且更常用的s平面到z平面映射是双线性变换，定义为：
$S = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}$
双线性变换是一个有理函数，它将s平面的左半部分映射到z平面单位圆的内部，是一个一对一的映射，将沿$j\Omega$轴从$-\infty$到$\infty$的整个频率仅映射到单位圆（$-\pi$到$\pi$）一次。模拟频率和数字频率之间的非线性关系表示为：
$\omega = 2\tan^{-1}(\frac{\Omega T}{2})$

以下是计算低通滤波器（LPF）系数并绘制响应的代码：

clc;
clear all;
close all;
rp = 0.2; % 通带波纹
rs = 50; % 阻带波纹
wp = 1700; % 通带频率
ws = 3400; % 阻带频率
fs = 8000; % 采样频率
w1 = 2*wp/fs;
w2 = 2*ws/fs;
[n,wn] = buttord(w1,w2,rp,rs,'s');
[b,a] = butter(n,wn,'low','s');
w = 0:.1:pi;
[h,om] = freqs(b,a,w);
m = 20*log10(abs(h));
an = angle(h);
figure
plot(om/pi,m);
title('IIR - LPF response ');
xlabel(' Normalized freq.');
ylabel('Gain (dB)');

图14.22显示了使用“butter”命令的LPF响应。

以下是计算高通滤波器（HPF）系数并绘制响应的代码：

clc;
clear all;
close all;
rp = 0.2; % 通带波纹
rs = 100; % 阻带波纹
wp = 1500; % 通带频率
ws = 3000; % 阻带频率
fs = 10000; % 采样频率
w1 = 2*wp/fs;
w2 = 2*ws/fs;
[n,wn] = buttord(w1,w2,rp,rs,'s');
[b,a] = butter(n,wn,'high','s');
w = 0:.01:pi;
[h,om] = freqs(b,a,w);
m = 20*log10(abs(h));
an = angle(h);
figure
plot(om/pi,m);
title('IIR-HPF response');
xlabel('Normalized freq.');
ylabel('Gain (dB)');

图14.23显示了使用“butter”命令的HPF响应。

1.7 相关问题

以下是一些相关问题及要求：
1. 给定二阶传递函数$H(z) = \frac{5(1 - z^{-1})}{1 + z^{-1} + 3.36z^{-2}}$，使用以下实现方式进行滤波器实现并写出差分方程：
- 直接形式I和直接形式II
- 通过一阶部分的级联形式
- 通过一阶部分的并行形式
2. 给定二阶传递函数$H(z) = \frac{0.5(1 - z^{-1})}{1 + 3z^{-1} + 2.5z^{-2}}$，使用以下实现方式进行滤波器实现并写出差分方程：
- 直接形式I和直接形式II
- 通过一阶部分的级联形式
- 通过一阶部分的并行形式
3. 求差分方程$y[n] = x[n] + x[n - 1]$的$H(z)$。
4. 求差分方程$y[n] = x[n] + x[n - 1]$的$H(z)$。
5. 求递归差分方程$y[n] = x[n] + x[n - 1] - b_1y[n - 1]$的$H(z)$。
6. 求递归差分方程$y[n] = 2a_0x[n] - a_1x[n - 1] + b_1y[n - 1]$的$H(z)$。
7. 给出二阶数字滤波器$H(z) = \frac{2 + 3z^{-1} + 4z^{-2}}{5 + 6z^{-1} + 7z^{-2}}$的信号流图。

2. FIR滤波器的实现与设计

2.1 有限脉冲响应滤波器表示

在数字信号处理中，有限脉冲响应（FIR）滤波器是一种滤波器，其脉冲响应或对任何有限长度输入的响应是有限持续时间的，因为它在有限时间内稳定到零。FIR滤波器的极点和零点数量相同，但所有极点都位于原点，由于所有极点都位于单位圆内，因此FIR滤波器显然是稳定的。

FIR滤波器完全由以下差分方程（输入 - 输出关系）指定：
$y(n) = b_0x(n) + b_1x(n - 1) + \cdots + b_Mx(n - M) = \sum_{i = 0}^{M}b_ix(n - i)$

理想低通滤波器的频率响应如图15.1a所示，数学表达式为：
$H(e^{j\Omega}) = \begin{cases}1, & 0 \leq \Omega \leq \Omega_c \ 0, & \Omega_c \leq \Omega \leq \pi\end{cases}$
其中$\Omega_c$是低通截止频率。对应的理想低通滤波器的脉冲响应如图15.1b所示，从$n = -M$到$n = M$的截断部分（实际响应在两侧延伸到无穷大）为：
$h(n) = \begin{cases}\frac{\Omega_c}{\pi}, & n = 0 \ \frac{\sin(\Omega_c n)}{\pi n}, & -M \leq n \leq M, n \neq 0\end{cases}$
脉冲响应关于$n = 0$对称。在这种情况下，脉冲响应的z变换为：
$H(z) = h(-M)z^M + h(-M + 1)z^{M - 1} + \cdots + h(0) + h(1)z^{-1} + \cdots + h(M)z^{-M}$
为了获得因果脉冲响应，将非因果脉冲响应向右移动$M$个样本，得到因果理想低通FIR滤波器的传递函数：
$H(z) = b_0 + b_1z^{-1} + b_2z^{-2} + \cdots + b_{2M}z^{-2M}$
其中$b_n = h(n - M)$（$n = 0, 1, 2, \ldots, 2M$）。

2.2 不同类型FIR滤波器的设计方程

还可以获得其他类型FIR滤波器（如高通、带通和带阻）的设计方程，表15.1给出了这些类型FIR滤波器的滤波器系数计算公式。

滤波器类型	理想脉冲响应$h(n)$（非因果FIR系数）
低通	$h(n) = \begin{cases}\frac{\Omega_c}{\pi}, & n = 0 \ \frac{\sin(\Omega_c n)}{\pi n}, & -M \leq n \leq M\end{cases}$
高通	$h(n) = \begin{cases}\frac{\pi - \Omega_c}{\pi}, & n = 0 \ -\frac{\sin(\Omega_c n)}{\pi n}, & -M \leq n \leq M\end{cases}$
带通	$h(n) = \begin{cases}\frac{\Omega_H - \Omega_L}{\pi}, & n = 0 \ \frac{\sin(\Omega_H n) - \sin(\Omega_L n)}{\pi n}, & -M \leq n \leq M\end{cases}$
带阻	$h(n) = \begin{cases}\frac{\pi - (\Omega_H + \Omega_L)}{\pi}, & n = 0 \ -\frac{\sin(\Omega_H n) + \sin(\Omega_L n)}{\pi n}, & -M \leq n \leq M\end{cases}$

2.3 FIR低通滤波器示例

2.3.1 计算滤波器系数

假设要设计一个3抽头的FIR低通滤波器，截止频率为1200Hz，采样率为12000Hz。首先确定归一化截止频率：
$\Omega_c = 2\pi f_cT_s = 2\pi\times\frac{1200}{12000} = 0.2\pi$ 弧度
由于$2M + 1 = 3$，所以$M = 1$。根据表7.1可得：
$h(0) = \frac{\Omega_c}{\pi} = \frac{0.2\pi}{\pi} = 0.2$
$h(1) = \frac{\sin(\Omega_c\times1)}{\pi\times1} = \frac{\sin(0.2\pi)}{\pi} = 0.1871$
利用对称性，$h(-1) = h(1) = 0.1871$。
延迟$h(n)$ 1个样本，得到滤波器系数：
$b_0 = h(0 - 1) = h(-1) = 0.1871$
$b_1 = h(1 - 1) = h(0) = 0.2$
$b_2 = h(2 - 1) = h(1) = 0.1871$

2.3.2 确定传递函数和差分方程

传递函数为：
$H(z) = b_0 + b_1z^{-1} + b_2z^{-2} = 0.1871 + 0.2z^{-1} + 0.1871z^{-2}$
差分方程为：
$y(n) = 0.1871x(n) + 0.2x(n - 1) + 0.1871x(n - 2)$

2.3.3 计算并绘制频率响应

频率响应为：
$H(e^{j\Omega}) = 0.1871 + 0.2e^{-j\Omega} + 0.1871e^{-2j\Omega}$
可改写为：
$H(e^{j\Omega}) = e^{-j\Omega}(0.1871e^{j\Omega} + 0.2 + 0.1871e^{-j\Omega}) = e^{-j\Omega}(0.2 + 0.1871\times2\cos(\Omega)) = e^{-j\Omega}(0.2 + 0.3742\cos(\Omega))$
幅度频率响应为：
$|H(e^{j\Omega})| = 0.2 + 0.3742\cos(\Omega)$
相位响应为：
$\angle H(e^{j\Omega}) = \begin{cases}-\Omega, & 0.2 + 0.3742\cos(\Omega) > 0 \ -\Omega + \pi, & 0.2 + 0.3742\cos(\Omega) < 0\end{cases}$

在具有对称系数的FIR滤波器中，相位响应在通带内具有线性行为，这意味着滤波器输入在通带内的所有频率分量在滤波器输出处都受到相同的时间延迟，这是音频和语音滤波应用的要求，因为需要避免相位失真。如果设计方法不能产生对称系数，则所得的FIR滤波器不具有线性相位特性，会导致滤波信号出现失真。

2.4 FIR带通滤波器示例

2.4.1 计算滤波器系数

设计一个5抽头的FIR带通滤波器，下截止频率为2000Hz，上截止频率为2400Hz，采样率为8000Hz。首先确定归一化截止频率：
$\Omega_L = 2\pi f_LT_s = 2\pi\times\frac{2000}{8000} = 0.5\pi$ 弧度
$\Omega_H = 2\pi f_HT_s = 2\pi\times\frac{2400}{8000} = 0.6\pi$ 弧度
由于$2M + 1 = 5$，所以$M = 2$。根据表7.1可得：
$h(0) = \frac{\Omega_H - \Omega_L}{\pi} = \frac{0.6\pi - 0.5\pi}{\pi} = 0.1$
$h(1) = \frac{\sin(\Omega_H\times1) - \sin(\Omega_L\times1)}{\pi\times1} = \frac{\sin(0.6\pi) - \sin(0.5\pi)}{\pi} = -0.01558$
$h(2) = \frac{\sin(\Omega_H\times2) - \sin(\Omega_L\times2)}{\pi\times2} = \frac{\sin(1.2\pi) - \sin(\pi)}{\pi\times2} = -0.09355$
利用对称性，$h(-1) = h(1) = -0.01558$，$h(-2) = h(2) = -0.09355$。
延迟$h(n)$ 2个样本，得到滤波器系数：
$b_0 = h(0 - 2) = h(-2) = -0.09355$
$b_1 = h(1 - 2) = h(-1) = -0.01558$
$b_2 = h(2 - 2) = h(0) = 0.1$
$b_3 = h(3 - 2) = h(1) = -0.01558$
$b_4 = h(4 - 2) = h(2) = -0.09355$

2.4.2 确定传递函数和差分方程

传递函数为：
$H(z) = b_0 + b_1z^{-1} + b_2z^{-2} + b_3z^{-3} + b_4z^{-4} = -0.09355 - 0.01558z^{-1} + 0.1z^{-2} - 0.01558z^{-3} - 0.09355z^{-4}$
差分方程为：
$y(n) = -0.09355x(n) - 0.01558x(n - 1) + 0.1x(n - 2) - 0.01558x(n - 3) - 0.09355x(n - 4)$

2.4.3 计算频率响应

频率响应为：
$H(e^{j\Omega}) = -0.09355 - 0.01558e^{-j\Omega} + 0.1e^{-2j\Omega} - 0.01558e^{-3j\Omega} - 0.09355e^{-4j\Omega}$

FIR滤波器具有良好的相位特性，但幅度频率响应可能不理想。在幅度频率响应的通带（主瓣）和阻带（旁瓣）中会出现振荡（波纹），这是由于低通滤波器的无限脉冲响应的突然截断引起的吉布斯振荡。可以通过加窗来避免这种行为。

2.5 FIR滤波器的优缺点总结

2.5.1 优点

• 稳定性 ：由于所有极点都位于原点，FIR滤波器天然稳定，不会出现因极点位置不当导致的不稳定问题，这在许多对稳定性要求较高的应用中非常重要。
• 线性相位特性 ：当滤波器系数对称时，FIR滤波器在通带内具有线性相位特性，这意味着所有频率分量在滤波器输出处受到相同的时间延迟，避免了相位失真，适用于音频和语音滤波等对相位敏感的应用。
• 易于设计 ：可以通过简单的数学方法（如傅里叶变换）来设计FIR滤波器，并且可以根据需要灵活调整滤波器的阶数和截止频率等参数。

2.5.2 缺点

• 幅度频率响应不理想 ：FIR滤波器的幅度频率响应可能会出现振荡（波纹），特别是在通带和阻带中，这是由于无限脉冲响应的突然截断引起的吉布斯振荡。
• 需要较高的阶数 ：为了达到较好的滤波效果，FIR滤波器通常需要较高的阶数，这会增加计算复杂度和存储需求。

2.6 FIR滤波器设计流程总结

下面是设计FIR滤波器的一般流程：
1. 确定滤波器类型和参数 ：根据具体应用需求，确定滤波器的类型（如低通、高通、带通、带阻）以及相关参数（如截止频率、采样率等）。
2. 计算归一化截止频率 ：将实际的截止频率转换为归一化频率，以便后续计算。
3. 计算滤波器系数 ：根据滤波器类型和归一化截止频率，使用相应的公式计算滤波器系数。
4. 确定传递函数和差分方程 ：根据滤波器系数，确定滤波器的传递函数和差分方程。
5. 计算并绘制频率响应 ：计算滤波器的频率响应（幅度和相位），并绘制响应曲线，以验证滤波器的性能。
6. 优化滤波器性能 ：如果幅度频率响应不理想，可以考虑使用加窗等方法来优化滤波器性能。

2.7 FIR滤波器设计示例流程图

graph TD
    A[确定滤波器类型和参数] --> B[计算归一化截止频率]
    B --> C[计算滤波器系数]
    C --> D[确定传递函数和差分方程]
    D --> E[计算并绘制频率响应]
    E --> F{性能是否满足要求}
    F -- 是 --> G[完成设计]
    F -- 否 --> H[优化滤波器性能]
    H --> C

3. IIR与FIR滤波器的比较

3.1 性能比较

比较项目	IIR滤波器	FIR滤波器
稳定性	可能不稳定，需要仔细设计极点位置	天然稳定，所有极点位于原点
相位特性	一般为非线性相位，可能导致相位失真	可以实现线性相位，避免相位失真
幅度频率响应	可以用较低的阶数实现较好的幅度频率响应	通常需要较高的阶数才能达到类似的幅度频率响应
计算复杂度	由于存在反馈环节，计算复杂度相对较低	没有反馈环节，计算复杂度相对较高
存储需求	相对较低	由于阶数较高，存储需求相对较高

3.2 应用场景比较

• IIR滤波器 ：适用于对相位要求不高，但对幅度频率响应要求较高，且希望使用较低阶数滤波器的应用场景，如音频和视频信号的滤波、通信系统中的信道均衡等。
• FIR滤波器 ：适用于对相位要求严格，需要避免相位失真的应用场景，如音频和语音处理、生物医学信号处理、雷达信号处理等。

3.3 选择建议

在选择使用IIR滤波器还是FIR滤波器时，可以考虑以下因素：
1. 相位要求 ：如果对相位失真敏感，应优先选择FIR滤波器；如果对相位要求不高，则可以考虑IIR滤波器。
2. 幅度频率响应要求 ：如果需要用较低的阶数实现较好的幅度频率响应，IIR滤波器可能更合适；如果对阶数没有严格限制，FIR滤波器也可以满足要求。
3. 计算资源和存储需求 ：如果计算资源有限或存储容量较小，IIR滤波器可能更适合；如果计算资源充足且对存储需求不敏感，FIR滤波器是一个不错的选择。

4. 总结

本文详细介绍了无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器的实现方法、设计流程以及它们的优缺点和应用场景。IIR滤波器通过反馈环节可以用较低的阶数实现较好的幅度频率响应，但可能存在稳定性和相位失真问题；FIR滤波器具有天然的稳定性和线性相位特性，但通常需要较高的阶数才能达到类似的幅度频率响应。在实际应用中，应根据具体需求和条件选择合适的滤波器类型。

同时，本文还给出了多个设计示例和相关代码，帮助读者更好地理解和掌握IIR和FIR滤波器的设计方法。通过合理设计和优化滤波器，可以满足不同应用场景对滤波性能的要求。希望本文能为从事数字信号处理的读者提供有益的参考。