/ concept-collection / timeseries-compressibility
Sign in
concept-collection / timeseries-compressibility
timeseries-compressibility / src / model / theory.ts
70 lines · 2.8 KBCodeBlameHistory
5bab85aRatio-first chart, quantization-floor theory formula, line-segment view, fixed latent dataJeremy Magland 2 * The theoretical bits/sample for the quantized filtered-Gaussian process.
5bab85aRatio-first chart, quantization-floor theory formula, line-segment view, fixed latent dataJeremy Magland 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):
5bab85aRatio-first chart, quantization-floor theory formula, line-segment view, fixed latent dataJeremy Magland 10 * S_z(f) = S(f) + σ_q², S(f) = σ² |H(f)|²
11 *
12 * The one-step Wiener prediction error of that process (Szegő/Kolmogorov)
13 *
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.
23const INTEGRATION_POINTS = 8192
26 * Exact entropy (bits) of round(N(0, s²)) on the unit lattice. Per-bin
27 * probabilities by Simpson integration of the density, so no erf is needed
28 * and the tail keeps relative accuracy.
29 */
30export function quantizedGaussianEntropy(s: number): number {
31 if (s <= 0.02) return 0
32 const zMax = Math.ceil(8 * s + 4)
33 // Enough points that a bin spans a few per standard deviation even when the
34 // bin is wide compared to the distribution.
35 const m = Math.min(401, Math.max(9, 2 * Math.ceil(3 / s) + 9)) | 1
36 const h = 1 / (m - 1)
37 const density = (u: number) => Math.exp((-u * u) / (2 * s * s)) / (Math.sqrt(2 * Math.PI) * s)
38 let sumH = 0
39 let total = 0
40 for (let z = -zMax; z <= zMax; z++) {
41 let acc = density(z - 0.5) + density(z + 0.5)
42 for (let i = 1; i < m - 1; i++) acc += (i % 2 === 1 ? 4 : 2) * density(z - 0.5 + i * h)
43 const p = (acc * h) / 3
44 if (p > 0) {
45 sumH -= p * Math.log2(p)
46 total += p
47 }
48 }
49 // Renormalize away the residual quadrature/truncation mass.
50 return sumH / total + Math.log2(total)
53export function theoreticalRateBits(kernel: Float64Array, sigma: number, dither: boolean): number {
54 const noiseVar = dither ? 1 / 6 : 1 / 12
56 let integral = 0
57 for (let k = 0; k < INTEGRATION_POINTS; k++) {
58 const f = (0.5 * (k + 0.5)) / INTEGRATION_POINTS
59 let re = 0
60 let im = 0
61 for (let i = 0; i < L; i++) {
62 re += kernel[i] * Math.cos(2 * Math.PI * f * i)
63 im -= kernel[i] * Math.sin(2 * Math.PI * f * i)
64 }
5bab85aRatio-first chart, quantization-floor theory formula, line-segment view, fixed latent dataJeremy Magland 65 integral += Math.log2(sigma * sigma * (re * re + im * im) + noiseVar)
67 integral *= 0.5 / INTEGRATION_POINTS
5bab85aRatio-first chart, quantization-floor theory formula, line-segment view, fixed latent dataJeremy Magland 68 // σ_e² = 2^(2·integral), so σ_e = 2^integral.
69 return quantizedGaussianEntropy(2 ** integral)
moveopenescclose