nx_lpc_synth.nx source
↩ module page · 181 lines · 6280 B
1// nx_lpc_synth.nx -- LPC synthesis filter (Arc 1 Phase A.3).
2// Counterpart to nx_lpc_autocorr + nx_lpc_levinson.
3//
4// Given LPC predictor coefs a[1..p] (Q30) and residual e[0..N-1]
5// (i16 PCM signed), reconstruct the speech signal:
6//
7// s[n] = e[n] + Σ_{k=1..p} a[k] * s[n-k]
8//
9// (Note: prediction was s_pred[n] = -Σ a[k] s[n-k] making residual
10// e[n] = s[n] - s_pred[n] = s[n] + Σ a[k] s[n-k]. Inverting that:
11// s[n] = e[n] - Σ a[k] s[n-k] ... but the standard convention used
12// in Levinson is a[k] satisfying s[n] + Σ a[k] s[n-k] = e[n], so
13// s[n] = e[n] - Σ a[k] s[n-k]. Match our sign convention.)
14//
15// Per [[feedback-bits-up-exceed-never-match]]: integer-only IIR pole
16// filter. Patent-clean foundation for sovereign nx_voice_codec.
17//
18// Inputs:
19// e_pcm: i16 LE residual signal (N * 2 bytes)
20// n: number of samples
21// a_q30: i64 LE LPC coefs (p * 8 bytes) from nx_lpc_levinson
22// p: LPC order
23// s_out: i16 LE reconstructed signal buffer (N * 2 bytes)
24//
25// Caller responsibility: s_out should be zeroed for n < 0 (we use
26// only n >= 0 history; initial taps treated as 0).
27
28import "nx_syscalls_x86_64.nx"
29const NX_MAGIC_32767: i64 = 32767
30const NX_MAGIC_32768: i64 = 32768
31const NX_MAGIC_2048: i64 = 2048
32const NX_MAGIC_4096: i64 = 4096
33const NX_MAGIC_6144: i64 = 6144
34const NX_MAGIC_536870912: i64 = 536870912
35const NX_MAGIC_268435456: i64 = 268435456
36const NX_MAGIC_134217728: i64 = 134217728
37const NX_MAGIC_67108864: i64 = 67108864
38
39const NX_LPC_Q: i64 = 30
40
41// Load i16 LE signed at byte offset
42func _i16_load_le(p: *u8, n: i64) -> i64 {
43 let lo: i64 = p[n * 2]
44 let hi: i64 = p[n * 2 + 1]
45 var v: i64 = lo | (hi << 8)
46 if v >= 0x8000 { v = v - 0x10000 }
47 return v
48}
49
50// Store i16 LE signed at byte offset, saturating to [-32768, 32767]
51func _i16_store_le_sat(p: *u8, n: i64, v: i64) -> i64 {
52 var x: i64 = v
53 if x > NX_MAGIC_32767 { x = NX_MAGIC_32767 }
54 if x < -NX_MAGIC_32768 { x = -NX_MAGIC_32768 }
55 if x < 0 { x = x + 0x10000 }
56 p[n * 2] = x & 0xff
57 p[n * 2 + 1] = (x >> 8) & 0xff
58 return 0
59}
60
61// Load i64 LE signed at byte offset
62func _i64_load_le(p: *u8, n: i64) -> i64 {
63 var v: i64 = 0
64 var i: i64 = 0
65 while i < 8 {
66 v = v | ((p[n * 8 + i]) << (i * 8))
67 i = i + 1
68 }
69 return v
70}
71
72// nx_lpc_synth: reconstruct s[n] = e[n] - Σ a[k] s[n-k]
73// a[k] in Q30; s in raw i16. Need s ≈ e in magnitude; clamp on output.
74func nx_lpc_synth(e_pcm: *u8, n_samples: i64,
75 a_q30: *u8, p_order: i64,
76 s_out: *u8) -> i64 {
77 if p_order < 1 { return -1 }
78 if p_order > 30 { return -1 }
79 var n: i64 = 0
80 while n < n_samples {
81 // pred = Σ_{k=1..p} a[k] * s[n-k] (a in Q30 ⇒ shift back >> 30)
82 var pred: i64 = 0
83 var k: i64 = 1
84 while k <= p_order {
85 let idx: i64 = n - k
86 var s_prev: i64 = 0
87 if idx >= 0 {
88 s_prev = _i16_load_le(s_out, idx)
89 }
90 let a_k: i64 = _i64_load_le(a_q30, k - 1)
91 pred = pred + ((a_k * s_prev) >> NX_LPC_Q)
92 k = k + 1
93 }
94 // s[n] = e[n] - pred
95 let e_n: i64 = _i16_load_le(e_pcm, n)
96 let s_n: i64 = e_n - pred
97 _i16_store_le_sat(s_out, n, s_n)
98 n = n + 1
99 }
100 return 0
101}
102
103// KAT smoke: forward-then-inverse — given a signal s, compute LPC
104// residual e = s + Σ a[k] s[n-k] (the prediction error), then
105// reconstruct s' via nx_lpc_synth, and compare. Should be bit-exact
106// modulo i16 saturation.
107//
108// Scratch layout (all 8-byte aligned):
109// scratch + 0 : s_in (n*2 bytes i16)
110// scratch + 2048 : e_residual (n*2 bytes i16)
111// scratch + 4096 : s_recon (n*2 bytes i16)
112// scratch + 6144 : a_q30 (p*8 bytes i64)
113//
114// Returns # samples that round-trip imperfectly.
115func nx_lpc_synth_roundtrip_test(scratch: *u8) -> i64 {
116 let s_in: *u8 = scratch
117 let e_resid: *u8 = (scratch as i64 + NX_MAGIC_2048) as *u8
118 let s_recon: *u8 = (scratch as i64 + NX_MAGIC_4096) as *u8
119 let a_q30: *u8 = (scratch as i64 + NX_MAGIC_6144) as *u8
120 let n: i64 = 64
121 let p: i64 = 4
122
123 // Generate a smooth signal: s[n] = 100 * sin-ish via integer
124 // approximation (we have no sin; use a polynomial: n*(64-n)*(n-32)/100)
125 var i: i64 = 0
126 while i < n {
127 let t: i64 = i
128 let v: i64 = (t * (64 - t) * (t - 32)) / 100
129 _i16_store_le_sat(s_in, i, v)
130 i = i + 1
131 }
132
133 // Use hand-picked LPC coefs a[1..4] = {0.5, -0.25, 0.125, -0.0625} in Q30
134 let HALF_Q30: i64 = NX_MAGIC_536870912 // 0.5 << 30
135 let Q_NEG_25: i64 = -NX_MAGIC_268435456 // -0.25 << 30
136 let Q_125: i64 = NX_MAGIC_134217728 // 0.125 << 30
137 let Q_NEG_0625: i64 = -NX_MAGIC_67108864 // -0.0625 << 30
138 var b: i64 = 0
139 while b < 8 {
140 a_q30[0 * 8 + b] = (HALF_Q30 >> (b * 8)) & 0xff
141 a_q30[1 * 8 + b] = (Q_NEG_25 >> (b * 8)) & 0xff
142 a_q30[2 * 8 + b] = (Q_125 >> (b * 8)) & 0xff
143 a_q30[3 * 8 + b] = (Q_NEG_0625 >> (b * 8)) & 0xff
144 b = b + 1
145 }
146
147 // Forward: e[n] = s[n] + Σ a[k] * s[n-k] (note + here, since
148 // prediction was s_pred[n] = -Σ a[k] s[n-k] ⇒ e = s + Σ a s).
149 // Wait, the synth formula above is s[n] = e[n] - pred where
150 // pred = Σ a[k] s[n-k]. Inverting: e[n] = s[n] + pred.
151 var nn: i64 = 0
152 while nn < n {
153 var pred: i64 = 0
154 var k: i64 = 1
155 while k <= p {
156 let idx: i64 = nn - k
157 var s_prev: i64 = 0
158 if idx >= 0 { s_prev = _i16_load_le(s_in, idx) }
159 let a_k: i64 = _i64_load_le(a_q30, k - 1)
160 pred = pred + ((a_k * s_prev) >> NX_LPC_Q)
161 k = k + 1
162 }
163 let e_n: i64 = _i16_load_le(s_in, nn) + pred
164 _i16_store_le_sat(e_resid, nn, e_n)
165 nn = nn + 1
166 }
167
168 // Reverse: synth from e using same coefs
169 nx_lpc_synth(e_resid, n, a_q30, p, s_recon)
170
171 // Compare s_in to s_recon
172 var diffs: i64 = 0
173 var j: i64 = 0
174 while j < n {
175 let a: i64 = _i16_load_le(s_in, j)
176 let r: i64 = _i16_load_le(s_recon, j)
177 if a != r { diffs = diffs + 1 }
178 j = j + 1
179 }
180 return diffs
181}