1/**
2 * Linear predictive coding, the generalisation of delta coding.
3 *
4 * Delta codes `x[n] - x[n-1]`. LPC fits an order-p predictor
5 *
6 * xhat[n] = sum_{k=1..p} a_k x[n-k]
7 *
8 * to the data and codes the residual `x[n] - xhat[n]` instead. Delta is the
9 * special case p = 1, a_1 = 1.
10 *
11 * For the residual to be *losslessly* invertible the decoder has to reproduce
12 * `xhat` bit-for-bit, so — as in FLAC and Shorten — the coefficients are
13 * quantised to integers with a shared right shift and the prediction is done
14 * in exact integer arithmetic:
15 *
16 * xhat[n] = floor( (sum_k q_k x[n-k]) / 2^shift )
17 *
18 * `sum_k q_k x[n-k]` reaches about 2^35 for the orders and widths used here,
19 * which is exact in a double, so no BigInt is needed.
20 */
22/** 15-bit signed coefficients, as FLAC's default. */
23const COEFF_PRECISION = 15
25export interface LpcModel {
26 order: number
27 shift: number
28 coeffs: Int32Array
29}
31/** Autocorrelation r[0..maxLag] of the signal. */
32function autocorrelation(x: Int16Array, maxLag: number): Float64Array {
33 const r = new Float64Array(maxLag + 1)
34 for (let lag = 0; lag <= maxLag; lag++) {
35 let sum = 0
36 for (let i = lag; i < x.length; i++) sum += x[i] * x[i - lag]
37 r[lag] = sum
38 }
39 return r
40}
42/**
43 * Levinson-Durbin: the order-p predictor minimising the mean squared
44 * prediction error, from the autocorrelation. Returns a[1..p] in a[0..p-1].
45 */
46function levinson(r: Float64Array, order: number): Float64Array | null {
47 if (r[0] <= 0) return null
48 const a = new Float64Array(order)
49 const tmp = new Float64Array(order)
50 let error = r[0]
52 for (let i = 0; i < order; i++) {
53 let acc = r[i + 1]
54 for (let j = 0; j < i; j++) acc -= a[j] * r[i - j]
55 const k = acc / error
56 if (!Number.isFinite(k)) return null
58 a[i] = k
59 for (let j = 0; j < i; j++) tmp[j] = a[j] - k * a[i - 1 - j]
60 for (let j = 0; j < i; j++) a[j] = tmp[j]
62 error *= 1 - k * k
63 if (error <= 0) return null
64 }
65 return a
66}
68/** Scale the real coefficients onto integers sharing one right shift. */
69function quantiseCoefficients(a: Float64Array): LpcModel | null {
70 let maxAbs = 0
71 for (const v of a) maxAbs = Math.max(maxAbs, Math.abs(v))
72 if (maxAbs === 0 || !Number.isFinite(maxAbs)) return null
74 let shift = COEFF_PRECISION - 1 - Math.floor(Math.log2(maxAbs)) - 1
75 shift = Math.max(0, Math.min(15, shift))
77 const limit = 2 ** (COEFF_PRECISION - 1)
78 const coeffs = new Int32Array(a.length)
79 for (let i = 0; i < a.length; i++) {
80 coeffs[i] = Math.max(-limit, Math.min(limit - 1, Math.round(a[i] * 2 ** shift)))
81 }
82 return { order: a.length, shift, coeffs }
83}
85export function fitLpc(x: Int16Array, order: number): LpcModel | null {
86 const a = levinson(autocorrelation(x, order), order)
87 return a ? quantiseCoefficients(a) : null
88}
90/** The unquantised order-p predictor, for coders that keep full precision. */
91export function fitRealCoeffs(x: Int16Array, order: number): Float64Array | null {
92 return levinson(autocorrelation(x, order), order)
93}
95/**
96 * Prediction residual. The first `order` entries are the samples themselves
97 * (the warm-up the decoder needs before it can predict). Differences are taken
98 * in int16 with wraparound, so the transform is exactly invertible.
99 */
100export function lpcResidual(x: Int16Array, model: LpcModel): Int16Array {
101 const { order, shift, coeffs } = model
102 const scale = 2 ** shift
103 const e = new Int16Array(x.length)
104 for (let n = 0; n < Math.min(order, x.length); n++) e[n] = x[n]
105 for (let n = order; n < x.length; n++) {
106 let sum = 0
107 for (let k = 0; k < order; k++) sum += coeffs[k] * x[n - 1 - k]
108 e[n] = ((x[n] - Math.floor(sum / scale)) << 16) >> 16
109 }
110 return e
111}
113/** Rebuild the signal from its residual — the exact inverse of `lpcResidual`. */
114export function lpcRestore(e: Int16Array, model: LpcModel): Int16Array {
115 const { order, shift, coeffs } = model
116 const scale = 2 ** shift
117 const x = new Int16Array(e.length)
118 for (let n = 0; n < Math.min(order, e.length); n++) x[n] = e[n]
119 for (let n = order; n < e.length; n++) {
120 let sum = 0
121 for (let k = 0; k < order; k++) sum += coeffs[k] * x[n - 1 - k]
122 x[n] = ((e[n] + Math.floor(sum / scale)) << 16) >> 16
123 }
124 return x
125}
127/** Bytes the decoder needs for the model: order, shift, and the coefficients. */
128export function modelSize(model: LpcModel): number {
129 return 2 + 2 * model.order
130}