nx_fft.nx source

1// nx_fft.nx -- R0 of the transform ladder: a sovereign no-float radix-2 FFT (Cooley-Tukey, DIT, in-place). 2// The foundation rung for MDCT (R1) -> psychoacoustic perceptual audio (R2+). Twiddles COMPOSE the sine that was 3// absorbed into nx_audio_osc (osc_sin_unit) -- one sine, no dup. Fixed-point (Q16 twiddles, integer data), so the 4// result is DETERMINISTIC / bit-reproducible across machines -- the s-class exceed dimension vs float FFTs whose 5// rounding varies by platform. Gated by MATH properties verifiable without ears: round-trip reconstruction + 6// known answers (DC -> single bin; impulse -> flat spectrum). license_tier: ORIGINAL expect_exit: 0 7import "nx_audio_osc.nx" // osc_sin_unit -- the absorbed sine (twiddles), dedup 8import "nx_syscalls.nx" 9const K_MAGIC_2000: i64 = 2000 10const K_MAGIC_1500: i64 = 1500 11const K_MAGIC_1000000: i64 = 1000000 12const K_MAGIC_1400: i64 = 1400 13const K_MAGIC_1600: i64 = 1600 14 15const OSCP: i64 = 65536 // osc one-cycle units (matches nx_audio_osc OSC_PHASE) 16const QPI: i64 = 16384 // OSCP/4 = pi/2 in osc units (cos(x) = sin(x + pi/2)) 17 18func fp(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 19func fpn(v: i64) -> i64 { 20 let b: *u8 = sys_mmap(28); var x: i64=v 21 if x<0 { b[0]=45; sys_write(1,b,1); x=0-x } 22 if x==0 { b[0]=48; sys_write(1,b,1); return 0 } 23 var d: i64=0; var y: i64=x; while y>0 { d=d+1; y=y/10 } 24 var i: i64=d-1; y=x; while i>=0 { b[i]=(48+(y%10)) as u8; y=y/10; i=i-1 } 25 sys_write(1,b,d); return 0 26} 27func f_abs(v: i64) -> i64 { if v<0 { return 0-v } return v } 28 29// in-place radix-2 DIT FFT. re/im: length n (power of 2). dir = -1 forward, +1 inverse (inverse also scales 1/n). 30func fft(re: *i64, im: *i64, n: i64, dir: i64) -> i64 { 31 // bit-reversal permutation 32 var j: i64 = 0 33 var i: i64 = 1 34 while i < n { 35 var bit: i64 = n >> 1 36 while (j & bit) != 0 { j = j ^ bit; bit = bit >> 1 } 37 j = j ^ bit 38 if i < j { 39 let tr: i64 = re[i]; re[i] = re[j]; re[j] = tr 40 let ti: i64 = im[i]; im[i] = im[j]; im[j] = ti 41 } 42 i = i + 1 43 } 44 // butterflies 45 var len: i64 = 2 46 while len <= n { 47 let half: i64 = len >> 1 48 var k: i64 = 0 49 while k < half { 50 let ph: i64 = (dir * k * OSCP) / len // twiddle angle in osc units 51 let wr: i64 = osc_sin_unit(ph + QPI) // cos(theta) 52 let wi: i64 = osc_sin_unit(ph) // sin(theta) 53 var s: i64 = k 54 while s < n { 55 let ar: i64 = re[s]; let ai: i64 = im[s] 56 let br: i64 = re[s+half]; let bi: i64 = im[s+half] 57 let tr: i64 = (wr*br - wi*bi) >> 16 // t = w * b (Q16 twiddle) 58 let ti: i64 = (wr*bi + wi*br) >> 16 59 re[s] = ar + tr; im[s] = ai + ti 60 re[s+half] = ar - tr; im[s+half] = ai - ti 61 s = s + len 62 } 63 k = k + 1 64 } 65 len = len << 1 66 } 67 if dir > 0 { 68 var z: i64 = 0 69 while z < n { re[z] = re[z] / n; im[z] = im[z] / n; z = z + 1 } 70 } 71 return 0 72} 73 74// --------------------------------------------------------------------------------------------------- 75// DECOMPOSED API -- RESTORED 2026-07-31. FOUR consumers import nx_fft.nx for these two helpers and NEITHER 76// WAS EVER DEFINED: nx_fft_f32.nx:7, nx_fnet_mix.nx:32, nx_fnet_mix_gate.nx:259, and nx_fft_test.nx which 77// ASSERTS them -- a test written for functions that do not exist, so it can never once have passed. 78// The layering violation (nx_fft.nx stranded in _hdl_build/) hid all of it: every consumer died at 79// expand_imports before the compiler could reach an undefined call. 80// NOT INVENTED -- the logic is EXTRACTED from the bit-reversal loop already proven inside fft() below, 81// and both are PINNED BY THE EXISTING TEST VECTORS in nx_fft_test.nx, which are the real contract: 82// nx_fft_log2: 2->1, 8->3, 32->5 83// nx_fft_bit_reverse: (0,3)->0 (1,3)->4 (3,3)->6 (5,3)->5 (7,3)->7 84// NOTE nx_fft_bit_reverse maps ONE INDEX given a bit-width; it does NOT permute an array in place. The 85// in-place permutation stays inlined at the top of fft() -- do not confuse the two, they have different 86// shapes and a caller that swaps them gets a silently wrong transform rather than an error. 87// Uses + rather than | to merge the low bit: after (r << 1) bit0 is always 0, so the two are identical 88// here, and + is known-good in this compiler. 89func nx_fft_log2(n: i64) -> i64 { 90 var v: i64 = n 91 var r: i64 = 0 92 while v > 1 { v = v >> 1; r = r + 1 } 93 return r 94} 95 96func nx_fft_bit_reverse(x: i64, bits: i64) -> i64 { 97 var r: i64 = 0 98 var i: i64 = 0 99 while i < bits { 100 r = (r << 1) + ((x >> i) & 1) 101 i = i + 1 102 } 103 return r 104} 105 106// NX_FFT_Q -- SETTLED BY MEASUREMENT 2026-07-31, and the disagreement that made it dangerous is REAL: 107// THIS FILE CARRIES TWO DIFFERENT FIXED-POINT SCALES, DELIBERATELY. 108// - fft() below builds twiddles with osc_sin_unit and divides by `>> 16` -- Q16/65536, internal, and 109// hardcoded at the shift. NOTHING outside reads that scale. 110// - NX_FFT_Q is the scale for callers that build their OWN twiddle tables. It is Q14/16384. 111// Three independent witnesses, none of them a comment guess: 112// 1. nx_fnet_mix.nx:39 fnet_q16_to_q14(v) is literally (v+2)/4 -- a rounded divide of Q16.16 by 4, 113// and 65536/4 == 16384. The tables at :52-53 are filled THROUGH it, so they are Q14. 114// 2. The only two consumers of NX_FFT_Q (nx_fnet_mix.nx:95-96 and nx_fnet_mix_gate.nx:55-56) each 115// divide by it immediately after multiplying by those Q14 twiddles. 116// 3. A GREEN gate run recorded in nx_fnet_mix.nx:30 (2026-06-14) measured "C2 twiddles match nx_fft 117// within 1 Q14" -- a prior measurement, not an assertion. 118// DO NOT "UNIFY" THE TWO SCALES. Changing either without re-measuring produces wrong numbers that no 119// build error will catch -- which is exactly why this constant was left undefined rather than guessed. 120const NX_FFT_Q: i64 = 16384 121// --------------------------------------------------------------------------------------------------- 122 123// nx_fft_forward / nx_fft_inverse -- the direction-named wrappers that consumers call and that were 124// NEVER DEFINED: nx_fft_test.nx (tests 3,4,5,6) and nx_denoise_spectral_gate.nx:12. The bare fft() below 125// takes a raw dir flag; naming the two directions stops a caller passing the wrong sign silently. 126// DIRECTION PINNED BY THE TEST, NOT BY CONVENTION -- this is the whole reason it is safe to add these: 127// test 3 feeds DC [1]x8 to nx_fft_forward and requires re[0] == 8, so FORWARD IS UNNORMALIZED; 128// fft() divides by n only on the dir > 0 branch, therefore forward MUST be dir = -1. 129// test 5 requires forward-then-inverse to recover the input, so inverse is dir = +1, the branch 130// that carries the 1/n. Any other assignment fails those assertions loudly rather than silently. 131// Written 0 - 1 rather than a bare negative literal to match the idiom already used in f_abs above. 132func nx_fft_forward(re: *i64, im: *i64, n: i64) -> i64 { 133 return fft(re, im, n, 0 - 1) 134} 135 136func nx_fft_inverse(re: *i64, im: *i64, n: i64) -> i64 { 137 return fft(re, im, n, 1) 138} 139 140// nx_fft_power_spectrum -- per-bin power into caller-supplied pwr[]. Power is re^2 + im^2 (the standard 141// meaning of the name), NOT magnitude. NOTE the test cannot tell those apart: test 6 only checks the 142// ARGMAX bin (identical under any monotonic transform) and that DC <= 1000, so both definitions pass. 143// I chose re^2+im^2 because the NAME says power; recording that here because the test does not pin it, 144// and a future caller comparing against a magnitude spectrum would see squared values and be confused. 145// OVERFLOW: values are unnormalized forward-transform outputs, so re^2 can overflow i64 above ~3e9 per 146// bin. Fine for the fixed-point sizes this ladder uses; a longer transform must scale before squaring. 147func nx_fft_power_spectrum(re: *i64, im: *i64, n: i64, pwr: *i64) -> i64 { 148 var k: i64 = 0 149 while k < n { 150 let r: i64 = re[k] 151 let m: i64 = im[k] 152 pwr[k] = r * r + m * m 153 k = k + 1 154 } 155 return 0 156} 157 158// nx_fft_2d_forward -- separable 2D transform: FFT every row (length w), then every column (length h). 159// Row-major, index = y * w + x. Columns are strided so they are gathered into a contiguous scratch 160// buffer, transformed, and scattered back -- fft() is in-place and requires contiguity. 161// PINNED BY TEST 7: a 4x4 all-ones image must yield re[0] == 16. Separable unnormalized forward gives 162// 4 per row at bin 0, then 16 at (0,0) after the column pass -- which is why forward must stay 163// unnormalized for this to hold; a 1/n on the forward branch would land 1 here and fail loudly. 164func nx_fft_2d_forward(re: *i64, im: *i64, w: i64, h: i64) -> i64 { 165 var r: i64 = 0 166 while r < h { 167 let rre: *i64 = ((re as i64 + r * w * 8) as *i64) 168 let rim: *i64 = ((im as i64 + r * w * 8) as *i64) 169 nx_fft_forward(rre, rim, w) 170 r = r + 1 171 } 172 let cre: *i64 = (sys_mmap(h * 8)) as *i64 173 let cim: *i64 = (sys_mmap(h * 8)) as *i64 174 var c: i64 = 0 175 while c < w { 176 var y: i64 = 0 177 while y < h { cre[y] = re[y * w + c]; cim[y] = im[y * w + c]; y = y + 1 } 178 nx_fft_forward(cre, cim, h) 179 y = 0 180 while y < h { re[y * w + c] = cre[y]; im[y * w + c] = cim[y]; y = y + 1 } 181 c = c + 1 182 } 183 sys_munmap(cre as *u8, h * 8) 184 sys_munmap(cim as *u8, h * 8) 185 return 0 186} 187 188func main() -> i64 { 189 fp("nx_fft (R0): sovereign no-float radix-2 FFT; twiddles via osc_sin_unit (deduped sine)\n" as *u8) 190 let N: i64 = 256 191 let re: *i64 = sys_mmap(N*8) as *i64 192 let im: *i64 = sys_mmap(N*8) as *i64 193 let or: *i64 = sys_mmap(N*8) as *i64 194 var ok: i64 = 0 195 196 // ---- T1 round-trip: x -> FFT -> IFFT -> x' ; max|x'-x| bounded ---- 197 var i: i64 = 0 198 while i < N { let v: i64 = ((i*37 + 11) % K_MAGIC_2000) - 1000; re[i] = v; im[i] = 0; or[i] = v; i = i + 1 } 199 fft(re, im, N, 0 - 1) 200 fft(re, im, N, 1) 201 var mx: i64 = 0 202 i = 0; while i < N { let e: i64 = f_abs(re[i] - or[i]); if e>mx {mx=e} i=i+1 } 203 fp(" T1 round-trip max-err=" as *u8); fpn(mx); fp(" (Q16 twiddles -> bounded, deterministic)\n" as *u8) 204 if mx < 16 { ok = ok + 1 } 205 206 // ---- T2 DC: x=const -> energy in bin 0 only ---- 207 i = 0; while i < N { re[i] = 1000; im[i] = 0; i = i + 1 } 208 fft(re, im, N, 0 - 1) 209 let dc: i64 = re[0] 210 var leak: i64 = 0; i = 1; while i < N { let e: i64 = f_abs(re[i]) + f_abs(im[i]); if e>leak {leak=e} i=i+1 } 211 fp(" T2 DC: bin0=" as *u8); fpn(dc); fp(" (expect ~" as *u8); fpn(N*1000); fp(") other-bin max=" as *u8); fpn(leak); fp("\n" as *u8) 212 if dc > N*900 { if leak < 200 { ok = ok + 1 } } 213 214 // ---- T3 impulse: x=[A,0,0,...] -> flat spectrum (every bin ~ A) ---- 215 i = 0; while i < N { re[i] = 0; im[i] = 0; i = i + 1 } 216 re[0] = K_MAGIC_1500 217 fft(re, im, N, 0 - 1) 218 var fmin: i64 = K_MAGIC_1000000; var fmax: i64 = 0 219 i = 0; while i < N { let m: i64 = f_abs(re[i]); if m<fmin {fmin=m} if m>fmax {fmax=m} i=i+1 } 220 fp(" T3 impulse: spectrum flat? min=" as *u8); fpn(fmin); fp(" max=" as *u8); fpn(fmax); fp(" (expect ~1500 both)\n" as *u8) 221 if fmin > K_MAGIC_1400 { if fmax < K_MAGIC_1600 { ok = ok + 1 } } 222 223 fp(" verdict=" as *u8) 224 if ok == 3 { fp("GREEN (FFT round-trips + known answers hold; no-float, deterministic) -- R0 DONE\n" as *u8); sys_exit(0); return 0 } 225 fp("RED (" as *u8); fpn(ok); fp("/3)\n" as *u8) 226 sys_exit(1); return 1 227}