1"""Semi-analytic cross-check of the entropy rate for
2x iid N(0, sigma^2) -> 31-tap bandpass 300-2000 Hz @ 30 kHz -> round.
4Model the rounding as additive iid U(-1/2, 1/2) noise (variance 1/12) that is
5independent of y. Then z = y + q is a stationary process with spectrum
6S_z(f) = sigma^2 |H(f)|^2 + 1/12, and
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).
16Also validated empirically on a long simulated realization with a
17real-coefficient Wiener predictor.
18"""
19import numpy as np
20from scipy.special import ndtr
22SIGMA = 5.0
23RATE = 30000.0
24LOW, HIGH, TAPS = 300.0, 2000.0, 31
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 windowed_sinc_bandpass(f_lo, f_hi, taps):
40 return windowed_sinc_lowpass(f_hi, taps) - windowed_sinc_lowpass(f_lo, taps)
43def 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)
53def 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())
63def 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')
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}')
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')
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')
90 # what pooling the fractional part costs
91 print(f' gap (frac-part information thrown away) = {Hmix - Hc:.4f} bits')
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')
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')
129if __name__ == '__main__':
130 main()