数字信号处理入门:数字滤波器设计方法与 IIR 滤波器设计

从数字滤波器的基本概念出发,系统整理技术指标、巴特沃斯模拟低通设计,以及脉冲响应不变法和双线性变换法两条 IIR 设计主线。

数字信号处理入门:数字滤波器设计方法与 IIR 滤波器设计

本文面向 0 基础读者,按照讲义顺序整理“数字滤波器设计方法与 IIR 滤波器设计”。正文已经去掉原 PPT 页码、模板页标题,改用博客式标题组织。为了便于发布,部分 PPT 中的示意图被重新绘制为配图,并集中放在 images/ 文件夹中。

写在前面:本讲到底解决什么问题

数字滤波器是数字信号处理里最典型、最常用的工程工具。前面学习过序列、卷积、傅里叶变换、Z 变换和系统频率响应之后,我们已经知道:一个线性时不变系统可以通过系统函数 H(z)H(z)、单位脉冲响应 h(n)h(n) 或频率响应 H(ejω)H(e^{j\omega}) 来描述。滤波器设计要做的事情,就是根据希望得到的频率选择效果,反过来求出一个可以实现的 H(z)H(z)h(n)h(n)

这一讲的主线是 IIR 数字滤波器设计。IIR 是 Infinite Impulse Response 的缩写,意思是“无限长单位脉冲响应”。这类滤波器通常具有反馈结构,单位脉冲响应理论上可以无限延续。IIR 滤波器的一个重要设计思路是:先设计成熟的模拟滤波器,再通过某种变换把模拟滤波器转化为数字滤波器。讲义中重点介绍了两种变换方法:脉冲响应不变法和双线性变换法。

学习路线

本讲可以分为四个部分。第一部分认识数字滤波器的基本概念,理解滤波的目的、分类、差分方程、IIR 与 FIR 的区别以及滤波器技术指标。第二部分学习模拟滤波器设计方法,重点是巴特沃斯低通滤波器,同时了解切比雪夫滤波器、椭圆滤波器以及从低通到高通、带通、带阻的频率变换思想。第三部分学习脉冲响应不变法,理解如何让数字滤波器的单位脉冲响应等于模拟滤波器冲激响应的采样值。第四部分学习双线性变换法,理解为什么它可以避免频谱混叠,以及为什么需要频率预校正。


数字滤波器的基本概念

滤波的直接目的,是抑制输入信号中的某些频率成分,从而改变信号频谱中各频率分量的相对比例。例如,一段语音信号中混入了高频噪声,可以用低通滤波器保留低频语音主体,削弱高频干扰;又如工频干扰集中在某个频率附近,可以用陷波器把这一小段频率压下去。

从更广义的角度看,滤波还包括信号检测与信号参量估计。信号检测关心的是:在干扰背景中,目标信号是否存在。参量估计关心的是:为了识别信号,需要估计它的某些参数,比如频率、幅值、相位、时延等。

滤波技术通常包括两件事。第一是滤波器设计,也就是根据给定的频率特性要求,求出满足要求的传递函数。第二是滤波过程实现,也就是拿到传递函数后,如何用程序、硬件或 DSP 芯片把这个系统真正运行起来。

什么是数字滤波器

数字滤波器是输入和输出都是数字信号的系统。它通过数值运算改变输入信号中各频率成分的比例,或者滤除某些不需要的频率成分。与模拟滤波器相比,数字滤波器具有处理精度高、稳定、重量轻、设计灵活、不存在阻抗匹配问题等优点,还可以实现某些模拟滤波器难以实现的特殊功能。

数字滤波器一般看作线性时不变系统。设计数字滤波器时,通常已知希望的频率响应 H(ejω)H(e^{j\omega}),然后要求系统函数 H(z)H(z) 或单位脉冲响应 h(n)h(n)。理想频率特性在数学上可以写出来,但实际系统不可能完全实现理想突变的频率响应,所以工程上常给出允许误差范围。只要滤波器满足指标,就可以认为设计合格。因此,在同一组指标下,满足要求的 H(z)H(z)h(n)h(n) 往往不是唯一的。

经典滤波器与现代滤波器

经典滤波器通常用于有用信号与干扰信号占据不同频带的情况。设输入信号为

x(n)=s(n)+u(n),x(n)=s(n)+u(n),

其中 s(n)s(n) 是有用成分,u(n)u(n) 是希望去除的噪声或干扰。如果它们在频域中分布在不同频段,就可以通过线性系统把 u(n)u(n) 有效削弱。

常见选频滤波器包括低通、高通、带通和带阻四类。

低通滤波器保留低频、抑制高频;高通滤波器保留高频、抑制低频;带通滤波器只保留某一段中间频率;带阻滤波器则抑制某一段频率而保留其他频率。陷波器可以看成带阻滤波器的一种特殊形式,它只压制很窄的一小段频率。

理想滤波器为什么不能直接实现

理想低通、高通、带通、带阻滤波器的幅频响应具有“突然跳变”的特点。例如理想低通滤波器在通带内幅度为 1,在截止频率之后立即变为 0。这样的频率响应没有过渡带,图像上看非常干净,但实际系统不能严格实现。

原因在于实际滤波器的频率响应必须由有限结构或稳定因果系统实现。过于陡峭的理想突变往往对应无限长冲激响应,或者带来不可实现的非因果特性。因此工程设计中一般允许通带有小波纹,阻带有小残留,并设置过渡带。换句话说,滤波器不是“完美切掉”,而是在允许误差内“足够好地逼近”。

由于离散时间频率响应具有 2π2\pi 周期性,只需考虑 02π0\sim 2\pi 范围。对于实序列,频谱又具有共轭对称性,实际常只讨论 0π0\sim \pi 的幅频特性。

选频滤波器与自适应滤波器

本课程讨论的数字滤波器主要是选频滤波器。它适合有用信号和干扰信号频谱位于不同频段的情形。用 IIR 或 FIR 系统实现时,滤波器系数是固定的。

如果有用信号与干扰信号共享同一个频段,普通选频滤波器就很难把两者分开。这时往往需要自适应滤波器。自适应滤波器与普通滤波器不同,它的系数会随着输入信号变化而调整,常用于回声抵消、噪声抵消、系统辨识等场景。

滤波前后信号的直观理解

设输入信号由低频正弦波和高频正弦波叠加而成。滤波前,在时域中可以看到信号波形较复杂,频域中可以看到低频和高频两个主要谱峰。经过低通滤波后,高频部分被削弱,时域波形更接近低频正弦波;频域中高频谱峰明显降低。这就是低通滤波器“保留低频、抑制高频”的直观效果。

滤波器的差分方程

数字滤波器的输入输出关系常用差分方程表示:

y(n)=i=1Naiy(ni)+i=0Mbix(ni).y(n)=\sum_{i=1}^{N}a_i y(n-i)+\sum_{i=0}^{M}b_i x(n-i).

这个式子说明,当前输出 y(n)y(n) 由两部分组成:过去若干个输出 y(ni)y(n-i) 的线性组合,以及当前和过去若干个输入 x(ni)x(n-i) 的线性组合。

如果输出还与过去输出有关,系统就具有反馈,也说明系统带有记忆。数字滤波器的基本操作其实只有三类:加法、乘法和延迟。复杂的数字滤波器最终都可以拆成这三种基本运算。

系统函数:IIR 与 FIR 的分界

数字滤波器的系统函数一般写成有理函数形式:

H(z)=Y(z)X(z)=i=0Mbizi1i=1Naizi.H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{i=0}^{M} b_i z^{-i}}{1-\sum_{i=1}^{N}a_i z^{-i}}.

如果反馈系数 aia_i 中至少有一个不为 0,那么分母中含有反馈项,这类系统通常是 IIR 滤波器。IIR 的单位脉冲响应一般是无限长的。

如果所有 ai=0a_i=0,则系统函数变为

H(z)=i=0Mbizi,H(z)=\sum_{i=0}^{M}b_i z^{-i},

这时系统没有递归反馈,只由有限个输入样本加权求和构成,称为 FIR 滤波器。FIR 是 Finite Impulse Response,即有限长单位脉冲响应。

IIR 系统的基本特点

IIR 滤波器有三个基本特点。

第一,单位脉冲响应 h(n)h(n) 是无限长的。即使输入只是一个单位冲激,输出也可能持续很久。

第二,系统函数 H(z)H(z) 在有限 zz 平面上通常存在极点。极点位置决定系统稳定性。如果极点落在单位圆外,系统可能不稳定。

第三,结构上是递归型的,即存在输出到输入计算过程的反馈。因此 IIR 滤波器通常用较低阶数就能实现较陡峭的幅频响应,但稳定性与相位特性需要特别注意。

数字滤波器的技术指标

数字滤波器的频率响应可以写成

H(ejω)=H(ejω)ejθ(ω).H(e^{j\omega})=|H(e^{j\omega})|e^{j\theta(\omega)}.

其中 H(ejω)|H(e^{j\omega})| 是幅频特性,表示不同频率成分通过滤波器后的幅度衰减情况;θ(ω)\theta(\omega) 是相频特性,反映不同频率成分通过滤波器后的时间延迟情况。

选频滤波器一般主要考虑幅频特性,对相频特性不作严格要求。如果输出波形形状很重要,例如音频、图像、通信波形等场景,就需要考虑线性相位问题。

以低通滤波器为例,技术指标通常包括通带截止频率 ωp\omega_p、阻带截止频率 ωs\omega_s、通带允许最大衰减 αp\alpha_p 和阻带应达到的最小衰减 αs\alpha_s

在通带中,希望幅度响应接近 1:

1δ1H(ejω)1,ω<ωp.1-\delta_1\le |H(e^{j\omega})|\le 1,\quad |\omega|<\omega_p.

在阻带中,希望幅度响应接近 0:

H(ejω)δ2,ωsωπ.|H(e^{j\omega})|\le \delta_2,\quad \omega_s\le |\omega|\le \pi.

通带和阻带之间是过渡带。过渡带越窄,滤波器越难设计,通常需要更高阶数。

用分贝表示滤波器指标

工程中常用分贝表示通带波纹和阻带衰减。通带最大衰减可写为

αp=20lgmaxH(ejω)minH(ejω),0ωωp.\alpha_p=20\lg\frac{\max |H(e^{j\omega})|}{\min |H(e^{j\omega})|},\quad 0\le |\omega|\le \omega_p.

ω=0\omega=0 处幅度归一化为 1,则可写成

αp=20lg(1δ1).\alpha_p=-20\lg(1-\delta_1).

阻带最小衰减为

αs=20lg通带中最大H(ejω)阻带中最大H(ejω)=20lgδ2.\alpha_s=20\lg\frac{\text{通带中最大}|H(e^{j\omega})|}{\text{阻带中最大}|H(e^{j\omega})|}=-20\lg\delta_2.

H(ejωc)=0.707=1/2|H(e^{j\omega_c})|=0.707=1/\sqrt{2} 时,对应功率下降一半,衰减约为 3 dB,因此 ωc\omega_c 称为 3 dB 通带截止频率。

需要注意:αp\alpha_p 越小,表示通带波纹越小;αs\alpha_s 越大,表示阻带衰减越强;ωp\omega_pωs\omega_s 越接近,过渡带越窄,设计难度越高。

IIR 数字滤波器设计方法概述

设计 IIR 数字滤波器一般有两类方法。

第一类是间接法。先设计满足要求的模拟滤波器,再把模拟滤波器转换成数字滤波器。这种方法适合设计幅频特性比较规则的滤波器,例如低通、高通、带通和带阻滤波器。讲义中重点讨论的脉冲响应不变法和双线性变换法都属于间接法。

第二类是直接法。直接在频域或时域中设计数字滤波器。由于往往需要联立方程或进行数值优化,实际设计常需要计算机辅助。


模拟滤波器设计方法

IIR 数字滤波器设计之所以常从模拟滤波器开始,是因为模拟滤波器理论已经非常成熟,尤其是低通原型滤波器有完整的设计公式和参数表。常用模拟滤波器包括巴特沃斯滤波器、切比雪夫滤波器和椭圆滤波器。

模拟低通滤波器的指标

模拟滤波器的单位冲激响应记作 ha(t)h_a(t),其频率响应为

Ha(jΩ)=+ha(t)ejΩtdt.H_a(j\Omega)=\int_{-\infty}^{+\infty}h_a(t)e^{-j\Omega t}\,dt.

系统函数为

Ha(s)=+ha(t)estdt.H_a(s)=\int_{-\infty}^{+\infty}h_a(t)e^{-st}\,dt.

设计模拟滤波器时,指标一般由幅频响应 Ha(jΩ)|H_a(j\Omega)| 给出,而设计任务是求系统函数 Ha(s)H_a(s)

工程中常用损耗函数,也称衰减函数,来描述模拟滤波器的幅频特性:

A(Ω)=20lgHa(jΩ)=10lgHa(jΩ)2dB.A(\Omega)=-20\lg |H_a(j\Omega)|=-10\lg |H_a(j\Omega)|^2\quad \text{dB}.

损耗函数把小幅度变化放大显示,因此可以同时观察通带和阻带的变化情况。例如当 Ha(jΩc)=1/2|H_a(j\Omega_c)|=1/\sqrt{2} 时,A(Ωc)=3 dBA(\Omega_c)=3\text{ dB}Ωc\Omega_c 就是 3 dB 截止频率。

模拟低通滤波器的四个关键指标是:通带截止频率 Ωp\Omega_p、阻带截止频率 Ωs\Omega_s、通带最大衰减 αp\alpha_p、阻带最小衰减 αs\alpha_s

由幅度平方函数确定系统函数

模拟滤波器幅度响应常用幅度平方函数表示:

Ha(jΩ)2=Ha(jΩ)Ha(jΩ).|H_a(j\Omega)|^2=H_a(j\Omega)H_a^*(j\Omega).

如果滤波器冲激响应 ha(t)h_a(t) 是实函数,则

Ha(jΩ)=Ha(jΩ).H_a^*(j\Omega)=H_a(-j\Omega).

因此

Ha(jΩ)2=Ha(jΩ)Ha(jΩ)=Ha(s)Ha(s)s=jΩ.|H_a(j\Omega)|^2=H_a(j\Omega)H_a(-j\Omega)=H_a(s)H_a(-s)\big|_{s=j\Omega}.

这说明,如果能由指标求出 Ha(jΩ)2|H_a(j\Omega)|^2,就可以得到 Ha(s)Ha(s)H_a(s)H_a(-s),再从中选出因果稳定的 Ha(s)H_a(s)。因果稳定的模拟滤波器要求极点在 ss 平面左半平面,因此在构造 Ha(s)H_a(s) 时只取左半平面的极点。

巴特沃斯低通滤波器

巴特沃斯逼近又称最平幅度逼近。它的特点是通带内幅度响应最大平坦,并且随着 Ω\Omega 增大,幅频响应单调下降。巴特沃斯低通滤波器的幅度平方函数为

Ha(jΩ)2=11+(ΩΩc)2N.|H_a(j\Omega)|^2=\frac{1}{1+\left(\frac{\Omega}{\Omega_c}\right)^{2N}}.

其中 NN 是滤波器阶数,Ωc\Omega_c 是 3 dB 截止频率。

从幅频曲线可以看出,当 Ω<Ωc\Omega<\Omega_c 时,阶数 NN 越大,通带越平坦;当 Ω>Ωc\Omega>\Omega_c 时,阶数 NN 越大,衰减越快,过渡带越陡。巴特沃斯滤波器的优点是幅度特性平坦,相位特性近似线性;缺点是为达到较陡过渡带,阶数往往较高。

Ω=0\Omega=0 时,Ha(jΩ)=1|H_a(j\Omega)|=1,表示直流分量无衰减。当 Ω=Ωc\Omega=\Omega_c 时,

Ha(jΩc)2=12,Ha(jΩc)=120.707,|H_a(j\Omega_c)|^2=\frac12,\quad |H_a(j\Omega_c)|=\frac{1}{\sqrt2}\approx 0.707,

即衰减 3 dB。

巴特沃斯圆与左半平面极点

将幅度平方函数写成 ss 的函数:

Ha(s)Ha(s)=11+(sjΩc)2N.H_a(s)H_a(-s)=\frac{1}{1+\left(\frac{s}{j\Omega_c}\right)^{2N}}.

这个函数有 2N2N 个极点,这些极点等间隔分布在半径为 Ωc\Omega_c 的圆上,这个圆称为巴特沃斯圆。极点可表示为

sk=Ωcejπ(12+2k+12N).s_k=\Omega_c e^{j\pi\left(\frac12+\frac{2k+1}{2N}\right)}.

为了形成因果稳定滤波器,2N2N 个极点中只取左半平面的 NN 个极点构成 Ha(s)H_a(s),右半平面的极点属于 Ha(s)H_a(-s)。因此

Ha(s)=ΩcNk=0N1(ssk).H_a(s)=\frac{\Omega_c^N}{\prod_{k=0}^{N-1}(s-s_k)}.

p=s/Ωcp=s/\Omega_c,则归一化低通原型为

Ga(p)=1k=0N1(ppk).G_a(p)=\frac{1}{\prod_{k=0}^{N-1}(p-p_k)}.

展开后可写成

Ga(p)=1pN+bN1pN1++b1p+b0.G_a(p)=\frac{1}{p^N+b_{N-1}p^{N-1}+\cdots+b_1p+b_0}.

实际设计时,归一化极点 pkp_k 和系数 bkb_k 可以查巴特沃斯归一化低通滤波器参数表得到。

巴特沃斯滤波器的阶数与截止频率

巴特沃斯设计的关键,是由技术指标求阶数 NN 和截止频率 Ωc\Omega_c。由

Ha(jΩ)2=11+(Ω/Ωc)2N|H_a(j\Omega)|^2=\frac{1}{1+(\Omega/\Omega_c)^{2N}}

可得阶数公式

N=log(100.1αs1100.1αp1)2log(Ωs/Ωp).N=\frac{\log\left(\frac{10^{0.1\alpha_s}-1}{10^{0.1\alpha_p}-1}\right)}{2\log(\Omega_s/\Omega_p)}.

实际阶数取大于或等于该数的最小整数。

截止频率可以由通带或阻带指标计算,例如

Ωc=Ωp(100.1αp1)1/(2N)\Omega_c=\frac{\Omega_p}{(10^{0.1\alpha_p}-1)^{1/(2N)}}

Ωc=Ωs(100.1αs1)1/(2N).\Omega_c=\frac{\Omega_s}{(10^{0.1\alpha_s}-1)^{1/(2N)}}.

巴特沃斯低通滤波器设计步骤

巴特沃斯低通设计可以概括为三步。

第一,根据 Ωp,αp,Ωs,αs\Omega_p,\alpha_p,\Omega_s,\alpha_s 求阶数 NN 和 3 dB 截止频率 Ωc\Omega_c。第二,求归一化极点 pkp_k,并得到归一化低通原型系统函数 Ga(p)G_a(p)。第三,将 p=s/Ωcp=s/\Omega_c 代入 Ga(p)G_a(p),完成去归一化,得到实际模拟滤波器系统函数

Ha(s)=Ga(sΩc).H_a(s)=G_a\left(\frac{s}{\Omega_c}\right).
解答题例 1:巴特沃斯低通设计的阶数与极点

讲义例题给出:通带截止频率 fp=5 kHzf_p=5\text{ kHz},通带最大衰减 αp=2 dB\alpha_p=2\text{ dB},阻带截止频率 fs=12 kHzf_s=12\text{ kHz},阻带最小衰减 αs=30 dB\alpha_s=30\text{ dB}。按公式算得 N=4.25N=4.25,取 N=5N=5。再求 Ωc\Omega_c,查表得到五阶巴特沃斯归一化极点为

0.3090±j0.9511,0.8090±j0.5878,1.0000.-0.3090\pm j0.9511,\quad -0.8090\pm j0.5878,\quad -1.0000.

归一化低通原型系统函数可写成五阶多项式分母形式。最后将 p=s/Ωcp=s/\Omega_c 代入即可得到实际滤波器。

解答

这道题的关键是两步:先由通阻带指标求最小阶数,再由阶数反推截止频率。

由巴特沃斯阶数公式算得

N=4.25,N=4.25,

由于滤波器阶数必须取整数,且必须满足指标,所以向上取整得到

N=5.N=5.

接着根据通带或阻带条件都可以求 Ωc\Omega_c。确定 Ωc\Omega_c 之后,再查五阶巴特沃斯归一化极点表,构造归一化原型系统函数,最后代入

p=sΩcp=\frac{s}{\Omega_c}

完成去归一化。

这类题目的主干流程固定为:指标 \to 阶数 NN \to 截止频率 Ωc\Omega_c \to 归一化原型 \to 实际滤波器。

巴特沃斯设计总结

巴特沃斯归一化表和公式都是相对于 3 dB 截止频率 Ωc\Omega_c 给出的。设计时最关键的两个参数是阶数 NN 和截止频率 Ωc\Omega_c。其中 NN 用来确定巴特沃斯多项式,Ωc\Omega_c 用来完成去归一化。

切比雪夫滤波器

巴特沃斯滤波器的缺点是阶数通常较高。原因是它的幅频特性在通带和阻带中都是单调变化的。如果某一边刚好满足指标,另一边可能存在较大富余量。

切比雪夫滤波器采用切比雪夫函数逼近给定指标,具有等波纹特性。它把误差更均匀地分布在通带或阻带内,因此在相同指标下通常可以降低阶数。

切比雪夫滤波器分为两类。切比雪夫 I 型为通带等波纹、阻带单调;切比雪夫 II 型为通带单调、阻带等波纹。切比雪夫 I 型幅度平方函数为

Ha(jΩ)2=11+ε2CN2(Ω/Ωp).|H_a(j\Omega)|^2=\frac{1}{1+\varepsilon^2 C_N^2(\Omega/\Omega_p)}.

其中 ε\varepsilon 是波纹参数,表示通带内幅度波动程度。ε\varepsilon 越大,通带波动越明显。CN(x)C_N(x)NN 阶切比雪夫多项式。

设计切比雪夫 I 型滤波器通常包括:根据指标确定 ε\varepsilon 和阶数 NN,求归一化传递函数 Ga(p)G_a(p),最后将 p=s/Ωpp=s/\Omega_p 代入完成去归一化。

椭圆滤波器

椭圆滤波器允许通带和阻带内都出现等波纹。因此在满足相同幅频指标时,椭圆滤波器通常可以进一步降低阶数。它的缺点是相位响应非线性更强,对波形的失真较大。

三类常用滤波器可以这样比较:巴特沃斯滤波器通带和阻带都单调,幅度最平坦,相位较好,但阶数较高;切比雪夫滤波器通带或阻带存在等波纹,阶数介于两者之间;椭圆滤波器通带和阻带都波动,阶数最低,但相位线性最差。

在只关心幅频指标的场合,椭圆滤波器往往性能价格比较高,应用也较广泛。

模拟高通、带通、带阻滤波器设计

模拟滤波器设计手册通常先给出低通滤波器公式,再通过频率变换从低通得到高通、带通和带阻滤波器。

一般过程为:先把希望设计的高通、带通或带阻指标转换为等效低通指标;再设计低通原型系统函数 Q(p)Q(p);最后对 Q(p)Q(p) 做频率变换,得到目标类型的模拟滤波器系统函数 Ha(s)H_a(s)


用模拟滤波器设计 IIR 数字滤波器

利用模拟滤波器设计数字滤波器,本质上是从已知模拟传递函数 Ha(s)H_a(s) 出发,求数字系统函数 H(z)H(z)。这是一个从 ss 平面到 zz 平面的映射问题。

这个映射必须遵循两个原则。第一,数字滤波器的频率响应要能模仿模拟滤波器的频率响应,因此 ss 平面的虚轴应映射到 zz 平面的单位圆。第二,模拟滤波器的因果稳定性映射后要保持不变,因此 ss 平面的左半平面应映射到 zz 平面的单位圆内。

一般设计流程是:先将数字滤波器指标转换为模拟滤波器指标;若目标不是低通,则先转换为低通原型指标;设计过渡模拟低通滤波器;再将模拟低通转换成目标类型的模拟滤波器;最后用脉冲响应不变法或双线性变换法得到数字滤波器。


脉冲响应不变法

脉冲响应不变法的核心思想是:让数字滤波器的单位脉冲响应 h(n)h(n) 等于模拟滤波器冲激响应 ha(t)h_a(t) 的采样值。

h(n)=ha(nT),h(n)=h_a(nT),

其中 TT 为采样周期。已知

Ha(s)=L[ha(t)],H(z)=Z[h(n)].H_a(s)=\mathcal{L}[h_a(t)],\quad H(z)=\mathcal{Z}[h(n)].

通过这个关系,就可以从 Ha(s)H_a(s) 推导出 H(z)H(z)

脉冲响应不变法的映射关系

脉冲响应不变法对应的平面映射关系为

z=esT.z=e^{sT}.

这说明 ss 平面中不同的虚轴频率会以指数形式映射到 zz 平面单位圆上。由于指数函数具有周期性,模拟频谱会周期延拓到数字频域中。

数字滤波器频率响应与模拟滤波器频率响应之间有如下关系:

H(ejω)=1Tk=+Ha[j(ωT2πkT)].H(e^{j\omega})=\frac{1}{T}\sum_{k=-\infty}^{+\infty}H_a\left[j\left(\frac{\omega}{T}-\frac{2\pi k}{T}\right)\right].

这个式子说明,数字滤波器频率响应是模拟滤波器频率响应以采样频率为周期进行周期延拓后叠加得到的。如果模拟滤波器频率响应不是严格带限,或者采样频率不满足采样定理,就会产生频谱混叠。

因此,脉冲响应不变法最适合低通和带通等在折叠频率以上衰减充分大的情况,不适合高通和带阻滤波器。

部分分式形式下的数字化

脉冲响应不变法特别适合模拟滤波器传递函数能写成部分分式形式的情况。设模拟滤波器只有单阶极点,且分母阶数高于分子阶数:

Ha(s)=i=1NAissi.H_a(s)=\sum_{i=1}^{N}\frac{A_i}{s-s_i}.

其中 sis_i 是单阶极点。其冲激响应为

ha(t)=i=1NAiesitu(t).h_a(t)=\sum_{i=1}^{N}A_i e^{s_i t}u(t).

采样后得到

h(n)=ha(nT)=i=1NAiesinTu(nT).h(n)=h_a(nT)=\sum_{i=1}^{N}A_i e^{s_i nT}u(nT).

h(n)h(n) 作 Z 变换:

H(z)=i=1NAi1esiTz1.H(z)=\sum_{i=1}^{N}\frac{A_i}{1-e^{s_iT}z^{-1}}.

若作幅度修正,则常写为

H(z)=i=1NTAi1esiTz1.H(z)=\sum_{i=1}^{N}\frac{T A_i}{1-e^{s_iT}z^{-1}}.

比较前后两式可知,系数 AiA_i 保持对应,而 ss 平面中的单极点 sis_i 映射到 zz 平面中的极点 esiTe^{s_iT}

解答题例 2:脉冲响应不变法的基本计算

设模拟滤波器系统函数为

Ha(s)=2(s+1)(s+3)=1s+11s+3.H_a(s)=\frac{2}{(s+1)(s+3)}=\frac{1}{s+1}-\frac{1}{s+3}.

极点为 s=1s=-1s=3s=-3,部分分式系数分别为 111-1。由脉冲响应不变法,有

H(z)=11eTz111e3Tz1.H(z)=\frac{1}{1-e^{-T}z^{-1}}-\frac{1}{1-e^{-3T}z^{-1}}.

若采用修正形式,还需要在系数前乘以 TT。讲义给出的具体化简结果中,当给定采样周期后,可以进一步写成 z1z^{-1} 的有理函数形式。

解答

先把模拟系统函数拆成部分分式:

Ha(s)=1s+11s+3.H_a(s)=\frac{1}{s+1}-\frac{1}{s+3}.

这样可以直接读出两个单极点

s1=1,s2=3.s_1=-1,\qquad s_2=-3.

脉冲响应不变法下,每个模拟极点 sis_i 都映射到数字极点 esiTe^{s_iT},因此得到

H(z)=11eTz111e3Tz1.H(z)=\frac{1}{1-e^{-T}z^{-1}}-\frac{1}{1-e^{-3T}z^{-1}}.

这道题最核心的记忆点是:脉冲响应不变法最适合先做部分分式展开,再逐项把

Aissi\frac{A_i}{s-s_i}

替换成

Ai1esiTz1.\frac{A_i}{1-e^{s_iT}z^{-1}}.

脉冲响应不变法的优缺点

脉冲响应不变法的优点是时域逼近好,因为数字滤波器的单位脉冲响应直接来自模拟滤波器冲激响应的采样值。同时,模拟频率与数字频率满足线性关系

ω=ΩT.\omega=\Omega T.

因此线性相位关系可以保持得较好。

它的缺点是存在频谱混叠。由于实际模拟滤波器不可能严格带限,频谱混叠通常不可完全避免。所以该方法适合低通和带通滤波器,而且要求模拟滤波器在折叠频率以上有足够衰减;不适合高通和带阻滤波器。


双线性变换法

双线性变换法的目标是建立 ss 平面和 zz 平面之间的一一对应关系,从而消除脉冲响应不变法中的多值映射和频谱混叠。

它可以理解为两步。第一步,把整个 ss 平面压缩到另一个 s1s_1 平面的一条横带中,横带宽度为 2π/T2\pi/T。第二步,再通过指数映射把这条横带映射到整个 zz 平面。

最终得到的双线性变换公式为

s=2T1z11+z1.s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}.

也可以写成

z=2T+s2Ts.z=\frac{\frac{2}{T}+s}{\frac{2}{T}-s}.

因此数字滤波器系统函数可由

H(z)=Ha(s)s=2T1z11+z1H(z)=H_a(s)\bigg|_{s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}

得到。

双线性变换的频率关系

双线性变换的频率映射为

Ω=2Ttanω2.\Omega=\frac{2}{T}\tan\frac{\omega}{2}.

它的优点是没有频谱混叠,因为 ss 平面和 zz 平面之间是一一映射。缺点是模拟频率 Ω\Omega 与数字频率 ω\omega 之间不是线性关系。零频率附近近似线性,但频率越高,非线性越明显,这会造成频率轴畸变。

因此,双线性变换会改变频率刻度。如果原模拟滤波器相位是线性的,变换后得到的数字滤波器也可能变成非线性相位。

频率预校正

由于双线性变换存在非线性频率映射,设计前需要把数字边界频率转换成对应的模拟边界频率,这叫频率预校正,公式为

Ω=2Ttanω2.\Omega=\frac{2}{T}\tan\frac{\omega}{2}.

例如设计数字低通滤波器时,给定 ωp\omega_pωs\omega_s,不能直接把它们当作模拟滤波器的 Ωp\Omega_pΩs\Omega_s。必须先用上式转换,再按模拟滤波器设计方法求 Ha(s)H_a(s),最后再用双线性变换得到 H(z)H(z)

数字高通、带通和带阻滤波器的设计

对于数字高通、带通和带阻滤波器,常用方法是双线性变换。基本思路是:先把数字滤波器边界频率通过预校正转换为相应模拟滤波器边界频率;再把目标类型的模拟滤波器指标转换为模拟低通原型指标;设计模拟低通原型;通过频率变换得到目标类型的模拟滤波器;最后用双线性变换得到所需数字滤波器。

这种方法适用于低通、高通、带通、带阻等选频滤波器,也是 IIR 数字滤波器设计中使用最普遍、最有效的工具之一。

解答题例 3:双线性变换的一阶系统替换

若模拟滤波器为一阶系统

Ha(s)=αs+α,H_a(s)=\frac{\alpha}{s+\alpha},

则用双线性变换

s=2T1z11+z1s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}

代入,得到

H(z)=Ha(s)s=2T1z11+z1.H(z)=H_a(s)\bigg|_{s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}.

这一步的本质是把模拟系统函数中的 ss 全部替换为 z1z^{-1} 的有理式,从而得到可由数字滤波器实现的 H(z)H(z)

解答

双线性变换法的做法很直接:不是先求冲激响应,而是把模拟系统函数里的 ss 统一替换成

s=2T1z11+z1.s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}.

因此一阶模拟低通

Ha(s)=αs+αH_a(s)=\frac{\alpha}{s+\alpha}

经过替换后,立即变成关于 z1z^{-1} 的有理式,这就是对应的数字滤波器系统函数 H(z)H(z)

这类题的重点不在展开代数细节,而在理解双线性变换保持了一一映射,所以不会出现脉冲响应不变法中的频谱混叠。

解答题例 4:双线性变换法设计数字低通滤波器

讲义还给出一个数字低通设计例题:给定 ωp=0.2π\omega_p=0.2\piαp=3 dB\alpha_p=3\text{ dB}ωs=0.3π\omega_s=0.3\piαs=15 dB\alpha_s=15\text{ dB}T=1 sT=1\text{ s},先用预畸变计算模拟低通指标,再设计巴特沃斯模拟低通滤波器,得到阶数 N=6N=6 和截止频率 Ωc=0.7663 rad/s\Omega_c=0.7663\text{ rad/s}。查表得归一化原型,再去归一化得到 Ha(s)H_a(s),最后用双线性变换得到数字滤波器

H(z)=0.0007378(1+z1)6(11.268z1+0.7051z2)(11.010z1+0.358z2)(10.9044z1+0.2155z2).H(z)=\frac{0.0007378(1+z^{-1})^6}{(1-1.268z^{-1}+0.7051z^{-2})(1-1.010z^{-1}+0.358z^{-2})(1-0.9044z^{-1}+0.2155z^{-2})}.

这个结果说明,高阶 IIR 滤波器常写成多个二阶节相乘的形式,这样更便于数值实现,也更稳定。

解答

这道题的完整思路分四步。

第一,先对数字边界频率做预校正,把 ωp,ωs\omega_p,\omega_s 转成模拟频率 Ωp,Ωs\Omega_p,\Omega_s

第二,用预校正后的模拟低通指标设计巴特沃斯模拟低通滤波器,得到阶数

N=6N=6

和截止频率

Ωc=0.7663 rad/s.\Omega_c=0.7663\text{ rad/s}.

第三,根据巴特沃斯归一化原型求出实际模拟系统函数 Ha(s)H_a(s)

第四,把 Ha(s)H_a(s) 中的 ss 用双线性变换公式替换,得到数字滤波器 H(z)H(z)

最后写成多个二阶节相乘的形式,不只是为了好看,而是因为高阶 IIR 直接写成一个高次分母时数值敏感性更强;拆成二阶节后实现更稳健。


课堂练习整理

将模拟滤波器 Ha(s)H_a(s) 转换为数字滤波器 H(z)H(z) 的常用方法有两种:脉冲响应不变法和双线性变换法。

双线性变换法的主要优点是避免频率响应的频谱混叠现象。

脉冲响应不变法的缺点是会产生频谱混叠,优点是数字频率与模拟角频率成线性关系,因此适合低通和带通滤波器设计。

数字滤波器的两个分支中,具有递归结构的是 IIR 滤波器;绝对稳定、无反馈结构的是 FIR 滤波器。

对于模拟滤波器

Ha(s)=12s2+3s+1,H_a(s)=\frac{1}{2s^2+3s+1},

可以分别用脉冲响应不变法和双线性变换法转换。先分解

Ha(s)=1s+1/21s+1.H_a(s)=\frac{1}{s+1/2}-\frac{1}{s+1}.

T=2T=2,脉冲响应不变法得到

H(z)=11e1z111e2z1.H(z)=\frac{1}{1-e^{-1}z^{-1}}-\frac{1}{1-e^{-2}z^{-1}}.

双线性变换法则使用

s=1z11+z1s=\frac{1-z^{-1}}{1+z^{-1}}

代入

H(z)=Ha(s)=12(1z11+z1)2+3(1z11+z1)+1.H(z)=H_a(s)=\frac{1}{2\left(\frac{1-z^{-1}}{1+z^{-1}}\right)^2+3\left(\frac{1-z^{-1}}{1+z^{-1}}\right)+1}.

两种方法的比较是:脉冲响应不变法时域逼近好,模拟频率与数字频率关系线性,但存在频谱混叠;双线性变换法避免频谱混叠,但频率映射非线性,需要频率预校正。

利用模拟滤波器设计 IIR 数字低通滤波器的步骤

完整步骤可以整理为四步。

第一,确定数字低通滤波器技术指标,包括通带边界频率 ωp\omega_p、通带最大衰减 αp\alpha_p、阻带截止频率 ωs\omega_s、阻带最小衰减 αs\alpha_s

第二,将数字低通滤波器技术指标转换成模拟低通滤波器技术指标。若使用双线性变换,需要用

Ω=2Ttanω2\Omega=\frac{2}{T}\tan\frac{\omega}{2}

进行频率预校正。

第三,按照模拟低通滤波器技术指标设计过渡模拟低通滤波器,例如巴特沃斯、切比雪夫或椭圆滤波器。

第四,用所选变换方法把模拟滤波器 Ha(s)H_a(s) 转换成数字低通滤波器系统函数 H(z)H(z)

本讲小结

这一讲的核心是:数字滤波器设计不是凭空写一个差分方程,而是先根据频率指标确定希望的幅频特性,再求能实现该特性的系统函数。

对于 IIR 滤波器,常用的工程路线是先设计模拟滤波器,再转换为数字滤波器。模拟滤波器设计中,巴特沃斯滤波器具有最大平坦特性,切比雪夫滤波器用等波纹降低阶数,椭圆滤波器在通带和阻带都允许波动,从而进一步降低阶数。

脉冲响应不变法保留时域冲激响应采样,频率关系线性,但会产生混叠;双线性变换法避免混叠,是最常用的 IIR 数字滤波器设计方法,但存在频率轴非线性畸变,因此要做频率预校正。

学完这一讲后,应该能回答三个问题:第一,数字滤波器的技术指标如何描述;第二,巴特沃斯模拟低通滤波器如何由指标求阶数和截止频率;第三,模拟滤波器如何通过脉冲响应不变法或双线性变换法转化为 IIR 数字滤波器。

Back to archive

Discussion

Comments

Post

Share questions, corrections, or extra notes about this post.