数字信号处理入门:FIR 滤波器线性相位、窗函数法与频率采样法

系统整理 FIR 滤波器的线性相位条件、四种对称类型、窗函数法设计步骤,以及频率采样法与频域采样定理之间的联系。

数字信号处理入门:FIR 滤波器线性相位、窗函数法与频率采样法

这篇文章面向刚入门数字信号处理的读者,主题是数字 FIR 滤波器设计方法。前面学习 IIR 滤波器时,重点往往是如何借助模拟滤波器设计数字滤波器;而在 FIR 滤波器中,最核心的关键词变成了“线性相位”“有限长冲激响应”“窗函数”和“频率采样”。

FIR 的英文是 Finite Impulse Response,意思是有限长单位冲激响应。通俗地说,FIR 滤波器的输出只由有限个输入样本的加权和决定,不需要把过去输出再反馈回来。因此,FIR 滤波器通常结构稳定,而且特别容易设计成严格线性相位,这一点对图像处理、数据传输、测量信号分析等场景非常重要。

IIR 与 FIR 的基本取舍

为什么还要学习 FIR 滤波器

IIR 数字滤波器有一个明显优点:设计方便。它可以利用成熟的模拟滤波器设计结果,例如巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器,然后通过脉冲响应不变法或双线性变换法得到数字滤波器。因此,在很多只关心幅频特性、不太关心波形形状的场合,IIR 滤波器非常经济。

但是 IIR 滤波器的缺点也很明显:相位一般不是线性的。如果一个系统只改变不同频率分量的幅度,而不希望改变它们之间的相对时间关系,那么非线性相位就可能造成波形失真。要让 IIR 具有线性相位,通常需要额外增加相位校正网络,系统会变复杂。

FIR 滤波器的优势正好在这里。它可以获得严格的线性相位,同时还能具有任意幅度特性。由于 FIR 通常没有反馈回路,所以它是无条件稳定的。进行长序列滤波时,还可以使用 FFT 加速计算。它的不足是,在同样幅频指标下,所需阶数通常比 IIR 高,计算量和延时也更大。

相位失真从哪里来

设输入信号的频谱为 X(ejω)X(e^{j\omega}),滤波器频率响应为 H(ejω)H(e^{j\omega}),输出频谱为 Y(ejω)Y(e^{j\omega})。线性时不变系统满足

Y(ejω)=X(ejω)H(ejω).Y(e^{j\omega})=X(e^{j\omega})H(e^{j\omega}).

如果把这三个量都写成幅度和相位的形式:

X(ejω)=X(ejω)ejargX(ejω),X(e^{j\omega})=|X(e^{j\omega})|e^{j\arg X(e^{j\omega})}, H(ejω)=H(ejω)ejθ(ω),H(e^{j\omega})=|H(e^{j\omega})|e^{j\theta(\omega)}, Y(ejω)=Y(ejω)ejargY(ejω).Y(e^{j\omega})=|Y(e^{j\omega})|e^{j\arg Y(e^{j\omega})}.

两边取模,得到

Y(ejω)=X(ejω)H(ejω).|Y(e^{j\omega})|=|X(e^{j\omega})|\cdot |H(e^{j\omega})|.

这说明幅频响应 H(ejω)|H(e^{j\omega})| 决定了不同频率分量被放大或衰减多少。两边取相位,得到

argY(ejω)=argX(ejω)+θ(ω).\arg Y(e^{j\omega})=\arg X(e^{j\omega})+\theta(\omega).

这说明相频特性 θ(ω)\theta(\omega) 决定了不同频率分量的相位变化。

如果一个信号由多个频率成分组成,而这些频率成分通过滤波器后获得的时间延迟不同,那么原来的波形结构就会被破坏。这就是相位失真的直观来源。

线性相位与群时延

为了不产生相位失真,希望不同频率分量通过滤波器后具有相同的时间延迟。这个时间延迟称为群时延,定义为

τ=dθ(ω)dω.\tau=-\frac{d\theta(\omega)}{d\omega}.

如果 τ\tau 是常数,那么对 ω\omega 积分可以得到线性相位。常见的准确线性相位有两类:

θ(ω)=τω,\theta(\omega)=-\tau\omega, θ(ω)=θ0τω,θ0=π2.\theta(\omega)=\theta_0-\tau\omega,\qquad \theta_0=-\frac{\pi}{2}.

第一类线性相位是一条过原点的直线;第二类线性相位是一条带有 π/2-\pi/2 截距的直线。二者的共同点是斜率都是 τ-\tau,因此群时延相同。

线性相位与恒定群时延

FIR 滤波器的基本形式

FIR 滤波器的单位冲激响应 h(n)h(n) 是有限长序列,通常规定

0nN1.0\leq n\leq N-1.

它的系统函数为

H(z)=n=0N1h(n)zn,H(z)=\sum_{n=0}^{N-1}h(n)z^{-n},

频率响应为

H(ejω)=n=0N1h(n)ejωn.H(e^{j\omega})=\sum_{n=0}^{N-1}h(n)e^{-j\omega n}.

h(n)h(n) 是实序列时,可以把频率响应写成

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

这里 Hg(ω)H_g(\omega) 称为幅度特性。注意,Hg(ω)H_g(\omega) 与普通幅度响应 H(ejω)|H(e^{j\omega})| 不完全一样,因为 Hg(ω)H_g(\omega) 可以为正,也可以为负;而 H(ejω)|H(e^{j\omega})| 总是非负。

第一类线性相位对应的时域条件

第一类线性相位满足

θ(ω)=ωτ.\theta(\omega)=-\omega\tau.

如果

H(ejω)=Hg(ω)ejωτ,H(e^{j\omega})=H_g(\omega)e^{-j\omega\tau},

另一方面又有

H(ejω)=n=0N1h(n)[cos(ωn)jsin(ωn)].H(e^{j\omega})=\sum_{n=0}^{N-1}h(n)[\cos(\omega n)-j\sin(\omega n)].

把实部和虚部分别比较,可以得到

Hg(ω)cos(ωτ)=n=0N1h(n)cos(ωn),H_g(\omega)\cos(\omega\tau)=\sum_{n=0}^{N-1}h(n)\cos(\omega n), Hg(ω)sin(ωτ)=n=0N1h(n)sin(ωn).H_g(\omega)\sin(\omega\tau)=\sum_{n=0}^{N-1}h(n)\sin(\omega n).

进一步整理后,得到

n=0N1h(n)sin[ω(nτ)]=0.\sum_{n=0}^{N-1}h(n)\sin[\omega(n-\tau)]=0.

要让这个式子对所有频率都成立,一个自然的选择是让

τ=N12,\tau=\frac{N-1}{2},

并要求 h(n)h(n) 关于 (N1)/2(N-1)/2 偶对称,即

h(n)=h(N1n),0nN1.h(n)=h(N-1-n),\qquad 0\leq n\leq N-1.

因此,第一类线性相位对应的时域条件可以概括为:冲激响应偶对称,平均延迟为 (N1)/2(N-1)/2

第二类线性相位对应的时域条件

第二类线性相位满足

θ(ω)=π2ωτ.\theta(\omega)=-\frac{\pi}{2}-\omega\tau.

类似地,把频率响应的实部和虚部分别比较,最终会得到

n=0N1h(n)cos[ω(nτ)]=0.\sum_{n=0}^{N-1}h(n)\cos[\omega(n-\tau)]=0.

这时仍取

τ=N12,\tau=\frac{N-1}{2},

但要求 h(n)h(n) 关于 (N1)/2(N-1)/2 奇对称,即

h(n)=h(N1n),0nN1.h(n)=-h(N-1-n),\qquad 0\leq n\leq N-1.

所以,第二类线性相位对应的时域条件可以概括为:冲激响应奇对称,平均延迟仍为 (N1)/2(N-1)/2,相位直线在零频处多出 π/2-\pi/2 的截距。

线性相位 FIR 的充要条件

FIR 滤波器具有线性相位的充要条件为

τ=N12,\tau=\frac{N-1}{2}, h(n)=±h(N1n),0nN1.h(n)=\pm h(N-1-n),\qquad 0\leq n\leq N-1.

其中正号表示偶对称,负号表示奇对称。由于 h(n)h(n) 可以是偶对称或奇对称,而长度 NN 又可以是奇数或偶数,因此线性相位 FIR 滤波器共有四种类型。

线性相位 FIR 的四种时域对称形式

四种线性相位 FIR 的幅度特性

偶对称、奇数长度

h(n)h(n) 为偶对称,且 NN 为奇数时,存在一个正好落在整数采样点上的对称中心。记

τ=N12.\tau=\frac{N-1}{2}.

频率响应可以写成

H(ejω)=Hg(ω)ejωτ,H(e^{j\omega})=H_g(\omega)e^{-j\omega\tau},

其中

Hg(ω)=h(τ)+n=0M12h(n)cos[ω(nτ)].H_g(\omega)=h(\tau)+\sum_{n=0}^{M-1}2h(n)\cos[\omega(n-\tau)].

由于余弦项关于 ω=0,π,2π\omega=0,\pi,2\pi 都具有偶对称性质,所以 Hg(ω)H_g(\omega) 也关于这些频点偶对称。这类 FIR 适合设计低通、高通、带通、带阻等常见选频滤波器。

偶对称、偶数长度

h(n)h(n) 为偶对称,且 NN 为偶数时,对称中心落在两个采样点之间,没有单独的中心项。这时

Hg(ω)=n=0M2h(n)cos[ω(nτ)].H_g(\omega)=\sum_{n=0}^{M}2h(n)\cos[\omega(n-\tau)].

由于 NN 为偶数,

τ=N212.\tau=\frac{N}{2}-\frac{1}{2}.

ω=π\omega=\pi 处,相关余弦项会导致 Hg(π)=0H_g(\pi)=0。因此,这类滤波器不适合设计在 ω=π\omega=\pi 处幅度不为零的滤波器,例如高通和带阻;它更适合低通和带通。

奇对称、奇数长度

h(n)h(n) 为奇对称,且 NN 为奇数时,中心点的值必须为 0。此时幅度函数可写成正弦项形式:

Hg(ω)=n=0M12h(n)sin[ω(nτ)].H_g(\omega)=\sum_{n=0}^{M-1}2h(n)\sin[\omega(n-\tau)].

ω=0,π,2π\omega=0,\pi,2\pi 时,正弦项为 0,所以 Hg(ω)H_g(\omega) 在这些点为 0。这类滤波器适合设计带通滤波器,也常用于需要奇对称特性的系统,例如 Hilbert 变换器、微分器等。

奇对称、偶数长度

h(n)h(n) 为奇对称,且 NN 为偶数时,同样有

Hg(ω)=n=0M2h(n)sin[ω(nτ)].H_g(\omega)=\sum_{n=0}^{M}2h(n)\sin[\omega(n-\tau)].

这类幅度函数在 ω=0\omega=02π2\pi 处为 0,但在 ω=π\omega=\pi 处不必为 0。因此,它适合设计高通和带通滤波器。

将四种情况记成一句话:偶对称适合常规选频,奇对称适合在零频处必须为零的响应;NN 的奇偶性决定了 ω=π\omega=\pi 处是否被强制为零。

线性相位 FIR 的零点分布

线性相位 FIR 满足

h(n)=±h(N1n).h(n)=\pm h(N-1-n).

将这个关系代入系统函数

H(z)=n=0N1h(n)zn,H(z)=\sum_{n=0}^{N-1}h(n)z^{-n},

可得到

H(z)=±z(N1)H(z1).H(z)=\pm z^{-(N-1)}H(z^{-1}).

这个式子说明,如果 ziz_iH(z)H(z) 的零点,那么 zi1z_i^{-1} 也是零点。又因为 h(n)h(n) 是实序列,所以零点还必须共轭成对。于是对于一般复零点 ziz_i,还会同时出现

zi,zi1,(zi1).z_i^*,\qquad z_i^{-1},\qquad (z_i^{-1})^*.

因此,线性相位 FIR 的零点通常呈“共轭 + 倒数”的成组分布。知道其中一个零点,其他相关零点也就基本确定了。

线性相位 FIR 零点的共轭倒数分布

FIR 设计的三类逼近方法

设计 FIR 滤波器时,通常先给出理想频率响应

Hd(ejω),H_d(e^{j\omega}),

然后希望找到一个 NN 点有限长序列 h(n)h(n),使得

H(ejω)=n=0N1h(n)ejωnH(e^{j\omega})=\sum_{n=0}^{N-1}h(n)e^{-j\omega n}

尽可能逼近 Hd(ejω)H_d(e^{j\omega})

常见逼近方法有三类:窗函数设计法、频率采样法和最优化设计。窗函数法是时域逼近,频率采样法是频域逼近,最优化设计则常用于等波纹逼近。

窗函数设计法的基本思想

理想单位冲激响应 hd(n)h_d(n) 可以由理想频率响应通过傅里叶反变换得到:

hd(n)=12πππHd(ejω)ejωndω.h_d(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}H_d(e^{j\omega})e^{j\omega n}d\omega.

对于理想低通、高通、带通、带阻等滤波器,理想频率响应往往是分段常数,在边界频率处有突变。因此,其反变换 hd(n)h_d(n) 往往是无限长、非因果序列。而实际可实现的 FIR 滤波器必须是有限长、因果序列。

怎样用有限长 h(n)h(n) 逼近无限长 hd(n)h_d(n)?最直接的方法就是截取一段。这个截取过程可以看作通过一个“窗口”观察 hd(n)h_d(n),所以称为窗函数法。

h(n)=w(n)hd(n).h(n)=w(n)h_d(n).

如果 w(n)w(n) 是矩形窗,就相当于直接截取;如果 w(n)w(n) 不是矩形窗,则相当于在截取范围内对 hd(n)h_d(n) 进行加权处理。

窗函数法设计流程

理想低通的单位冲激响应

以截止频率为 ωc\omega_c 的理想低通滤波器为例。若希望它具有线性相位,理想频率响应可写成

Hd(ejω)={ejωα,ωωc,0,ωc<ωπ.H_d(e^{j\omega})= \begin{cases} e^{-j\omega\alpha}, & |\omega|\leq \omega_c,\\ 0, & \omega_c<|\omega|\leq \pi. \end{cases}

这里 α\alpha 是滤波器的延迟时间。由傅里叶反变换可得

hd(n)=12πωcωcejωαejωndω.h_d(n)=\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c}e^{-j\omega\alpha}e^{j\omega n}d\omega.

整理后得到

hd(n)=sin[ωc(nα)]π(nα).h_d(n)=\frac{\sin[\omega_c(n-\alpha)]}{\pi(n-\alpha)}.

也可以写为

hd(n)=ωcπSa[ωc(nα)].h_d(n)=\frac{\omega_c}{\pi}\operatorname{Sa}[\omega_c(n-\alpha)].

这个序列是以 α\alpha 为中心的 sinc 型无限长序列。为了让截取后的 h(n)h(n) 仍具有线性相位,需要令

α=N12.\alpha=\frac{N-1}{2}.

然后取

h(n)={hd(n),0nN1,0,其他.h(n)= \begin{cases} h_d(n), & 0\leq n\leq N-1,\\ 0, & \text{其他}. \end{cases}

理想低通冲激响应与有限长度截取

矩形窗的频率响应

矩形窗可写成

wR(n)=RN(n).w_R(n)=R_N(n).

它的频率响应为

WR(ejω)=ejω(N1)/2sin(Nω/2)sin(ω/2).W_R(e^{j\omega}) = e^{-j\omega (N-1)/2} \frac{\sin(N\omega/2)}{\sin(\omega/2)}.

其中

WRg(ω)=sin(Nω/2)sin(ω/2)W_{Rg}(\omega)=\frac{\sin(N\omega/2)}{\sin(\omega/2)}

是起主要作用的幅度函数。矩形窗的主瓣宽度决定了过渡带的宽度,旁瓣高度决定了通带和阻带的波动程度。

窗函数截取本质上是时域相乘,因此在频域中会变成卷积。理想低通的矩形频率响应和窗函数频率响应卷积后,原来突变的边沿会被抹宽,形成过渡带;通带和阻带会出现波动。

吉布斯效应与窗函数取舍

加窗处理对理想矩形频率响应的影响主要有三点。

第一,窗函数的频谱会使滤波器频率响应出现通带波动和阻带波动,在过渡带两侧波动最大。矩形窗对应的最大波动约为 8.95%,并且与 NN 无关,这称为吉布斯效应。

第二,理想频率响应的不连续边沿会被加宽,形成过渡带。过渡带宽度大致由窗函数频率响应的主瓣宽度决定。

第三,增大窗口长度 NN 可以减小主瓣宽度,因此能减小过渡带宽度;但如果窗函数形状不变,主瓣与旁瓣的相对比例不会改变,所以不能从根本上降低波动。

因此,调整 NN 主要控制过渡带宽度;想减少通带和阻带波动、提高阻带衰减,则需要改变窗函数形状。一般来说,旁瓣越低,阻带衰减越大;但代价是主瓣往往变宽,过渡带变宽。

常见窗函数形状

窗函数主瓣与旁瓣的折中关系

常见窗函数

矩形窗

矩形窗最简单:

w(n)=RN(n).w(n)=R_N(n).

它的主瓣较窄,过渡带较窄,但旁瓣较高,阻带衰减较差。

Bartlett 窗

Bartlett 窗也叫三角窗,它的形状从两端逐渐上升到中间,再逐渐下降。相比矩形窗,它降低了旁瓣,但主瓣变宽。

Hanning 窗

Hanning 窗也称汉宁窗或升余弦窗,形式为

w(n)=0.5[1cos(2πnN1)]RN(n).w(n)=0.5\left[1-\cos\left(\frac{2\pi n}{N-1}\right)\right]R_N(n).

它比矩形窗具有更低旁瓣,因此阻带衰减更好,但过渡带比矩形窗更宽。

Hamming 窗

Hamming 窗也叫哈明窗,是改进的升余弦窗:

w(n)=[0.540.46cos(2πnN1)]RN(n).w(n)=\left[0.54-0.46\cos\left(\frac{2\pi n}{N-1}\right)\right]R_N(n).

它通常能提供比 Hanning 窗更大的阻带衰减。

Blackman 窗与 Kaiser 窗

Blackman 窗进一步降低旁瓣,阻带衰减更大,但过渡带更宽。Kaiser 窗比较灵活,可以通过参数 β\beta 调节主瓣宽度和旁瓣衰减,是工程中常见的可调窗函数。

常见窗函数参数可概括如下:

窗函数旁瓣峰值/dB过渡带近似宽度过渡带精确宽度阻带最小衰减/dB
矩形窗-134π/N4\pi/N1.8π/N1.8\pi/N-21
三角窗-258π/N8\pi/N6.1π/N6.1\pi/N-25
汉宁窗-318π/N8\pi/N6.2π/N6.2\pi/N-44
哈明窗-418π/N8\pi/N6.6π/N6.6\pi/N-53
布莱克曼窗-5712π/N12\pi/N11π/N11\pi/N-74
Kaiser 窗 (β=7.865)(\beta=7.865)-57视参数而定10π/N10\pi/N-80

这张表的使用方法很简单:如果题目给出阻带最小衰减,就先选一个能满足该衰减的窗函数;如果题目给出过渡带宽度,就用该窗函数的过渡带公式估计 NN

窗函数法设计 FIR 的步骤

窗函数法的基本步骤如下。

第一,根据过渡带宽度和阻带衰减指标选择窗函数类型,并估计窗口长度 NN

第二,构造希望逼近的理想频率响应 Hd(ejω)H_d(e^{j\omega})

第三,计算理想单位冲激响应

hd(n)=12πππHd(ejω)ejωndω.h_d(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}H_d(e^{j\omega})e^{j\omega n}d\omega.

第四,加窗得到

h(n)=hd(n)w(n).h(n)=h_d(n)w(n).

第五,检查实际频率响应是否满足误差指标。如果不满足,就重新选择窗函数类型或窗口长度。

解答题例题:用窗函数法设计高通 FIR 滤波器

要求设计线性相位高通 FIR 数字滤波器,指标为

ωp=π2,ωs=π4,\omega_p=\frac{\pi}{2},\qquad \omega_s=\frac{\pi}{4}, αp=1 dB,αs=40 dB.\alpha_p=1\text{ dB},\qquad \alpha_s=40\text{ dB}.

首先根据阻带最小衰减 αs=40 dB\alpha_s=40\text{ dB},可选择 Hanning 窗或 Hamming 窗。这里选择 Hanning 窗。过渡带宽度为

Btωpωs=π4.B_t\leq \omega_p-\omega_s=\frac{\pi}{4}.

Hanning 窗的精确过渡带宽度为

Bt=6.2πN.B_t=\frac{6.2\pi}{N}.

为了满足要求,需要

6.2πNπ4,\frac{6.2\pi}{N}\leq \frac{\pi}{4},

所以

N24.8.N\geq 24.8.

高通滤波器通常要求选择奇数长度,因此取

N=25.N=25.

Hanning 窗为

w(n)=0.5[1cos(πn12)]R25(n).w(n)=0.5\left[1-\cos\left(\frac{\pi n}{12}\right)\right]R_{25}(n).

由于

τ=N12=12,\tau=\frac{N-1}{2}=12,

理想高通的截止频率取通带和阻带边界的中点:

ωc=ωp+ωs2=3π8.\omega_c=\frac{\omega_p+\omega_s}{2}=\frac{3\pi}{8}.

理想高通可以看作全通减去低通,因此

hd(n)=δ(n12)sin[3π(n12)/8]π(n12).h_d(n)=\delta(n-12)-\frac{\sin[3\pi(n-12)/8]}{\pi(n-12)}.

最后加窗得到

h(n)=hd(n)w(n).h(n)=h_d(n)w(n).

也就是

h(n)={δ(n12)sin[3π(n12)/8]π(n12)}0.5[1cos(πn12)]R25(n).h(n)= \left\{ \delta(n-12)-\frac{\sin[3\pi(n-12)/8]}{\pi(n-12)} \right\} \cdot 0.5\left[1-\cos\left(\frac{\pi n}{12}\right)\right]R_{25}(n).

这个例题体现了窗函数法的典型套路:先由阻带衰减选窗,再由过渡带宽度估算 NN,再写出理想响应,最后加窗。

解答

这道题的核心步骤是四段式。

第一,先由阻带衰减

αs=40 dB\alpha_s=40\text{ dB}

选择合适窗函数。这里选 Hanning 窗,是因为它能满足题目要求的阻带衰减水平。

第二,用过渡带宽度估算长度 NN。题目给出

ωp=π2,ωs=π4,\omega_p=\frac{\pi}{2},\qquad \omega_s=\frac{\pi}{4},

所以允许的过渡带宽度为

Btωpωs=π4.B_t\le \omega_p-\omega_s=\frac{\pi}{4}.

再代入 Hanning 窗的经验公式

Bt=6.2πN,B_t=\frac{6.2\pi}{N},

即可求得

N24.8,N\ge 24.8,

因此取奇数长度

N=25.N=25.

第三,写出理想高通滤波器的冲激响应。先取延迟

τ=N12=12,\tau=\frac{N-1}{2}=12,

再取截止频率中点

ωc=ωp+ωs2=3π8,\omega_c=\frac{\omega_p+\omega_s}{2}=\frac{3\pi}{8},

从而得到理想高通响应 hd(n)h_d(n)

第四,用窗函数加权:

h(n)=hd(n)w(n).h(n)=h_d(n)w(n).

这一题真正要记住的是设计顺序:先按阻带选窗,再按过渡带宽度定长度,最后写理想响应并加窗。

课堂练习整理

FIR 数字滤波器与 IIR 数字滤波器相比,最大的优点是可以保证系统具有线性相位特性。

关于 IIR 和 FIR 的说法中,“IIR 滤波器可以得到严格线性相位”通常是错误的。IIR 具有递归结构,可用较低阶数获得较高选择性,但严格线性相位不是它的天然优势。

关于 FIR 的说法中,“FIR 滤波器容易设计成线性相位特性”是正确的;“FIR 的冲激响应长度无限”是错误的,因为 FIR 本身就表示有限长冲激响应;“FIR 可直接利用模拟滤波器间接设计”也不是它的主要方法;相同幅频要求下,FIR 阶数一般比 IIR 高。

如果 FIR 滤波器满足

h(n)=h(N1n),h(n)=h(N-1-n),

NN 为奇数,那么它属于偶对称、奇数长度类型,可以实现低通、高通、带通、带阻滤波器。

窗函数法设计 FIR 的步骤可以概括为:先选窗函数类型,再确定窗口长度,再构造理想频率响应,求 hd(n)h_d(n),最后加窗得到 h(n)h(n)。选择窗函数类型的依据主要是阻带最小衰减,选择窗函数长度的依据主要是过渡带宽度。

频率采样法为什么更直接

窗函数法从时域出发,先求理想冲激响应 hd(n)h_d(n),再截断成有限长 h(n)h(n)。但是实际滤波器指标通常是在频域给出的,例如要求某些频段通过,某些频段衰减。因此,从频域出发的频率采样法更加直接。

频率采样法特别适合两种情况:第一,理想频率响应 Hd(ejω)H_d(e^{j\omega}) 公式比较复杂;第二,理想频率响应无法用封闭公式表示,只能用若干离散频率点上的数值表示。

频率采样法的基本原理

频率采样法从频域出发,把给定的理想频率响应等间隔采样:

Hd(k)=Hd(ejω)ω=2πNk,k=0,1,2,,N1.H_d(k)=H_d(e^{j\omega})\bigg|_{\omega=\frac{2\pi}{N}k}, \qquad k=0,1,2,\cdots,N-1.

再对 Hd(k)H_d(k)NN 点 IDFT,得到

h(n)=1Nk=0N1Hd(k)ej2πNkn.h(n)=\frac{1}{N}\sum_{k=0}^{N-1}H_d(k)e^{j\frac{2\pi}{N}kn}.

这样得到的 h(n)h(n) 就是一个有限长 FIR 序列。

从系统函数角度看,也可以利用频率采样值直接构造

H(z)=1zNNk=0N1Hd(k)1ej2πNkz1.H(z)=\frac{1-z^{-N}}{N}\sum_{k=0}^{N-1} \frac{H_d(k)}{1-e^{j\frac{2\pi}{N}k}z^{-1}}.

这个形式适合频率采样结构型 FIR。直接由 h(n)h(n) 实现,则适合直接型 FIR。

频率采样法示意图

窗函数法与频率采样法的区别

窗函数法从时域出发。它先得到理想的 hd(n)h_d(n),再用窗函数截成有限长的 h(n)h(n),然后用这个 h(n)h(n) 逼近理想频率响应。

频率采样法从频域出发。它先对 Hd(ejω)H_d(e^{j\omega}) 等间隔采样,得到 H(k)H(k),再由这些频域样本通过 IDFT 唯一确定 h(n)h(n)

两种方法的核心区别是:窗函数法控制的是时域截取方式,频率采样法控制的是频域采样点。前者更适合规则理想响应,后者更适合频域指标以离散数据给出的情况。

频率采样法中的线性相位约束

频率采样法若要设计线性相位 FIR,也必须满足线性相位的时域条件:

h(n)=h(N1n)h(n)=h(N-1-n)

h(n)=h(N1n).h(n)=-h(N-1-n).

对于偶对称 FIR,如果 NN 为奇数,幅度函数 Hg(k)H_g(k) 必须满足偶对称:

Hg(k)=Hg(Nk).H_g(k)=H_g(N-k).

这对应于幅度函数关于 N/2N/2 偶对称。

如果 NN 为偶数,则幅度函数应满足

Hg(k)=Hg(Nk),H_g(k)=-H_g(N-k),

并且

Hg(N/2)=0.H_g(N/2)=0.

也就是说,偶数长度情况下,ω=π\omega=\pi 对应的采样点幅度受到约束。

对于理想低通,若采样点数 NN 为奇数,可以先根据截止频率确定通过区间的最大整数 kck_c,再设置低频采样点为 1,高频采样点为 0,并按对称关系补齐另一半。

对于高通和带阻滤波器,通常要求 NN 取奇数,因为这类滤波器往往要求 ω=π\omega=\pi 处不为 0,而偶数长度偶对称结构可能会强制 ω=π\omega=\pi 处为 0。

IIR 与 FIR 的比较

IIR 与 FIR 滤波器对比

从性能上比较

IIR 滤波器可以用较低阶数获得较高选择性,所需存储单元少,计算量小,因此经济高效。但是这种效率往往以相位非线性为代价。选择性越好,相位非线性可能越严重。

FIR 滤波器可以得到严格线性相位,但如果要求高选择性,通常需要较高阶数。对于同样指标,FIR 阶数可能比 IIR 高 5 到 10 倍,因此成本较高,信号延时也较大。

从结构上比较

IIR 必须采用递归结构,极点位置必须在单位圆内,否则系统不稳定。由于存在反馈,有限字长效应也需要格外注意。

FIR 主要采用非递归结构。有限长冲激响应决定了它本身稳定,不论从理论上还是在实际有限精度运算中,稳定性问题都更容易控制。FIR 还可以采用 FFT 算法提高运算速度。

从设计工具上比较

IIR 可以借助模拟滤波器设计成果,一般有封闭形式的公式,计算较直接,对计算工具要求不高。

FIR 的设计通常没有完全封闭的统一公式。窗函数法虽然能给出窗函数公式,但通带、阻带衰减等指标通常不能用简单显式公式完全准确表达,往往需要数值验证。

从应用灵活性比较

IIR 更适合低通、高通、带通、带阻等幅频特性比较规则的选频滤波器。

FIR 更灵活,尤其适合对波形保真度要求高的场合,也适合构成微分器、Hilbert 变换器、特殊幅频响应滤波器等。

实际选择时可以按一句话判断:如果主要追求效率、相位要求不高,可以选 IIR;如果波形携带重要信息、要求线性相位,优先考虑 FIR。

补充复习:频域采样理论

讲义最后补充了频域采样理论,它是理解 DFT、频率采样法和 FIR 设计的重要基础。

时域抽样定理告诉我们,若想从抽样序列无失真恢复连续信号,抽样频率必须满足

fs2fh.f_s\geq 2f_h.

恢复连续信号时,可用抽样内插公式

xa(t)=m=xa(mT)sin[π(tmT)/T]π(tmT)/T.x_a(t)=\sum_{m=-\infty}^{\infty}x_a(mT) \frac{\sin[\pi(t-mT)/T]}{\pi(t-mT)/T}.

这说明时域采样会导致频域周期延拓;如果采样频率不足,就会发生频域混叠。

而频域采样理论讨论的是相反问题:如果在频域对 X(z)X(z)X(ejω)X(e^{j\omega}) 等间隔采样,时域会发生什么?

对于有限长序列 x(n)x(n),其 DFT 样本为

X(k)=X(z)z=ej2πNk.X(k)=X(z)\bigg|_{z=e^{j\frac{2\pi}{N}k}}.

也可以理解为

X(k)=X(ejω)ω=2πNk.X(k)=X(e^{j\omega})\bigg|_{\omega=\frac{2\pi}{N}k}.

也就是说,DFT 是对 X(z)X(z) 在单位圆上的 NN 个等分点采样,或者说是对 DTFT 在 [0,2π][0,2\pi] 上作 NN 点等间隔采样。

频域采样定理

设任意序列 x(n)x(n) 的 Z 变换为

X(z)=n=x(n)zn,X(z)=\sum_{n=-\infty}^{\infty}x(n)z^{-n},

并且收敛域包含单位圆。若在单位圆上对 X(z)X(z) 等间隔采样 NN 点:

X(k)=X(z)z=WNk=n=x(n)WNnk,X(k)=X(z)\bigg|_{z=W_N^{-k}} = \sum_{n=-\infty}^{\infty}x(n)W_N^{nk},

再把 X(k)X(k) 看成长度为 NN 的序列并作 IDFT,得到 xN(n)x_N(n)。这时 xN(n)x_N(n) 与原序列的关系为

xN(n)=r=x(n+rN).x_N(n)=\sum_{r=-\infty}^{\infty}x(n+rN).

这句话非常重要:频域抽样,会导致时域周期延拓并相加。

如果原序列长度为 MM,只有当

NMN\geq M

时,周期延拓的各段不会在主值区间发生重叠,因此可以无失真恢复原序列。若 N<MN<M,延拓后的各段会相互叠加,产生时域混叠。

频域采样点数与时域混叠

频域采样定理可以总结为:长度为 MM 的有限长序列,频域采样不失真的条件是频域采样点数 NMN\geq M。此时,长度不超过 NN 的有限长序列可以由它的 Z 变换在单位圆上 NN 个均分点处的采样值精确表示。

频域内插公式

NMN\geq M 时,可以用 NN 个频率采样值 X(k)X(k) 表示整个 X(z)X(z)。内插公式为

X(z)=k=0N1X(k)Φk(z),X(z)=\sum_{k=0}^{N-1}X(k)\Phi_k(z),

其中

Φk(z)=1N1zN1WNkz1.\Phi_k(z)=\frac{1}{N}\frac{1-z^{-N}}{1-W_N^{-k}z^{-1}}.

在单位圆上令 z=ejωz=e^{j\omega},得到

X(ejω)=k=0N1X(k)Φk(ejω).X(e^{j\omega})=\sum_{k=0}^{N-1}X(k)\Phi_k(e^{j\omega}).

Φk(ejω)\Phi_k(e^{j\omega}) 可以理解为频域内插函数,它把离散采样点之间的频谱连续地补出来。换句话说,只要频域采样点足够多,就可以从这些采样点恢复出原来的连续频谱形状。

本讲总结

FIR 滤波器设计的主线可以概括为四句话。

第一,FIR 的最大优势是可以实现严格线性相位,并且由于无反馈结构,稳定性好。

第二,线性相位 FIR 的关键在于冲激响应的对称性:

h(n)=±h(N1n).h(n)=\pm h(N-1-n).

第三,窗函数法从时域出发,用有限长窗函数截取理想无限长冲激响应;窗函数类型主要由阻带衰减决定,窗函数长度主要由过渡带宽度决定。

第四,频率采样法从频域出发,对理想频率响应等间隔采样,再通过 IDFT 得到 FIR 系数;它与频域采样理论密切相关,频域采样点数不足会导致时域混叠。

对于 0 基础学习者,建议先掌握“线性相位等于不改变波形相对时间结构”这一核心观念,再理解“偶对称或奇对称的 h(n)h(n) 会带来线性相位”。随后再学习窗函数法和频率采样法,就会自然很多。

Back to archive

Discussion

Comments

Post

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