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 */
32const TWO_PI_E = 2 * Math.PI * Math.E
33const INTEGRATION_POINTS = 8192
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 */
40export 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}
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 */
75export 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}
113export 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}