滤波器:一条频响曲线管什么
全课题骨架:采样不是"取点",是"用一串离散值 + 一个重建核,去猜那个连续的东西";
核决定它像不像;而"先滤再采"是唯一的治本手段。
上一课反复说"先滤再采",可滤波器到底是什么?它凭什么"知道"该把 700 Hz 压掉?
这一课把它拆开:一台滤波器就是一个差分方程 —— 五项系数就是那条频响曲线;
而"设计滤波器"就是从"我要一条什么形状"反解出这五项系数(RBJ 的公式)。
唯一的坑是:从连续域搬过来的时候,频率轴被 tan 扭了一下,所以要预畸变。
四段路:五个系数就是一台机器(17.1)→ 它到底把哪个频率放了、哪个压了(17.2)→ 反解系数(17.3)→ 唯一的坑:频率轴被扭了一下(17.4)。
★ 本页同样不发声(判据必须能在这里逐样本复算),并且不拿浏览器原生的 BiquadFilterNode 当判据
—— 但规范给了完整系数公式,它就是本页 17.3 的一手源(见页末)。
保真度:真机可证(可测) —— 滤波器全部逐样本算,四张表、十条判据都能离线复算。
一手源:Web Audio 规范(BiquadFilterNode 的系数公式)、RBJ Audio EQ Cookbook、CCRMA(双线性变换 / 频率扭曲 / 预畸变 / Butterworth)、DAFx 论文(预畸变的局限)。
本页的原理与四张数值表照常可读(画面只是把同一份状态画出来)。
(脚本启动中…)
定点观察:(脚本启动中…)
本帧的读数
(脚本启动中…)
(点「跑一遍读数」或「扫一遍设计频率」就在这里出数)
17.1 五个系数就是一台机器
本段只增加一样东西:反馈项(a 系数)(旋钮是「反馈系数 a₁」;一阶、只有五项)。
上一段的缺口:上一课反复说"先滤再采",可滤波器到底是什么?它凭什么"知道"该压掉 700 Hz?
答案朴素得有点意外:一台滤波器就是一个差分方程 —— 每个输出样本 = 现在与过去几个输入、加上过去几个输出的加权和:
y[n] = b₀x[n] + b₁x[n−1] + b₂x[n−2] − a₁y[n−1] − a₂y[n−2]。
五个系数(b₀b₁b₂a₁a₂)就是全部。本段先只看最简的一阶:y[n] = b₀x[n] + a₁y[n−1]。
★ 能拿笔核:喂一个单位脉冲(x[0] = 1,其余 0),输出就是冲激响应,而它恰好是
h[n] = a₁n(实测与它逐位一致,最大差 1.1e−16)。于是整台机器的"总账"是
Σh[n] = 1/(1−a₁) —— 也就是直流增益(CCRMA 的原话:一阶节的 `H(z) = 1/(1−pz⁻¹)`,逆变换是 `pⁿu(n)`;
Stanford EE152 讲义:一阶 IIR `yₙ = αy_{n−1} + βx` 的直流增益是 `β/(1−α)`)。
★ 表① 里有两条反直觉的读数:
① a₁ = 0.9、b₀ = 1 时,这台"看起来只是平均一下"的机器直流增益是 10 倍;
真实设计里写 b₀ = 1 − a₁ 正是为了把它拉回 1。
② a₁ = 0 与 a₁ ≠ 0 是两种东西:前者冲激响应有限长(FIR,a₁ = 0 时 h[200] 恰好 0)、
后者无限长(IIR,a₁ = 0.9 时 h[200] 还有 7.06e−10)—— 因为它们一个没有反馈、一个有。
★ 还有一条精确恒等式:归一化后截断到前 N 项的和恰好是 1 − a₁N
(a₁ = 0.99、N = 300 时实测 0.950959,手算 1 − 0.99³⁰⁰ = 0.950959)——
它顺便说明"Σh 要取到足够长才等于 1/(1−a₁)":极点越靠近 1,冲激响应越长。
| a₁ | 与闭式 a₁ⁿ 的最大差 | Σh(前 300 项) | 闭式 1/(1−a₁) | 归一化后 Σh / 手算 1 − a₁³⁰⁰ | h[200] |
|---|
17.2 它到底把哪个频率放了、哪个压了
本段只增加一样东西:频响(把复响在单位圆上求值)(旋钮是「探针频率」,类型固定为低通)。
上一段的缺口:17.1 只看时间侧(冲激响应),可说清"它把 700 Hz 压掉多少"还没着落。
做法是把同一个差分方程换个问法:喂一个 cos(ωn) 进去,稳态输出的幅度就是"这个频率被放大了多少倍"。
而这个倍数有闭式(KAIST 数字滤波器讲义,以及一切教科书都这么写):
|H(ejω)| = |b₀ + b₁e−jω + b₂e−2jω| / |1 + a₁e−jω + a₂e−2jω|。
★ 于是就有了本课最硬的一条:两条完全独立的路线 ——
① 时域:真的把 cos 喂进去、等暂态过去(跳过 4800 个样本)、在 2400 个整周期上量幅度;
② 频域:把 e^{jω} 代进那条闭式。
实测(表②):三种滤波器 × 六个频率,最大差 ≤ 2e−14 —— 它们算的是同一个数。
★ 一阶的闭式也能拿笔核:|H(ω)| = 1/√(1 + a₁² − 2a₁cos ω)。
在 ω = 0 处它给出 1/(1−a₁) = 10(a₁ = 0.9)—— 与 17.1 的 Σh[n] 是同一条账的两条路。
★ 有一条口径要单独说清(与上一课的判定同源):直流(f = 0)与 Nyquist(f = f_s/2)是两个退化格 ——
正负频率折进同一格,时域幅度读数恰好是闭式的 2 倍(实测 20.000000 vs 10.000000)。
所以本页的"两条路线一致"这句话,只在 0 < f < f_s/2 上成立(判据④ 专门把那个 2 倍断言下来)。
| 滤波器 | 频率 | 时域测 | 频域代 | 差 |
|---|
17.3 从"我要一条什么形状"反解出系数
本段只增加一样东西:设计式(RBJ 的系数公式)(旋钮是「滤波器档」四档;探针频率不动)。
上一段的缺口:17.2 会算了,可"我想要一条低通"该给什么系数?没人愿意从五个系数正着猜。
RBJ 的 Audio EQ Cookbook(Web Audio 规范给了完整公式,并注明 "inspired by formulas found in the Audio EQ Cookbook")
把这件事反过来了:你给 f₀(工作频率)、Q(形状)、G(增益 dB),
它直接算出五个系数。中间量只有两个:ω₀ = 2πf₀/f_s、α = sin ω₀/(2Q)。
低通是 b₀ = (1−cos ω₀)/2、b₁ = 1−cos ω₀、b₂ = (1−cos ω₀)/2、a₀ = 1+α、a₁ = −2cos ω₀、a₂ = 1−α;
高通把 cos 号翻一下;带通用 b₀ = α、b₁ = 0、b₂ = −α(就是规范里那个 "constant 0 dB peak gain" 版);
峰值用 b₀ = 1+αA、b₂ = 1−αA、a₀ = 1+α/A、a₂ = 1−α/A,其中 A = 10^(G/40)。
★ 于是有三条可以拿笔核的特值恒等式(表③ 逐档实测):
① 低通与高通在 f₀ 处恰好等于 Q(Q = 1/√2 ⇒ 0.707107,实测 0.707107);
② 带通在 f₀ 处恰好是 1.000000(这就是"0 dB 峰值增益"这个名字的来处);
③ 峰值 EQ 在 f₀ 处恰好是 A² = 10^(G/20)(G = +6.02 dB ⇒ 2.000000,实测 2.000000,差 0)。
★ 顺便说清一件很容易踩的事:同一个"低通",Q 的口径有两种 ——
RBJ cookbook 用线性 Q(f₀ 处增益 = Q),而 Web Audio 规范把低通/高通的 Q 当 dB 解释
(α = sin ω₀/(2·10^(Q/20)),原话:"For lowpass and highpass filters the Q value is interpreted to be in dB")。
所以 Web Audio 里写 Q = −3.0103 得到的正是 Butterworth(f₀ 处 0.707107)——
这条换算(Q_dB = 20·log₁₀ Q)实测逐档吻合。
★ 还有一件上一课留下的尾巴要在这一段收掉:4 阶 Butterworth 不是"把 2 阶的节级联两次" ——
它要两个 Q 不同的节:0.5412 与 1.3066(6 阶是 0.5176 / 0.7071 / 1.9319)。
实测:同一个 4 阶 fc=450 @ 8000 Hz,"两个相同节"的 |H(700)| 是 0.138482,而"真 4 阶"是 0.158705,差 14.60%
—— 上一课表③ 那条 21.6% 里,主要就是这 14.6 个点(另外约 6 个点才是双线性变换的扭曲,17.4 会算清)。
★ 附送一条好看的恒等式:偶数阶 Butterworth 的各级 Q 之积恰好 = 1/√2
(0.5412 × 1.3066 = 0.707107)—— 所以级联后在截止点恰好 −3 dB。
| 档 | f₀ 处实测 | f₀ 处闭式 | 手算 | 差 | Web Audio Q_dB |
|---|
17.4 唯一的坑:双线性变换把频率轴扭了一下
本段只增加一样东西:预畸变(旋钮是「设计频率 f₀」;画面同时画"预畸变"与"不预畸变"两条曲线)。
上一段的缺口:17.3 的系数是"从 f₀ 算来的",可数字滤波器的 −3 dB 点到底落在哪个频率上?
要回答它得先说清系数是怎么来的:RBJ 的公式来自双线性变换(CCRMA 原话)
s = c·(1−z⁻¹)/(1+z⁻¹),它把整条模拟频率轴压进单位圆:
s = 0 映到 z = 1,s = ∞ 映到 z = −1,而频率轴的对应关系是
jω_a = c·tan(ω_dT/2) —— 也就是被 tan 扭了一下。
CCRMA 把它说得更直白:"Some kind of frequency warping is obviously unavoidable in any one-to-one map
because the analog frequency axis is infinite while the digital frequency axis is finite."
★ 于是预畸变就是把这个扭曲在"一个点"上扭回来:想让数字截止恰好落在 ω₀,就取
c = 1/tan(ω₀/2)(CCRMA:"we set c = cot(ω_cT/2) to obtain a digital cut-off frequency at ω_c")。
RBJ 的公式里用 cos ω₀、sin ω₀ 而不是 ω₀ 的小角形式,这就是预畸变。
★ 能拿笔核:如果不预畸变(把 tan(ω₀/2) 当成 ω₀/2),实际 −3 dB 点会跑到
f_actual = (f_s/π)·atan(π f₀/f_s)。表④ 五档实测逐行命中:
f₀ = 0.05·f_s 时只偏 −0.81%,0.1·f_s 时偏 −3.11%,0.2 时 −10.72%,0.3 时 −19.81%,
0.45 时 −32.44% —— 越靠近 Nyquist 越糟,而预畸变之后五档的偏差都是 0.0000%。
★ 这一段顺手收掉上一课留下的那条尾巴:上一课表③ 里"两个口径差 21.6%",拆开就是两笔账 ——
形状差 14.60%(两个相同 Q 的节 ≠ 真 4 阶 Butterworth,17.3 已量)+ 扭曲约 6.1%(本条);
两个因子相乘 1.1460 × 1.0608 = 1.2156,正好是那个 21.6%。
★ 最后一条要如实标注(DAFx 2021 那篇原话):预畸变"allows for some correction of the effects of
frequency warping by directly mapping one continuous frequency to one discrete frequency,
but it will not completely correct for frequency warping near the Nyquist limit" ——
它只能精确对齐一个频率,靠近 Nyquist 的部分仍然被压扁。
| 设计 f₀ | f₀ / f_s | Nyquist 距离 | 预畸变 −3 dB 点 | 偏差 | 不预畸变 −3 dB 点 | 偏差 | 手算 (f_s/π)·atan(π f₀/f_s) |
|---|
这一课带走了什么
三行:① 一台滤波器就是一个差分方程,五项系数是全部;② 那条频响曲线有闭式,
时域测与频域代算的是同一个数(实测差 2e−14);③ "设计"就是从形状反解系数 ——
三条特值恒等式把这件事变得能拿笔核(低通/高通 @f₀ = Q、带通 @f₀ = 1、峰值 @f₀ = A²)。
而唯一的坑是频率轴被 tan 扭了一下:不预畸变时 −3 dB 点会整体下移(0.45·f_s 处偏 −32.44%),
预畸变能把它精确钉在一个频率上,但靠近 Nyquist 的部分仍然被压扁(这是 IIR 的固有代价,一手源原话如此)。
下一课指针:到这一课为止,"滤波器"都是现成的解析设计(给 Q 给 f₀ 出系数)。
下一课换一种做法:把一根弦的物理写下来,让滤波器自己长出来 —— 一个延迟线加一个两点的平均,
几十行就能拨响一根弦(Karplus-Strong),而它的音高与衰减也能拿笔算。
本页用到的做法,出处在哪
- 五种滤波器类型的完整系数公式 —— W3C Web Audio API 规范 §Filters Characteristics: 给出 lowpass/highpass/bandpass/notch/allpass/peaking 的 b₀…a₂ 全部表达式(含 `ω0 = 2πf0/Fs`、`α_Q = sin ω0/(2Q)`、`α_QdB = sin ω0/(2·10^(Q/20))`、`A = 10^(G/40)`), 并声明 implementations "MUST implement",同时注明 "They are inspired by formulas found in the Audio EQ Cookbook";传递函数写作 `H(z) = (b0/a0 + b1/a0 z⁻¹ + b2/a0 z⁻²)/(1 + a1/a0 z⁻¹ + a2/a0 z⁻²)`: w3.org/TR/webaudio/
- "低通/高通的 Q 以 dB 计" —— 同规范
Q属性原文: "The Q factor of the filter. For lowpass and highpass filters the Q value is interpreted to be in dB. … For the bandpass, notch, allpass, and peaking filters, this value is a linear value. … This is not used for the lowshelf and highshelf filters." ★ 本页表③ 的 Q_dB 那一列就是这条的换算。 - RBJ Audio EQ Cookbook —— Robert Bristow-Johnson,经授权由 W3C 转载 ("Adapted from Audio-EQ-Cookbook.txt, by Robert Bristow-Johnson, with permission"): 本页 17.3 的低通/高通/带通/峰值公式与 `α = sin ω0/(2Q)` 均出自这里; 带通的两个版本原文分别标为 "peak gain = Q" 与 "constant 0 dB peak gain"(本页用的是后者): webaudio.github.io/Audio-EQ-Cookbook (W3C Note 版 w3.org/TR/2021/NOTE-audio-eq-cookbook)
- 双线性变换、频率扭曲、预畸变 —— Julius O. Smith III(CCRMA / Stanford): `s = c·(1−z⁻¹)/(1+z⁻¹)`; "analog dc (s=0) maps to digital dc (z=1) and the highest analog frequency (s=∞) maps to the highest digital frequency (z=−1)";`jω_a = c·tan(ω_dT/2)`; "the amplitude response takes on exactly the same values over both axes, with the only defect being a frequency warping…"; "Some kind of frequency warping is obviously unavoidable in any one-to-one map because the analog frequency axis is infinite while the digital frequency axis is finite."; 预畸变 "we set c = cot(ω_cT/2) to obtain a digital cut-off frequency at ω_c": ccrma.stanford.edu/~jos/filters/Bilinear_Transformation.html、 …/Frequency_Warping.html、 …/Example_Second_Order_Butterworth_Lowpass.html
- 预畸变也治不了 Nyquist 附近 —— Darabundit & Abel(CCRMA), DAFx20in21: "a byproduct of mapping infinite continuous frequencies to a finite range of digital frequencies is frequency warping. This distorts the original response and is particularly severe near half the sampling rate: the Nyquist limit. A technique known as "prewarping" allows for some correction … but it will not completely correct for frequency warping near the Nyquist limit." dafx.de/paper-archive/2021/…
- 一阶 IIR 的冲激响应与直流增益 —— CCRMA:`H(z) = 1/(1−pz⁻¹) = 1 + pz⁻¹ + p²z⁻² + …`, 逆变换 "H(z) ↔ h(n) = u(n)·pⁿ"; Stanford EE152 实验讲义:一阶 IIR `yₙ = αy_{n−1} + βx` 的 "DC gain is y/x = β/(1−α)": …/One_Pole_Transfer_Functions.html、 web.stanford.edu/class/ee152/labs/lab2.pdf
- 二阶节频响的闭式 —— KAIST GCT535 讲义:`H(ω) = (b₀ + b₁e^{−jω} + b₂e^{−2jω}) / (1 + a₁e^{−jω} + a₂e^{−2jω})`: mac.kaist.ac.kr/~juhan/gct535/…
- Butterworth:最大平坦、截止处 −3 dB —— CCRMA("maximally flat amplitude response"、 原型 `1/(s²+√2s+1)`、"the amplitude response is G_a(j)=1/√2"、`20log₁₀|H|= −3 dB`); TI SLOA049D:"Butterworth filters are termed maximally-flat-magnitude-response filters… The attenuation is –3 dB at the cutoff frequency." …/Butterworth_Lowpass_Design.html、 ti.com/lit/an/sloa049d
- ★ 四处如实标注:① "低通/高通在 f₀ 处等于 Q"是推导结论,不是 cookbook 的原话 (cookbook 只给归一化原型 `H(s) = 1/(s²+s/Q+1)`;带通那两个 "peak gain = Q" / "constant 0 dB peak gain" 标签才是原文);② "Q = 1/√2 ≈ 0.7071" 的逐字表述 我们只在课程材料里找到,CCRMA/TI 的一手页只给"最大平坦"与"−3 dB 截止"(CCRMA 用 `s²+√2s+1` 隐含它); ③ Butterworth 各级 Q 的分配(0.5412 / 1.3066 / …)在课本上是标准结果,但我们没有抓到可引用的 原页,本页只用"各级 Q 之积 = 1/√2"这条自己实测过的性质来说明;④ 一阶 IIR 的直流增益 `β/(1−α)` 那一句出自 Stanford 课程讲义(一手性弱于 CCRMA 的书),已如实标出来源。