Clear the analytic rate: R now comes from the timeseries-entropy estimator
Remove the spectral/sub-threshold formula (theory.ts), the math section,
the dashed R line, the S(f)=1 threshold, dither, and true_rate.py. The
reference rate R is a placeholder; the UI shows the uvx command that
runs the unbiased Monte-Carlo estimator from the companion
timeseries-entropy package, and a short note explains the method.
In-browser estimation is planned.
15 changed files+121−870
README.mdmodified+25−42View file
@@ -3,41 +3,27 @@
33 Interactive exploration of how compressible quantized time series are.
44
55 The generating model is: i.i.d. Gaussian noise (std σ, measured in quantization
6-steps) → FIR filter → optional additive uniform dither on [-½, ½) → round to
7-integers. The app shows the filter (convolution kernel and frequency response,
8-with cutoffs in Hz against a chosen sample rate), a window of the generated
9-integer signal (stationary by default, with a play toggle to let it stream
10-endlessly), and the measured compression of a block of the generated integers
11-under nine methods — zlib, zstd, and an rANS entropy coder, each raw,
12-delta-coded, and LPC-residual-coded — as bits per sample and as ratio against
13-raw int16 storage. The predictor order and the block size are controls, so the
14-measurement can be pushed from 10k to a million samples and LPC from order 1
15-to 128. Each prefilter group also carries a hollow bar: the order-0 entropy of
16-the stream being coded, the limit a per-sample entropy coder cannot beat, which
17-ANS misses by 1–2% (its symbol table plus its own arithmetic loss).
18-
19-Alongside the measurements it plots a theoretical bits/sample R — the smaller
20-of a spectral estimate and a rigorous one-sample ceiling:
21-
22-```
23-R = min(Rspec, Rsamp)
24-Rspec = ∫₀¹ ½ log₂( 2πe (S(f) + ν) ) df S(f) = σ²|H(f)|², ν = 1/12 (1/6 dithered)
25-Rsamp = H( round(N(0, v) [+ U(-½,½) with dither]) ) v = σ² Σ h²
26-```
27-
28-Rspec is the Zamir–Feder rate of the dithered quantizer counted per Fourier
29-mode: the signal modes are independent Gaussians of variance S(f), and the
30-i.i.d. roundoff(+dither) noise is Gaussianized per mode by the CLT, so it
31-enters at its full variance ν — not at the entropy power 1/(2πe) an aligned
32-scalar quantizer would charge (the lattice lives in the sample basis; a dead
33-band inside a live process costs ≈0.25 bits/mode, not zero). At high SNR it
34-reduces to the Kolmogorov rate ½log₂(2πe σ²) + ∫log₂|H| df. Where the whole
35-process sits below the quantization step, Rspec bottoms out while the true
36-rate collapses; subadditivity H(z) ≤ Σ H(zₙ) makes Rsamp — the exact marginal
37-entropy of one stored sample — a true upper bound with the right collapse,
38-and the min selects it exactly there. Monte-Carlo puts R within ~0.01–0.02
39-bits/sample for v ≳ 0.25 (worst ~+0.03 at the branch crossover). LPC + ANS
40-should approach R; probing where the approximation holds is the point.
6+steps) → FIR filter → round to integers. The app shows the filter (convolution
7+kernel and frequency response, with cutoffs in Hz against a chosen sample
8+rate), a window of the generated integer signal (stationary by default, with a
9+play toggle to let it stream endlessly), and the measured compression of a
10+block of the generated integers under nine methods — zlib, zstd, and an rANS
11+entropy coder, each raw, delta-coded, and LPC-residual-coded — as bits per
12+sample and as ratio against raw int16 storage. The predictor order and the
13+block size are controls, so the measurement can be pushed from 10k to a million
14+samples and LPC from order 1 to 128. Each prefilter group also carries a hollow
15+bar: the order-0 entropy of the stream being coded, the limit a per-sample
16+entropy coder cannot beat, which ANS misses by 1–2% (its symbol table plus its
17+own arithmetic loss).
18+
19+The reference rate R — the entropy rate of the process, the bits/sample limit
20+no lossless method can beat — comes from the companion
21+[timeseries-entropy](https://github.com/concept-collection/timeseries-entropy)
22+package: an unbiased Monte-Carlo estimator of H(z_next | a long past), by Gibbs
23+sampling the latent Gaussian under the rounding constraints and applying
24+Rhee–Glynn randomized telescoping to the sampled chain. The app shows the exact
25+command to run it at the current settings; estimating R in the browser is
26+planned, and until then the UI shows a placeholder for it.
4127
4228 ## Run it
4329
@@ -50,17 +36,14 @@ npm run dev
5036
5137 ```
5238 src/model/ the latent source (fixed seeded randomness indexed by sample
53- position, convolved zero-phase with the kernel on demand),
54- FIR presets, and the theoretical-rate formula
39+ position, convolved zero-phase with the kernel on demand)
40+ and the FIR presets
5541 src/compress/ lossless codecs run in the browser: zlib (fflate), zstd (wasm),
5642 ans.ts (a bit-identical port of simple_ans), and FLAC-style
5743 integer LPC; borrowed from entropy-quantized-linear-transform
5844 src/worker/ the codecs run off the main thread on a debounced parameter set
59-src/components/ controls, filter plots, signal canvas, compression chart
60-scripts/ true_rate.py — Monte-Carlo ground truth for R (Gibbs over the
61- latent Gaussian given the observed integers, Rao-Blackwellized
62- next-sample pmf); the app prints the exact command to run.
63- Requires numpy only.
45+src/components/ controls, filter plots, signal canvas, compression chart,
46+ and the reference-rate method note
6447 ```
6548
6649 Every reported size round-trips through the decoder and includes whatever the
scripts/true_rate.pydeleted+0−357View file
@@ -1,357 +0,0 @@
1-#!/usr/bin/env python3
2-"""Monte-Carlo ground truth for the theoretical rate R shown in the app.
3-
4-Model (matching the app): x iid N(0, sigma^2) -> y = h * x -> optional
5-additive uniform dither on [-1/2, 1/2) -> z = round(.). The entropy rate
6-H = E[-log2 P(z_next | past)] is the true lossless limit in bits/sample.
7-(The app applies the kernel zero-phase; a time shift does not change the
8-law of the process, so causal convolution is used here.)
9-
10-Method
11- 1. Draw a past: sample x (and dither) from the prior and push it through
12- the pipeline to get z_1..z_M.
13- 2. The generating x is itself an exact draw from p(x | z_1..z_M), so a
14- Gibbs chain started there is already in stationarity - no burn-in
15- bias, only autocorrelation.
16- 3. Gibbs-sample x | z: this posterior is a box-truncated multivariate
17- normal, and each conditional x_i | rest is N(0, sigma^2) truncated to
18- an interval read off the <= L constraint boxes x_i appears in. With
19- dither on, the dither values are extra latents with uniform
20- conditionals. Coordinates a multiple of L apart share no constraint,
21- so each of the L "colors" is updated as one vectorized block.
22- 4. Rao-Blackwellization: given a chain state, the next sample is
23- z_next = round(c + h_0 * x_free (+ d)) with x_free ~ N(0, sigma^2)
24- still unconstrained, so P(z_next = j | state) has a closed form.
25- Averaging these pmfs over the chain gives P(z_next | past) exactly in
26- the limit; its entropy is the conditional entropy for that past.
27- 5. Average over independent pasts; report mean +/- standard error.
28-
29-The estimate targets H(z_{M+1} | z_1..z_M), which is an upper bound on the
30-rate and decreases toward it as --past grows beyond the memory of the
31-process. More --sweeps reduces the (downward) plug-in bias from noise in
32-the averaged pmf; narrowband filters and large sigma mix more slowly and
33-deserve more sweeps.
34-
35-Requires numpy only.
36-"""
37-
38-import argparse
39-import math
40-
41-import numpy as np
42-
43-SQRT2PI = math.sqrt(2 * math.pi)
44-
45-
46-# ---------------------------------------------------------------- kernels
47-# Ported from src/model/filters.ts; must stay in step with it.
48-
49-def windowed_sinc_lowpass(fc, taps):
50- n = taps | 1
51- mid = (n - 1) / 2
52- i = np.arange(n)
53- t = i - mid
54- sinc = np.where(t == 0, 2 * fc, np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
55- w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
56- h = sinc * w
57- return h / h.sum()
58-
59-
60-def design_kernel(args):
61- if args.filter == 'none':
62- return np.array([1.0])
63- if args.filter == 'moving-average':
64- return np.full(args.width, 1.0 / args.width)
65- if args.filter == 'lowpass':
66- return windowed_sinc_lowpass(args.high / args.rate, args.taps)
67- if args.filter == 'bandpass':
68- lo = windowed_sinc_lowpass(args.low / args.rate, args.taps)
69- hi = windowed_sinc_lowpass(args.high / args.rate, args.taps)
70- return hi - lo
71- if args.filter == 'first-difference':
72- return np.array([1.0, -1.0])
73- raise ValueError(args.filter)
74-
75-
76-# ------------------------------------------------- the app's formula for R
77-# Ported from src/model/theory.ts. Phi via math.erf (machine precision).
78-
79-TWO_PI_E = 2 * math.pi * math.e
80-
81-
82-def quantized_gaussian_entropy(s):
83- if s <= 0.02:
84- return 0.0
85- # The discrete entropy approaches 1/2 log2(2 pi e s^2) from above like
86- # log2(e)/(24 s^2) (the delta^2/24 Fisher-information correction); at
87- # s >= 6 the corrected asymptote is within 2e-6 bits.
88- if s >= 6:
89- return 0.5 * math.log2(TWO_PI_E * s * s) + math.log2(math.e) / (24 * s * s)
90- zmax = int(math.ceil(8 * s + 4))
91- H = 0.0
92- prev = 0.5 * (1 + math.erf((-zmax - 0.5) / (s * math.sqrt(2))))
93- for z in range(-zmax, zmax + 1):
94- cur = 0.5 * (1 + math.erf((z + 0.5) / (s * math.sqrt(2))))
95- p = cur - prev
96- prev = cur
97- if p > 0:
98- H -= p * math.log2(p)
99- return H
100-
101-
102-def dithered_quantized_gaussian_entropy(s):
103- """Exact entropy of round(N(0, s^2) + U[-1/2, 1/2)) — the marginal of a
104- stored sample with dither on. The pmf has the closed form
105- p_j = s * (G((j+1)/s) - 2 G(j/s) + G((j-1)/s)) with G(t) = t Phi(t) + phi(t)
106- the antiderivative of Phi; machine-exact at every s via math.erf."""
107- if s <= 0:
108- return 0.0
109-
110- def G(t):
111- return t * 0.5 * (1 + math.erf(t / math.sqrt(2))) + math.exp(-0.5 * t * t) / SQRT2PI
112-
113- jmax = int(math.ceil(8 * s + 2))
114- H = 0.0
115- for j in range(-jmax, jmax + 1):
116- p = s * (G((j + 1) / s) - 2 * G(j / s) + G((j - 1) / s))
117- if p > 0:
118- H -= p * math.log2(p)
119- return H
120-
121-
122-def formula_rates(kernel, sigma, dither, points=8192):
123- """The app's R = min(Rspec, Rsamp). Rspec charges the rounding(+dither)
124- noise at its full variance nu per Fourier mode; Rsamp is the exact
125- marginal entropy of one stored sample, a subadditivity upper bound that
126- takes over when the whole process is sub-threshold. Midpoint grid on
127- [0, 1/2]; |H| is symmetric, so the grid mean equals the unit-circle
128- integral. Returns (rspec, rsamp)."""
129- nu = 1 / 6 if dither else 1 / 12
130- f = (0.5 * (np.arange(points) + 0.5)) / points
131- w = -2j * np.pi * np.outer(f, np.arange(len(kernel)))
132- S = sigma * sigma * np.abs(np.exp(w) @ kernel) ** 2
133- rspec = float(np.mean(0.5 * np.log2(TWO_PI_E * (S + nu))))
134- v = sigma * sigma * float(np.sum(np.asarray(kernel) ** 2))
135- rsamp = (dithered_quantized_gaussian_entropy(math.sqrt(v)) if dither
136- else quantized_gaussian_entropy(math.sqrt(v)))
137- return rspec, rsamp
138-
139-
140-# ------------------------------------------------ vectorized normal helpers
141-
142-def norm_pdf(t):
143- return np.exp(-0.5 * t * t) / SQRT2PI
144-
145-
146-def norm_cdf(t):
147- """Abramowitz-Stegun 26.2.17; |error| < 7.5e-8, plenty for sampling and
148- for pmf weights whose entropy is wanted to ~1e-4 bits."""
149- t = np.asarray(t, dtype=float)
150- z = np.abs(t)
151- k = 1.0 / (1.0 + 0.2316419 * z)
152- poly = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + k * 1.330274429))))
153- tail = norm_pdf(z) * poly
154- return np.where(t >= 0, 1.0 - tail, tail)
155-
156-
157-def norm_ppf(p):
158- """Acklam's rational approximation to the standard normal quantile."""
159- a = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
160- 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00)
161- b = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
162- 6.680131188771972e+01, -1.328068155288572e+01)
163- c = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
164- -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00)
165- d = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
166- 3.754408661907416e+00)
167- p = np.asarray(p, dtype=float)
168- x = np.empty_like(p)
169- plow, phigh = 0.02425, 1 - 0.02425
170-
171- lo = p < plow
172- hi = p > phigh
173- mid = ~(lo | hi)
174-
175- if mid.any():
176- q = p[mid] - 0.5
177- r = q * q
178- x[mid] = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / \
179- (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
180- if lo.any():
181- q = np.sqrt(-2 * np.log(p[lo]))
182- x[lo] = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
183- ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
184- if hi.any():
185- q = np.sqrt(-2 * np.log(1 - p[hi]))
186- x[hi] = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
187- ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
188- return x
189-
190-
191-def trunc_std_normal(lo, hi, rng):
192- """Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into
193- the lower tail so the CDF differences keep precision."""
194- flip = (lo + hi) > 0
195- a = np.where(flip, -hi, lo)
196- b = np.where(flip, -lo, hi)
197- Fa = norm_cdf(a)
198- Fb = norm_cdf(b)
199- u = Fa + (Fb - Fa) * rng.random(a.shape)
200- x = norm_ppf(np.clip(u, 1e-300, 1 - 1e-16))
201- x = np.where(flip, -x, x)
202- return np.clip(x, lo, hi)
203-
204-
205-# --------------------------------------------------------- the RB next-pmf
206-
207-def big_g(t):
208- """G(t) = t Phi(t) + phi(t), the antiderivative of Phi."""
209- return t * norm_cdf(t) + norm_pdf(t)
210-
211-
212-def next_pmf(c, s0, dither):
213- """P(z_next = j | chain state): round(c + N(0, s0^2) (+ U(-1/2,1/2)))."""
214- s = max(s0, 1e-12)
215- half = 8 * s + (1.0 if dither else 0.0) + 1.0
216- js = np.arange(math.floor(c - half), math.ceil(c + half) + 1)
217- if dither:
218- # Integrating the Gaussian bin probability over the dither gives a
219- # second difference of G; as s -> 0 it degrades gracefully to the
220- # uniform-overlap width.
221- p = s * (big_g((js + 1 - c) / s) - 2 * big_g((js - c) / s) + big_g((js - 1 - c) / s))
222- else:
223- edges = norm_cdf((np.append(js, js[-1] + 1) - 0.5 - c) / s)
224- p = np.diff(edges)
225- return js, np.maximum(p, 0.0)
226-
227-
228-# ---------------------------------------------------------------- one past
229-
230-def conditional_entropy_of_one_past(kernel, sigma, dither, M, sweeps, rng):
231- h = np.asarray(kernel, dtype=float)
232- L = len(h)
233- hr = h[::-1]
234- N = M + L - 1 # latents covering the windows of z_1..z_M
235-
236- # The past, with its true latents as the (stationary) chain start.
237- x = sigma * rng.standard_normal(N)
238- y = np.convolve(x, h, mode='valid')
239- d = (rng.random(M) - 0.5) if dither else None
240- z = np.floor(y + (d if dither else 0.0) + 0.5)
241-
242- # Boxes and y live in padded arrays so that every coordinate x_i sees
243- # exactly L constraint rows (rows outside the data are unconstrained).
244- P = L - 1
245- ypad = np.zeros(M + 2 * P)
246- lo = np.full(M + 2 * P, -np.inf)
247- hi = np.full(M + 2 * P, np.inf)
248-
249- def set_boxes():
250- dd = d if dither else 0.0
251- lo[P:P + M] = z - 0.5 - dd
252- hi[P:P + M] = z + 0.5 - dd
253-
254- set_boxes()
255-
256- # Color classes: coordinates L apart share no constraint row, so a class
257- # updates as one vectorized block. Row i+j (padded) carries coefficient
258- # h[j] for coordinate i.
259- classes = [np.arange(c0, N, L) for c0 in range(L)]
260- rowmats = [idx[:, None] + np.arange(L)[None, :] for idx in classes]
261- nonzero = h != 0
262-
263- s0 = sigma * abs(h[0])
264- pmf = {}
265-
266- for _ in range(sweeps):
267- ypad[P:P + M] = np.convolve(x, h, mode='valid') # kill fp drift
268- if dither:
269- ycur = ypad[P:P + M]
270- dlo = np.maximum(-0.5, z - 0.5 - ycur)
271- dhi = np.minimum(0.5, z + 0.5 - ycur)
272- d = dlo + np.maximum(dhi - dlo, 0.0) * rng.random(M)
273- set_boxes()
274- for idx, rows in zip(classes, rowmats):
275- r = ypad[rows] - np.outer(x[idx], h)
276- with np.errstate(divide='ignore', invalid='ignore'):
277- b1 = (lo[rows] - r) / h[None, :]
278- b2 = (hi[rows] - r) / h[None, :]
279- xlo = np.where(h[None, :] > 0, b1, b2)
280- xhi = np.where(h[None, :] > 0, b2, b1)
281- xlo[:, ~nonzero] = -np.inf
282- xhi[:, ~nonzero] = np.inf
283- xlo = xlo.max(axis=1)
284- xhi = xhi.min(axis=1)
285- xnew = trunc_std_normal(xlo / sigma, xhi / sigma, rng) * sigma
286- delta = xnew - x[idx]
287- x[idx] = xnew
288- ypad[rows] += delta[:, None] * h[None, :]
289-
290- c = float(hr[:-1] @ x[M:]) if L > 1 else 0.0
291- js, p = next_pmf(c, s0, dither)
292- for j, pj in zip(js, p):
293- if pj > 1e-15:
294- pmf[int(j)] = pmf.get(int(j), 0.0) + pj
295-
296- total = sum(pmf.values())
297- return -sum((p / total) * math.log2(p / total) for p in pmf.values() if p > 0)
298-
299-
300-# --------------------------------------------------------------------- main
301-
302-def main():
303- ap = argparse.ArgumentParser(
304- description='Monte-Carlo estimate of the true entropy rate, for '
305- 'checking the analytic rate R shown in the app.')
306- ap.add_argument('--sigma', type=float, required=True, help='input std, in quantization steps')
307- ap.add_argument('--filter', required=True,
308- choices=['none', 'moving-average', 'lowpass', 'bandpass', 'first-difference'])
309- ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
310- ap.add_argument('--high', type=float, help='lowpass cutoff / bandpass high edge, Hz')
311- ap.add_argument('--taps', type=int, default=101, help='windowed-sinc kernel length')
312- ap.add_argument('--width', type=int, default=8, help='moving-average width')
313- ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
314- ap.add_argument('--dither', action='store_true')
315- ap.add_argument('--past', type=int, help='conditioning window M (default max(512, 4·taps))')
316- ap.add_argument('--pasts', type=int, default=24, help='independent pasts to average')
317- ap.add_argument('--sweeps', type=int, default=600, help='Gibbs sweeps per past')
318- ap.add_argument('--seed', type=int, default=0)
319- args = ap.parse_args()
320-
321- if args.filter == 'bandpass' and (args.low is None or args.high is None):
322- ap.error('bandpass needs --low and --high')
323- if args.filter == 'lowpass' and args.high is None:
324- ap.error('lowpass needs --high')
325-
326- kernel = design_kernel(args)
327- L = len(kernel)
328- M = args.past if args.past is not None else max(512, 4 * L)
329-
330- R_spec, R_samp = formula_rates(kernel, args.sigma, args.dither)
331- R = min(R_spec, R_samp)
332- print(f'model: sigma={args.sigma} filter={args.filter} L={L} dither={args.dither}')
333- ratio = f' (ratio vs int16: {16 / R:.3f}x)' if R > 0 else ''
334- print(f'formula R = min(spec {R_spec:.4f}, samp {R_samp:.4f}) = {R:.4f} bits/sample{ratio}')
335- print(f'MC: {args.pasts} pasts x {args.sweeps} sweeps, conditioning on M={M} samples')
336-
337- rng = np.random.default_rng(args.seed)
338- Hs = []
339- for i in range(args.pasts):
340- Hs.append(conditional_entropy_of_one_past(
341- kernel, args.sigma, args.dither, M, args.sweeps, rng))
342- mean = float(np.mean(Hs))
343- se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs))) if len(Hs) > 1 else float('nan')
344- print(f' past {i + 1:3d}/{args.pasts}: H = {Hs[-1]:.4f} running mean {mean:.4f} +/- {se:.4f}')
345-
346- mean = float(np.mean(Hs))
347- se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs)))
348- print(f'\nMC entropy rate: {mean:.4f} +/- {se:.4f} bits/sample'
349- f' (ratio vs int16: {16 / mean:.3f}x)')
350- print(f'formula R: {R:.4f} bits/sample (formula - MC = {R - mean:+.4f}; '
351- f'spec {R_spec:.4f}, samp {R_samp:.4f})')
352- print('note: the MC value estimates H(z_next | M past samples), an upper bound '
353- 'on the rate that tightens as --past grows.')
354-
355-
356-if __name__ == '__main__':
357- main()
src/App.tsxmodified+33−41View file
@@ -4,9 +4,8 @@ import Controls from './components/Controls'
44 import FilterViz from './components/FilterViz'
55 import ScrollingView from './components/ScrollingView'
66 import CompressionChart from './components/CompressionChart'
7-import MathSection from './components/MathSection'
7+import MethodNote from './components/MethodNote'
88 import { DEFAULT_SPEC, clampSpec, designKernel, kernelNorm } from './model/filters'
9-import { theoreticalRateBits } from './model/theory'
109 import { LATENT_SEED } from './model/latent'
1110 import { DEFAULT_LPC_ORDER, LPC_ORDERS } from './compress/codecs'
1211 import type { BoundResult, CodecResult } from './compress/codecs'
@@ -27,7 +26,6 @@ interface CompressionState {
2726 function useCompression(
2827 kernel: Float64Array,
2928 sigma: number,
30- dither: boolean,
3129 lpcOrder: number,
3230 blockSize: number,
3331 ): CompressionState {
@@ -70,7 +68,6 @@ function useCompression(
7068 id,
7169 kernel,
7270 sigma,
73- dither,
7471 blockSize,
7572 lpcOrder,
7673 // The same seed the signal view draws from, so the block really is
@@ -80,14 +77,22 @@ function useCompression(
8077 workerRef.current?.postMessage(request)
8178 }, 250)
8279 return () => clearTimeout(timer)
83- }, [kernel, sigma, dither, lpcOrder, blockSize])
80+ }, [kernel, sigma, lpcOrder, blockSize])
8481
8582 return state
8683 }
8784
88-/** The terminal command for scripts/true_rate.py at the current settings. */
89-function mcCommand(sigma: number, spec: ReturnType<typeof clampSpec>, rate: number, dither: boolean): string {
90- const parts = ['python scripts/true_rate.py', `--sigma ${sigma}`]
85+/**
86+ * The terminal command that estimates the reference rate R at the current
87+ * settings, using the unbiased Monte-Carlo estimator from the companion
88+ * timeseries-entropy package.
89+ */
90+function mcCommand(sigma: number, spec: ReturnType<typeof clampSpec>, rate: number): string {
91+ const parts = [
92+ 'uvx --from git+https://github.com/concept-collection/timeseries-entropy',
93+ 'timeseries-entropy',
94+ `--sigma ${sigma}`,
95+ ]
9196 switch (spec.family) {
9297 case 'none':
9398 parts.push('--filter none')
@@ -111,7 +116,6 @@ function mcCommand(sigma: number, spec: ReturnType<typeof clampSpec>, rate: numb
111116 parts.push('--filter first-difference')
112117 break
113118 }
114- if (dither) parts.push('--dither')
115119 return parts.join(' ')
116120 }
117121
@@ -119,26 +123,21 @@ export default function App() {
119123 const [sigma, setSigma] = useState(5)
120124 const [sampleRateHz, setSampleRateHz] = useState(30000)
121125 const [spec, setSpec] = useState(DEFAULT_SPEC)
122- const [dither, setDither] = useState(false)
123126 const [lpcOrder, setLpcOrder] = useState(DEFAULT_LPC_ORDER)
124127 const [blockSize, setBlockSize] = useState(DEFAULT_BLOCK_SIZE)
125128
126129 const kernel = useMemo(() => designKernel(spec, sampleRateHz), [spec, sampleRateHz])
127- const sigmaY = useMemo(() => {
128- const filtered = sigma * kernelNorm(kernel)
129- return dither ? Math.sqrt(filtered * filtered + 1 / 12) : filtered
130- }, [kernel, sigma, dither])
131- const theoryBits = useMemo(() => theoreticalRateBits(kernel, sigma, dither), [kernel, sigma, dither])
132- const compression = useCompression(kernel, sigma, dither, lpcOrder, blockSize)
130+ const sigmaY = useMemo(() => sigma * kernelNorm(kernel), [kernel, sigma])
131+ const compression = useCompression(kernel, sigma, lpcOrder, blockSize)
133132
134133 return (
135134 <div className="app">
136135 <header className="app-header">
137136 <h1>Time-series compressibility</h1>
138137 <p>
139- Gaussian noise → FIR filter → optional dither → round to integers. How well can the
140- integer stream be losslessly compressed, and does the spectral entropy-rate formula
141- predict the limit?
138+ Gaussian noise → FIR filter → round to integers. How well can the integer stream be
139+ losslessly compressed, and how close do practical codecs get to the entropy rate of
140+ the process?
142141 </p>
143142 </header>
144143
@@ -157,8 +156,6 @@ export default function App() {
157156 }}
158157 spec={spec}
159158 setSpec={setSpec}
160- dither={dither}
161- setDither={setDither}
162159 />
163160 </section>
164161
@@ -178,16 +175,14 @@ export default function App() {
178175 <small>steps</small>
179176 </span>
180177 </div>
178+ {/* Placeholder: R will be estimated in the browser by the unbiased
179+ estimator; until then the command below produces it locally. */}
181180 <div className="stat">
182- <span className="label">theoretical rate R</span>
181+ <span className="label">reference rate R</span>
183182 <span className="value">
184- {theoryBits.toFixed(2)} <small>bits/sample</small>
183+ — <small>bits/sample</small>
185184 </span>
186185 </div>
187- <div className="stat">
188- <span className="label">implied best ratio</span>
189- <span className="value">{theoryBits > 0 ? `${(16 / theoryBits).toFixed(2)}×` : '—'}</span>
190- </div>
191186 </div>
192187 {/* Settings of the measurement, not of the model — so they live with
193188 the chart they change rather than in the model bar. */}
@@ -219,7 +214,6 @@ export default function App() {
219214 <CompressionChart
220215 results={compression.results}
221216 bounds={compression.bounds}
222- theoryBits={theoryBits}
223217 computing={compression.computing}
224218 />
225219 )}
@@ -229,36 +223,34 @@ export default function App() {
229223 coefficients). Baseline is raw int16 (16 bits/sample). The hollow bar in each group is
230224 that group's entropy limit — the order-0 entropy of the stream being coded, which no
231225 per-sample entropy coder can beat and ANS falls short of by its symbol table plus its
232- own arithmetic loss. The dashed line is the theoretical rate R from the math section —
233- a spectral estimate with the rounding charged at its full variance, capped by the
234- exact one-sample entropy when the whole process sits below the quantization step
235- (see the S(f) = 1 threshold on the response plot); it is least certain at the
236- crossover between those two regimes.
226+ own arithmetic loss. The reference rate R — the entropy rate of the process itself,
227+ the limit no lossless method can beat — is not yet computed in the browser; the
228+ command below estimates it locally (see the method section at the bottom).
237229 </p>
238230 <CopyableCommand
239- label="check R against a Monte-Carlo ground truth (runs locally, ~2 min):"
240- command={mcCommand(sigma, spec, sampleRateHz, dither)}
231+ label="reference rate R by unbiased Monte-Carlo (runs locally):"
232+ command={mcCommand(sigma, spec, sampleRateHz)}
241233 />
242234 </section>
243235
244236 <section className="card">
245237 <h2>Quantized signal z</h2>
246- <ScrollingView kernel={kernel} sigma={sigma} dither={dither} sigmaY={sigmaY} />
238+ <ScrollingView kernel={kernel} sigma={sigma} sigmaY={sigmaY} />
247239 <p className="card-note">
248240 A window of samples from the model, drawn from a fixed latent noise sequence — changing
249- σ, the filter, or dither transforms the same underlying data, so the trace morphs
250- rather than resampling. Press play to advance through the sequence.
241+ σ or the filter transforms the same underlying data, so the trace morphs rather than
242+ resampling. Press play to advance through the sequence.
251243 </p>
252244 </section>
253245
254246 <section className="card">
255247 <h2>Filter</h2>
256- <FilterViz kernel={kernel} sampleRateHz={sampleRateHz} sigma={sigma} />
248+ <FilterViz kernel={kernel} sampleRateHz={sampleRateHz} />
257249 </section>
258250
259251 <section className="card">
260- <h2>The math</h2>
261- <MathSection />
252+ <h2>The reference rate</h2>
253+ <MethodNote />
262254 </section>
263255 </div>
264256 )
src/app.cssmodified+0−47View file
@@ -165,26 +165,6 @@ body {
165165 width: 170px;
166166 }
167167
168-.control-toggle {
169- flex-direction: row;
170- align-items: center;
171- gap: 8px;
172- min-width: 0;
173- padding-bottom: 6px;
174-}
175-
176-.control-toggle input {
177- accent-color: var(--series-1);
178- width: 16px;
179- height: 16px;
180- margin: 0;
181-}
182-
183-.control-toggle label {
184- font-size: 13px;
185- color: var(--ink);
186-}
187-
188168 /* ---- plots ---- */
189169
190170 .filter-panels {
@@ -507,33 +487,6 @@ body {
507487 max-width: 76ch;
508488 }
509489
510-.math-section h3 {
511- font-size: 12px;
512- font-weight: 650;
513- text-transform: uppercase;
514- letter-spacing: 0.05em;
515- color: var(--ink-2);
516- margin: 26px 0 10px;
517-}
518-
519490 .math-section .katex-display {
520491 margin: 14px 0;
521492 }
522-
523-.defs {
524- display: grid;
525- grid-template-columns: max-content 1fr;
526- gap: 8px 18px;
527- margin: 18px 0 0;
528- max-width: 76ch;
529-}
530-
531-.defs dt {
532- text-align: right;
533- color: var(--ink);
534-}
535-
536-.defs dd {
537- margin: 0;
538- color: var(--ink-2);
539-}
src/components/CompressionChart.tsxmodified+2−41View file
@@ -89,7 +89,6 @@ function buildRows(results: CodecResult[], bounds: BoundResult[]): Row[] {
8989 export default function CompressionChart(props: {
9090 results: CodecResult[]
9191 bounds: BoundResult[]
92- theoryBits: number
9392 computing: boolean
9493 }) {
9594 const ref = useRef<HTMLDivElement>(null)
@@ -98,7 +97,7 @@ export default function CompressionChart(props: {
9897 const [tip, setTip] = useState<Tip | null>(null)
9998 const [hovered, setHovered] = useState<string | null>(null)
10099
101- const { results, bounds, theoryBits } = props
100+ const { results, bounds } = props
102101 if (results.length === 0) {
103102 return <p className="card-note">Computing compression on the first block…</p>
104103 }
@@ -106,12 +105,8 @@ export default function CompressionChart(props: {
106105 const rows = buildRows(results, bounds)
107106 const value = (r: { bitsPerSample: number; ratio: number }) =>
108107 metric === 'bits' ? r.bitsPerSample : r.ratio
109- const theoryValue = metric === 'bits' ? theoryBits : theoryBits > 0 ? 16 / theoryBits : NaN
110- const theoryVisible = Number.isFinite(theoryValue) && theoryValue > 0
111108 const xMax =
112- metric === 'bits'
113- ? Math.max(16, ...rows.map(value), theoryVisible ? theoryValue : 0) * 1.02
114- : Math.max(...rows.map(value), theoryVisible ? theoryValue : 0) * 1.1
109+ metric === 'bits' ? Math.max(16, ...rows.map(value)) * 1.02 : Math.max(...rows.map(value)) * 1.1
115110
116111 const plotW = width - LABEL_W - RIGHT_PAD
117112 const height = AXIS_H + GROUPS.length * (GROUP_H + ROWS_PER_GROUP * ROW_H) + 6
@@ -129,10 +124,6 @@ export default function CompressionChart(props: {
129124 setTip({ x: e.clientX - box.left, y: e.clientY - box.top, row })
130125 }
131126
132- const theoryX = theoryVisible ? xOf(theoryValue) : 0
133- const theoryLabel =
134- metric === 'bits' ? `R = ${theoryBits.toFixed(2)}` : `R ⇒ ${(16 / theoryBits).toFixed(2)}×`
135-
136127 return (
137128 <div>
138129 <div className="chart-header">
@@ -233,28 +224,6 @@ export default function CompressionChart(props: {
233224 </g>
234225 )
235226 })}
236- {theoryVisible && (
237- <g>
238- <line
239- x1={theoryX}
240- x2={theoryX}
241- y1={AXIS_H - 2}
242- y2={height - 4}
243- stroke="var(--ink-2)"
244- strokeWidth={1.5}
245- strokeDasharray="5 4"
246- />
247- <text
248- x={theoryX + (theoryX > width - 150 ? -6 : 6)}
249- y={AXIS_H + 10}
250- textAnchor={theoryX > width - 150 ? 'end' : 'start'}
251- className="bar-value"
252- fill="var(--ink)"
253- >
254- {theoryLabel}
255- </text>
256- </g>
257- )}
258227 </svg>
259228 {tip && (
260229 <div className="viz-tooltip" style={{ left: tip.x + 14, top: tip.y - 8 }}>
@@ -291,14 +260,6 @@ export default function CompressionChart(props: {
291260 <td>{r.ratio.toFixed(3)}</td>
292261 </tr>
293262 ))}
294- {theoryBits > 0 && (
295- <tr>
296- <td>theory: rate R</td>
297- <td>—</td>
298- <td>{theoryBits.toFixed(3)}</td>
299- <td>{(16 / theoryBits).toFixed(3)}</td>
300- </tr>
301- )}
302263 </tbody>
303264 </table>
304265 </details>
src/components/Controls.tsxmodified+0−11View file
@@ -17,8 +17,6 @@ export interface ControlsProps {
1717 setSampleRateHz: (v: number) => void
1818 spec: FilterSpec
1919 setSpec: (v: FilterSpec) => void
20- dither: boolean
21- setDither: (v: boolean) => void
2220 }
2321
2422 function formatHz(hz: number): string {
@@ -125,15 +123,6 @@ export default function Controls(p: ControlsProps) {
125123 ))}
126124 </select>
127125 </div>
128- <div className="control control-toggle">
129- <input
130- id="dither"
131- type="checkbox"
132- checked={p.dither}
133- onChange={e => p.setDither(e.target.checked)}
134- />
135- <label htmlFor="dither">dither (uniform ±½ before rounding)</label>
136- </div>
137126 </div>
138127 )
139128 }
src/components/FilterViz.tsxmodified+4−33View file
@@ -118,16 +118,8 @@ function KernelPanel({ kernel, sampleRateHz }: { kernel: Float64Array; sampleRat
118118 )
119119 }
120120
121-/** |H(f)| in dB up to Nyquist, with the S(f) = 1 step² threshold. */
122-function ResponsePanel({
123- kernel,
124- sampleRateHz,
125- sigma,
126-}: {
127- kernel: Float64Array
128- sampleRateHz: number
129- sigma: number
130-}) {
121+/** |H(f)| in dB up to Nyquist. */
122+function ResponsePanel({ kernel, sampleRateHz }: { kernel: Float64Array; sampleRateHz: number }) {
131123 const ref = useRef<HTMLDivElement>(null)
132124 const width = useWidth(ref)
133125 const [tip, setTip] = useState<Tip | null>(null)
@@ -146,11 +138,6 @@ function ResponsePanel({
146138
147139 const path = db.map((v, i) => `${i ? 'L' : 'M'}${xOf(i).toFixed(1)},${yOf(v).toFixed(1)}`).join('')
148140
149- // S(f) = σ²|H|² = 1 (one step² of spectral power) sits at |H| = 1/σ. Below
150- // this line the high-resolution formula is on thin ice.
151- const stepDb = -20 * Math.log10(sigma)
152- const showStep = stepDb < dbMax && stepDb > DB_FLOOR
153-
154141 const onMove = (e: React.PointerEvent<SVGSVGElement>) => {
155142 const box = e.currentTarget.getBoundingClientRect()
156143 const px = e.clientX - box.left
@@ -189,22 +176,6 @@ function ResponsePanel({
189176 <text x={width - MARGIN.right} y={HEIGHT - 8} textAnchor="end" className="axis-tick">
190177 Hz
191178 </text>
192- {showStep && (
193- <g>
194- <line
195- x1={MARGIN.left}
196- x2={width - MARGIN.right}
197- y1={yOf(stepDb)}
198- y2={yOf(stepDb)}
199- stroke="var(--muted)"
200- strokeWidth={1}
201- strokeDasharray="4 3"
202- />
203- <text x={width - MARGIN.right} y={yOf(stepDb) - 4} textAnchor="end" className="axis-tick">
204- S(f) = 1 step²
205- </text>
206- </g>
207- )}
208179 <path d={path} fill="none" stroke="var(--series-1)" strokeWidth={2} strokeLinejoin="round" strokeLinecap="round" />
209180 </svg>
210181 <Tooltip tip={tip} />
@@ -212,11 +183,11 @@ function ResponsePanel({
212183 )
213184 }
214185
215-export default function FilterViz(props: { kernel: Float64Array; sampleRateHz: number; sigma: number }) {
186+export default function FilterViz(props: { kernel: Float64Array; sampleRateHz: number }) {
216187 return (
217188 <div className="filter-panels">
218189 <KernelPanel kernel={props.kernel} sampleRateHz={props.sampleRateHz} />
219- <ResponsePanel kernel={props.kernel} sampleRateHz={props.sampleRateHz} sigma={props.sigma} />
190+ <ResponsePanel kernel={props.kernel} sampleRateHz={props.sampleRateHz} />
220191 </div>
221192 )
222193 }
src/components/MathSection.tsxdeleted+0−127View file
@@ -1,127 +0,0 @@
1-import { useMemo } from 'react'
2-import katex from 'katex'
3-
4-function Tex({ tex, display }: { tex: string; display?: boolean }) {
5- const html = useMemo(
6- () => katex.renderToString(tex, { displayMode: !!display, throwOnError: false }),
7- [tex, display],
8- )
9- return <span dangerouslySetInnerHTML={{ __html: html }} />
10-}
11-
12-/** One term of the formula: the symbol, then what it is. */
13-function Def({ tex, children }: { tex: string; children: React.ReactNode }) {
14- return (
15- <>
16- <dt>
17- <Tex tex={tex} />
18- </dt>
19- <dd>{children}</dd>
20- </>
21- )
22-}
23-
24-/**
25- * The theoretical rate R exactly as `model/theory.ts` computes it, with every
26- * symbol defined, followed by a sketch of the derivation: the dither identity,
27- * the per-mode count with the noise at full variance, and the subadditivity
28- * ceiling that takes over below threshold.
29- */
30-export default function MathSection() {
31- return (
32- <div className="math-section">
33- <p>
34- The dashed line on the compression chart is R, the predicted bits per sample — the
35- smaller of a spectral estimate and a rigorous one-sample ceiling:
36- </p>
37-
38- <Tex display tex="R \;=\; \min\big(R_{\mathrm{spec}},\, R_{\mathrm{samp}}\big)" />
39- <Tex display tex="R_{\mathrm{spec}} \;=\; \int_0^1 \tfrac{1}{2}\log_2\!\big(2\pi e\,(S(f) + \nu)\big)\, df, \qquad S(f) \;=\; \sigma^2\,\big|H(f)\big|^2" />
40- <Tex display tex="R_{\mathrm{samp}} \;=\; H\big(\operatorname{round}(\mathcal N(0, v) \,[+\, U(-\tfrac12,\tfrac12)\ \text{with dither}])\big), \qquad v \;=\; \sigma^2 \textstyle\sum_m h_m^2" />
41-
42- <dl className="defs">
43- <Def tex="\sigma">
44- standard deviation of the i.i.d. Gaussian input, in quantization steps (the step is the
45- unit, so rounding is to the nearest integer)
46- </Def>
47- <Def tex="h_0,\dots,h_{L-1}">
48- the FIR kernel the input is convolved with — the taps drawn in the kernel plot, L of
49- them
50- </Def>
51- <Def tex="H(f)">
52- the kernel's frequency response, the quantity plotted in dB as |H(f)|; f is in cycles
53- per sample, symmetric about ½ (Nyquist), and the plots label the same axis in Hz
54- </Def>
55- <Def tex="S(f)">
56- power spectrum of the filtered signal alone, in steps² per unit frequency — the dither
57- is <em>not</em> folded in here
58- </Def>
59- <Def tex="\nu">
60- variance charged to the rounding: 1/12 for the roundoff, 1/6 with dither (the dither's
61- own 1/12 is stored in the integers and adds)
62- </Def>
63- <Def tex="v">
64- variance of a single output sample, the integral of S(f)
65- </Def>
66- <Def tex="R">
67- bits per sample; the compression ratio the chart marks is 16/R, against 16-bit integer
68- storage
69- </Def>
70- </dl>
71-
72- <h3>The spectral branch</h3>
73- <p>
74- The <em>dither identity</em> starts it off: for z = round(y) and u an independent
75- uniform on [-½, ½)<sup>N</sup>, the discrete entropy of z equals the differential
76- entropy of z + u, exactly. When the process is live on the unit-cell scale (v ≳ ¼),
77- z + u has nearly the law of y + u, so R is the entropy rate of the signal plus a white
78- unit-cell noise — the Zamir–Feder universal-quantization rate; with physical dither the
79- smoothing noise is d + u and ν doubles to 1/6.
80- </p>
81- <p>
82- Counting that entropy per Fourier mode, the signal modes are independent Gaussians of
83- variance S(f), and each mode of the i.i.d. cube noise mixes all N samples'
84- contributions — so the central limit theorem Gaussianizes it, and it enters at its{' '}
85- <em>full variance</em> ν. It does not enter at the entropy power 1/(2πe) that a scalar
86- quantizer aligned with the mode would charge: the quantization lattice lives in the
87- sample basis, and only for the trivial kernel do modes and quantizers align. (An earlier
88- version of this app charged entropy power — additive constant 1 instead of 2πe·ν ≈ 1.42
89- — and systematically underestimated the measured rate by up to ~0.23 bits/sample. Its
90- "exact per-mode" refinement was worse still: it modeled the wrong physics more
91- faithfully.) One consequence worth naming: a dead band inside a live process contributes
92- ½log₂(2πe/12) ≈ 0.25 bits per mode, not zero.
93- </p>
94-
95- <h3>The sub-threshold ceiling</h3>
96- <p>
97- When v ≪ 1 nearly every sample rounds to zero and the true rate collapses
98- exponentially, while R<sub>spec</sub> bottoms out at ½log₂(2πe ν) > 0. Subadditivity
99- rescues the estimate rigorously: H(z) ≤ Σ<sub>n</sub> H(z<sub>n</sub>), and each stored
100- sample is exactly round(N(0, v)) — plus the uniform dither first when it is on — so
101- R<sub>samp</sub> is a true upper bound on the rate with exactly the right collapse. The
102- min selects it precisely where the spectral branch fails.
103- </p>
104-
105- <h3>Checks and accuracy</h3>
106- <p>
107- For the identity kernel the two branches agree with the exact i.i.d. entropy at every σ
108- (both carry the Fisher correction log₂e/(24σ²) at large σ; below one step the min
109- switches to the exact R<sub>samp</sub>). At high SNR, R<sub>spec</sub> →
110- ½log₂(2πe σ²) + ∫log₂|H| df — the Kolmogorov formula. What the spectral branch ignores
111- is the cross-mode dependence of the cube noise, at most ½log₂(2πe/12) ≈ 0.2546
112- bits/sample and recoverable only when nearly the whole spectrum is noise-dominated;
113- Monte-Carlo puts the estimate within ~0.01–0.02 bits/sample for v ≳ 0.25, with the
114- worst observed error ~+0.03 near the crossover between branches, slightly positive
115- everywhere — as befits a formula whose sample branch is a genuine bound.
116- </p>
117-
118- <p className="card-note">
119- The integral is evaluated by the midpoint rule on 8192 points over [0, ½] (symmetry
120- supplies the other half). R<sub>samp</sub> sums the exact bin probabilities of the
121- rounded Gaussian — integrated against the triangular dither-overlap window when dither
122- is on. The Monte-Carlo command under the chart estimates the true entropy rate of the
123- same process, for checking R where the approximations are in doubt.
124- </p>
125- </div>
126- )
127-}
src/components/MethodNote.tsxadded+41−0View file
@@ -0,0 +1,41 @@
1+import { useMemo } from 'react'
2+import katex from 'katex'
3+
4+function Tex({ tex, display }: { tex: string; display?: boolean }) {
5+ const html = useMemo(
6+ () => katex.renderToString(tex, { displayMode: !!display, throwOnError: false }),
7+ [tex, display],
8+ )
9+ return <span dangerouslySetInnerHTML={{ __html: html }} />
10+}
11+
12+/** Where the reference rate R comes from, in brief. */
13+export default function MethodNote() {
14+ return (
15+ <div className="math-section">
16+ <p>
17+ R is the entropy rate of the quantized process z — the bits per sample that no lossless
18+ code can beat. It has no usable closed form here, so it is estimated by an unbiased
19+ Monte-Carlo method from the companion{' '}
20+ <a href="https://github.com/concept-collection/timeseries-entropy">timeseries-entropy</a>{' '}
21+ package. The estimand is the conditional entropy of the next sample given a long past,
22+ </p>
23+ <Tex display tex="H\big(z_{M+1} \,\big|\, z_1, \dots, z_M\big) \;\searrow\; R \qquad (M \to \infty)," />
24+ <p>
25+ which reaches R once M exceeds the memory of the process. A past is drawn from the
26+ model, and the latent Gaussian input is Gibbs-sampled under the rounding constraints —
27+ the latents that generated the past are an exact draw from the conditional, so the chain
28+ starts in stationarity with no burn-in bias — emitting exact draws of z<sub>M+1</sub>.
29+ Rhee–Glynn randomized telescoping with antithetic half-block corrections then turns the
30+ plug-in entropies of that chain into an estimate whose expectation is exactly the
31+ conditional entropy, despite the finite-sample bias of every plug-in estimate and the
32+ autocorrelation of the Gibbs draws. Averaging over independent pasts gives R with an
33+ honest standard error.
34+ </p>
35+ <p className="card-note">
36+ Running the estimator in the browser is planned; until then, the command under the
37+ compression chart runs it locally at the current settings.
38+ </p>
39+ </div>
40+ )
41+}
src/components/ScrollingView.tsxmodified+2−3View file
@@ -23,7 +23,6 @@ function niceStep(span: number): number {
2323 export default function ScrollingView(props: {
2424 kernel: Float64Array
2525 sigma: number
26- dither: boolean
2726 /** Predicted std of the quantized signal, for a stable y-scale. */
2827 sigmaY: number
2928 }) {
@@ -84,7 +83,7 @@ export default function ScrollingView(props: {
8483 // The window ends at posRef and never reaches before index 0, so the
8584 // first thing shown is the start of the compression block.
8685 if (posRef.current < visible) posRef.current = visible
87- const win = latent.window(posRef.current - visible, visible, props.kernel, props.sigma, props.dither)
86+ const win = latent.window(posRef.current - visible, visible, props.kernel, props.sigma)
8887
8988 const surface = styles.getPropertyValue('--surface')
9089 ctx.fillStyle = surface
@@ -140,7 +139,7 @@ export default function ScrollingView(props: {
140139 cancelAnimationFrame(raf)
141140 scheme.removeEventListener('change', refreshStyles)
142141 }
143- }, [props.kernel, props.sigma, props.dither, props.sigmaY])
142+ }, [props.kernel, props.sigma, props.sigmaY])
144143
145144 return (
146145 <div className="scroll-wrap">
src/model/filters.tsmodified+4−4View file
@@ -2,10 +2,10 @@
22 * FIR filter presets and their frequency response.
33 *
44 * The filter is always realized as an explicit convolution kernel, so the
5- * pipeline is exactly x → h*x → (+dither) → round, and the theory can work
6- * from |H(f)| of the same taps the data actually went through. Cutoffs are
7- * given in Hz against a user-set sample rate; internally everything is in
8- * normalized frequency (cycles/sample, Nyquist = 0.5).
5+ * pipeline is exactly x → h*x → round, and the |H(f)| plotted is the response
6+ * of the same taps the data actually went through. Cutoffs are given in Hz
7+ * against a user-set sample rate; internally everything is in normalized
8+ * frequency (cycles/sample, Nyquist = 0.5).
99 */
1010
1111 export type FilterFamily = 'none' | 'movingAverage' | 'lowpass' | 'bandpass' | 'firstDifference'
src/model/latent.tsmodified+8−19View file
@@ -1,10 +1,10 @@
11 /**
2- * A fixed latent randomness underlying everything: standard normals (and
3- * dither uniforms) indexed by absolute sample position. The pipeline
4- * x → h*x → (+dither) → round is evaluated on demand against these indices,
5- * so changing σ, the filter, or dither transforms the *same* underlying data
6- * — the display morphs smoothly instead of resampling — and the compression
7- * block (indices 0…N) shares its randomness with the displayed window.
2+ * A fixed latent randomness underlying everything: standard normals indexed
3+ * by absolute sample position. The pipeline x → h*x → round is evaluated on
4+ * demand against these indices, so changing σ or the filter transforms the
5+ * *same* underlying data — the display morphs smoothly instead of resampling
6+ * — and the compression block (indices 0…N) shares its randomness with the
7+ * displayed window.
88 */
99 import { GaussianStream } from './random'
1010
@@ -12,19 +12,15 @@ export const LATENT_SEED = 20260729
1212
1313 export class LatentSource {
1414 private xs: number[] = []
15- private ds: number[] = []
1615 private xStream: GaussianStream
17- private dStream: GaussianStream
1816
1917 constructor(seed: number) {
2018 this.xStream = new GaussianStream(seed)
21- this.dStream = new GaussianStream((seed ^ 0x9e3779b9) >>> 0)
2219 }
2320
2421 private ensure(n: number) {
2522 while (this.xs.length <= n) {
2623 this.xs.push(this.xStream.normal())
27- this.ds.push(this.dStream.uniformCentered())
2824 }
2925 }
3026
@@ -34,17 +30,11 @@ export class LatentSource {
3430 * does not shift features along the time axis. Latent indices before 0 read
3531 * as zero input.
3632 */
37- window(
38- start: number,
39- count: number,
40- kernel: Float64Array,
41- sigma: number,
42- dither: boolean,
43- ): Int16Array {
33+ window(start: number, count: number, kernel: Float64Array, sigma: number): Int16Array {
4434 const L = kernel.length
4535 const mid = (L - 1) >> 1
4636 this.ensure(start + count - 1 + mid)
47- const { xs, ds } = this
37+ const { xs } = this
4838 const out = new Int16Array(count)
4939 for (let j = 0; j < count; j++) {
5040 const n = start + j
@@ -54,7 +44,6 @@ export class LatentSource {
5444 if (idx >= 0) y += kernel[k] * xs[idx]
5545 }
5646 y *= sigma
57- if (dither) y += ds[n]
5847 out[j] = Math.max(-32768, Math.min(32767, Math.round(y)))
5948 }
6049 return out
src/model/random.tsmodified+0−5View file
@@ -21,11 +21,6 @@ export class GaussianStream {
2121 return ((t ^ (t >>> 14)) >>> 0) / 4294967296
2222 }
2323
24- /** Uniform on [-1/2, 1/2), for dither. */
25- uniformCentered(): number {
26- return this.uniform() - 0.5
27- }
28-
2924 /** Standard normal. */
3025 normal(): number {
3126 if (this.spare !== null) {
src/model/theory.tsdeleted+0−137View file
@@ -1,137 +0,0 @@
1-/**
2- * The theoretical bits/sample for the quantized filtered-Gaussian process:
3- *
4- * R̂ = min(R_spec, R_samp)
5- *
6- * R_spec = ∫₀¹ ½ log₂(2πe (S(f) + ν)) df, S(f) = σ²|H(f)|²
7- * ν = 1/12 (1/6 with dither)
8- * R_samp = exact entropy of one stored sample, N(0, v) rounded
9- * (+ uniform dither first when it is on), v = σ² Σ h²
10- *
11- * R_spec is the Zamir–Feder rate of the dithered quantizer, counted per
12- * Fourier mode: the signal modes are independent Gaussians of variance S(f),
13- * and each mode of the i.i.d. roundoff-plus-dither noise mixes all N samples'
14- * contributions, so it is Gaussianized by the CLT and enters at its full
15- * variance ν — not at the entropy power 1/(2πe) an aligned scalar quantizer
16- * would charge (the quantization lattice lives in the sample basis, not the
17- * Fourier basis; charging entropy power is what made earlier versions of this
18- * estimate underestimate the rate). A consequence worth naming: a dead band
19- * inside a live process contributes ½log₂(2πe ν) ≈ 0.25 bits per mode, not
20- * zero. What R_spec ignores is the cross-mode dependence of the cube noise —
21- * at most ½log₂(2πe/12) ≈ 0.2546 bits/sample, in practice ≲ 0.02 unless
22- * nearly the whole spectrum is noise-dominated.
23- *
24- * That failure mode is exactly the globally sub-threshold process, and there
25- * subadditivity gives a rigorous ceiling with the right collapse: H(z) ≤
26- * Σ H(zₙ) = N·R_samp, the marginal entropy of a single output sample. The
27- * min selects it precisely where the spectral branch fails. Monte Carlo puts
28- * R̂ within ~0.01–0.02 bits/sample for v ≳ 0.25, worst ~+0.03 near the
29- * crossover; testing that against LPC+ANS is the app's point.
30- */
31-
32-const TWO_PI_E = 2 * Math.PI * Math.E
33-const INTEGRATION_POINTS = 8192
34-
35-/**
36- * Exact entropy (bits) of round(N(0, s²)) on the unit lattice. Per-bin
37- * probabilities by Simpson integration of the density, so no erf is needed
38- * and the tail keeps relative accuracy.
39- */
40-export function quantizedGaussianEntropy(s: number): number {
41- if (s <= 0.02) return 0
42- // The discrete entropy approaches the differential entropy ½log₂(2πe s²)
43- // from above like log₂e/(24 s²) — the Δ²/24 Fisher-information correction,
44- // with the next term O(1/s⁴). At s ≥ 6 the corrected asymptote is within
45- // 2·10⁻⁶ bits, so the sum is only ever taken over a handful of bins.
46- if (s >= 6) return 0.5 * Math.log2(TWO_PI_E * s * s) + Math.LOG2E / (24 * s * s)
47- const zMax = Math.ceil(8 * s + 4)
48- // Enough points that a bin spans a few per standard deviation even when the
49- // bin is wide compared to the distribution.
50- const m = Math.min(401, Math.max(9, 2 * Math.ceil(3 / s) + 9)) | 1
51- const h = 1 / (m - 1)
52- const density = (u: number) => Math.exp((-u * u) / (2 * s * s)) / (Math.sqrt(2 * Math.PI) * s)
53- let sumH = 0
54- let total = 0
55- for (let z = -zMax; z <= zMax; z++) {
56- let acc = density(z - 0.5) + density(z + 0.5)
57- for (let i = 1; i < m - 1; i++) acc += (i % 2 === 1 ? 4 : 2) * density(z - 0.5 + i * h)
58- const p = (acc * h) / 3
59- if (p > 0) {
60- sumH -= p * Math.log2(p)
61- total += p
62- }
63- }
64- // Renormalize away the residual quadrature/truncation mass.
65- return sumH / total + Math.log2(total)
66-}
67-
68-/**
69- * Exact entropy (bits) of round(N(0, s²) + U[-½,½)) — the marginal of a
70- * stored sample when dither is on. Conditioned on the Gaussian landing at t,
71- * bin j is hit with probability equal to the overlap of the dither interval
72- * with the bin, the triangular hat Λ(j−t) = max(0, 1−|j−t|); so p_j is the
73- * density integrated against Λ, done by Simpson on each side of the kink.
74- */
75-export function ditheredQuantizedGaussianEntropy(s: number): number {
76- if (s <= 0) return 0
77- // Approaches ½log₂(2πe s²) from above like log₂e/(12 s²) — the dither's
78- // 1/12 of variance plus the Δ²/24 quantization correction, each worth
79- // log₂e/(24 s²). Within 10⁻⁵ bits at s ≥ 6.
80- if (s >= 6) return 0.5 * Math.log2(TWO_PI_E * s * s) + Math.LOG2E / (12 * s * s)
81- if (s <= 0.1) {
82- // Only the neighbors of zero are reachable, through the tip of the hat:
83- // p±1 = ∫₀^∞ t φ_s(t) dt = s/√(2π), machine-exact in this range.
84- const p1 = s / Math.sqrt(2 * Math.PI)
85- const p0 = 1 - 2 * p1
86- return -p0 * Math.log2(p0) - 2 * p1 * Math.log2(p1)
87- }
88- const zMax = Math.ceil(8 * s + 2)
89- const m = Math.min(401, Math.max(9, 2 * Math.ceil(3 / s) + 9)) | 1
90- const h = 1 / (m - 1)
91- const density = (u: number) => Math.exp((-u * u) / (2 * s * s)) / (Math.sqrt(2 * Math.PI) * s)
92- // Simpson of φ_s(t)·w(t) over [a, a+1] with w linear from w0 to w1.
93- const half = (a: number, w0: number, w1: number) => {
94- let acc = density(a) * w0 + density(a + 1) * w1
95- for (let i = 1; i < m - 1; i++) {
96- const t = i * h
97- acc += (i % 2 === 1 ? 4 : 2) * density(a + t) * (w0 + (w1 - w0) * t)
98- }
99- return (acc * h) / 3
100- }
101- let sumH = 0
102- let total = 0
103- for (let j = -zMax; j <= zMax; j++) {
104- const p = half(j - 1, 0, 1) + half(j, 1, 0)
105- if (p > 0) {
106- sumH -= p * Math.log2(p)
107- total += p
108- }
109- }
110- return sumH / total + Math.log2(total)
111-}
112-
113-export function theoreticalRateBits(kernel: Float64Array, sigma: number, dither: boolean): number {
114- const nu = dither ? 1 / 6 : 1 / 12
115- const L = kernel.length
116- let rspec = 0
117- // Midpoints on [0, ½]; |H| is symmetric about ½ for a real kernel, so the
118- // grid average equals the integral over the full frequency circle.
119- for (let k = 0; k < INTEGRATION_POINTS; k++) {
120- const f = (0.5 * (k + 0.5)) / INTEGRATION_POINTS
121- let re = 0
122- let im = 0
123- for (let i = 0; i < L; i++) {
124- re += kernel[i] * Math.cos(2 * Math.PI * f * i)
125- im -= kernel[i] * Math.sin(2 * Math.PI * f * i)
126- }
127- rspec += 0.5 * Math.log2(TWO_PI_E * (sigma * sigma * (re * re + im * im) + nu))
128- }
129- rspec /= INTEGRATION_POINTS
130- let v = 0
131- for (let i = 0; i < L; i++) v += kernel[i] * kernel[i]
132- v *= sigma * sigma
133- const rsamp = dither
134- ? ditheredQuantizedGaussianEntropy(Math.sqrt(v))
135- : quantizedGaussianEntropy(Math.sqrt(v))
136- return Math.min(rspec, rsamp)
137-}
src/worker/compressWorker.tsmodified+2−3View file
@@ -23,7 +23,6 @@ export interface CompressRequest {
2323 id: number
2424 kernel: Float64Array
2525 sigma: number
26- dither: boolean
2726 blockSize: number
2827 lpcOrder: number
2928 seed: number
@@ -44,10 +43,10 @@ const PLAIN_CODECS = [ZLIB, ZSTD, ANS, DELTA_ZLIB, DELTA_ZSTD, DELTA_ANS]
4443 const post = self.postMessage as (message: CompressResponse) => void
4544
4645 self.onmessage = async (e: MessageEvent<CompressRequest>) => {
47- const { id, kernel, sigma, dither, blockSize, lpcOrder, seed } = e.data
46+ const { id, kernel, sigma, blockSize, lpcOrder, seed } = e.data
4847 try {
4948 await initCodecs()
50- const samples = new LatentSource(seed).window(0, blockSize, kernel, sigma, dither)
49+ const samples = new LatentSource(seed).window(0, blockSize, kernel, sigma)
5150 let sum = 0
5251 let sumSq = 0
5352 for (let i = 0; i < samples.length; i++) {