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

采样:为什么"多采几个点"不是答案

全课题骨架:采样不是"取点",是"用一串离散值 + 一个重建核,去猜那个连续的东西"。 核决定它像不像;而"先滤再采"是唯一的治本手段 —— 图像里的摩尔纹与声音里的混叠,是同一道数学。
上一课的输入是"已经离好的东西"(顶点、像素、瓦片、样例图)。可现实里的信号是连续的 —— 把它变成一串数的那一步,偷走了什么?答案是 频率:欠采样时,一串 700 Hz 的样本与一串 300 Hz 的样本 逐位相同(不是"像",是同一个序列),这就是混叠。
四段路:采样偷偷改掉了频率(16.1)→ 把点连回去,核决定一切(16.2)→ 先滤再采(16.3)→ 图像里的同一道数学(16.4)。
★ 两条如实说明:① 本页不发声 —— 判据必须能在这里逐样本复算,而发声要过一次浏览器自己的重采样与设备 (W3C 规范原话:解码出来的 AudioBuffer 会被重采样到上下文的采样率),读到的数与听到的数会对不上; ② 浏览器原生的 OscillatorNode 可能自己带限(规范说"实际产生的波形可能为了防混叠而不同", 且"必须等价于用一个带相应傅里叶级数的 PeriodicWave"),但本页没有核验具体实现,所以不拿它当判据。 保真度:真机可证(可测) —— DSP 全部逐样本算,四张表、十条判据都能离线复算。 一手源见页末(Shannon 1949 / J.O. Smith / PBRT / Crow 1977 / Whitted 1980 / W3C)。

分段:
采样率: Hz ← 16.1 的旋钮(信号 700 Hz 固定)
重建核: ← 16.2 的旋钮
采前低通: ← 16.3 的旋钮(管线:8000 Hz 造信号 → 低通 → 降到目标采样率)
图像宽度: 像素 ← 16.4 的旋钮(条纹数恒为 W 的 11/16)
(四张表与十条判据在后台跑)

(脚本启动中…)

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

本帧的读数

(脚本启动中…)

(点「扫一遍采样率」或「跑一遍读数」就在这里出数)

16.1 采样就是"偷看":它偷偷改掉了频率

本段只增加一样东西:采样(在连续信号上按等距取点)(旋钮是「采样率」,信号固定 700 Hz)。
上一段的缺口:前八课的输入都是已经离好的东西 —— 顶点、像素、瓦片、样例图。可现实里的信号是连续的, 把它变成一串数的那一步偷走了什么?
取点这件事看起来无害:y[n] = x(n / f_s)。可它一旦取少了,频率就被改了 —— 而且改法是折叠(J. O. Smith 的原话:"the new aliased frequency is f_alias = |f − k·Fs|… This phenomenon is also called 'folding', since… 'mirror image'"): `f_alias = |f − k·f_s|`,`k = round(f / f_s)`,把原频率折进 `[0, f_s/2]`。
能拿笔核:700 Hz 用 1000 Hz 采 ⇒ `k = round(0.7) = 1` ⇒ `|700 − 1000| = 300 Hz`; 用 800 Hz 采 ⇒ `|700 − 800| = 100 Hz`;用 400 Hz 采 ⇒ `k = 2` ⇒ `|700 − 800| = 100 Hz`。
★ 而"混叠"这两个字最严格的形态不是"听起来像",是逐位相同: 采样出来的那一串数与"直接采一个 300 Hz 的音"一个比特都不差(实测最大差 5.0e−11, 纯粹是 `cos` 在 ±3.8 万弧度上的舍入)。表① 把这个比较逐档列出来。
★ 还有一档要如实说:`f = f_s/2`(Nyquist 上)是退化档 —— 正负频率折进同一格, 谱读数是 2 而不是 1,折叠闭式在这一档不唯一(表① 注脚)。
★ 顺带说清一件容易误读的事:对一个采样序列,"在 700 Hz 处的能量"这句话没有定义 —— 它只能表达 `[0, f_s/2]` 以内的频率。所以判据写成"两个序列逐位相同",不写成"原频率处为零"。

输入 这个阶段在执行的规则(照抄件二真源码) 输出 连续信号 x(t) f = 700 Hz(固定) 旋钮:采样率 f_s 4000 / … / 600 Hz 样本数 N = 15360 (240 × 64,各档都整除) ★ 本段只加"采样" // 取点:连续 → 离散 const y = new Float64Array(N); for (let n = 0; n < N; n++) y[n] = Math.cos(2 * Math.PI * 700 * n / fs); // 折叠:把原频率折进 [0, fs/2] const k = Math.round(700 / fs); const fAlias = Math.abs(700 - k * fs); // 判据:与"直接采 fAlias"逐位比较 for (let n = 0; n < N; n++) d = Math.max(d, Math.abs(y[n] - Math.cos(2*Math.PI*fAlias*n/fs))); fs = 1000: f_alias = 300 Hz 逐位最大差 2.2e−11 fs = 800: f_alias = 100 Hz fs = 4000: f_alias = 700 Hz (f_alias 就是 f) ★ 样本把频率改了
定点观察:把「采样率」从 4000 一路切到 600 → 看画面里那条从样本连出来的亮线 与读数行里的"折叠频率" → 采样率降到 1000 以下,700 Hz 的真值被采成了 300 / 100 Hz 的慢波; 而它与"直接采那个慢频率"的序列逐位相同(不是像,是同一个序列)。
表① 七档采样率(信号 700 Hz,N = 15360)。「逐位最大差」= 这串样本与"直接采 `f_alias`"的样本逐位比较的最大差: ★ 这不是相似,是同一个序列(残下的 1e−11 是 `cos` 在大幅角上的舍入)。 「投影@f_alias」= 在这一串样本上做单频投影(整周期采样,所以读数精确)。
★ 最后两列是"能量全部在折叠那一个频率上"的精确证据:取一段整周期的短记录(长度 = `f_s / gcd(700, f_s)` 个样本, 于是 700 Hz 与 `f_alias` 都正好落在 DFT 的格点上),在那一小段上做 DFT —— 折叠那一格的幅度恰好 1.000000,其余每一格恰好 0(正交性,不是近似)。
★ 注脚:`f = f_s/2`(Nyquist 上)是退化档 —— 正负频率折进同一格,投影读数变成 2.000000, 折叠闭式在这一档不唯一。本表不含该档(它单列在页面正文与判据③)。
f_s (Hz)Nyquistk = round(700/f_s)手算 f_alias逐位最大差投影 @ f_alias整周期记录长其余格最大

16.2 把点连回去:重建核决定一切

本段只增加一样东西:重建核(旋钮在「零阶保持 / 线性 / 加窗 sinc」之间切;采样率不动)。
上一段的缺口:16.1 只说清了"采样的那一步偷走了什么",可"怎么把点连回去"还没有说 —— 而连法决定了它像不像原来那个连续的东西。
参考实现的定义只有一句(J. O. Smith 的 Theory of Ideal Bandlimited Interpolation): x̂(t) = Σn x(nTs)·hs(t − nTs), 而理想核是 hs(t) = sinc(Fs·t) = sin(πFst)/(πFst), 它的原话是"peak is at amplitude 1, and zero-crossings occur at all nonzero integers… Thus, the interpolation goes precisely through the existing samples, as it should."
★ 于是无论用哪个核,采样点上都必须精确过点:h(0) = 1、h(整数 ≠ 0) = 0。 这条可以拿笔逐点核:盒式核 `h(u) = 1 (|u| < 0.5)` 在 u = ±1 处是 0 ✓;线性核 `h(u) = max(0, 1−|u|)` 在 u = ±1 处是 0 ✓;sinc 在非零整数处是 0 ✓。实测:盒式与线性恰好 0,sinc 是 5.6e−17。
★ 差别全在两点之间(表②):一个 300 Hz 的带限信号,两点之间的最大误差 —— 盒式 0.951、线性 0.392、加窗 sinc M=8 是 2.9e−3、M=32 是 3.6e−5。
★ 而核的代价要一起说(PBRT 原话):盒式"just about the worst filter possible… allows high-frequency sample data to leak into the reconstructed values. This causes postaliasing.", 而 sinc"prone to ringing when the underlying function has frequencies beyond the Nyquist limit (Gibbs phenomenon)… the sinc filter has infinite support" —— 本页的加窗 sinc 就是拿这个代价换精度的: M 越大越准,但支撑越宽、振铃越明显(表② 的"支撑"列)。
★ 还有一条必须先说清,否则 16.3 就没有理由:如果信号本身越过了 Nyquist,换核救不回来 —— 同一个 700 Hz 信号(Nyquist = 500),四个核在两点之间的误差是 1.670 / 1.530 / 1.910 / 1.913, 全都在"完全重建不出来"的量级。核只能把带内的东西拼得更好。

输入 这个阶段在执行的规则(照抄件二真源码) 输出 样本序列 x[n] (16.1 采出来的那串) 旋钮:重建核 盒式 / 线性 / 加窗 sinc 位移以样本为单位 u = t·f_s − n ★ 本段只加"核" // 理想核(J.O. Smith 的定义式) const sinc = u => Math.sin(Math.PI*u)/(Math.PI*u); const hBox = u => Math.abs(u) < 0.5 ? 1 : 0; const hLin = u => Math.max(0, 1 - Math.abs(u)); const hSinc = (u,M) => sinc(u) * hann(u,M); // 重建:每个样本各出一个核,叠起来 let v = 0; for (let m = n - M; m <= n + M; m++) v += x[m] * h(t * fs - m, M); // 过点:t = m/fs 时只有 m 自己贡献 ⇒ v = x[m] 300 Hz 两点之间: 盒式 0.951 线性 0.392 sinc M=8 2.9e−3 sinc M=32 3.6e−5 700 Hz(越 Nyquist): 1.670 / 1.530 1.910 / 1.913 ★ 核救不回来
定点观察:把「重建核」在四档之间切 → 看采样点之间那几段(采样点本身四个核都过)→ 300 Hz 的带限信号:盒式两点之间能差 0.951、线性 0.392,加窗 sinc M=32 只差 3.6e−5; 而把信号换成 700 Hz(越过 Nyquist)之后,四个核全都差 1.5 以上 —— 核只能把带内的东西拼得更好。
表② 四个重建核(测试信号 300 Hz 与 700 Hz @ 1000 Hz 采样,Nyquist = 500 Hz)。
★ 「过点最大误差」在所有采样点上取最大:盒式与线性恰好 0,加窗 sinc 是 5.6e−17 / 2.2e−16 (浮点噪声)—— 这就是 h(0)=1、h(整数≠0)=0 这条条件,可以逐点拿笔核。
★ 「两点之间」在内部区间取最大(刻意避开两端:那里核的支撑伸到序列外,是截断误差不是核的误差)。
★ 最后两列并排看:带内(300 Hz)误差随支撑单调下降,越过 Nyquist(700 Hz)之后四个核全部失效。
核支撑h(0) / h(±1)过点最大误差两点之间(300 Hz)两点之间(700 Hz)

16.3 先滤再采:唯一的治本手段

本段只增加一样东西:采之前的低通(旋钮是五档低通;采样率回到 1000 Hz 不动)。
上一段的缺口:16.2 有了最好的核,可如果一开始就采错了,再好的核也救不回来(表② 最后两列)。 于是唯一治本的做法只有一个:在采样之前,先把高于 Nyquist 的东西滤掉。 Whitted 1980 的原话:"Since aliasing is the result of undersampling during the display process, the most straightforward cure is to low-pass filter the entire image before sampling for display." PBRT 把它写成三步里的第二步:"Prefilter the function L̃ to remove any frequencies past the Nyquist limit for the pixel spacing."
★ 这里有一条必须写死的工程事实:低通要加在"还没采样"的那一侧。 本页的管线是 8000 Hz 上造信号 → 低通 → 降到 1000 Hz。 (本轮量测时我第一版把低通加在已经欠采样的那串数上 —— 那时信息早就丢了, 量出来是"滤波对混叠毫无用处"这个错结论。)
★ 而代价可以手算:低通把 700 Hz 压掉多少,混叠分量就只剩多少 —— 混叠幅度 = `|H(f)|`。本页用的二阶节级联(Butterworth,Q = 1/√2)的双线性变换精确式是 |H| = (1 + q⁴)−n/4,其中 q = tan(πf/f_s) / tan(πf_c/f_s)。
★★ 表③ 里并排摆了三个口径,这是本段最值钱的一列: 同一个 4 阶 fc = 450 的滤波器,双线性精确式给 0.138482,教科书上的模拟原型 `1/√(1+(f/f_c)^(2n))` 给 0.168351 —— 差 21.6%。
★ 这 21.6% 当初被本文写成"双线性变换扭了一下频率轴"一件事 —— 第 2 课(17.3 / 17.4)把它量清楚了, 其实是两笔账:形状 14.60%(`(1+q⁴)^(−n/4)` 对应的是"两个相同 Q = 1/√2 的节", 而名副其实的 4 阶 Butterworth 要两个 Q 不同的节:0.5412 与 1.3066 —— 同一个 700 Hz 上两者给 0.138482 与 0.158705)× 扭曲 6.08%(同一台数字机 vs 教科书模拟原型)= 1.1460 × 1.0608 = 1.2157。 所以手算要用带 `tan` 的那个式子;实测与它的差 ≤ 1e−4。(判据⑧ 现在同时断言这条乘法恒等式。)
★ 最后把账单摊开:混叠分量从 1.000000(不滤)一路压到 0.000947(6 阶 fc=225), 可700 Hz 那一条本身也被压到 0.000947 —— 这就是"先滤再采"的代价: 想在 Nyquist 以下采样,就得接受"高于 Nyquist 的内容保不住"。另一条路是抬高采样率(表① 的 fs = 4000 档: 完全不混叠,但同样 2048 个样本只覆盖 1/4 的时间)。

输入 这个阶段在执行的规则(照抄件二真源码) 输出 连续信号(8000 Hz 上) → 低通 → 降到 1000 Hz 旋钮:低通五档 关 / 2 阶 450 / 4 阶 450 / 4 阶 225 / 6 阶 225 Q = 1 / √2 (Butterworth 二阶节) ★ 本段只加"采前低通" // 二阶节系数(双线性变换,带预畸变) const K = Math.tan(Math.PI * fc / fs), Q = 1/Math.SQRT2; const d = 1 + K/Q + K*K; a1 = 2*(K*K-1)/d; a2 = (1-K/Q+K*K)/d; b0 = K*K/d; b1 = 2*b0; b2 = b0; // 级联 order/2 级,然后每 R 个取一个 const y = new Float64Array(N); for (let m = 0; m < N; m++) y[m] = hi[m * R]; // 手算:混叠幅度 = |H(700)| = (1+q⁴)^(−n/4) 混叠分量: 关 1.000000 2阶450 0.372144 4阶450 0.138491 6阶225 0.000947 4 阶 450 的两个口径: 双线性 0.138482 模拟原型 0.168351 差 21.6% = 形状 14.6% × 扭曲 6.1%
定点观察:把「采前低通」从「关」一路切到「6 阶 225」→ 看读数行里的"混叠分量"与 表③ 第一列 → 1.000000 → 0.372144 → 0.138491 → 0.000947; 再对比表③ 的第二、三列:手算要用带 tan 的双线性精确式(差 ≤1e−4),教科书那条模拟原型在 4 阶 450 上差 21.6%。
表③ 采前低通五档(信号 700 Hz,管线:8000 Hz 造信号 → 低通 → 降到 1000 Hz 采样,测 300 Hz 处的混叠分量)。
★ 「手算(双线性精确式)」= (1 + q⁴)−n/4,q = tan(πf/f_s)/tan(πf_c/f_s); 「手算(模拟原型)」= 1/√(1+(f/f_c)^(2n)),就是教科书里那条。
★ 最右两列是本段的重点:两个口径在 4 阶 fc=450 上差 21.6%,因为双线性变换把频率轴扭了一下, 所以"手算"要用带 `tan` 的那一条。实测与它逐档差 ≤ 1e−4。
★ 这 21.6% 不是一笔账(第 2 课 17.3 / 17.4 拆清):形状 14.60%("两个相同 Q = 1/√2 的节" ≠ 真 4 阶 —— 真 4 阶要 Q = 0.5412 与 1.3066)× 扭曲 6.08%(同一台数字机 vs 模拟原型)= 1.2157。 判据⑧ 现在断言这条乘法恒等式(`形状 × 扭曲 === 总`,差 0)。
★ 代价是同一件事:低通把 700 Hz 压到多少,混叠就只剩多少(两列本质是同一个数)。
采前低通实测混叠分量手算(双线性精确式)差手算(模拟原型)两口径之差

16.4 图像里的同一道数学

本段只增加一样东西:维度(把同一条式子搬到 2D 图像上)(旋钮是「图像宽度」,条纹数恒为宽度的 11/16)。
上一段的缺口:前三段都是一条波形(1D)。可图像是 2D 的 —— "采样率不足"在图像里长什么样?
答案就在 16.1 的那条折叠闭式里,只要把"每秒采几个点"换成"整幅画几个像素": 一个宽 `W` 像素、画了 `c` 条条纹的图,pixel[i] = cos(2π c i / W) —— 而它逐位等于"画 c mod W 条"(以及它的镜像 W − (c mod W) 条)的图。 所以可见条数 = `min(c mod W, W − (c mod W))`。
能拿笔核:W = 16 画 11 条 ⇒ `11 mod 16 = 11`、镜像 `16 − 11 = 5` ⇒ 看见 5 条; W = 64 画 44 条 ⇒ `44 mod 64 = 44`、镜像 `64 − 44 = 20` ⇒ 看见 20 条。 实测:两种画法的像素序列逐位最大差 1.0e−14,且折叠那条的投影恰好是 1.000000、其它整数频率 ≤ 4.6e−15。
★ 这就是图形学里那个老故事的同一件事 —— Crow 1977 的原话: "The consequence of failing to filter the signal properly before sampling is known as 'aliasing.'… highly periodic images of scenes involving, for example, picket fences or venetian blinds may appear, when sampled, to be made up of a few broad strips rather than many fine lines." 画 11 条、只有 16 个像素,看到的就是 5 条粗纹 —— 摩尔纹不是画错了,是采样。
★ 而 16.2 与 16.3 的结论在这里原封不动地落地:重建核就是"放大时怎么填那格" (本页画面下半条带用的就是零阶保持:每个像素涂满它自己那一格); 先滤再采就是 Mipmap 与多采样(Williams 1983:"a set of prefiltered source images may be created… minimizes aliasing effects";GPU Gems 2:"each final pixel is a straight average of the samples it covers. This is equivalent to a box (that is, constant) filter of width 1.")。

输入 这个阶段在执行的规则(照抄件二真源码) 输出 条纹信号 连续:画 c 条 旋钮:宽度 W c = round(11/16 · W) 一维、只有 W 个格子 (整幅=周期采样) ★ 本段只加"维度" // 画 W 个像素、c 条条纹 for (let i = 0; i < W; i++) px[i] = Math.cos(2 * Math.PI * c * i / W); // 折叠:与 16.1 同一个式子,只是单位换成"像素" const cm = ((c % W) + W) % W; const rep = Math.min(cm, W - cm); // 判据:与"画 rep 条"逐位比较 d = Math.max(d, Math.abs(px[i] - Math.cos(2*Math.PI*rep*i/W))); W=16 画 11 条: 看见 5 条 逐位差 1.0e−14 W=64 画 44 条: 看见 20 条 分辨率越低: 看不见的细节越多 ★ 摩尔纹 = 折叠
定点观察:把「图像宽度」从 12 切到 64 → 看画面里从上到下那三条带 与读数行里的"看见几条"+"三条带振幅" → 宽度越小、看见的条纹越少(12→4 条、16→5、24→7、32→10、64→20); 而且它与"直接画折叠后那个条数"的图逐位相同。 三条带是"预滤波到底做了什么"的全部答案:A 直接采(满对比度假条纹)、 B 像素内积分平均(假条纹还在,但振幅只剩 |sinc(c/W)| ≈ 0.385 —— 盒式核是"最差的滤波器",PBRT 原话)、 C 按 Nyquist 限带(那一整个频率成分被滤掉,振幅 0 ⇒ 一片均匀灰)。这就是 Mipmap 在做的事。
表④ 五档图像宽度(条纹数恒为宽度的 11/16 ⇒ `c = round(11/16·W)`)。
★ 「折叠代表」= `min(c mod W, W − (c mod W))`;「与代表逐位差」= 这两种画法的像素序列逐位比较的最大差 —— 这不是相似,是同一张图(残下的 1e−14 是 `cos` 的舍入)。
★ 「投影 @ 代表」在折叠代表那个整数条数上是 1.000000,「其它条数」是其余所有整数频率上的最大投影(≤ 4.6e−15)。
★ 最后两列是"采前滤波"的账单,也是画面三条带的振幅:B = 像素内积分平均,振幅恰好 `|sinc(c/W)|` (画面里那条带是像素内 64 点中点法数值积分算出来的,与这条闭式逐像素对上,残差 ≤ 1e−4 —— 那 1e−4 就是这条数值积分自己的分辨率);C = 按 Nyquist 限带, 因为 `c/W ≈ 0.69 > 0.5`,那一个频率成分整个被滤掉 ⇒ 振幅恰好 0、条纹变成一片均匀灰。 两列几乎不随 W 变,是因为比值恒定 —— 这本身也说明"账单与分辨率无关,只与 c/W 有关"。
宽度 W画的条数 cc mod W镜像 W − (c mod W)折叠代表与代表逐位差投影 @ 代表其它条数B 振幅 |sinc(c/W)|C 振幅

这一课带走了什么

三行式子,全部可手算:① 折叠律 `f_alias = |f − k·f_s|`;② 重建核过点 `h(0)=1、h(整数≠0)=0`; ③ 混叠幅度 = 采前低通的 `|H(f)|`。
而最该记住的一句是:混叠不是"像",是"相同" —— 700 Hz 采成 1000 Hz 与直接采 300 Hz, 是两个逐位相同的序列;画 11 条、只有 16 个像素,与画 5 条是同一张图。 信息在那一步就没了,任何事后修补都补不回来(表② 最后两列:越过 Nyquist 之后四个核全废)。
所以只有两条路:要么把采样率抬到 Nyquist 以上(表① fs = 4000 档,完全不混叠,代价是样本数);要么先滤再采(表③, 代价是高频内容本身也被压掉)。这就是 Mipmap、MSAA、纹理过滤在做的事。
下一课指针:这一课反复用到"低通",可低通到底是什么形状的一条曲线? 下一课把二阶节的系数摊开,看系数怎么决定频响曲线 —— 顺便会看到双线性变换把频率轴扭了一下这件事的真身。

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

运行时自检 等待脚本…