nx_biquad.nx source
↩ module page · 240 lines · 9825 B
1// nx_biquad.nx -- direct-form-II transposed biquad filter in Q14.
2//
3// Cross-cutting Substrate B primitive per
4// nxc2/docs/NISHI_GAME_ENGINE_ROADMAP.md.
5//
6// Biquad = second-order recursive filter section. Building block for
7// every parametric EQ, low-pass, high-pass, band-pass, band-stop,
8// peak, and shelf. Coefficients computed per Bristow-Johnson's
9// "Audio EQ Cookbook" (W3C public-domain reference).
10//
11// This primitive ships the EVALUATION (the per-sample state-update
12// step + block apply); coefficient DESIGN helpers (computing b0..b2,
13// a1..a2 from filter_kind, cutoff, Q) ship as a follow-up primitive
14// once the in-house Q14 sin/cos helpers are factored out. Caller
15// supplies pre-normalized coefficients (a0 = 1 implicit, all others
16// already divided by a0).
17//
18// Transposed Direct Form II:
19//
20// y[n] = b0 * x[n] + s1
21// s1' = b1 * x[n] + s2 - a1 * y[n]
22// s2' = b2 * x[n] - a2 * y[n]
23//
24// Two state vars instead of four. Numerically better-behaved for
25// IIR than Direct Form I; standard DSP-textbook recommendation.
26//
27// Fixed-point: Q14 throughout (matches nx_sound's NX_SOUND_Q).
28// Coefficient * signal products are Q28 internally; divided down to
29// Q14 immediately to keep state vars in-range.
30//
31// Loss audit: integer division (/Q) truncates 14 fractional bits per
32// multiply. For stable filters with |poles| << 1 the error is tiny;
33// for filters near the unit circle the error can accumulate (the same
34// limitation any fixed-point biquad has). Documented at the top of
35// every filter design helper that ships against this primitive.
36//
37// genealogy_id: bristow_johnson_eq_cookbook + smith_julius_iir_filters_textbook
38// lineage_id: nx_biquad_tdf2_q14
39
40// nx_safety_envelope:
41// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
42// sil_target: SIL1
43// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
44// verdict: NOT_YET_EVALUATED
45
46import "nx_syscalls.nx"
47import "nx_tier.nx"
48const NX_MAGIC_12345: i64 = 12345
49const NX_MAGIC_67890: i64 = 67890
50const NX_MAGIC_8000: i64 = 8000
51const NX_MAGIC_4000: i64 = 4000
52const NX_MAGIC_8192: i64 = 8192
53const NX_MAGIC_4096: i64 = 4096
54const NX_MAGIC_2048: i64 = 2048
55
56// ===== Q14 constants ================================================
57// Audio convention -- 1.0 = 16384. Matches nx_sound for downstream
58// composition.
59const NX_BIQUAD_Q: nx_int = 16384
60
61// ===== Sealed-enum filter kinds =====================================
62// Used by future design helpers (nx_biquad_design_lp_q14, etc.).
63// The evaluation primitive itself is kind-agnostic -- caller selects
64// the kind by which coefficient set they pass.
65const NX_BIQUAD_KIND_LPF: nx_int = 0
66const NX_BIQUAD_KIND_HPF: nx_int = 1
67const NX_BIQUAD_KIND_BPF: nx_int = 2
68const NX_BIQUAD_KIND_BSF: nx_int = 3
69const NX_BIQUAD_KIND_PEAK: nx_int = 4
70const NX_BIQUAD_KIND_LSHELF: nx_int = 5
71const NX_BIQUAD_KIND_HSHELF: nx_int = 6
72
73// ===== Filter state struct =========================================
74// Two state variables for transposed direct form II.
75struct BiquadState {
76 s1: nx_int,
77 s2: nx_int,
78}
79
80// ===== Validity predicate ==========================================
81func nx_biquad_kind_is_valid(k: nx_int) -> nx_int {
82 if k == NX_BIQUAD_KIND_LPF { return 1 }
83 if k == NX_BIQUAD_KIND_HPF { return 1 }
84 if k == NX_BIQUAD_KIND_BPF { return 1 }
85 if k == NX_BIQUAD_KIND_BSF { return 1 }
86 if k == NX_BIQUAD_KIND_PEAK { return 1 }
87 if k == NX_BIQUAD_KIND_LSHELF { return 1 }
88 if k == NX_BIQUAD_KIND_HSHELF { return 1 }
89 return 0
90}
91
92// ===== State init ==================================================
93func nx_biquad_init(s: *BiquadState) {
94 s.s1 = 0
95 s.s2 = 0
96}
97
98// ===== Per-sample step =============================================
99// s: state pointer (carries s1, s2 across calls)
100// x_q14: input sample (audio convention: 16384 = 1.0)
101// b0, b1, b2: feedforward coefficients in Q14, pre-normalized by a0
102// a1, a2: feedback coefficients in Q14, pre-normalized by a0
103// (and POSITIVE -- the formula already negates them)
104// returns: output sample y[n] in Q14
105//
106// Q14 arithmetic: each Q14 * Q14 product is Q28; divide by Q14 to
107// return to Q14 amplitude scale.
108func nx_biquad_step(
109 s: *BiquadState,
110 x_q14: nx_int,
111 b0: nx_int,
112 b1: nx_int,
113 b2: nx_int,
114 a1: nx_int,
115 a2: nx_int
116) -> nx_int {
117 // y[n] = b0 * x[n] + s1
118 let b0x: nx_int = (b0 * x_q14) / NX_BIQUAD_Q
119 let y: nx_int = b0x + s.s1
120
121 // s1' = b1 * x[n] + s2 - a1 * y[n]
122 let b1x: nx_int = (b1 * x_q14) / NX_BIQUAD_Q
123 let a1y: nx_int = (a1 * y) / NX_BIQUAD_Q
124 let s1_new: nx_int = b1x + s.s2 - a1y
125
126 // s2' = b2 * x[n] - a2 * y[n]
127 let b2x: nx_int = (b2 * x_q14) / NX_BIQUAD_Q
128 let a2y: nx_int = (a2 * y) / NX_BIQUAD_Q
129 let s2_new: nx_int = b2x - a2y
130
131 s.s1 = s1_new
132 s.s2 = s2_new
133 return y
134}
135
136// ===== Block apply ==================================================
137// Apply the filter to N samples. Updates state across the block.
138// in_buf and out_buf may alias (in-place filtering supported).
139func nx_biquad_block(
140 s: *BiquadState,
141 in_buf: *i64,
142 out_buf: *i64,
143 n: nx_int,
144 b0: nx_int,
145 b1: nx_int,
146 b2: nx_int,
147 a1: nx_int,
148 a2: nx_int
149) {
150 var i: nx_int = 0
151 while i < n {
152 let x: nx_int = in_buf[i]
153 let y: nx_int = nx_biquad_step(s, x, b0, b1, b2, a1, a2)
154 out_buf[i] = y
155 i = i + 1
156 }
157}
158
159// ===== Self-test ====================================================
160func main() -> i64 {
161 let state: *BiquadState = (sys_mmap(2 * NX_SIZEOF_NX_INT)) as *BiquadState
162
163 // T1: nx_biquad_init zeros the state.
164 state.s1 = NX_MAGIC_12345
165 state.s2 = NX_MAGIC_67890
166 nx_biquad_init(state)
167 if state.s1 != 0 { return __syscall(93, 1, 0, 0, 0, 0, 0) }
168 if state.s2 != 0 { return __syscall(93, 2, 0, 0, 0, 0, 0) }
169
170 // T2: Identity filter (b0=Q, others=0). Output equals input
171 // sample-for-sample.
172 let y_id_a: nx_int = nx_biquad_step(state, NX_MAGIC_8000, NX_BIQUAD_Q, 0, 0, 0, 0)
173 if y_id_a != NX_MAGIC_8000 { return __syscall(93, 3, 0, 0, 0, 0, 0) }
174 let y_id_b: nx_int = nx_biquad_step(state, 0 - NX_MAGIC_4000, NX_BIQUAD_Q, 0, 0, 0, 0)
175 if y_id_b != (0 - NX_MAGIC_4000) { return __syscall(93, 4, 0, 0, 0, 0, 0) }
176 let y_id_c: nx_int = nx_biquad_step(state, 0, NX_BIQUAD_Q, 0, 0, 0, 0)
177 if y_id_c != 0 { return __syscall(93, 5, 0, 0, 0, 0, 0) }
178
179 // T3: Pure-FIR impulse response. b0=Q/2, b1=Q/4, b2=Q/8, a1=a2=0.
180 // For input x=[Q, 0, 0, 0, ...] the output should be the tap
181 // sequence: y = [Q/2, Q/4, Q/8, 0, 0, ...] -- the (b0, b1, b2)
182 // coefficients themselves.
183 nx_biquad_init(state)
184 let b0: nx_int = NX_BIQUAD_Q / 2 // NX_MAGIC_8192
185 let b1: nx_int = NX_BIQUAD_Q / 4 // NX_MAGIC_4096
186 let b2: nx_int = NX_BIQUAD_Q / 8 // NX_MAGIC_2048
187
188 let y0: nx_int = nx_biquad_step(state, NX_BIQUAD_Q, b0, b1, b2, 0, 0)
189 if y0 != NX_MAGIC_8192 { return __syscall(93, 6, 0, 0, 0, 0, 0) }
190 let y1: nx_int = nx_biquad_step(state, 0, b0, b1, b2, 0, 0)
191 if y1 != NX_MAGIC_4096 { return __syscall(93, 7, 0, 0, 0, 0, 0) }
192 let y2: nx_int = nx_biquad_step(state, 0, b0, b1, b2, 0, 0)
193 if y2 != NX_MAGIC_2048 { return __syscall(93, 8, 0, 0, 0, 0, 0) }
194 let y3: nx_int = nx_biquad_step(state, 0, b0, b1, b2, 0, 0)
195 if y3 != 0 { return __syscall(93, 9, 0, 0, 0, 0, 0) }
196 let y4: nx_int = nx_biquad_step(state, 0, b0, b1, b2, 0, 0)
197 if y4 != 0 { return __syscall(93, 10, 0, 0, 0, 0, 0) }
198
199 // T4: Zero-input + zero-state -> zero-output (forever).
200 nx_biquad_init(state)
201 var k: nx_int = 0
202 while k < 8 {
203 let y: nx_int = nx_biquad_step(state, 0, NX_BIQUAD_Q / 2, NX_BIQUAD_Q / 4, NX_BIQUAD_Q / 8, 0, 0)
204 if y != 0 { return __syscall(93, 20, 0, 0, 0, 0, 0) }
205 k = k + 1
206 }
207
208 // T5: Block apply -- equivalent to N individual steps.
209 nx_biquad_init(state)
210 let in_buf: *i64 = (sys_mmap(8 * NX_SIZEOF_NX_INT)) as *i64
211 let out_buf: *i64 = (sys_mmap(8 * NX_SIZEOF_NX_INT)) as *i64
212 in_buf[0] = NX_BIQUAD_Q
213 in_buf[1] = 0
214 in_buf[2] = 0
215 in_buf[3] = 0
216 in_buf[4] = 0
217 in_buf[5] = 0
218 in_buf[6] = 0
219 in_buf[7] = 0
220 nx_biquad_block(state, in_buf, out_buf, 8, b0, b1, b2, 0, 0)
221 if out_buf[0] != NX_MAGIC_8192 { return __syscall(93, 30, 0, 0, 0, 0, 0) }
222 if out_buf[1] != NX_MAGIC_4096 { return __syscall(93, 31, 0, 0, 0, 0, 0) }
223 if out_buf[2] != NX_MAGIC_2048 { return __syscall(93, 32, 0, 0, 0, 0, 0) }
224 if out_buf[3] != 0 { return __syscall(93, 33, 0, 0, 0, 0, 0) }
225 if out_buf[7] != 0 { return __syscall(93, 34, 0, 0, 0, 0, 0) }
226
227 // T6: Filter kind validity predicate.
228 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_LPF) != 1 { return __syscall(93, 40, 0, 0, 0, 0, 0) }
229 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_HPF) != 1 { return __syscall(93, 41, 0, 0, 0, 0, 0) }
230 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_BPF) != 1 { return __syscall(93, 42, 0, 0, 0, 0, 0) }
231 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_BSF) != 1 { return __syscall(93, 43, 0, 0, 0, 0, 0) }
232 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_PEAK) != 1 { return __syscall(93, 44, 0, 0, 0, 0, 0) }
233 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_LSHELF) != 1 { return __syscall(93, 45, 0, 0, 0, 0, 0) }
234 if nx_biquad_kind_is_valid(NX_BIQUAD_KIND_HSHELF) != 1 { return __syscall(93, 46, 0, 0, 0, 0, 0) }
235 if nx_biquad_kind_is_valid(0 - 1) != 0 { return __syscall(93, 47, 0, 0, 0, 0, 0) }
236 if nx_biquad_kind_is_valid(7) != 0 { return __syscall(93, 48, 0, 0, 0, 0, 0) }
237 if nx_biquad_kind_is_valid(99) != 0 { return __syscall(93, 49, 0, 0, 0, 0, 0) }
238
239 return 0
240}