采样:为什么"多采几个点"不是答案
全课题骨架:采样不是"取点",是"用一串离散值 + 一个重建核,去猜那个连续的东西"。
核决定它像不像;而"先滤再采"是唯一的治本手段 —— 图像里的摩尔纹与声音里的混叠,是同一道数学。
上一课的输入是"已经离好的东西"(顶点、像素、瓦片、样例图)。可现实里的信号是连续的 ——
把它变成一串数的那一步,偷走了什么?答案是 频率:欠采样时,一串 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)。
本页的原理与四张数值表照常可读(画面只是把同一份状态画出来)。
(脚本启动中…)
定点观察:(脚本启动中…)
本帧的读数
(脚本启动中…)
(点「扫一遍采样率」或「跑一遍读数」就在这里出数)
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]` 以内的频率。所以判据写成"两个序列逐位相同",不写成"原频率处为零"。
| f_s (Hz) | Nyquist | k = 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,
全都在"完全重建不出来"的量级。核只能把带内的东西拼得更好。
| 核 | 支撑 | 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 的时间)。
| 采前低通 | 实测混叠分量 | 手算(双线性精确式) | 差 | 手算(模拟原型) | 两口径之差 |
|---|
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.")。
| 宽度 W | 画的条数 c | c 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、纹理过滤在做的事。
下一课指针:这一课反复用到"低通",可低通到底是什么形状的一条曲线?
下一课把二阶节的系数摊开,看系数怎么决定频响曲线 —— 顺便会看到双线性变换把频率轴扭了一下这件事的真身。
本页用到的做法,出处在哪
- 采样定理 —— C. E. Shannon, Communication in the Presence of Noise(Proc. IRE, 1949): "Theorem 1: If a function f(t) contains no frequencies higher than W cps, it is completely determined by giving its ordinates at a series of points spaced 1/2W seconds apart." ★ 如实标注:折叠闭式 `f_alias = |f − k·F_s|` 不在这篇论文里,它出自下面的 J. O. Smith。
- 折叠闭式与"镜像" —— Julius O. Smith III, Mathematics of the Discrete Fourier Transform, §D.1: "…sampling at Fs samples per second causes that energy to alias to a lower frequency. If we write the original frequency as f, then the new aliased frequency is f_alias = |f − k·Fs|… This phenomenon is also called 'folding', since… 'mirror image' of… about Fs/2." ccrma.stanford.edu/~jos/mdft/Aliasing_Sampled_Signals.html
- 带限可精确重建 —— 同书 §D.3("x(t) can be exactly reconstructed from its samples x(nT) if…"): ccrma.stanford.edu/~jos/mdft/Sampling_Theorem.html
- 理想重建核与"必然过点" —— J. O. Smith, Theory of Ideal Bandlimited Interpolation:
"
h_s(t) ≜ sinc(Fs·t) ≜ sin(π·Fs·t)/(π·Fs·t)… 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." ccrma.stanford.edu/~jos/resample/Theory_Ideal_Bandlimited_Interpolation.html ★ 如实标注:插值核的一般过点条件 `h(0)=1、h(n≠0)=0` 的教科书原页我们没抓到 (Keys 1981 / Mitchell–Netravali 1988 正文被墙),所以页面只引 sinc 特例,并把条件本身写成可逐点核的断言。 - 先滤再采、以及三个核的代价 —— PBRT, Image Reconstruction: 盒式 "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"; 三步 "1. Reconstruct… 2. Prefilter the function L̃ to remove any frequencies past the Nyquist limit… 3. Sample": pbr-book.org/3ed-2018/…/Image_Reconstruction
- 低通是治本手段 —— Turner Whitted, An Improved Illumination Model for Shaded Display(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." eng.utah.edu/~cs6958/papers/p343-whitted.pdf
- "栅栏变宽条" —— Franklin Crow, The Aliasing Problem in Computer-Generated Shaded Images(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." dl.acm.org/doi/pdf/10.1145/280811.280976
- Mipmap 与多采样 —— Lance Williams, Pyramidal Parametrics(SIGGRAPH 1983): "a set of prefiltered source images may be created… advances a 'pyramidal parametric' prefiltering and sampling geometry which minimizes aliasing effects"(★ 正文 PDF 未获取,引 ACM DL 官方摘要); GPU Gems 2, Ch.21(NVIDIA):"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." developer.nvidia.com/gpugems/gpugems2/…/chapter-21
- 像素该是"滤波后的样本" —— Marc Levoy, The early history of point-based graphics: "To avoid aliasing artifacts in computer-generated images, each pixel should be assigned not a point sample of the scene but instead a sample of the convolution of the scene by a filter function." scroll.stanford.edu/papers/points/levoy-pointbook-ch2.pdf
- 规范里的抗混叠原话 —— W3C Web Audio API: "…care MUST be taken to discard (filter out) the high-frequency information higher than the Nyquist frequency before converting the waveform to a digital form. If this is not done, then aliasing of higher frequencies… will fold back as mirror images into frequencies lower than the Nyquist frequency."; 以及振荡器那条(★ 本页据此不声称原生振荡器一定带限): "The actual waveforms produced by the oscillator may differ to prevent aliasing affects." w3.org/TR/webaudio/
- ★ 四处如实标注:① 折叠闭式的一手出处是 J. O. Smith,不是 Shannon; ② 我们没有找到任何一手文献用一句话并列断言"摩尔纹与音频混叠是同一个现象" —— 可确证的是权威课程把它俩归入同一个 "aliasing"(UMD CMSC427:"Spatial: Jaggies, Moire"; UC Berkeley CS184 L26:"sampling a fence (periodic structure) at a sampling rate similar to the lattice period will cause Moiré patterns"),所以本页写的是"同一道数学(同一个折叠律)",不声称某篇论文这样说过; ③ 插值核一般过点条件的教科书原页未获取;④ Williams 1983 只能引摘要。