通常一个 FIR 的时域响应序列是纯实数,也就是 purely real 的 h(n),那这意味着其频域需要,那么必然频域的响应是中心对称的
H(ω)=H∗(−ω)
这是因为从离散时间傅里叶变换(DTFT)的公式出发
H(ω)=n=∞∑∞h(n)e−jωn
这里已知 h(n) 是纯实数,所以 h(n)=h∗(n),再用到 (AB)∗=(A∗)(B∗),(eA)∗=eA∗ 这两个小公式,就可以推导到这样的结果
H(−ω)H∗(−ω)=n=∞∑∞h(n)ejωn=n=∞∑∞h∗(n)(ejωn)∗=n=∞∑∞h(n)e−jωn=H(ω)
接着,我们还是观察这个 IDTFT 的表达式,不考虑平移的问题
h(n)=2π1∫−ππH(ω)cos(ωn)dω+j⋅2π1∫−ππH(ω)sin(ωn)dω
如果 H(ω)=R(ω)+jI(ω)=R(ω) 是 Purely Real,第二项必须等于 0,才能保证纯实数的 h(n),也就是说 H(ω) 的实数部分是偶对称,这样 偶对称 × 奇对称 = 奇对称,第二项积分结果为零
j⋅2π1∫−ππR(ω)sin(ωn)dω=0 ⟶ R(−ω)=R(ω)
此时我们代入 h(−n), 发现 h(−n)=h(n),意味着 Purely Realy 的 H(ω) 在频域是偶对称的,时域也是偶对称的
同理,如果 H(ω)=R(ω)+jI(ω)=jI(ω) 是 Puly Imaginary,第一项必须等于 0,才能保证纯实数的 h(n),也就是说 H(ω) 的实数部分是偶对称,这样 奇对称 × 偶对称 = 奇对称,第一项积分结果为零
2π1∫−ππjI(ω)cos(ωn)dω=0 ⟶ I(−ω)=−I(ω)
接着我们代入 h(−n) ,发现 h(−n)=−h(n),意味着 Purely Imaginary 的 H(ω) 在频域是奇对称的,时域也是奇对称的
首先,解释相位线性的定义:在频域里,任何系统的频率响应 H(ω) 都是一个复数
H(ω)=∣H(ω)∣ejθ(ω)
其中,∣H(ω)∣ 是幅度响应(决定声音变大还是变小),θ(ω) 是相位响应(决定不同频率的正弦波要延迟多久)。如果一个系统的相位响应 θ(ω) 和频率 ω 满足以下公式
θ(ω)=β−αω
其中 α 是一个常数,我们就称这个系统具有严格线性相位。直观理解 相位 = 角频率 × 时间,其实就是 路程 = 速度 × 时间。
如果 θ(ω)=β−αω,这意味着系统对所有频率的信号分量,在时间上的相对延迟都是一样的,全都是:
τg=−dωdθ(ω)=α (个采样点)
需要注意的是,这里的 ω 是数学上归一化的 ω,其含义是 “弧度/采样点”,如果要代入物理世界的采样频率,那么 ω=2πf/fs
这里的 β,是一个固定的相位差,可以理解为对信号的一种变形,比如将 cos 信号变成了 sin 信号
之前的讨论中,我们发现 h(n) 在 n<0 有值,这是一个物理不可实现的,着意对于一个冲击函数,n=0 才有值的 δ(n),响应 h(−3) 已经出现了,这是本来应该在 h(3) 时刻才出现的。物理系统没有“未卜先知”的能力,这就是所谓的物理不可实现系统。在学术上,我们称这种必须依赖未来输入才能计算当前输出的系统为非因果系统(Non-causal System)
真实世界中,我们需要给这样的系统一个延迟,才能物理可实现。现在我们来观察这个延迟会对响应造成什么样的影响。离散时间傅里叶反变换 (IDTFT) 的公式如下
h(n)=2π1∫−ππH(jω)ejωndω
根据傅里叶变换时移性质,时域上延迟 t0/Ts=α 时间,那么频域上就需要乘以 e−jωα,也就是说
hnew(n)=h(n−α)=2π1∫−ππH(jω)ejω(n−α)dω=2π1∫−ππ(H(jω)e−jωα)ejωndω
此时,我们引入欧拉公式 ejωn=cos(ωn)+jsin(ωn),再代入θ(ω)=−αω,求其时域的的冲击响应,可以看到由于我们加了一个延迟后,IDTFT 变得很麻烦。
h(n)=2π1∫−ππH(ω)cos(ωn)dω+j⋅2π1∫−ππH(ω)sin(ωn)dω
h(n)=2π1∫−ππ∣H(jω)∣cos(ω(n−α))+j⋅2π1∫−ππ∣H(jω)∣sin(ω(n−α))
那么故事就是这样的:首先,由于关于 n=0 对称的 h(n) 物理不可实现,对时域进行了延迟,整体向右平移了 α 个相位把,这样 hnew(n)=h(n−α),虽然此时频域发生了改变,多了一个 e−jωα,但是这个并不影响频域的 manitude,只是影响了相位,而且这个相位对所有频率都是线性的(时域上的固定延迟),所以对整体信号不会导致改变,所以设计滤波器可以从更方便的非因果开始。
这里额外补充,虽然线性相位的广义定义是 θ(ω)=β−αω,多了一个固定的 β 相位,对于信号是有涌现的。不同频率的分量,由于都需要走一个 β 的相位,但是由于频率不一样,其时域上的相对位移时间是不一致的。这个复合信号的内部成分发生了相对位移,波形的“长相”被彻底改变了。原来可能是个方波,现在可能变成了一堆乱七八糟的尖刺。这在工程上叫做相位失真(Phase Distortion)或色散(Dispersion)。
更神奇的是,早在 1843 年,物理学家格奥尔格·欧姆(没错,就是提出欧姆定律、研究电阻的那位大神)提出了欧姆声学定律(Ohm's acoustic law)。他指出,人类的听觉系统(耳蜗内的基底膜和毛细胞)本质上是在对声音做频谱分析。我们的耳朵极其敏感于不同频率的幅度(能量大小,∣H(ω)∣),这构成了我们对“音色”的认知。但是,我们的耳朵对不同频率分量之间的相对相位(θ)极其迟钝,甚至可以说是“相位耳聋”的。但是像底鼓(Kick)或军鼓的“啪”的一声。这些声音的能量高度集中在短短的几毫秒内(所有频率瞬间对齐)。如果 β 把高频和低频扯开了,这个干脆的“啪”声就会变成软绵绵的“噗”声,失去了冲击力(Smearing)
时遇上,根据 h(n) 的长度的奇偶性以及对称的奇偶性,可以将 FIR 分为以下四大类型,主要就是观察其在 DC 和 fs/2 处的频率特性。这个可以通过一个非常简单的示例给出
- Type1 - 奇数长度偶对称:低通/高通/带通/带阻全能
- Type2 - 偶数长度偶对称:fs/2 必然是 0 ,偶对称 FIR = {+1, +1},AC SIG = {+1 -1} 卷积后就是 0,DC SIG = {+1, +1} 卷积后结果是 2 总是有值的;
- Type3 - 奇数长度奇对称:DC 和 fs/2 必然是 0,奇对称决定中点值必然是 0,奇对称 FIR = {+1, 0, -1},AC SIG = {+1,-1...} 代入后都是0
- Type4 - 偶数长度奇对称:DC 必然是 0
更一般的,认为一个输入信号 x(n) 滤波器的 h(n) 长度是 M,从 0 到 (M−1) ,则响应 y(n) 的表达式为
y(n)=k=0∑M−1h(k)x(n−k)
对称无论是奇对称或者偶对称,0 和 M−1 一定是满足 h(0)=±h(M−1), h(1)=±h(M−2) ,因此有
h(n)=±h(M−1−n)
Group Delay 定义的其实是波形额能量重心的平移,Group Delay = (M−1)/2
在伯德图中,横坐标是频率 ω,纵坐标是 manitude;在 Z 平面中,绘制的不是 幅值vs频率,也不是 相位vs频率,绘制的是零点/极点 的位置。
reciprocal pair:首先我们证明零点 z0 和 z0−1 一定是成对出现的,一定都是零点,从 DTFT 出发
H(ω)=n=∞∑∞h(n)e−jωn
考虑有限长度的 Z 域表达式
H(z)=n=0∑M−1h(n)z−n
那么根据 h(n) 对称性,h(n)=h(M−1−n),代入就会有
H(z)=n=0∑M−1h(M−1−n)z−n
换元法,令 k=M−1−n,当 n=0 时,k=M−1;当 n=M−1 时,k=0,所以求和范围只是倒叙来一遍,还是 k=[0,M−1] 这个范围内,同时替换掉表达式中的 n,可以得到
H(z)=k=0∑M−1h(k)z−(M−1−k)=z−(M−1)k=0∑M−1h(k)(z−1)−k=z−(M−1)H(z−1)
假设 z0 时 H(z) 的零点,满足 H(z0) = 0,z−(M−1) 一定时一个非零的复数,那这样的话 H(z0−1) 必须是零了。所以只要时域 h(n) 是对称的,z0 是零点,那么 z0−1 也一定是零点。
conjugate pair : 接着证明, z0 和 z0∗ 一定是成对出现的,一定也都是零点,同理
H(z0)=n=0∑M−1h(n)z0−n=0
那么 0 的共轭复数复述仍然是 0,利用共轭的乘法性质,以及由于 h(n) 一定是实数的等式 h(n)=h∗(n),化简就可以得到如下结果
(H(z0))∗=0 ⟶ n=0∑M−1h∗(n)⋅(z0−n)∗=n=0∑M−1h(n)(z0∗)−n=0
这样,就得到了如下这张 四点完美对称星座图 (Quadruplet)

还可以看到下面这个表达式,
Hquad(z)=(1−z1z−1)(1−z1∗z−1)(1−z11z−1)(1−z1∗1z−1)=(1−(z1+z1∗)z−1+(z1z1∗)z−2)(1−(z11+z1∗1)z−1+(z1z1∗1)z−2)
从这里可以看到,共轭复述的零点,合并成了实数,消除了虚部。
任何一个 H(z) 都可以进行应式分解,得到
H(z)=1+a1z−1+a2z−2+⋯+aNz−Nb0+b1z−1+b2z−2+⋯+bMz−M
H(z)=C⋅zK(z−p1)(z−p2)…(z−pN)(z−z1)(z−z2)…(z−zM)
其中 C 是一个常数增益,zK 是原点处的极点或零点。现在,我们将观察点拉回到现实的频域,也就是 z=ejω,这里没有 zσ+jω 中的 σ=0,这样就会得到
H(ejω)=C⋅ejωK(ejω−p1)(ejω−p2)…(ejω−z1)(ejω−z2)…
我们要看幅度响应,就是对等式两边取绝对值,在复数运算中,乘积/商的绝对值,等于绝对值的乘积/商
∣∣∣∣∣CA⋅B∣∣∣∣∣=∣C∣∣A∣⋅∣B∣
而且 ∣ejωK∣=1,C 是一个常数认为是1,那么把绝对值分配进去,就可以得到
∣H(ejω)∣=∣ejω−p1∣⋅∣ejω−p2∣…∣ejω−z1∣⋅∣ejω−z2∣…
所以我们根据 Z 平面观察频率响应,就是想象一个点,绕着单位圆,从角度 0 出发,转一圈转到 π,剩下的 [π,2π] 没必要转,那是负频率,是为了抵消虚部的,是实信号的必然数学要求。复数的减法绝对值 ∣ejω−z1∣ 就是图纸上两点之间的距离,从一通过距离的相乘,就可以得到 幅度vs频率 响应了。
我们对零点的 Z 平面图,可以有如下理解:
- 对于一个 M 抽头的 FIR 滤波器,那么其零点一定有 M-1 个,因为 M 抽头对应的多项式阶数就是 M-1。这个准则,可以用来反推 FIR 的抽头数。对于奇数抽头的滤波器,其零点一定是偶数个
- 在实轴上的零点,在 DC 频率处或者 fs/2 这个频率,对应 ω=0 或者 ω=π,只能是成对出现,买1得2;
- 如果不在两端频率上的零点,一定是买1得4
- 对于偶数个抽头的滤波器,其零点一定是奇数个,那这个多余出来的奇数,必然只能出现在单位圆上,不在 ω=0, 就一定会在 ω=π,这个取决于其系数是偶对称还是奇对称
一个理想的低通滤波器,其频域是一个矩形窗,意味着他有无限的低的阻带,和无限垂直的过渡区间,以及非常平坦的通带。把这个完美的矩形做离散时间傅里叶反变换(IDTFT),我们在时域上会得到一个以 n=0 为中心的 sinc 函数。

它有两个致命缺点:
- 非因果: 在 n<0 的地方有无穷无尽的值。
- 无限长: 它的尾巴向两端无限延伸,计算机根本无法存储和计算无限多的滤波器系数(Tap)
为了让计算机能算出来,我们必须把这个无限长的序列“一刀切断”,只保留中间有限的 M 个点。然后再把这保留下来的部分向右平移 (M−)/2 个点,使其变成一个合法的、因果的系统。 这个“一刀切断”的动作,在数学上就等同于给原始信号乘以了一个“矩形窗(Rectangular Window),这就会导致一个实际的滤波是如下样子

这里的故事是这样的,频域的矩形,意味着时域的无限长度;为了解决这个无限长度,对这个时域的无限长度,进行截断,而阶段相当于在时域乘以了一个矩形窗(左图),时域的矩形意味着频域的无限 SINC 函数(右图),

矩形窗的瑕疵有如下 2 条:
- 过渡带变宽(由主瓣引起):原本垂直切断的频率边缘,现在变成了一个有坡度的滑梯。
- 频域起伏震荡(由旁瓣引起):这是最要命的!在截断边缘,通带和阻带会出现剧烈的波纹震荡。在阻带依然有巨大的能量残留),并且边缘有强烈的震荡 。无论你把滤波器长度 M 增加到多长,这个震荡的幅度都不会减小,只会变得更密集,这就是著名的吉布斯现象。如下图显示

既然矩形窗会带来各种各样的问题,那么是否可以通过其它窗来权衡,这就是 windows based design 的初衷。
采用该技术设计滤波器的步骤为:首先确定阻带最大衰减指标,该指标由最终实际应用场景决定。举例来说,若某系统的本底噪声为 - 50 分贝,汉明窗便是最合适的选择,因其阻带衰减可低于系统本底噪声。布莱克曼窗虽也能满足该项指标,但需要更多抽头系数,这会无谓地提升算法复杂度,并带来更大的处理时延。选定窗函数后,抽头数量M由过渡带的带宽确定。
| Window |
Max Stopband |
Transition Bandwith |
| Rectangular |
-21 dB |
1.8 π/M |
| Hanning |
-44 dB |
6.2 π/M |
| Hamming |
-53dB |
6.6 π/M |
| Blackman |
-74dB |
11 π/M |
如果我们要设计
- cut-off frequency = 0.1 π
- stopband attenuation = -40dB
- Transition bandwidth = 0.02 π
第一步,我们认为 Hanning 窗满足了阻带要求,选择 hanning 窗。
w(n)=0.5−0.5cos(M−12πn)
第二步,根据阻带的要求,设计滤波器的抽头数
M6.2π=0.02π → M=310
第三步,我们设计一个奇数抽头的滤波器,选择 M=311,这样它的 Group Delay τg=(M−1)/2=155;接着计算理想低通滤波器(矩形)的系数,我们先根据设计要求(比如截止频率 ωc)写出理想低通滤波器的冲激响应。在时域上,它是一个无限延伸的 sinc 函数:
H_{RD}(ω)=\cases{1;\quad|\omega|<0.1\pi \\0;\quad\text{otherwise} }
hd(n)=∫−ωcωce−jτgωejωndω=π(n−2M−1)sin(ωc(n−2M−1))=π(n−155)sin(0.1π(n−155))
第四步,我们要取其中 n∈[0,1,…,310] 一共 311 个系数,然后乘以Hanning 窗。它的公式是
w(n)=0.5−0.5cos(M−12πn)
得到的结果如下

其它窗函数的公式
矩形窗(Rectangular),主瓣最窄(过渡带最陡峭),但旁瓣最大(阻带衰减仅 −21 dB)。适用于对过渡带要求极高、完全不在乎噪音残留的糙汉场景
ω(n)=1
汉宁窗(Hanning),滑的升余弦窗。主瓣变宽一倍(过渡带变缓),但阻带衰减提升至 −44 dB
w(n)=0.5−0.5cos(M−12πn)
哈明窗(Hamming),微调了 Hanning 窗的系数,利用边缘微小的突变产生一个相反相位的旁瓣,与其第一旁瓣发生相消干涉。阻带起始衰减可达 −53 dB,是工程界最常用的万金油
w(n)=0.54−0.46cos(M−12πn)
布莱克曼窗(Blackman),引入了二次谐波进一步抹平边缘。过渡带最宽,但阻带衰减达到了恐怖的 −74 dB,适合追求极致静音的精密医疗/高保真音频滤波
w(n)=0.42−0.5cos(M−12πn)+0.08sin(M−14πn)