Charge rounding at full variance: R = min(spectral rate, one-sample ceiling)
The previous spectral estimate charged the quantizer at its entropy power
(additive constant 1), which systematically underestimated the MC ground
truth. Per Fourier mode the roundoff(+dither) noise is CLT-Gaussianized and
enters at full variance nu = 1/12 (1/6 dithered), and a subadditivity ceiling
- the exact marginal entropy of one stored sample - covers the sub-threshold
collapse. Remove the closed-form/exact-per-mode toggle, both variants being
superseded; update the math section, README, and the true_rate.py port, and
stop tracking scripts/__pycache__.
8 changed files+238−73
.gitignoremodified+2−0View file
@@ -1,2 +1,4 @@
11 node_modules
22 dist
3+__pycache__
4+tmp
README.mdmodified+18−13View file
@@ -16,23 +16,28 @@ to 128. Each prefilter group also carries a hollow bar: the order-0 entropy of
1616 the stream being coded, the limit a per-sample entropy coder cannot beat, which
1717 ANS misses by 1–2% (its symbol table plus its own arithmetic loss).
1818
19-Alongside the measurements it plots a theoretical bits/sample R: quantization
20-is modeled as an additive white noise floor on the spectrum, the one-step
21-Wiener prediction error of the resulting process comes from the
22-Szegő–Kolmogorov formula, and R is the exact entropy of that innovation
23-quantized at unit step:
19+Alongside the measurements it plots a theoretical bits/sample R — the smaller
20+of a spectral estimate and a rigorous one-sample ceiling:
2421
2522 ```
26-S_z(f) = σ²|H(f)|² + σ_q² σ_q² = 1/12 (1/6 with dither)
27-σ_e² = exp( 2 ∫₀^½ ln S_z(f) df )
28-R = H_Δ(σ_e) (exact quantized-Gaussian entropy)
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²
2926 ```
3027
31-Where the spectrum sits well above one step² this reduces to the classical
32-Gaussian entropy rate ½log₂(2πe) + ∫log₂S df; the noise floor keeps it finite
33-and positive where a deep stopband would send that integral to −∞. LPC + ANS
34-should approach R; probing where the approximation holds is the point. The
35-math section is a stub for the full derivation.
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.
3641
3742 ## Run it
3843
scripts/__pycache__/true_rate.cpython-313.pycdeleted+0−0View file
Binary file not shown.
scripts/true_rate.pymodified+47−8View file
@@ -76,9 +76,17 @@ def design_kernel(args):
7676 # ------------------------------------------------- the app's formula for R
7777 # Ported from src/model/theory.ts. Phi via math.erf (machine precision).
7878
79+TWO_PI_E = 2 * math.pi * math.e
80+
81+
7982 def quantized_gaussian_entropy(s):
8083 if s <= 0.02:
8184 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)
8290 zmax = int(math.ceil(8 * s + 4))
8391 H = 0.0
8492 prev = 0.5 * (1 + math.erf((-zmax - 0.5) / (s * math.sqrt(2))))
@@ -91,13 +99,42 @@ def quantized_gaussian_entropy(s):
9199 return H
92100
93101
94-def formula_rate(kernel, sigma, dither, points=8192):
95- noise_var = 1 / 6 if dither else 1 / 12
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
96130 f = (0.5 * (np.arange(points) + 0.5)) / points
97131 w = -2j * np.pi * np.outer(f, np.arange(len(kernel)))
98132 S = sigma * sigma * np.abs(np.exp(w) @ kernel) ** 2
99- integral = float(np.log2(S + noise_var).mean()) * 0.5
100- return quantized_gaussian_entropy(2.0 ** integral)
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
101138
102139
103140 # ------------------------------------------------ vectorized normal helpers
@@ -290,10 +327,11 @@ def main():
290327 L = len(kernel)
291328 M = args.past if args.past is not None else max(512, 4 * L)
292329
293- R = formula_rate(kernel, args.sigma, args.dither)
330+ R_spec, R_samp = formula_rates(kernel, args.sigma, args.dither)
331+ R = min(R_spec, R_samp)
294332 print(f'model: sigma={args.sigma} filter={args.filter} L={L} dither={args.dither}')
295- print(f'formula R (as shown in the app): {R:.4f} bits/sample'
296- f' (ratio vs int16: {16 / R:.3f}x)' if R > 0 else f'formula R: {R:.4f} bits/sample')
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}')
297335 print(f'MC: {args.pasts} pasts x {args.sweeps} sweeps, conditioning on M={M} samples')
298336
299337 rng = np.random.default_rng(args.seed)
@@ -309,7 +347,8 @@ def main():
309347 se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs)))
310348 print(f'\nMC entropy rate: {mean:.4f} +/- {se:.4f} bits/sample'
311349 f' (ratio vs int16: {16 / mean:.3f}x)')
312- print(f'formula R: {R:.4f} bits/sample (formula - MC = {R - mean:+.4f})')
350+ print(f'formula R: {R:.4f} bits/sample (formula - MC = {R - mean:+.4f}; '
351+ f'spec {R_spec:.4f}, samp {R_samp:.4f})')
313352 print('note: the MC value estimates H(z_next | M past samples), an upper bound '
314353 'on the rate that tightens as --past grows.')
315354
src/App.tsxmodified+5−3View file
@@ -229,9 +229,11 @@ export default function App() {
229229 coefficients). Baseline is raw int16 (16 bits/sample). The hollow bar in each group is
230230 that group's entropy limit — the order-0 entropy of the stream being coded, which no
231231 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 spectral
233- formula in the math section — approximate where quantization dominates the spectrum
234- (see the S(f) = 1 threshold on the response plot).
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.
235237 </p>
236238 <CopyableCommand
237239 label="check R against a Monte-Carlo ground truth (runs locally, ~2 min):"
src/app.cssmodified+9−0View file
@@ -507,6 +507,15 @@ body {
507507 max-width: 76ch;
508508 }
509509
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+
510519 .math-section .katex-display {
511520 margin: 14px 0;
512521 }
src/components/MathSection.tsxmodified+69−28View file
@@ -23,20 +23,21 @@ function Def({ tex, children }: { tex: string; children: React.ReactNode }) {
2323
2424 /**
2525 * The theoretical rate R exactly as `model/theory.ts` computes it, with every
26- * symbol defined. The derivation that justifies it is still to be written.
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.
2729 */
2830 export default function MathSection() {
2931 return (
3032 <div className="math-section">
3133 <p>
32- The dashed line on the compression chart is R, the predicted bits per sample. It is
33- computed in three steps: the spectrum of the stored signal, the residual an ideal
34- predictor leaves, and the entropy of that residual on the integer grid.
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:
3536 </p>
3637
37- <Tex display tex="S_z(f) \;=\; \sigma^2\,\big|H(f)\big|^2 \;+\; \sigma_q^2, \qquad H(f) \;=\; \sum_{n=0}^{L-1} h_n\, e^{-2\pi i f n}" />
38- <Tex display tex="\sigma_e \;=\; 2^{\,\int_0^{1/2} \log_2 S_z(f)\,df}" />
39- <Tex display tex="R \;=\; -\sum_{z \in \mathbb{Z}} p_z \log_2 p_z, \qquad p_z \;=\; \Phi\!\left(\frac{z + \tfrac12}{\sigma_e}\right) - \Phi\!\left(\frac{z - \tfrac12}{\sigma_e}\right)" />
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" />
4041
4142 <dl className="defs">
4243 <Def tex="\sigma">
@@ -48,27 +49,19 @@ export default function MathSection() {
4849 them
4950 </Def>
5051 <Def tex="H(f)">
51- the kernel's frequency response, the quantity plotted in dB as |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
5254 </Def>
53- <Def tex="f">
54- frequency in cycles per sample, running from 0 to ½ (Nyquist); the plots label the same
55- axis in Hz, as f times the sample rate
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
5658 </Def>
57- <Def tex="\sigma_q^2">
58- variance charged to rounding, treated as additive white noise: 1/12 for the roundoff
59- alone, 1/6 when dither is on (the dither is stored in the integers, so its 1/12 adds)
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)
6062 </Def>
61- <Def tex="S_z(f)">
62- power spectrum of the stored integer signal, in steps² per unit frequency
63- </Def>
64- <Def tex="\sigma_e">
65- standard deviation of the innovation — what an ideal linear predictor still cannot
66- predict from all earlier samples. The exponent is the Szegő–Kolmogorov formula for the
67- one-step prediction error, the geometric mean of the spectrum.
68- </Def>
69- <Def tex="\Phi">standard normal cumulative distribution function</Def>
70- <Def tex="p_z">
71- probability that the innovation, rounded to the integer grid, lands on z
63+ <Def tex="v">
64+ variance of a single output sample, the integral of S(f)
7265 </Def>
7366 <Def tex="R">
7467 bits per sample; the compression ratio the chart marks is 16/R, against 16-bit integer
@@ -76,10 +69,58 @@ export default function MathSection() {
7669 </Def>
7770 </dl>
7871
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+
79118 <p className="card-note">
80- The integral is evaluated by the midpoint rule on 8192 points and the sum over z is taken
81- out to where the remaining mass is negligible. A derivation — and an account of where
82- modeling the roundoff as white noise stops being fair — is still to be written.
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.
83124 </p>
84125 </div>
85126 )
src/model/theory.tsmodified+88−21View file
@@ -1,25 +1,35 @@
11 /**
2- * The theoretical bits/sample for the quantized filtered-Gaussian process.
2+ * The theoretical bits/sample for the quantized filtered-Gaussian process:
33 *
4- * The pure high-resolution entropy rate ½log₂(2πe) + ∫log₂ S(f) df diverges
5- * to -∞ wherever the spectrum falls far below the quantization step, so it is
6- * useless for filters with deep stopbands. Instead, model the roundoff as an
7- * additive white noise floor σ_q² (1/12 without dither; 1/6 with, since the
8- * dither itself is carried into the stored integers):
4+ * R̂ = min(R_spec, R_samp)
95 *
10- * S_z(f) = S(f) + σ_q², S(f) = σ² |H(f)|²
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²
1110 *
12- * The one-step Wiener prediction error of that process (Szegő/Kolmogorov)
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.
1323 *
14- * σ_e² = exp( 2 ∫₀^{1/2} ln S_z(f) df )
15- *
16- * is what an ideal predictor leaves behind; the rate is the exact entropy of
17- * that innovation quantized at unit step, R = H_Δ(σ_e). Where S ≫ 1 this
18- * reduces to the classical ½log₂(2πe σ_e²); in the coarse regime it stays
19- * positive and finite. It remains an approximation — roundoff is not truly
20- * white or independent — and testing it against LPC+ANS is the app's point.
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.
2130 */
2231
32+const TWO_PI_E = 2 * Math.PI * Math.E
2333 const INTEGRATION_POINTS = 8192
2434
2535 /**
@@ -29,6 +39,11 @@ const INTEGRATION_POINTS = 8192
2939 */
3040 export function quantizedGaussianEntropy(s: number): number {
3141 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)
3247 const zMax = Math.ceil(8 * s + 4)
3348 // Enough points that a bin spans a few per standard deviation even when the
3449 // bin is wide compared to the distribution.
@@ -50,10 +65,57 @@ export function quantizedGaussianEntropy(s: number): number {
5065 return sumH / total + Math.log2(total)
5166 }
5267
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+
53113 export function theoreticalRateBits(kernel: Float64Array, sigma: number, dither: boolean): number {
54- const noiseVar = dither ? 1 / 6 : 1 / 12
114+ const nu = dither ? 1 / 6 : 1 / 12
55115 const L = kernel.length
56- let integral = 0
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.
57119 for (let k = 0; k < INTEGRATION_POINTS; k++) {
58120 const f = (0.5 * (k + 0.5)) / INTEGRATION_POINTS
59121 let re = 0
@@ -62,9 +124,14 @@ export function theoreticalRateBits(kernel: Float64Array, sigma: number, dither:
62124 re += kernel[i] * Math.cos(2 * Math.PI * f * i)
63125 im -= kernel[i] * Math.sin(2 * Math.PI * f * i)
64126 }
65- integral += Math.log2(sigma * sigma * (re * re + im * im) + noiseVar)
127+ rspec += 0.5 * Math.log2(TWO_PI_E * (sigma * sigma * (re * re + im * im) + nu))
66128 }
67- integral *= 0.5 / INTEGRATION_POINTS
68- // σ_e² = 2^(2·integral), so σ_e = 2^integral.
69- return quantizedGaussianEntropy(2 ** integral)
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)
70137 }