nx_lpc.nx source
↩ module page · 180 lines · 6463 B
1// nx_lpc.nx -- Linear Predictive Coding (LPC) analysis via the
2// Levinson-Durbin recursion. Pure integer Q15 fixed-point; no
3// floats; tier-0 ready (runs on MCUs without an FPU).
4//
5// Phase 4a of the Nishi voice codec (docs/NISHI_COMMS_ROADMAP.md).
6// LPC is the spectral-envelope estimator the codec wraps around the
7// excitation signal. Order 10 covers 8 kHz speech; order 16 covers
8// 16 kHz; order 20 for 24 kHz.
9//
10// Q15 convention:
11// - signed coefficients in [-1.0, +1.0) map to int16 in
12// [-32768, +32767]
13// - multiply: (a * b) >> 15 keeps the Q15 invariant
14// - the LPC reflection coefficients k[i] are bounded |k| < 1 in
15// a stable filter; the recursion's stability check is k^2 < 1
16//
17// We store everything as i64 to give headroom for autocorrelation
18// sums (which can be O(N) where N is samples in a frame, ~320 at 16
19// kHz / 20 ms). Output coefficients are bounded to Q15 range.
20//
21// Algorithm (standard Levinson-Durbin):
22// R[k] = autocorr(s, k) for k in 0..M
23// e[0] = R[0]
24// for i in 1..M:
25// k_i = -(R[i] + sum(a[j] * R[i-j] for j in 1..i-1)) / e[i-1]
26// a[i] = k_i
27// for j in 1..i-1: a'[j] = a[j] + k_i * a[i-j]
28// a = a'
29// e[i] = (1 - k_i^2) * e[i-1]
30//
31// genealogy_id: levinson_1947 + durbin_1960 + rfc_3951_ilbc +
32// opus_silk_lpc_tables
33// lineage_id: nishi_voice_lpc_q15
34
35// nx_safety_envelope:
36// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
37// sil_target: SIL1
38// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
39// verdict: NOT_YET_EVALUATED
40
41import "nx_syscalls.nx"
42
43// Sealed verdict per analysis run.
44const NX_LPC_VERDICT_UNKNOWN: i64 = 0
45const NX_LPC_VERDICT_OK: i64 = 1
46const NX_LPC_VERDICT_BAD_ORDER: i64 = 2
47const NX_LPC_VERDICT_ZERO_ENERGY: i64 = 3
48const NX_LPC_VERDICT_UNSTABLE: i64 = 4
49const NX_LPC_VERDICT_N: i64 = 5
50
51const NX_LPC_Q15_ONE: i64 = 32768 // 1.0 in Q15 (clamped to 32767 in storage)
52const NX_LPC_Q15_MASK: i64 = 0x7fff
53
54// Q15 multiply with rounding: (a * b + 0x4000) >> 15.
55func _q15_mul(a: i64, b: i64) -> i64 {
56 return (a * b + 16384) >> 15
57}
58
59// Compute autocorrelation R[0..order] of samples[0..n). samples
60// are i64 (caller decides sign/scaling). Result is stored in
61// R[0..order]. R[k] is sum over i in [0..n-k) of samples[i] *
62// samples[i+k].
63//
64// Returns R[0] (= signal energy) so caller can short-circuit on
65// zero-energy frames.
66func nx_lpc_autocorrelate(samples: *i64, n: i64, R: *i64, order: i64) -> i64 {
67 var k: i64 = 0
68 while k <= order {
69 var acc: i64 = 0
70 var i: i64 = 0
71 while i + k < n {
72 acc = acc + samples[i] * samples[i + k]
73 i = i + 1
74 }
75 R[k] = acc
76 k = k + 1
77 }
78 return R[0]
79}
80
81// Solve the LPC system via Levinson-Durbin recursion. Returns the
82// sealed verdict. On NX_LPC_VERDICT_OK, `a_q15[1..order]` holds the
83// LPC coefficients in Q15 (a_q15[0] is by convention 1.0 = 32767;
84// we don't write it). `reflection_q15[1..order]` holds the
85// reflection coefficients k_i in Q15.
86//
87// Caller pre-allocates:
88// a_q15: (order+1) i64 slots
89// reflection_q15: (order+1) i64 slots
90// a_tmp: (order+1) i64 slots (scratch)
91// R: (order+1) i64 slots (the autocorrelation, filled
92// by nx_lpc_autocorrelate)
93func nx_lpc_levinson(
94 R: *i64, order: i64,
95 a_q15: *i64, reflection_q15: *i64,
96 a_tmp: *i64
97) -> i64 {
98 if order < 1 { return NX_LPC_VERDICT_BAD_ORDER }
99 if order > 32 { return NX_LPC_VERDICT_BAD_ORDER }
100 if R[0] == 0 { return NX_LPC_VERDICT_ZERO_ENERGY }
101
102 // Working in "scaled Q15". R values are i64 (no fixed point);
103 // we normalise on the fly when computing k_i = -num / e.
104 a_q15[0] = NX_LPC_Q15_MASK // 1.0 in Q15 (clamped to 32767)
105 var ei: i64 = R[0] // e[0]
106
107 var i: i64 = 1
108 while i <= order {
109 // num = R[i] + sum_{j=1..i-1} a[j] * R[i-j]
110 var num: i64 = R[i]
111 var j: i64 = 1
112 while j < i {
113 // a[j] is in Q15 -- divide by 32768 at multiply time
114 num = num + (a_q15[j] * R[i - j] >> 15)
115 j = j + 1
116 }
117 // k_i = -num / e[i-1]. Scale to Q15 by shifting num left.
118 // Stability requires |k_i| < 1; here we clamp.
119 var k_q15: i64 = 0
120 if ei != 0 {
121 // Compute -num * 32768 / ei carefully to avoid overflow.
122 // num and ei are O(R[0]) which is O(N * sample^2).
123 // Keep as i64; final result clipped to int16 range.
124 k_q15 = -((num << 15) / ei)
125 }
126 if k_q15 > NX_LPC_Q15_MASK { k_q15 = NX_LPC_Q15_MASK }
127 if k_q15 < -NX_LPC_Q15_MASK { k_q15 = -NX_LPC_Q15_MASK }
128 reflection_q15[i] = k_q15
129
130 // Update a[]: a_new[j] = a[j] + k_i * a[i-j] for j in 1..i-1;
131 // a_new[i] = k_i.
132 var j2: i64 = 1
133 while j2 < i {
134 a_tmp[j2] = a_q15[j2] + _q15_mul(k_q15, a_q15[i - j2])
135 j2 = j2 + 1
136 }
137 a_tmp[i] = k_q15
138
139 // Commit a_tmp -> a.
140 var j3: i64 = 1
141 while j3 <= i {
142 a_q15[j3] = a_tmp[j3]
143 j3 = j3 + 1
144 }
145
146 // e[i] = (1 - k_i^2) * e[i-1]. Q15 * Q15 -> Q15 via >> 15.
147 let k_sq: i64 = _q15_mul(k_q15, k_q15)
148 let one_minus_k_sq_q15: i64 = NX_LPC_Q15_ONE - k_sq
149 // Update e: ei * (1 - k^2) in Q15 means (ei * one_minus) >> 15.
150 ei = (ei * one_minus_k_sq_q15) >> 15
151 if ei <= 0 {
152 // Energy went negative or zero -- numerical instability,
153 // probably a degenerate input. Caller can fall back to
154 // a flat spectrum or higher precision.
155 return NX_LPC_VERDICT_UNSTABLE
156 }
157 i = i + 1
158 }
159 return NX_LPC_VERDICT_OK
160}
161
162// Convenience: do both autocorrelate + Levinson in one call.
163// Same out-buffers as nx_lpc_levinson. Caller also provides a R
164// scratch buffer (order+1 slots).
165func nx_lpc_analyze(
166 samples: *i64, n: i64, order: i64,
167 R_scratch: *i64,
168 a_q15: *i64, reflection_q15: *i64,
169 a_tmp: *i64
170) -> i64 {
171 nx_lpc_autocorrelate(samples, n, R_scratch, order)
172 return nx_lpc_levinson(R_scratch, order, a_q15, reflection_q15, a_tmp)
173}
174
175// Sealed-enum validity gate.
176func nx_lpc_verdict_is_valid(v: i64) -> i64 {
177 if v < 0 { return 0 }
178 if v >= NX_LPC_VERDICT_N { return 0 }
179 return 1
180}