← 图形原理 · 交互课 课题九 · 采样与重建 · 第 2 课

滤波器:一条频响曲线管什么

全课题骨架:采样不是"取点",是"用一串离散值 + 一个重建核,去猜那个连续的东西"; 核决定它像不像;而"先滤再采"是唯一的治本手段。
上一课反复说"先滤再采",可滤波器到底是什么?它凭什么"知道"该把 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 论文(预畸变的局限)。

分段:
反馈系数 a₁: ← 17.1 的旋钮
探针频率: 1000 Hz ← 17.2 的旋钮(在频响曲线上滑一条竖线,同时报"时域测"与"闭式算"两个数)
滤波器档: ← 17.3 的旋钮
设计频率 f₀: 1600 Hz(0.200·f_s) ← 17.4 的旋钮(画面同时画"预畸变 / 不预畸变"两条曲线)
(四张表与十条判据在后台跑)

(脚本启动中…)

定点观察:(脚本启动中…)

本帧的读数

(脚本启动中…)

(点「跑一遍读数」或「扫一遍设计频率」就在这里出数)

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,冲激响应越长。

输入 这个阶段在执行的规则(照抄件二真源码) 输出 单位脉冲 x[n] x[0] = 1,其余 0 旋钮:a₁ 0 / 0.5 / 0.9 / 0.99 b₀ = 1(未归一化) 或 b₀ = 1 − a₁ ★ 本段只加"反馈项" // 一阶差分方程(反馈只有一项) let z = 0; for (let n = 0; n < N; n++) { z = b0 * x[n] + a1 * z; // ← a1 就是反馈 y[n] = z; } // 闭式:h[n] = a1^n,Σh = 1/(1−a1) const dc = 1 / (1 - a1); // 归一化后截断到 N 项:Σ = 1 − a1^N const cut = 1 - Math.pow(a1, N); a₁ = 0: Σh = 1.000000 h[200] = 0(FIR) a₁ = 0.9: Σh = 10.000000 h[200] = 7.06e−10 a₁ = 0.99: 前 300 项和 = 0.950959 (手算 1 − 0.99³⁰⁰) ★ 越慢的极点在 1 附近
定点观察:把「反馈系数 a₁」从 0 切到 0.99 → 看画面里那串竖条(冲激响应) 与读数行的 Σh → a₁ = 0 时一格之后恰好归零(FIR);a₁ ≠ 0 时永远不归零(IIR), 而 Σh 一路从 1 涨到 10(a₁ = 0.9)、再涨到 100(a₁ = 0.99 的极限) —— 一个反馈系数就把"平均一下"变成了"放大十倍"。
表① 一阶 IIR 的四档(b₀ = 1;N = 300 个样本)。
★ 第 2 列把冲激响应与闭式 a₁ⁿ 逐点比较 —— 它们就是同一个东西(残下的 1e−16 是浮点舍入)。
★ 第 3~5 列是同一个量的三种说法:500 项的实测和、闭式 1/(1−a₁)、 归一化(b₀ = 1−a₁)后截断到 N 项的精确闭式 1 − a₁ᴺ。 a₁ = 0.99 那一行两者对不上,是截断而不是算错:a₁³⁰⁰ = 0.0490 还不小。
★ 最后一列是"它有没有尽头":a₁ = 0(无反馈)时 h[200] 恰好 0,其余永远不为 0。
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.1 的那一台) 旋钮:探针频率 20 … 3990 Hz 两条路线: 时域量 ≈ 频域代 ★ 本段只加"频响" // 时域:喂 cos,等暂态过去,量稳态幅度 const x = new Float64Array(SKIP + WIN); for (let n = 0; n < x.length; n++) x[n] = Math.cos(2*Math.PI*f*n/fs); const tm = proj(biquad(x, b, a), fs, f, SKIP, WIN); // 频域:把 e^{jω} 代进闭式(同一台机器的另一个问法) const w = 2 * Math.PI * f / fs; const cl = |b0 + b1·e^{−jw} + b2·e^{−2jw}| / |1 + a1·e^{−jw} + a2·e^{−2jw}|; // 判据:|tm − cl| 应 ≈ 0 低通 @1000 Hz: f=250 → 0.998405 f=1000 → 0.707107 f=3000 → 0.029425 两条路线最大差: 2e−14(三类型六点) 退化格(DC/Nyquist): 读数恰好是 2 倍 ★ 口径要单独说
定点观察:拖「探针频率」从 20 到 3990 → 看频响曲线上那条竖线 与读数行里的两个数 → "时域测出来的幅度"与"闭式代出来的幅度"始终相等(最大差 2e−14), 而它们扫过的就是那条频响曲线:低通在 250 Hz 是 0.998、在 1000 Hz 掉到 0.707、在 3000 Hz 只剩 0.029。
表② 两条独立路线(三种滤波器,f₀ = 1000 Hz、Q = 1/√2、f_s = 8000 Hz)。
★ 「时域测」= 喂 cos 进去、跳过 4800 个样本的暂态、在 2400 个整周期上做单频投影; 「频域代」= 把 e^{jω} 代进 `|Σb e^{−jkω}| / |1+Σa e^{−jkω}|`。 两者最大差 2e−14 —— 同一条算式,两条互不相干的算法。
★ 「时域测/频域代」两列并排写出来,正是为了让读者自己核对;探针频率那一档由滑块给(表里默认报 1000 Hz 那一列)。
滤波器频率时域测频域代差

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₀、Q、G(dB) 旋钮:四档 低通/高通/带通/峰值 Q = 1/√2 (Butterworth 那档) ★ 本段只加"设计式" // 两个中间量(照抄 RBJ / Web Audio) const w0 = 2 * Math.PI * f0 / fs; const alp = Math.sin(w0) / (2 * Q); // 低通(含 BLT 的预畸变:用 cos w0 / sin w0 而不是 w0) b0 = (1 - cos w0)/2; b1 = 1 - cos w0; b2 = b0; a0 = 1 + alp; a1 = -2 cos w0; a2 = 1 - alp; // 带通:b0 = alp, b1 = 0, b2 = -alp(0 dB 峰值增益版) // 峰值:b0 = 1+alp*A, b2 = 1-alp*A, a2 = 1-alp/A const A = Math.pow(10, dBgain / 40); // 特值:低通/高通 @f0 = Q、带通 @f0 = 1、峰值 @f0 = A² f₀ 处的实测: 低通 0.707107 高通 0.707107 带通 1.000000 峰值 2.000000 4 阶的正确 Q: 0.5412 与 1.3066 (积 = 0.707107) 两个相同节 vs 真 4 阶 差 14.60%
定点观察:把「滤波器档」在四档之间切 → 看读数行里"f₀ 处的 |H|"与表③ 的手算列 → 低通与高通都恰好 0.707107(= Q)、带通恰好 1.000000、峰值恰好 2.000000(= A² = 10^(6.02/20)) —— 三个都能拿笔算出来,而且都不随 f₀ 变。
表③ 四档(外加一档负增益的峰值)在 f₀ 处的读数(f₀ = 1000 Hz、Q = 1/√2、G = ±6.02 dB、f_s = 8000 Hz)。
★ 「手算」那一列就是三条特值恒等式:低通/高通 = Q、带通 = 1、峰值 = A² = 10^(G/20)。 第三列「Web Audio Q_dB」= 20·log₁₀(Q),也就是"同一个滤波器在 Web Audio 里该写多少"(低通/高通的 Q 以 dB 计)。
★ 最后一行是 4 阶 Butterworth 的两级 Q(0.5412 与 1.3066)—— 它们的积恰好是 0.707107 = 1/√2, 所以级联后在截止点恰好 −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 的部分仍然被压扁。

输入 这个阶段在执行的规则(照抄件二真源码) 输出 模拟原型(∞ 长) H(s) = 1/(s²+s/Q+1) 旋钮:设计频率 f₀ 0.05 … 0.45 · f_s 两档一起画 预畸变 / 不预畸变 ★ 本段只加"预畸变" // 双线性变换 + 预畸变常数(CCRMA: c = cot(w0 T/2)) const K = Math.tan(Math.PI * f0 / fs); // = tan(w0 / 2) const cw = (1 - K*K) / (1 + K*K); // cos w0 const sw = 2 * K / (1 + K*K); // sin w0 // 「不预畸变」= 把 tan(w0/2) 当成 w0/2(小角近似) const Kn = Math.PI * f0 / fs; const cwn = (1 - Kn*Kn) / (1 + Kn*Kn); // 手算:实际 −3dB 点 = (fs/π)·atan(π f0/fs) const hand = (fs / Math.PI) * Math.atan(Math.PI * f0 / fs); // 预畸变后偏差应为 0(实测 0.0000%) 预畸变(−3dB 点): 五档偏差都是 0.0000% 不预畸变: 0.05·fs → −0.81% 0.20·fs → −10.72% 0.45·fs → −32.44% 上一课的 21.6%: 形状 14.6% × 扭曲 6.1% = 21.6% ★ 只能对齐一个频率
定点观察:拖「设计频率 f₀」从 200 到 3600 Hz → 看画面里那两条频响曲线与 −3 dB 水平线的两个交点 与读数行 → f₀ 低的时候两条曲线几乎重合(0.05·f_s 时截止点只差 0.81%),越往 Nyquist 拖分得越开 (0.45·f_s 时差了 32.44%);而"预畸变"那条的 −3 dB 交点始终钉在 f₀ 上。 ★ 这就是上一课那 21.6% 的真身。
表④ 五档设计频率(低通、Q = 1/√2、f_s = 8000 Hz)。「−3 dB 点」由二分法找出 (|H| 降到 1/√2 的频率)。
★ 「预畸变」那一列的偏差全是 0.0000% —— 这就是预畸变在做的事(把一个频率精确钉住)。
★ 「不预畸变」那两列:实测与手算 (f_s/π)·atan(π f₀/f_s) 逐行吻合,偏差随 f₀ 单调恶化。
★ 「Nyquist 距离」= 1 − f₀/(f_s/2):越贴近 Nyquist,扭曲越厉害(DAFx 原话:预畸变也不能完全校正 Nyquist 附近)。
设计 f₀f₀ / f_sNyquist 距离预畸变 −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),而它的音高与衰减也能拿笔算。

本页用到的做法,出处在哪

运行时自检 等待脚本…