9d25feaKeep the Python evidence for the entropy-rate gapJeremy Magland 1"""How long must the conditioning window M be?
3The latent y = h*x is exactly (L-1)-dependent, so a naive reading says the
4past only needs to be ~L samples. That is wrong for prediction: y is an
5MA(L-1) process, and the optimal *predictor* of an MA process from its own
6past is AR(infinity) — the whitening filter is 1/H(z), whose impulse response
7decays at a rate set by how close H's zeros sit to the unit circle. For a
8narrowband filter those zeros are very close, so the memory that matters runs
9many multiples of L.
11What keeps it finite is the rounding: z = round(y) has spectrum
12S_z = sigma^2 |H|^2 + 1/12, which never reaches zero, so its whitening filter
13always converges. Deep nulls (large sigma, narrow band, many taps) mean slow
14decay; the 1/12 floor sets where the decay finally bites.
16This computes, from that spectral model, the prediction order (= the past
17length that matters) needed to get within 0.01 bits/sample of the infinite-
18past limit — the M the estimator needs.
19"""
20import numpy as np
21from scipy.special import ndtr
23RATE = 30000.0
24LOW, HIGH = 300.0, 2000.0
27def 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()
39def bandpass(taps):
40 return windowed_sinc_lowpass(HIGH / RATE, taps) - windowed_sinc_lowpass(LOW / RATE, taps)
43def 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
64def 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())
75def 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
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
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
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
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}')
116if __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()