9d25feaKeep the Python evidence for the entropy-rate gapJeremy Magland 1"""Reproduce the app's LPC(order) + ANS rate in Python.
3Exact port of src/compress/lpc.ts semantics: Levinson fit on the block's
4autocorrelation, coefficients quantized to 15-bit signed ints with a shared
5right shift, integer prediction with floor division, residual wrapped to
6int16. The coded size is the order-0 entropy of the residual (the ANS bars in
7the app sit 1-2% above this, for the symbol table plus arithmetic loss).
8"""
9import numpy as np
10from scipy.linalg import solve_toeplitz
12SIGMA = 5.0
13RATE = 30000.0
14LOW, HIGH, TAPS = 300.0, 2000.0, 31
15COEFF_PRECISION = 15
18def 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()
30def 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)
38def 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
51def 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
64def 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())
70def 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)')
84if __name__ == '__main__':
85 main()