Keep the Python evidence for the entropy-rate gap
The Python work behind the new coder, so the claim that the estimator is right
and the gap is real can be re-run rather than taken on faith. Every number is
reproducible from these scripts, with the raw run logs alongside.
analysis_spectral models rounding as additive U(-1/2,1/2) noise, takes the
Kolmogorov innovation variance of sigma^2|H|^2 + 1/12, and
splits the ideal conditional entropy (0.723 bits) from the
mixture entropy an integer residual is stuck with (1.067)
lpc_ans_repro the app's integer LPC, ported exactly -- 15-bit coeffs,
shared shift, floor prediction, int16 wrap -- reaching the
same 14.3x, so the plateau is not a browser artifact
codec_gaussian the conditional-Gaussian coder in Python: 0.7184 bits
(22.3x) on 1M samples, round-trip exact. This is the proof
-- real bytes, not a bound
gibbs_predictive_rate the nonlinear optimum: both sides Gibbs-sample the latents
given the decoded past and code against the
Rao-Blackwellized posterior predictive, 0.674 +/- 0.065
(23.7x), landing on the estimator's own number
memory_length why the conditioning window M is sized as it is: y is
MA(L-1), but its optimal predictor is AR(infinity), and it
is the 1/12 rounding floor that makes the whitening filter
converge at all. 100 taps needs M >~ 53-96; the default
max(512, 4L) clears it several times over
11 changed files+974−0
exploration/RESULTS.mdadded+102−0View file
@@ -0,0 +1,102 @@
1+# Why LPC+ANS gets 14x when the entropy rate says 23x
2+
3+Setting: `--sigma 5 --filter bandpass --low 300 --high 2000 --taps 31 --rate 30000`.
4+All numbers are bits/sample on the quantized signal z (ratio = 16 / bits).
5+
6+## Verdict
7+
8+**The entropy estimate is right, and the gap is real.** LPC+ANS is not
9+suboptimal by accident — integer-residual prediction followed by a memoryless
10+entropy coder structurally cannot reach the conditional entropy of this
11+process. A codec that codes each sample against a discretized Gaussian
12+centred at the *real-valued* prediction reaches 22.3x with a verified
13+lossless round-trip; a Monte-Carlo model-based codec reaches ~23.7x.
14+
15+## The mechanism
16+
17+After optimal linear prediction from the quantized past, the prediction error
18+of y_next has std s ~= 0.28 quantization steps. The conditional distribution
19+of z_next = round(y_next) is therefore a Gaussian bump covering ~1-3 integer
20+bins, and its shape depends strongly on the *fractional part* of the
21+prediction mu: mu near a bin centre gives ~0.4 bits, mu near a bin edge gives
22+~1.1 bits.
23+
24+- An ideal conditional coder uses mu per sample: mean entropy ~= **0.72 bits**.
25+- LPC+ANS codes the integer residual z - floor-prediction with one pooled
26+ histogram: that histogram is the *mixture* over all fractional parts,
27+ entropy ~= **1.07 bits** (semi-analytic; 1.10 measured with the app's
28+ 15-bit integer coefficients).
29+
30+The ~0.35-0.38 bit gap **is** the 14x-vs-22x discrepancy. Raising the LPC
31+order does nothing (14.6x plateau by order 64): the loss is in the residual
32+coding, not the prediction.
33+
34+## Measurements (4M-sample realizations; codecs charged for all side info)
35+
36+| method | bits/sample | ratio |
37+|---|---|---|
38+| raw int16 | 16 | 1x |
39+| order-0 entropy of z | 1.954 | 8.2x |
40+| LPC(32) + ANS (app's method, reproduced) | 1.101 (+~1.5% ANS) | 14.3x |
41+| LPC(128) + ANS ceiling | 1.098 | 14.6x |
42+| **conditional-Gaussian arithmetic codec (real bytes, exact round-trip)** | **0.7184** | **22.27x** |
43+| Gibbs posterior-predictive codec (achievable rate, 60 sweeps) | 0.674 +/- 0.065 | 23.7x |
44+| timeseries-entropy estimator, M=512 (default settings) | 0.615 +/- 0.065 | 26x |
45+| timeseries-entropy estimator, M=1024, thin=2 | 0.658 +/- 0.046 | 24.3x |
46+| semi-analytic conditional entropy (rounding-as-noise model) | 0.723 | 22.1x |
47+
48+The true entropy rate is bracketed: <= 0.718 is *proven* by real compressed
49+bytes that round-trip; ~0.674 is achievable by the model-based codec; the
50+Monte-Carlo estimator's 0.62-0.68 is consistent with both. The estimator's
51+Rhee-Glynn pilot (`--pilot 6`) shows RMS Delta_m <= 0.016 decaying at
52+exponent ~2, so its r=1.5 truncation is sound at these settings.
53+
54+## The codec that proves it (`codec_gaussian.py`)
55+
56+1. Fit an order-64 real-coefficient predictor on the block (Levinson);
57+ send coefficients as float32 in the header (270 B total).
58+2. Fit a single residual scale s (grid search on the block); send it too.
59+3. For each sample: mu_t = dot(coeffs, previous 64 decoded samples);
60+ code z_t with an arithmetic coder under P(k) proportional to
61+ Phi((k+1/2-mu_t)/s) - Phi((k-1/2-mu_t)/s), tails folded into edge bins,
62+ 16-bit frequencies. Decoder computes the identical mu_t from decoded
63+ samples, so tables agree bit for bit.
64+4. First 64 samples are coded under the marginal N(0, std_z) table.
65+
66+1,048,576 samples -> 94,157 bytes = 0.7184 bits/sample = **22.27x**,
67+`np.array_equal(z, decoded) == True`. Encode/decode ~7 s each in pure
68+Python+numpy.
69+
70+`gibbs_predictive_rate.py` pushes further: both sides could run the
71+timeseries-entropy Gibbs sampler on the decoded past (shared seed) and code
72+against the Rao-Blackwellized posterior predictive of the generative model —
73+the nonlinear optimum. Its measured cross-entropy is 0.674 +/- 0.065
74+bits/sample (23.7x). The -log2-of-a-mean is Jensen-biased upward, so more
75+sweeps would only lower it toward the true conditional entropy.
76+
77+## Files
78+
79+- `analysis_spectral.py` — semi-analytic model: Kolmogorov innovation
80+ variance of S_z = sigma^2|H|^2 + 1/12; conditional vs mixture entropy;
81+ empirical Wiener-predictor rates by order.
82+- `lpc_ans_repro.py` — exact port of the app's integer LPC (15-bit coeffs,
83+ shared shift, floor prediction, int16 wrap); reproduces 14.3x.
84+- `codec_gaussian.py` — the round-tripping conditional-Gaussian codec.
85+- `gibbs_predictive_rate.py` — achievable rate of the Monte-Carlo
86+ model-based codec, reusing timeseries_entropy.ConditionalChain.
87+- `entropy_cli_run.log`, `entropy_pilot.log`, `codec_gaussian_1M.log`,
88+ `gibbs_predictive.log` — raw outputs.
89+
90+## Notes for the app / estimator
91+
92+- The app's hollow "order-0 entropy of the residual" bar is the LPC+ANS
93+ *ceiling*, and the measured 14x sits right on it — the coder is fine; the
94+ prefilter family is the limit.
95+- The tenth method now exists in the app: `src/compress/conditionalGaussian.ts`,
96+ computed with the other nine in the compression worker and drawn as its own
97+ "conditional" chart group. In-browser at these settings, 100k block:
98+ LPC(32)+ANS 1.119 bits (14.30x) vs cond. Gaussian AC 0.738 bits (21.68x).
99+- With prediction-error std s, the fractional-phase loss of integer-residual
100+ coding is ~E[H(mix)] - E[H(cond)]; it vanishes for s >> 1 (white/large
101+ sigma) and grows as s drops below ~0.5 — worst exactly in the
102+ narrowband, moderate-sigma regime this parameter set sits in.
exploration/analysis_spectral.pyadded+130−0View file
@@ -0,0 +1,130 @@
1+"""Semi-analytic cross-check of the entropy rate for
2+x iid N(0, sigma^2) -> 31-tap bandpass 300-2000 Hz @ 30 kHz -> round.
3+
4+Model the rounding as additive iid U(-1/2, 1/2) noise (variance 1/12) that is
5+independent of y. Then z = y + q is a stationary process with spectrum
6+S_z(f) = sigma^2 |H(f)|^2 + 1/12, and
7+
8+ - one-step innovation variance of z: P_z = exp( mean_f ln S_z(f) ) (Kolmogorov)
9+ - Var(y_next | past z) = s^2 = P_z - 1/12 (q_next independent of past)
10+ - H(z_next | past z) ~= E_frac [ entropy of round-binned N(mu, s^2) ]
11+ averaged over the fractional part of mu (approximately uniform)
12+ - LPC + memoryless coder ~= entropy of the MIXTURE distribution of
13+ round(y_next) - floor(mu): the integer-residual histogram pools all
14+ fractional parts, so it cannot use frac(mu).
15+
16+Also validated empirically on a long simulated realization with a
17+real-coefficient Wiener predictor.
18+"""
19+import numpy as np
20+from scipy.special import ndtr
21+
22+SIGMA = 5.0
23+RATE = 30000.0
24+LOW, HIGH, TAPS = 300.0, 2000.0, 31
25+
26+
27+def windowed_sinc_lowpass(fc, taps):
28+ n = taps | 1
29+ mid = (n - 1) / 2
30+ i = np.arange(n)
31+ t = i - mid
32+ sinc = np.where(t == 0, 2 * fc,
33+ np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
34+ w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
35+ h = sinc * w
36+ return h / h.sum()
37+
38+
39+def windowed_sinc_bandpass(f_lo, f_hi, taps):
40+ return windowed_sinc_lowpass(f_hi, taps) - windowed_sinc_lowpass(f_lo, taps)
41+
42+
43+def binned_gaussian_entropy(mu_frac, s, kmax=60):
44+ """Entropy (bits) of round(N(mu_frac, s^2)) for scalar/array mu_frac."""
45+ mu_frac = np.atleast_1d(mu_frac)[:, None]
46+ k = np.arange(-kmax, kmax + 1)[None, :]
47+ p = ndtr((k + 0.5 - mu_frac) / s) - ndtr((k - 0.5 - mu_frac) / s)
48+ with np.errstate(divide='ignore', invalid='ignore'):
49+ terms = np.where(p > 0, -p * np.log2(p), 0.0)
50+ return terms.sum(axis=1)
51+
52+
53+def mixture_residual_entropy(s, kmax=60, nfrac=4001):
54+ """Entropy (bits) of round(y) - floor(mu), y ~ N(mu, s^2), frac(mu) uniform."""
55+ u = (np.arange(nfrac) + 0.5) / nfrac # frac part of mu in (0,1)
56+ k = np.arange(-kmax, kmax + 1)[None, :]
57+ p = ndtr((k + 0.5 - u[:, None]) / s) - ndtr((k - 0.5 - u[:, None]) / s)
58+ pmix = p.mean(axis=0)
59+ pmix = pmix[pmix > 0]
60+ return float(-(pmix * np.log2(pmix)).sum())
61+
62+
63+def main():
64+ h = windowed_sinc_bandpass(LOW / RATE, HIGH / RATE, TAPS)
65+ var_y = SIGMA**2 * (h**2).sum()
66+ print(f'kernel taps={len(h)} sum h^2 = {(h**2).sum():.6f}')
67+ print(f'Var(y) = {var_y:.4f} std(y) = {np.sqrt(var_y):.4f} steps')
68+
69+ nfft = 1 << 18
70+ Hf2 = np.abs(np.fft.fft(h, nfft))**2
71+ S_y = SIGMA**2 * Hf2
72+ S_z = S_y + 1.0 / 12.0
73+ P_z = np.exp(np.mean(np.log(S_z)))
74+ s2 = P_z - 1.0 / 12.0
75+ s = np.sqrt(s2)
76+ print(f'\nadditive-noise model:')
77+ print(f' innovation var of z: P_z = {P_z:.5f}')
78+ print(f' Var(y_next | past z): s^2 = {s2:.5f} s = {s:.4f}')
79+
80+ # conditional entropy: average over fractional part of the predictive mean
81+ fracs = (np.arange(4001) + 0.5) / 4001 - 0.5
82+ Hc = binned_gaussian_entropy(fracs, s).mean()
83+ print(f'\n H(z_next | past z) ~= {Hc:.4f} bits/sample '
84+ f'(ratio {16/Hc:.1f}x) <- what an ideal conditional coder gets')
85+
86+ Hmix = mixture_residual_entropy(s)
87+ print(f' order-0 entropy of integer LPC residual ~= {Hmix:.4f} bits/sample '
88+ f'(ratio {16/Hmix:.1f}x) <- LPC + ANS ceiling')
89+
90+ # what pooling the fractional part costs
91+ print(f' gap (frac-part information thrown away) = {Hmix - Hc:.4f} bits')
92+
93+ # ---- empirical check on a long realization -------------------------------
94+ rng = np.random.default_rng(1)
95+ N = 1 << 22
96+ x = SIGMA * rng.standard_normal(N + len(h) - 1)
97+ y = np.convolve(x, h, mode='valid')
98+ z = np.floor(y + 0.5)
99+ counts = np.unique(z, return_counts=True)[1]
100+ p = counts / counts.sum()
101+ H0 = -(p * np.log2(p)).sum()
102+ print(f'\nempirical ({N} samples): Var(z) = {z.var():.4f} '
103+ f'order-0 H(z) = {H0:.4f} bits')
104+
105+ # real-coefficient Wiener predictor of z_next from past z, various orders
106+ from scipy.linalg import solve_toeplitz
107+ from scipy.signal import lfilter
108+ zc = z - z.mean()
109+ maxlag = 256
110+ r = np.array([zc @ zc if lag == 0 else zc[lag:] @ zc[:-lag]
111+ for lag in range(maxlag + 1)]) / len(zc)
112+ for order in (8, 16, 32, 64, 128, 256):
113+ a = solve_toeplitz(r[:order], r[1:order + 1])
114+ pred = lfilter(np.concatenate(([0.0], a)), [1.0], zc) # strictly causal
115+ e = zc - pred
116+ ev = e[order:].var()
117+ s2e = ev - 1.0 / 12.0
118+ # ideal conditional-Gaussian coding rate at this order
119+ mu = pred[order:] + z.mean()
120+ zt = z[order:]
121+ se = np.sqrt(max(s2e, 1e-6))
122+ pz = (ndtr((zt + 0.5 - mu) / se) - ndtr((zt - 0.5 - mu) / se))
123+ pz = np.maximum(pz, 1e-30)
124+ rate = float(-np.log2(pz).mean())
125+ print(f' order {order:4d}: pred-err var {ev:.5f} '
126+ f'(s^2={s2e:.5f}) conditional-Gaussian rate {rate:.4f} bits')
127+
128+
129+if __name__ == '__main__':
130+ main()
exploration/codec_gaussian.pyadded+315−0View file
@@ -0,0 +1,315 @@
1+"""A codec that beats LPC+ANS on the quantized filtered-Gaussian process.
2+
3+The insight: with prediction-error std ~0.28 quantization steps, the
4+conditional distribution of z_t is a Gaussian bump over 1-3 integer bins whose
5+placement depends on the *fractional part* of the real-valued prediction.
6+Integer-residual LPC + a memoryless entropy coder pools all fractional parts
7+into one histogram and pays the mixture entropy (~1.10 bits). Coding z_t
8+against a discretized Gaussian centred at the real-valued prediction pays the
9+conditional entropy (~0.72 bits) instead.
10+
11+Codec = order-p linear predictor (coefficients fitted on the block, sent as
12+float32 in the header) + adaptive binary arithmetic coding of z_t under
13+N(mu_t, s^2) discretized to integer bins, where mu_t is the float prediction
14+from already-decoded samples and s is a single fitted scale sent in the
15+header. Encoder and decoder compute mu_t with the identical np.dot on
16+identical float64 data, so the frequency tables agree bit for bit.
17+
18+Everything the decoder needs is charged: order, coefficients, s, marginal
19+std for the warm-up samples, and the sample count.
20+"""
21+import math
22+import numpy as np
23+from scipy.linalg import solve_toeplitz
24+from scipy.signal import lfilter
25+from scipy.special import ndtr
26+
27+SIGMA = 5.0
28+RATE = 30000.0
29+LOW, HIGH, TAPS = 300.0, 2000.0, 31
30+
31+TOTAL_BITS = 16
32+TOTAL = 1 << TOTAL_BITS
33+KWIN = 32 # symbols are residuals in [-KWIN, KWIN-1] around round(mu)
34+SQRT2 = math.sqrt(2.0)
35+
36+
37+# ---------------------------------------------------------------- the process
38+def windowed_sinc_lowpass(fc, taps):
39+ n = taps | 1
40+ mid = (n - 1) / 2
41+ i = np.arange(n)
42+ t = i - mid
43+ sinc = np.where(t == 0, 2 * fc,
44+ np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
45+ w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
46+ h = sinc * w
47+ return h / h.sum()
48+
49+
50+def make_signal(n, seed=1):
51+ h = windowed_sinc_lowpass(HIGH / RATE, TAPS) - windowed_sinc_lowpass(LOW / RATE, TAPS)
52+ rng = np.random.default_rng(seed)
53+ x = SIGMA * rng.standard_normal(n + len(h) - 1)
54+ y = np.convolve(x, h, mode='valid')
55+ return np.floor(y + 0.5).astype(np.int64)
56+
57+
58+# ------------------------------------------------- arithmetic coder (WNC-style)
59+class BitWriter:
60+ def __init__(self):
61+ self.bytes = bytearray()
62+ self.acc = 0
63+ self.nbits = 0
64+
65+ def write(self, bit):
66+ self.acc = (self.acc << 1) | bit
67+ self.nbits += 1
68+ if self.nbits == 8:
69+ self.bytes.append(self.acc)
70+ self.acc = 0
71+ self.nbits = 0
72+
73+ def flush(self):
74+ while self.nbits:
75+ self.write(0)
76+
77+
78+class BitReader:
79+ def __init__(self, data):
80+ self.data = data
81+ self.pos = 0
82+
83+ def read(self):
84+ byte = self.data[self.pos >> 3] if (self.pos >> 3) < len(self.data) else 0
85+ bit = (byte >> (7 - (self.pos & 7))) & 1
86+ self.pos += 1
87+ return bit
88+
89+
90+class ArithEncoder:
91+ FULL = (1 << 32) - 1
92+ HALF = 1 << 31
93+ QUARTER = 1 << 30
94+
95+ def __init__(self):
96+ self.low = 0
97+ self.high = self.FULL
98+ self.pending = 0
99+ self.out = BitWriter()
100+
101+ def _emit(self, bit):
102+ self.out.write(bit)
103+ while self.pending:
104+ self.out.write(1 - bit)
105+ self.pending -= 1
106+
107+ def encode(self, cum, freq, tot):
108+ span = self.high - self.low + 1
109+ self.high = self.low + span * (cum + freq) // tot - 1
110+ self.low = self.low + span * cum // tot
111+ while True:
112+ if self.high < self.HALF:
113+ self._emit(0)
114+ elif self.low >= self.HALF:
115+ self._emit(1)
116+ self.low -= self.HALF
117+ self.high -= self.HALF
118+ elif self.low >= self.QUARTER and self.high < self.HALF + self.QUARTER:
119+ self.pending += 1
120+ self.low -= self.QUARTER
121+ self.high -= self.QUARTER
122+ else:
123+ break
124+ self.low <<= 1
125+ self.high = (self.high << 1) | 1
126+
127+ def finish(self):
128+ self.pending += 1
129+ self._emit(0 if self.low < self.QUARTER else 1)
130+ self.out.flush()
131+ return bytes(self.out.bytes)
132+
133+
134+class ArithDecoder:
135+ FULL = (1 << 32) - 1
136+ HALF = 1 << 31
137+ QUARTER = 1 << 30
138+
139+ def __init__(self, data):
140+ self.low = 0
141+ self.high = self.FULL
142+ self.inp = BitReader(data)
143+ self.value = 0
144+ for _ in range(32):
145+ self.value = (self.value << 1) | self.inp.read()
146+
147+ def target(self, tot):
148+ span = self.high - self.low + 1
149+ return ((self.value - self.low + 1) * tot - 1) // span
150+
151+ def consume(self, cum, freq, tot):
152+ span = self.high - self.low + 1
153+ self.high = self.low + span * (cum + freq) // tot - 1
154+ self.low = self.low + span * cum // tot
155+ while True:
156+ if self.high < self.HALF:
157+ pass
158+ elif self.low >= self.HALF:
159+ self.low -= self.HALF
160+ self.high -= self.HALF
161+ self.value -= self.HALF
162+ elif self.low >= self.QUARTER and self.high < self.HALF + self.QUARTER:
163+ self.low -= self.QUARTER
164+ self.high -= self.QUARTER
165+ self.value -= self.QUARTER
166+ else:
167+ break
168+ self.low <<= 1
169+ self.high = (self.high << 1) | 1
170+ self.value = (self.value << 1) | self.inp.read()
171+
172+
173+# --------------------------------------------- discretized-Gaussian frequencies
174+def freq_table(d, s):
175+ """Integer frequencies (sum TOTAL) for residual symbols -KWIN..KWIN-1:
176+ bin k has probability P(round(N(d, s^2)) = k), tails folded into the edge
177+ bins, every bin floored at 1. Deterministic in (d, s)."""
178+ inv = 1.0 / (s * SQRT2)
179+ m = min(KWIN - 1, int(8.0 * s) + 2)
180+ freqs = [1] * (2 * KWIN)
181+ spare = TOTAL - 2 * KWIN
182+ lo_cdf = 0.0
183+ probs = []
184+ for k in range(-m, m + 1):
185+ hi_cdf = 1.0 if k == m else 0.5 * math.erfc(-(k + 0.5 - d) * inv)
186+ probs.append(hi_cdf - lo_cdf)
187+ lo_cdf = hi_cdf
188+ scaled = [int(p * spare) for p in probs]
189+ deficit = spare - sum(scaled)
190+ scaled[max(range(len(scaled)), key=scaled.__getitem__)] += deficit
191+ for i, k in enumerate(range(-m, m + 1)):
192+ freqs[k + KWIN] += scaled[i]
193+ return freqs
194+
195+
196+def cumulative(freqs):
197+ cum = [0] * (len(freqs) + 1)
198+ for i, f in enumerate(freqs):
199+ cum[i + 1] = cum[i] + f
200+ return cum
201+
202+
203+# ------------------------------------------------------------------- the codec
204+def fit_model(z, order):
205+ zf = z.astype(np.float64)
206+ r = np.array([zf @ zf if lag == 0 else zf[lag:] @ zf[:-lag]
207+ for lag in range(order + 1)]) / len(zf)
208+ a = solve_toeplitz(r[:order], r[1:order + 1]).astype(np.float32)
209+ # residual scale: fit s by minimizing the ideal code length on the block
210+ pred = lfilter(np.concatenate(([0.0], a.astype(np.float64))), [1.0], zf)
211+ e = zf[order:] - pred[order:]
212+ s0 = math.sqrt(max(e.var() - 1.0 / 12.0, 1e-6))
213+ best = (None, np.inf)
214+ zt, mu = zf[order:], pred[order:]
215+ for s in s0 * np.linspace(0.85, 1.15, 13):
216+ p = ndtr((zt + 0.5 - mu) / s) - ndtr((zt - 0.5 - mu) / s)
217+ bits = float(-np.log2(np.maximum(p, 1e-12)).mean())
218+ if bits < best[1]:
219+ best = (s, bits)
220+ s = np.float32(best[0])
221+ std_z = np.float32(max(zf.std(), 1e-3))
222+ return a, s, std_z, best[1]
223+
224+
225+HEADER_BYTES = 4 + 2 + 4 + 4 # n, order, s, std_z (+ coefficients, counted below)
226+
227+
228+def encode(z, order):
229+ a, s, std_z, ideal_bits = fit_model(z, order)
230+ arev = a[::-1].astype(np.float64)
231+ zf = z.astype(np.float64)
232+ enc = ArithEncoder()
233+
234+ warm = freq_table(0.0, float(std_z))
235+ warm_cum = cumulative(warm)
236+ s_f = float(s)
237+ for t in range(len(z)):
238+ if t < order:
239+ table, cum = warm, warm_cum
240+ c = 0
241+ else:
242+ mu = float(np.dot(arev, zf[t - order:t]))
243+ c = math.floor(mu + 0.5)
244+ table = freq_table(mu - c, s_f)
245+ cum = cumulative(table)
246+ sym = int(z[t]) - c + KWIN
247+ if not 0 <= sym < 2 * KWIN:
248+ raise ValueError(f'residual out of range at {t}') # production: escape code
249+ enc.encode(cum[sym], table[sym], TOTAL)
250+ payload = enc.finish()
251+ header = HEADER_BYTES + 4 * order
252+ return payload, header, a, s, std_z, ideal_bits
253+
254+
255+def decode(payload, n, order, a, s, std_z):
256+ arev = a[::-1].astype(np.float64)
257+ zf = np.zeros(n, dtype=np.float64)
258+ z = np.zeros(n, dtype=np.int64)
259+ dec = ArithDecoder(payload)
260+ warm = freq_table(0.0, float(std_z))
261+ warm_cum = cumulative(warm)
262+ s_f = float(s)
263+ for t in range(n):
264+ if t < order:
265+ table, cum = warm, warm_cum
266+ c = 0
267+ else:
268+ mu = float(np.dot(arev, zf[t - order:t]))
269+ c = math.floor(mu + 0.5)
270+ table = freq_table(mu - c, s_f)
271+ cum = cumulative(table)
272+ tgt = dec.target(TOTAL)
273+ # find symbol: cum[sym] <= tgt < cum[sym+1]
274+ lo, hi = 0, 2 * KWIN
275+ while hi - lo > 1:
276+ mid = (lo + hi) // 2
277+ if cum[mid] <= tgt:
278+ lo = mid
279+ else:
280+ hi = mid
281+ sym = lo
282+ dec.consume(cum[sym], table[sym], TOTAL)
283+ z[t] = sym - KWIN + c
284+ zf[t] = float(z[t])
285+ return z
286+
287+
288+def main():
289+ import time
290+ n = 1 << 20
291+ order = 64
292+ z = make_signal(n)
293+
294+ t0 = time.time()
295+ payload, header, a, s, std_z, ideal_bits = encode(z, order)
296+ t1 = time.time()
297+ total_bytes = len(payload) + header
298+ bps = 8.0 * total_bytes / n
299+ print(f'encoded {n} samples in {t1 - t0:.1f}s')
300+ print(f' ideal model cross-entropy (block fit): {ideal_bits:.4f} bits/sample')
301+ print(f' payload {len(payload)} B + header {header} B = {total_bytes} B')
302+ print(f' -> {bps:.4f} bits/sample ratio vs int16: {16 / bps:.2f}x')
303+
304+ t0 = time.time()
305+ zdec = decode(payload, n, order, a, s, std_z)
306+ t1 = time.time()
307+ ok = bool(np.array_equal(z, zdec))
308+ print(f'decoded in {t1 - t0:.1f}s round-trip exact: {ok}')
309+ if not ok:
310+ bad = np.nonzero(z != zdec)[0][:5]
311+ print(f' first mismatches at {bad}')
312+
313+
314+if __name__ == '__main__':
315+ main()
exploration/codec_gaussian_1M.logadded+5−0View file
@@ -0,0 +1,5 @@
1+encoded 1048576 samples in 7.1s
2+ ideal model cross-entropy (block fit): 0.7149 bits/sample
3+ payload 93887 B + header 270 B = 94157 B
4+ -> 0.7184 bits/sample ratio vs int16: 22.27x
5+decoded in 7.7s round-trip exact: True
exploration/entropy_cli_M1024.logadded+53−0View file
@@ -0,0 +1,53 @@
1+model: sigma=5.0 filter=bandpass L=31 M=1024
2+48 pasts x 8 reps, n0=128 r=1.5 thin=2
3+ past 1/48: H = 0.7113 running mean 0.7113 +/- nan
4+ past 2/48: H = 0.9976 running mean 0.8544 +/- 0.1431
5+ past 3/48: H = 0.5336 running mean 0.7475 +/- 0.1352
6+ past 4/48: H = 1.0042 running mean 0.8117 +/- 0.1151
7+ past 5/48: H = 0.4127 running mean 0.7319 +/- 0.1196
8+ past 6/48: H = 0.7989 running mean 0.7430 +/- 0.0983
9+ past 7/48: H = 0.7610 running mean 0.7456 +/- 0.0831
10+ past 8/48: H = 0.5111 running mean 0.7163 +/- 0.0777
11+ past 9/48: H = 0.9950 running mean 0.7473 +/- 0.0752
12+ past 10/48: H = 0.9849 running mean 0.7710 +/- 0.0714
13+ past 11/48: H = 0.6413 running mean 0.7592 +/- 0.0656
14+ past 12/48: H = 0.3640 running mean 0.7263 +/- 0.0684
15+ past 13/48: H = 0.0534 running mean 0.6745 +/- 0.0814
16+ past 14/48: H = 0.8967 running mean 0.6904 +/- 0.0771
17+ past 15/48: H = 0.0997 running mean 0.6510 +/- 0.0818
18+ past 16/48: H = 0.2239 running mean 0.6243 +/- 0.0811
19+ past 17/48: H = 0.2020 running mean 0.5995 +/- 0.0801
20+ past 18/48: H = 0.0822 running mean 0.5707 +/- 0.0808
21+ past 19/48: H = 0.7783 running mean 0.5817 +/- 0.0772
22+ past 20/48: H = 0.9495 running mean 0.6001 +/- 0.0755
23+ past 21/48: H = 0.7712 running mean 0.6082 +/- 0.0723
24+ past 22/48: H = 0.2220 running mean 0.5907 +/- 0.0711
25+ past 23/48: H = 0.9836 running mean 0.6077 +/- 0.0701
26+ past 24/48: H = 0.7568 running mean 0.6140 +/- 0.0674
27+ past 25/48: H = 0.7307 running mean 0.6186 +/- 0.0648
28+ past 26/48: H = 0.9518 running mean 0.6314 +/- 0.0636
29+ past 27/48: H = 0.8170 running mean 0.6383 +/- 0.0616
30+ past 28/48: H = 0.8275 running mean 0.6451 +/- 0.0597
31+ past 29/48: H = 0.7904 running mean 0.6501 +/- 0.0578
32+ past 30/48: H = 0.6843 running mean 0.6512 +/- 0.0559
33+ past 31/48: H = 0.4105 running mean 0.6435 +/- 0.0546
34+ past 32/48: H = 0.9649 running mean 0.6535 +/- 0.0538
35+ past 33/48: H = 0.8406 running mean 0.6592 +/- 0.0525
36+ past 34/48: H = 0.1328 running mean 0.6437 +/- 0.0532
37+ past 35/48: H = 0.8578 running mean 0.6498 +/- 0.0520
38+ past 36/48: H = 0.0236 running mean 0.6324 +/- 0.0535
39+ past 37/48: H = 0.9863 running mean 0.6420 +/- 0.0529
40+ past 38/48: H = 0.8351 running mean 0.6471 +/- 0.0517
41+ past 39/48: H = 0.9460 running mean 0.6547 +/- 0.0509
42+ past 40/48: H = 0.8007 running mean 0.6584 +/- 0.0498
43+ past 41/48: H = 0.1739 running mean 0.6466 +/- 0.0500
44+ past 42/48: H = 0.5200 running mean 0.6435 +/- 0.0489
45+ past 43/48: H = 0.2239 running mean 0.6338 +/- 0.0487
46+ past 44/48: H = 0.8977 running mean 0.6398 +/- 0.0480
47+ past 45/48: H = 0.3833 running mean 0.6341 +/- 0.0472
48+ past 46/48: H = 0.9892 running mean 0.6418 +/- 0.0468
49+ past 47/48: H = 1.0626 running mean 0.6508 +/- 0.0467
50+ past 48/48: H = 0.9816 running mean 0.6576 +/- 0.0462
51+
52+H(z_next | 1024 past samples) = 0.6576 +/- 0.0462 bits/sample (ratio vs int16: 24.329x)
53+note: an upper bound on the entropy rate that tightens as --past grows.
exploration/entropy_cli_run.logadded+29−0View file
@@ -0,0 +1,29 @@
1+model: sigma=5.0 filter=bandpass L=31 M=512
2+24 pasts x 8 reps, n0=128 r=1.5 thin=1
3+ past 1/24: H = 0.7268 running mean 0.7268 +/- nan
4+ past 2/24: H = 0.6755 running mean 0.7012 +/- 0.0257
5+ past 3/24: H = 0.7716 running mean 0.7247 +/- 0.0278
6+ past 4/24: H = 0.2857 running mean 0.6149 +/- 0.1115
7+ past 5/24: H = 0.6559 running mean 0.6231 +/- 0.0867
8+ past 6/24: H = 0.4595 running mean 0.5958 +/- 0.0759
9+ past 7/24: H = 0.9782 running mean 0.6505 +/- 0.0843
10+ past 8/24: H = 0.6822 running mean 0.6544 +/- 0.0731
11+ past 9/24: H = 0.3457 running mean 0.6201 +/- 0.0730
12+ past 10/24: H = 0.7768 running mean 0.6358 +/- 0.0672
13+ past 11/24: H = 0.0000 running mean 0.5780 +/- 0.0838
14+ past 12/24: H = 0.3871 running mean 0.5621 +/- 0.0782
15+ past 13/24: H = 0.9726 running mean 0.5937 +/- 0.0785
16+ past 14/24: H = 0.0791 running mean 0.5569 +/- 0.0815
17+ past 15/24: H = 0.1500 running mean 0.5298 +/- 0.0806
18+ past 16/24: H = 0.8008 running mean 0.5467 +/- 0.0772
19+ past 17/24: H = 0.2154 running mean 0.5272 +/- 0.0751
20+ past 18/24: H = 1.1716 running mean 0.5630 +/- 0.0794
21+ past 19/24: H = 0.5296 running mean 0.5613 +/- 0.0751
22+ past 20/24: H = 1.0329 running mean 0.5849 +/- 0.0750
23+ past 21/24: H = 0.9684 running mean 0.6031 +/- 0.0737
24+ past 22/24: H = 0.7480 running mean 0.6097 +/- 0.0706
25+ past 23/24: H = 0.5959 running mean 0.6091 +/- 0.0674
26+ past 24/24: H = 0.7638 running mean 0.6155 +/- 0.0649
27+
28+H(z_next | 512 past samples) = 0.6155 +/- 0.0649 bits/sample (ratio vs int16: 25.993x)
29+note: an upper bound on the entropy rate that tightens as --past grows.
exploration/entropy_pilot.logadded+9−0View file
@@ -0,0 +1,9 @@
1+model: sigma=5.0 filter=bandpass L=31 M=512
2+pilot: 8 pasts, levels 1..6, n0=128
3+ m=1: rms Delta = 0.00777
4+ m=2: rms Delta = 0.01553 decay exponent -2.00
5+ m=3: rms Delta = 0.00699 decay exponent 2.30
6+ m=4: rms Delta = 0.00233 decay exponent 3.17
7+ m=5: rms Delta = 0.00241 decay exponent -0.10
8+ m=6: rms Delta = 0.00051 decay exponent 4.48
9+choose r safely below the E[Delta^2] decay exponent (and > 1); r=1.5 suits decay near 2.
exploration/gibbs_predictive.logadded+9−0View file
@@ -0,0 +1,9 @@
1+sigma*|h0| (innovation scale given full latent past) = 0.02161
2+ 50/300: rate = 1.0758 +/- 0.2135 bits/sample (floor hits: 0)
3+ 100/300: rate = 0.8102 +/- 0.1261 bits/sample (floor hits: 0)
4+ 150/300: rate = 0.7488 +/- 0.0916 bits/sample (floor hits: 0)
5+ 200/300: rate = 0.6868 +/- 0.0726 bits/sample (floor hits: 0)
6+ 250/300: rate = 0.6161 +/- 0.0598 bits/sample (floor hits: 0)
7+ 300/300: rate = 0.6742 +/- 0.0646 bits/sample (floor hits: 1)
8+
9+MC-codec achievable rate (M=512, 60 sweeps): 0.6742 +/- 0.0646 bits/sample -> ratio 23.7x
exploration/gibbs_predictive_rate.pyadded+115−0View file
@@ -0,0 +1,115 @@
1+"""Achievable rate of a model-based Monte-Carlo codec (the nonlinear optimum).
2+
3+The Gaussian-conditional codec (codec_gaussian.py) is linear: it conditions on
4+the past only through a linear prediction and charges a fixed Gaussian. The
5+true conditional law P(z_t | z_{t-M..t-1}) is the posterior predictive of the
6+generative model — computable by Gibbs-sampling the latent x under the box
7+constraints round(h*x) = z (exactly the machinery of timeseries-entropy) and
8+Rao-Blackwellizing:
9+
10+ P_hat(k) = mean over sweeps of Phi((k+1/2-c)/(sigma h0)) - Phi((k-1/2-c)/(sigma h0)),
11+ c = hr_head . x[M:]
12+
13+An encoder/decoder pair sharing the RNG seed can both compute P_hat from
14+already-decoded samples, so rate = E[-log2 P_tilde(z_t)] is achievable
15+(P_tilde mixes in a wide floor distribution so no symbol gets probability 0).
16+The -log2 of a finite-sweep average is on average >= the -log2 of the true
17+predictive (Jensen), so this measures an upper bound that tightens with more
18+sweeps: whatever number comes out IS achievable by this codec family.
19+
20+Chain start: the generating latents of the window, an exact draw from
21+p(x | all z) — near the target p(x | window z); BURN sweeps then wash the
22+difference. Reuses ConditionalChain's sweep by overwriting its state.
23+"""
24+import numpy as np
25+from scipy.special import ndtr
26+
27+import sys
28+from pathlib import Path
29+
30+# The companion package, cloned alongside this repo.
31+sys.path.insert(0, str(Path(__file__).resolve().parents[2] / 'timeseries-entropy' / 'src'))
32+from timeseries_entropy.model import ConditionalChain
33+
34+SIGMA = 5.0
35+RATE = 30000.0
36+LOW, HIGH, TAPS = 300.0, 2000.0, 31
37+
38+
39+def windowed_sinc_lowpass(fc, taps):
40+ n = taps | 1
41+ mid = (n - 1) / 2
42+ i = np.arange(n)
43+ t = i - mid
44+ sinc = np.where(t == 0, 2 * fc,
45+ np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
46+ w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
47+ h = sinc * w
48+ return h / h.sum()
49+
50+
51+def make_kernel():
52+ return windowed_sinc_lowpass(HIGH / RATE, TAPS) - windowed_sinc_lowpass(LOW / RATE, TAPS)
53+
54+
55+def posterior_predictive_rate(M=512, burn=30, sweeps=60, positions=300,
56+ spacing=64, seed=7, eps=1e-3):
57+ h = make_kernel()
58+ L = len(h)
59+ rng = np.random.default_rng(seed)
60+
61+ n = positions * spacing + M + 4 * L
62+ x = SIGMA * rng.standard_normal(n + L - 1)
63+ y = np.convolve(x, h, mode='valid')
64+ z = np.floor(y + 0.5)
65+
66+ h0 = abs(h[0])
67+ print(f'sigma*|h0| (innovation scale given full latent past) = {SIGMA * h0:.5f}')
68+
69+ # wide fallback so P_tilde is never 0 (what a real codec would also do)
70+ std_z = z.std()
71+ kmax = 32
72+ ks = np.arange(-kmax, kmax + 1)
73+ fallback = ndtr((ks + 0.5) / std_z) - ndtr((ks - 0.5) / std_z)
74+ fallback /= fallback.sum()
75+
76+ proto = ConditionalChain(h, SIGMA, past=M, rng=rng) # template; state overwritten
77+ rates = []
78+ hit_floor = 0
79+ for i in range(positions):
80+ t = M + i * spacing # predict z[t] from z[t-M:t]
81+ i0 = t - M
82+ proto.z = z[i0:t].copy()
83+ proto.x = x[i0:t + L - 1].copy() # exact posterior draw (given all z)
84+ P = L - 1
85+ proto.lo[P:P + M] = proto.z - 0.5
86+ proto.hi[P:P + M] = proto.z + 0.5
87+
88+ for _ in range(burn):
89+ proto._sweep()
90+ acc = 0.0
91+ zt = int(z[t])
92+ for _ in range(sweeps):
93+ proto._sweep()
94+ c = float(proto.hr_head @ proto.x[M:])
95+ acc += (ndtr((zt + 0.5 - c) / (SIGMA * h0))
96+ - ndtr((zt - 0.5 - c) / (SIGMA * h0)))
97+ p_hat = acc / sweeps
98+ p_tilde = (1 - eps) * p_hat + eps * float(fallback[np.clip(zt, -kmax, kmax) + kmax])
99+ if p_hat < eps:
100+ hit_floor += 1
101+ rates.append(-np.log2(p_tilde))
102+ if (i + 1) % 50 == 0:
103+ r = np.array(rates)
104+ print(f' {i + 1:4d}/{positions}: rate = {r.mean():.4f} '
105+ f'+/- {r.std(ddof=1) / np.sqrt(len(r)):.4f} bits/sample '
106+ f'(floor hits: {hit_floor})')
107+ r = np.array(rates)
108+ print(f'\nMC-codec achievable rate (M={M}, {sweeps} sweeps): '
109+ f'{r.mean():.4f} +/- {r.std(ddof=1) / np.sqrt(len(r)):.4f} bits/sample '
110+ f'-> ratio {16 / r.mean():.1f}x')
111+ return r
112+
113+
114+if __name__ == '__main__':
115+ posterior_predictive_rate()
exploration/lpc_ans_repro.pyadded+85−0View file
@@ -0,0 +1,85 @@
1+"""Reproduce the app's LPC(order) + ANS rate in Python.
2+
3+Exact port of src/compress/lpc.ts semantics: Levinson fit on the block's
4+autocorrelation, coefficients quantized to 15-bit signed ints with a shared
5+right shift, integer prediction with floor division, residual wrapped to
6+int16. The coded size is the order-0 entropy of the residual (the ANS bars in
7+the app sit 1-2% above this, for the symbol table plus arithmetic loss).
8+"""
9+import numpy as np
10+from scipy.linalg import solve_toeplitz
11+
12+SIGMA = 5.0
13+RATE = 30000.0
14+LOW, HIGH, TAPS = 300.0, 2000.0, 31
15+COEFF_PRECISION = 15
16+
17+
18+def windowed_sinc_lowpass(fc, taps):
19+ n = taps | 1
20+ mid = (n - 1) / 2
21+ i = np.arange(n)
22+ t = i - mid
23+ sinc = np.where(t == 0, 2 * fc,
24+ np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
25+ w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
26+ h = sinc * w
27+ return h / h.sum()
28+
29+
30+def make_signal(n, seed=1):
31+ h = windowed_sinc_lowpass(HIGH / RATE, TAPS) - windowed_sinc_lowpass(LOW / RATE, TAPS)
32+ rng = np.random.default_rng(seed)
33+ x = SIGMA * rng.standard_normal(n + len(h) - 1)
34+ y = np.convolve(x, h, mode='valid')
35+ return np.floor(y + 0.5).astype(np.int64)
36+
37+
38+def fit_lpc_quantized(z, order):
39+ zf = z.astype(np.float64)
40+ r = np.array([zf @ zf if lag == 0 else zf[lag:] @ zf[:-lag]
41+ for lag in range(order + 1)])
42+ a = solve_toeplitz(r[:order], r[1:order + 1])
43+ max_abs = np.abs(a).max()
44+ shift = COEFF_PRECISION - 1 - int(np.floor(np.log2(max_abs))) - 1
45+ shift = max(0, min(15, shift))
46+ limit = 2 ** (COEFF_PRECISION - 1)
47+ q = np.clip(np.round(a * 2.0**shift), -limit, limit - 1).astype(np.int64)
48+ return q, shift
49+
50+
51+def lpc_residual(z, q, shift):
52+ order = len(q)
53+ # pred[n] = floor( sum_k q[k] z[n-1-k] / 2^shift ) for n >= order
54+ acc = np.convolve(z, q, mode='full')[:len(z)] # acc[n] = sum_k q[k] z[n-k]
55+ pred = np.empty_like(z)
56+ pred[:] = 0
57+ pred[1:] = acc[:-1] // (1 << shift) # floor division, exact int64
58+ e = z.copy()
59+ e[order:] = z[order:] - pred[order:]
60+ e = ((e + (1 << 15)) & 0xFFFF) - (1 << 15) # int16 wraparound
61+ return e
62+
63+
64+def entropy_bits(v):
65+ counts = np.unique(v, return_counts=True)[1]
66+ p = counts / counts.sum()
67+ return float(-(p * np.log2(p)).sum())
68+
69+
70+def main():
71+ n = 1 << 22
72+ z = make_signal(n)
73+ print(f'{n} samples var(z) = {z.var():.4f} order-0 H(z) = {entropy_bits(z):.4f} bits '
74+ f'({16/entropy_bits(z):.1f}x)')
75+ for order in (8, 16, 32, 64, 128):
76+ q, shift = fit_lpc_quantized(z, order)
77+ e = lpc_residual(z, q, shift)
78+ H = entropy_bits(e[order:])
79+ coeff_bits = (2 + 2 * order) * 8 / n
80+ print(f' LPC({order:3d}): shift={shift:2d} residual H0 = {H:.4f} bits '
81+ f'-> ratio {16/H:.2f}x (with ~1.5% ANS overhead: {16/(H*1.015):.2f}x)')
82+
83+
84+if __name__ == '__main__':
85+ main()
exploration/memory_length.pyadded+122−0View file
@@ -0,0 +1,122 @@
1+"""How long must the conditioning window M be?
2+
3+The latent y = h*x is exactly (L-1)-dependent, so a naive reading says the
4+past only needs to be ~L samples. That is wrong for prediction: y is an
5+MA(L-1) process, and the optimal *predictor* of an MA process from its own
6+past is AR(infinity) — the whitening filter is 1/H(z), whose impulse response
7+decays at a rate set by how close H's zeros sit to the unit circle. For a
8+narrowband filter those zeros are very close, so the memory that matters runs
9+many multiples of L.
10+
11+What keeps it finite is the rounding: z = round(y) has spectrum
12+S_z = sigma^2 |H|^2 + 1/12, which never reaches zero, so its whitening filter
13+always converges. Deep nulls (large sigma, narrow band, many taps) mean slow
14+decay; the 1/12 floor sets where the decay finally bites.
15+
16+This computes, from that spectral model, the prediction order (= the past
17+length that matters) needed to get within 0.01 bits/sample of the infinite-
18+past limit — the M the estimator needs.
19+"""
20+import numpy as np
21+from scipy.special import ndtr
22+
23+RATE = 30000.0
24+LOW, HIGH = 300.0, 2000.0
25+
26+
27+def windowed_sinc_lowpass(fc, taps):
28+ n = taps | 1
29+ mid = (n - 1) / 2
30+ i = np.arange(n)
31+ t = i - mid
32+ sinc = np.where(t == 0, 2 * fc,
33+ np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
34+ w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
35+ h = sinc * w
36+ return h / h.sum()
37+
38+
39+def bandpass(taps):
40+ return windowed_sinc_lowpass(HIGH / RATE, taps) - windowed_sinc_lowpass(LOW / RATE, taps)
41+
42+
43+def levinson_all(r, P):
44+ """Prediction-error variance for every order 0..P, one O(P^2) pass.
45+ Also returns the order-P coefficients (the whitening filter's AR part)."""
46+ a = np.zeros(P)
47+ err = np.empty(P + 1)
48+ err[0] = r[0]
49+ e = r[0]
50+ for i in range(P):
51+ acc = r[i + 1] - (a[:i] @ r[i:0:-1] if i else 0.0)
52+ k = acc / e
53+ if i:
54+ a[:i] = a[:i] - k * a[i - 1::-1]
55+ a[i] = k
56+ e *= 1 - k * k
57+ err[i + 1] = e
58+ if e <= 0:
59+ err[i + 1:] = e
60+ break
61+ return err, a
62+
63+
64+def entropy_at(var_pred, nfrac=1001, kmax=60):
65+ """Ideal conditional-coding rate (bits) given prediction-error variance."""
66+ s = np.sqrt(max(var_pred - 1.0 / 12.0, 1e-9))
67+ d = ((np.arange(nfrac) + 0.5) / nfrac - 0.5)[:, None]
68+ k = np.arange(-kmax, kmax + 1)[None, :]
69+ p = ndtr((k + 0.5 - d) / s) - ndtr((k - 0.5 - d) / s)
70+ with np.errstate(divide='ignore', invalid='ignore'):
71+ h = np.where(p > 0, -p * np.log2(p), 0.0).sum(axis=1)
72+ return float(h.mean())
73+
74+
75+def report(sigma, taps, P=6000, tol_bits=0.01):
76+ h = bandpass(taps)
77+ L = len(h)
78+ # r_z[k] = sigma^2 (h corr h)[k] + (1/12) delta[k]; zero past lag L-1.
79+ rr = np.correlate(h, h, 'full')[L - 1:]
80+ r = np.zeros(P + 2)
81+ r[:L] = sigma**2 * rr
82+ r[0] += 1.0 / 12.0
83+
84+ err, a = levinson_all(r, P)
85+ nfft = 1 << 16
86+ S = sigma**2 * np.abs(np.fft.fft(h, nfft))**2 + 1.0 / 12.0
87+ P_inf = float(np.exp(np.mean(np.log(S)))) # Kolmogorov limit
88+
89+ H_inf = entropy_at(P_inf)
90+ # The rate is monotone in the prediction-error variance, and err[] is
91+ # monotone in the order — so bisect for the variance that costs tol bits,
92+ # then read off the first order that reaches it.
93+ vlo, vhi = P_inf, P_inf * 16
94+ for _ in range(60):
95+ vmid = 0.5 * (vlo + vhi)
96+ if entropy_at(vmid) - H_inf < tol_bits:
97+ vlo = vmid
98+ else:
99+ vhi = vmid
100+ reached = np.nonzero(err[1:] <= vlo)[0]
101+ need = int(reached[0]) + 1 if reached.size else None
102+
103+ # where the whitening filter's tail falls below 1e-3 of its peak
104+ tail = np.abs(a) / np.abs(a).max()
105+ decay = int(np.max(np.nonzero(tail > 1e-3)[0])) + 1 if (tail > 1e-3).any() else P
106+
107+ cli_default = max(512, 4 * L)
108+ flag = 'OK' if cli_default >= (need or P) else '** TOO SHORT **'
109+ print(f' sigma={sigma:<4g} taps={taps:<4d} L={L:<4d} '
110+ f'H_inf={H_inf:.3f} bits ({16/H_inf:5.1f}x) '
111+ f'need M>~{need if need else ">"+str(P):<5} '
112+ f'(={(need or P)/L:4.1f} L) AR tail {decay:<5d} '
113+ f'CLI default M={cli_default} -> {flag}')
114+
115+
116+if __name__ == '__main__':
117+ print('bandpass 300-2000 Hz @ 30 kHz; "need M" = order within 0.01 bits '
118+ 'of the infinite-past limit\n')
119+ for sigma in (5, 20):
120+ for taps in (15, 31, 101, 301):
121+ report(sigma, taps)
122+ print()