本文是“离散时间信号处理”系列的第 07 章，主题为“滤波器设计方法”。

上一篇：离散时间系统结构 · 下一篇：离散傅里叶变换

滤波器设计方法

IIR滤波器的设计

脉冲响应不变法

离散时间滤波器与连续时间滤波器的时域关系
\[ h[n]=T_dh_c(nT_d)\quad\text{(}T_d\text{表示采样间隔)} \]

(频域设计) 离散时间滤波器和连续时间滤波器的频域关系
\[ H(e^{j\omega})=\sum_{k=-\infty}^{\infty}H_c\left(j\frac{\omega}{T_d}+j\frac{2\pi}{T_d}k\right)\quad\text{(注意无}\frac{1}{T_d}\text{加权)} \]

如果连续时间滤波器是带限的, 即
\[ H_c(j\Omega)=0,|\Omega|\geq\frac{\pi}{T_d} \]

则离散时间滤波器在一个周期内可表示为
\[ H(e^{j\omega})=H_c(j\Omega)\big|_{\Omega=\omega/T_d}=H_c\left(j\frac{\omega}{T_d}\right),|\omega|\leq\pi \]

(时域设计) 连续时间滤波器的系统函数
\[ H_c(s)=\sum_{k=1}^{N}\frac{A_k}{s-s_k} \]

对应的脉冲响应
\[ h_c(t)=\begin{cases} \displaystyle\sum\limits_{k=1}^{N}A_ke^{s_kt} & t\geq0\\[2pt] 0 & t<0 \end{cases} \]

采样后得到离散时间滤波器的脉冲响应
\[ h[n]=T_dh_c(nT_d)=\sum_{k=1}^{N}T_dA_ke^{s_knT_d}u[n] \]

经过\(z\)变换后得到离散时间滤波器的系统函数
\[ H(z)=\sum_{k=1}^{N}\frac{T_dA_k}{1-e^{s_kT_d}z^{-1}} \]

双线性变换法

复变量\(s\)和\(z\)的关系
\[ s=\frac{2}{T_d}\left(\frac{1-z^{-1}}{1+z^{-1}}\right)\quad\text{(}\frac{2}{T_d}\text{的作用是预畸变)} \]

由此可得系统函数的关系为
\[ H(z)=H_c\left[\frac{2}{T_d}\left(\frac{1-z^{-1}}{1+z^{-1}}\right)\right] \]

离散时间频率\(\omega\)和连续时间频率\(\Omega\)的关系
\[ \Omega=\frac{2}{T_d}\tan\frac{\omega}{2} \]
\[ \begin{aligned} \omega&=2\arctan\frac{\Omega T_d}{2}\quad\text{(非线性压缩)}\\ &=\Omega_1T\quad\text{(}\Omega_1\text{为压缩后的连续时间频率)} \end{aligned} \]

复变量关系和频率关系的推导

(1) 由频率关系导出复变量关系, 已知\(\Omega=\dfrac{2}{T_d}\tan\dfrac{\omega}{2}\)
\[ \begin{aligned} \Omega&=\frac{2}{T_d}\left(\frac{\sin\dfrac{\omega}{2}}{\cos\dfrac{\omega}{2}}\right)\\ &=\frac{2}{jT_d}\left(\frac{e^{j\omega/2}-e^{-j\omega/2}}{e^{j\omega/2}+e^{-j\omega/2}}\right)\\ &=\frac{2}{jT_d}\left(\frac{1-e^{-j\omega}}{1+e^{-j\omega}}\right)\Rightarrow s=\frac{2}{T_d}\left(\frac{1-z^{-1}}{1+z^{-1}}\right) \end{aligned} \]

(2) 由复变量关系导出频率关系, 已知\(s=\dfrac{2}{T_d}\left(\dfrac{1-z^{-1}}{1+z^{-1}}\right)\)
\[ \begin{aligned} j\Omega&=\frac{2}{T_d}\left(\frac{1-e^{-j\omega}}{1+e^{-j\omega}}\right)\\ &=\frac{2}{T_d}\left(\frac{e^{-j\omega/2}}{e^{-j\omega/2}}\cdot\frac{e^{j\omega/2}-e^{-j\omega/2}}{e^{j\omega/2}+e^{-j\omega/2}}\right)\\ &=\frac{2}{T_d}\left(\frac{2j\sin\dfrac{\omega}{2}}{2\cos\dfrac{\omega}{2}}\right)\Rightarrow\Omega=\frac{2}{T_d}\tan\frac{\omega}{2} \end{aligned} \]

巴特沃兹低通滤波器的设计

离散时间滤波器的指标

离散时间滤波器的指标(容限图)

(1) 通带截止频率\(\omega_p\)、阻带截止频率\(\omega_s\)

截止频率\(\omega_c=\dfrac{\omega_p+\omega_s}{2}\)

(2) 通带逼近误差\(\delta_1\)、阻带逼近误差\(\delta_2\)

最大通带增益(分贝)\(=20\log_{10}(1+\delta_1)\)、最大阻带增益(分贝)\(=20\log_{10}(\delta_2)\)

最大通带衰减(分贝)\(=-20\log_{10}(1-\delta_1)\)、最小阻带衰减(分贝)\(=-20\log_{10}(\delta_2)\)

\(N\)阶连续时间巴特沃兹低通滤波器的幅度平方函数
\[ |H_c(j\Omega)|^2=\frac{1}{1+\left(\dfrac{j\Omega}{j\Omega_c}\right)^{2N}}=\frac{1}{1+\left(\dfrac{\Omega}{\Omega_c}\right)^{2N}} \]

连续时间巴特沃兹滤波器的幅度平方曲线

巴特沃兹滤波器的幅度特性曲线与阶次\(N\)的关系

(1) 随着\(N\)的增加, 滤波器的特性曲线变得更加尖锐, 更加趋近于理想滤波器.

(2) 在截止频率\(\Omega_c\)处, 幅度平方函数始终为\(\dfrac{1}{2}\), 与滤波器的阶次无关.

巴特沃兹低通滤波器的零极点特性

将\(j\Omega=s\)代入幅度平方函数可得
\[ H_c(s)H_c(-s)=\frac{1}{1+\left(\dfrac{s}{j\Omega_c}\right)^{2N}} \]

幅度平方函数的极点\(s_k\)满足\(1+\left(\dfrac{s_k}{j\Omega_c}\right)^{2N}=0\), 则
\[ \begin{aligned} s_k&=j\Omega_c\cdot(-1)^{1/2N}\\ &=\Omega_ce^{j\pi/2}\cdot(e^{-j\pi+j2k\pi})^{1/2N}\\ &=\Omega_ce^{(j\pi/2N)(N-1+2k)},k=0,1,\cdots,2N-1 \end{aligned} \]

(1) 滤波器的幅度平方函数在\(s\)平面中半径为\(\Omega_c\)的圆周上有\(2N\)个等分排列的极点, 相邻极点间的角度为\(\dfrac{\pi}{N}\)弧度, 在无穷处有\(2N\)阶零点(滤波器的系统函数为\(N\)阶).

(2) 幅度平方函数的极点关于虚轴(原点)对称分布, 但没有极点落在虚轴上, 并且当\(N\)为奇数时有一个极点位于实轴上, 而当\(N\)为偶数时则没有.

(3) 幅度平方函数的极点总是成对出现, 因此为了得到一个稳定和因果的滤波器, 应当选取位于\(s\)平面左半平面上的极点.

一个三阶巴特沃兹滤波器在\(s\)平面上极点的位置

窗函数法设计FIR滤波器

使用矩形窗得到因果FIR滤波器的脉冲响应
\[ \begin{aligned} h[n]&=\begin{cases} h_d[n] & 0\leq n\leq M\\ 0 & \text{其他} \end{cases}\\ &=h_d[n]\cdot w[n] \end{aligned} \]

其中\(h_d[n]\)为理想滤波器的脉冲响应, 则FIR滤波器的频率响应
\[ H(e^{j\omega})=\frac{1}{2\pi}\int_{-\pi}^{\pi}H_d(e^{j\theta})W(e^{j(\omega-\theta)})\mathrm{d}\theta \]

矩形窗的脉冲响应
\[ w[n]= \begin{cases} 1 & 0\leq n\leq M\\ 0 & \text{其他} \end{cases} \]

矩形窗的频率响应
\[ W(e^{j\omega})=\sum_{n=0}^{M}e^{-j\omega n}=\frac{1-e^{-j\omega(M+1)}}{1-e^{-j\omega}}=e^{-j\omega M/2}\cdot\frac{\sin\dfrac{\omega(M+1)}{2}}{\sin\dfrac{\omega}{2}} \]

吉布斯现象的形成

(1) 矩形窗\(W(e^{j\omega})\)的主瓣宽度为\(\Delta\omega_m=\dfrac{4\pi}{M+1}\), 大于\(H(e^{j\omega})\)的过渡带宽度.

(2) 矩形窗的主瓣幅度为\((M+1)\), 每个瓣的幅度随\(M\)的增大而增大, 宽度随\(M\)的增大而减小, 但是每个瓣的面积却是常量.

(3) 随着\(M\)的增大, \(H(e^{j\omega})\)出现的振荡将加快, 但是其幅度不随\(M\)的增大而减小.

(a) 截断理想脉冲响应所包含的卷积过程; (b) 对理想脉冲响应加窗后得到的典型逼近

矩形窗傅里叶变换的幅度(\(M=7\))

常用窗函数的性质

常用的窗函数

\(M=50\)时各种窗函数的傅里叶变换(对数幅度)；(a) 矩形窗; (b) Bartlett窗; (c) Hanning窗; (d) Hamming窗; (e) Blackman窗

窗的类型	最大旁瓣幅度 （相对值）	主瓣近似宽度	最大逼近误差 \(20\log_{10}\delta\,(\mathrm{dB})\)	等效 Kaiser 窗 \(\beta\)	等效 Kaiser 窗 过渡带宽度
矩形	\(-13\)	\(4\pi/(M+1)\)	\(-21\)	\(0\)	\(1.81\pi/M\)
Bartlett	\(-25\)	\(8\pi/M\)	\(-25\)	\(1.33\)	\(2.37\pi/M\)
Hanning	\(-31\)	\(8\pi/M\)	\(-44\)	\(3.86\)	\(5.01\pi/M\)
Hamming	\(-41\)	\(8\pi/M\)	\(-53\)	\(4.86\)	\(6.27\pi/M\)
Blackman	\(-57\)	\(12\pi/M\)	\(-74\)	\(7.04\)	\(9.19\pi/M\)

广义线性相位的合并

全部窗函数均关于点\(\dfrac{M}{2}\)偶对称, 即
\[ w[n]= \begin{cases} w[M-n] & 0\leq n\leq M\\ 0 & \text{其他} \end{cases} \]

对应的傅里叶变换
\[ W(e^{j\omega})=W_e(e^{j\omega})e^{-j\omega M/2} \]

其中\(W_e(e^{j\omega})\)是\(\omega\)的实偶函数.

(1) 若理想滤波器的脉冲响应关于点\(\dfrac{M}{2}\)偶对称, 即\(h_d[M-n]=h_d[n]\), 则加窗后的脉冲响应也是关于点\(\dfrac{M}{2}\)偶对称的, 且所得出的频率响应将有广义线性相位.

(2) 若理想滤波器的脉冲响应关于点\(\dfrac{M}{2}\)奇对称, 即\(h_d[M-n]=-h_d[n]\), 则加窗后的脉冲响应也是关于点\(\dfrac{M}{2}\)奇对称的, 且所得出的频率响应将有附带\(90^\circ\)常数相移的广义线性相位.

Kaiser窗的定义
\[ w[n]= \begin{cases} \displaystyle\frac{I_0[\beta(1-[(n-\alpha)/\alpha]^2)^{1/2}]}{I_0(\beta)} & 0\leq n\leq M\\[2pt] 0 & \text{其他} \end{cases} \]

其中
\[ \alpha=\frac{M}{2} \]

对于低通滤波器逼近, 其过渡区的宽度为
\[ \Delta\omega=\omega_s-\omega_p \]

定义
\[ A=-20\log_{10}\delta \]

Kaiser窗的长度参数
\[ M=\frac{A-8}{2.285\Delta\omega} \]

Kaiser窗的形状参数
\[ \beta= \begin{cases} 0.1102(A-8.7) & A>50\\ 0.5842(A-21)^{0.4}+0.07886(A-21) & 21\leq A\leq50\\ 0 & A<21 \end{cases} \]

(a) \(\beta=0,3,6\)以及\(M=20\)时的Kaiser窗; (b) 各窗函数的傅里叶变换; ；(c) 取\(\beta=6\)以及\(M=10,20,40\)时Kaiser窗的傅里叶变换

(1) 窗的两端越尖, 其傅里叶变换的旁瓣就越低, 但是主瓣也就越宽.

(2) 若增大\(M\)而同时保持\(\beta\)不变, 可使主瓣宽度减小, 且不影响旁瓣的幅度.

设计广义的多频带滤波器

计算窗参数的 Kaiser 公式可以用于预估多频带滤波器的逼近误差和过渡带宽度。应当注意，逼近误差的大小将与产生这些误差的跳变幅度成正比：如果幅度为 1 的间断点能产生 \(\delta\) 的峰值逼近误差，则幅度为 2 的间断点将产生 \(2\delta\) 的峰值逼近误差。

补充知识点

IIR滤波器

优点: 采用递归结构, 能以较低阶数实现.

缺点: 不能设计成严格线性相位.

FIR滤波器

优点: 很容易设计成严格线性相位.

缺点: 要达到同样性能指标情况, 比IIR阶数高, 且设计过程需反复多次.

首先根据要求选择滤波器类型: 如果要求线性相位则必须选择FIR滤波器, 否则最好选择IIR.

IIR滤波器的设计(模拟滤波器法)

(1) 等效的模拟系统指标: \(\omega_p=\Omega_pT_d\), \(\omega_s=\Omega_sT_d\)

(2) 数字指标: \(\Omega_p=\dfrac{\omega_p}{T_d}\), \(\Omega_s=\dfrac{\omega_s}{T_s}\); \(\Omega_p=\dfrac{2}{T_d}\tan\dfrac{\omega_p}{2}\), \(\Omega_s=\dfrac{2}{T_s}\tan\dfrac{\omega_s}{2}\)

(3) 原型模拟滤波器的指标

i. 脉冲响应不变法
\[ H_c(s)=\sum_{k=1}^{N}\frac{A_k}{s-s_k}\Rightarrow H(z)=\sum_{k=1}^{N}\frac{T_dA_k}{1-e^{s_kT_d}z^{-1}} \]

ii. 双线性变换法
\[ s=\frac{2}{T_d}\left(\frac{1-z^{-1}}{1+z^{-1}}\right) \]

(4) 特点

i. 脉冲响应不变法: 无频率畸变, 有频响混叠, 不适用非带限的模拟滤波器.

ii. 双线性变换法: 有频率畸变(注意预畸变), 无频响混叠, 不适用微分器/积分器.

线性相位FIR滤波器的设计步骤(窗函数法)

设计之前分析通带阻带等波动频响特性, 明确窗形状和窗长对通带阻带误差以及过渡带宽的影响, 注意四类线性相位FIR各自适用的滤波器类型.

(1) 利用IDTFT求出\(h_d[n]=\mathscr{F}^{-1}[H_d(e^{j\omega})]\).

(2) 根据最大逼近误差\(\delta=\min\{\delta_p,\delta_s\}\)确定窗形状.

(3) 根据主瓣宽度\(\Delta\omega_{m}\approx2|\omega_p-\omega_s|\)确定窗长\(L\).

(4) 求\(h[n]=h_d[n]\cdot w[n]\), 检验\(H(e^{j\omega})=\mathscr{F}[h(n)]\)是否满足所给滤波器的性能要求, 如不满足, 则需考虑改变窗形状或改变窗长\(L\), 到满足要求为止.

给定一个连续时间全通系统
\[ H(s)=K\frac{s-a}{s+a} \]

若利用双线性变换法将其变换成一个离散时间系统, 得到的也是一个全通系统.
\[ H(z)=-K\frac{z^{-1}-b}{1-bz^{-1}} \]

其中
\[ b=\frac{\dfrac{2}{T_d}-a}{\dfrac{2}{T_d}+a} \]

假设\(H_c(s)\)在\(s=s_0\)处有一个\(r\)阶极点, 使得\(H_c(s)\)可以表示成
\[ H_c(s)=\sum_{k=1}^{r}\frac{A_k}{(s-s_0)^k}+G_c(s) \]

式中\(G_c(s)\)只有一阶极点. 设\(H_c(s)\)是因果的, 则有
\[ A_k=\frac{1}{(r-k)!}\cdot\frac{\mathrm{d}^{r-k}}{\mathrm{d}s^{r-k}}[(s-s_0)^rH_c(s)]\big|_{s=s_0} \]

若\(g_c(t)\)为\(G_c(s)\)的\(s\)反变换, 则\(H_c(s)\)的脉冲响应可表示为
\[ h_c(t)=\mathscr{L}^{-1}[H_c(s)]=\sum_{k=1}^{r}\frac{A_kt^{k-1}}{(k-1)!}e^{s_0t}\varepsilon(t)+g_c(t) \]

阶跃响应不变法设计IIR滤波器

连续时间单位脉冲响应对应的\(s\)变换
\[ H_c(s)=\sum_{k=1}^{N}\frac{A_k}{s-s_k} \]

连续时间单位阶跃响应对应的\(s\)变换
\[ S_c(s)=\frac{1}{s}H_c(s)=\sum_{k=1}^{N}\frac{A'_k}{s-s_k}+\frac{B}{s} \]

离散时间单位阶跃响应对应的\(z\)变换
\[ S(z)=\sum_{k=1}^{N}\frac{T_dA'_k}{1-e^{s_kT_d}z^{-1}}+\frac{T_dB}{1-z^{-1}} \]

离散时间滤波器的系统函数
\[ H(z)=(1-z^{-1})S(z)=(1-z^{-1})\sum_{k=1}^{N}\frac{T_dA'_k}{1-e^{s_kT_d}z^{-1}}+T_dB \]

自相关函数不变法设计IIR滤波器

(1) 考虑一个具有脉冲响应\(h_c(t)\)和系统函数\(H_c(s)\)的稳定连续时间系统.该系统的自相关函数定义为
\[ \phi_c(\tau)=\int_{-\infty}^{\infty}h_c(t)h_c(t+\tau)\mathrm{d}t \]

对于实脉冲响应, 其拉氏变换为
\[ \begin{aligned} \Phi_c(s)&=\int_{-\infty}^{\infty}\phi_c(\tau)e^{-s\tau}\mathrm{d}\tau\\ &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}h_c(t)h_c(t+\tau)e^{-s\tau}\mathrm{d}\tau\mathrm{d}t\\ &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}h_c(t)h_c(t+\tau)e^{-s(\tau+t)}e^{st}\mathrm{d}\tau\mathrm{d}t\\ &=\int_{-\infty}^{\infty}h_c(t+\tau)e^{-s(\tau+t)}\mathrm{d}(t+\tau)\cdot\int_{-\infty}^{\infty}h_c(t)e^{st}\mathrm{d}t\\ &=H_c(s)H_c(-s) \end{aligned} \]

考虑一个具有脉冲响应\(h[n]\)和系统函数\(H(z)\)的稳定离散时间系统.该系统的自相关函数定义为
\[ \phi[m]=\sum_{n=-\infty}^{\infty}h[n]h[n+m] \]

对于实脉冲响应, 其\(z\)变换为
\[ \begin{aligned} \Phi(z)&=\sum_{m=-\infty}^{\infty}\phi[m]z^{-m}\\ &=\sum_{m=-\infty}^{\infty}\sum_{n=-\infty}^{\infty}h[n]h[n+m]z^{-m}\\ &=\sum_{m=-\infty}^{\infty}\sum_{n=-\infty}^{\infty}h[n]h[n+m]z^{-(m+n)}z^{n}\\ &=\sum_{m=-\infty}^{\infty}h[n+m]z^{-(m+n)}\cdot\sum_{n=-\infty}^{\infty}h[n]z^{n}\\ &=H(z)H(z^{-1}) \end{aligned} \]

(2) 自相关函数不变法的时域关系
\[ \phi[m]=T_d\phi_c(mT_d),\infty<m<\infty \]

(3) 当\(H_c(s)\)是一个有理函数, 在\(s_k(k=1,2,\cdots,N)\)处有\(N\)个一阶极点并有\(M<N\)个零点时, 自相关函数不变法的设计步骤

i. 根据\(\Phi_c(s)=H_c(s)H_c(-s)\)求出\(\phi_c(\tau)\)的\(s\)变换, 然后进行部分分式展开, 形式为
\[ \Phi_c(s)=\sum_{k=1}^{N}\left(\frac{A_k}{s-s_k}+\frac{B_k}{s+s_k}\right) \]

ii. 根据时域关系求出\(\phi[m]\)的\(z\)变换, 其具有和脉冲响应不变法相同的形式, 即
\[ \Phi(z)=\sum_{k=1}^{N}\left(\frac{T_dA_k}{1-e^{s_kT_d}z^{-1}}+\frac{T_dB_k}{1-e^{-s_kT_d}z^{-1}}\right) \]

iii. 根据\(\Phi(z)=H(z)H(z^{-1})\)得到\(H(z)\). 求出\(\Phi(z)\)的极点和零点, \(H(z)\)是由\(\Phi(z)\)在单位圆内的极点和零点构成的最小相位系统函数.

(4) 自相关函数不变法的频域关系
\[ \Phi(e^{j\omega})=\sum_{k=-\infty}^{\infty}\Phi_c\left(j\frac{\omega}{T_d}+j\frac{2\pi}{T_d}k\right) \]

因此, 当频域不发生混叠, 即
\[ \Phi_c(j\Omega)\simeq0,|\Omega|\geq\dfrac{\pi}{T_d} \]

此时\(\Phi(e^{j\omega})\simeq\Phi(j\dfrac{\omega}{T_d})\), 因此\(|H(e^{j\omega})|^2\simeq|H(j\dfrac{\omega}{T_d})|^2\).

(5) \(H(z)\)与全通系统相乘后的新系统依然满足\(\Phi(z)=H(z)H(z^{-1})\), 因此, 使用自相关函数不变法得出的系统函数不唯一.

---

上一篇：离散时间系统结构 · 下一篇：离散傅里叶变换