_fft256_q15_throughput_bench.nx source

1// _fft256_q15_throughput_bench.nx -- 256-point Q15 FFT throughput. 2// 3// THE primary W1 (automotive sensor / leak detection) workload 4// baseline per WORKLOAD_TARGETS.md. Acoustic leak detection, brake 5// pressure analysis, battery thermal classification, motor vibration 6// FFT-fingerprinting -- every one of these starts with a real-input 7// FFT in the 256..2048-point range. 8// 9// Q15 fixed-point (16-bit signed in -1.0..+1.0 range): the industry 10// default for low-power DSP per ARM CMSIS-DSP, TI C5000, and 11// every automotive sensor MCU. NishiLang has no native fp, so Q15 12// is also a natural fit for the substrate. 13// 14// Algorithm: Cooley-Tukey radix-2 decimation-in-time (DIT). In-place 15// 256-point complex FFT. Twiddle factors generated iteratively at 16// startup via the angle-addition formula: 17// cos(t+dt) = cos(t)*cos(dt) - sin(t)*sin(dt) 18// sin(t+dt) = sin(t)*cos(dt) + cos(t)*sin(dt) 19// so we hard-code ONE step (cos/sin of 2*pi/256) and grow the table. 20// 21// Bench: 1000 iterations of 256-point FFT. Reports FFTs/sec + 22// per-FFT latency in microseconds. Reference targets: 23// ARM CMSIS-DSP (Cortex-M4 @ 100 MHz, SIMD) ~50K FFTs/sec 24// FFTW double precision (modern x86_64) ~100K FFTs/sec 25// Pure-software (this bench) ~1-10K FFTs/sec est 26// CUDA cuFFT batched ~10M FFTs/sec 27// 28// Silicon-feedback consumption: rv64im_min_hot_report.nx attribution 29// over this bench identifies which functions dominate (likely the 30// butterfly inner loop's multiply + accumulate). PROP_MAC_FUSED 31// silicon candidate from rv64im_min_hot_report.nx (3000 gates, 32// 2x speedup) is the prime addition target -- THIS bench is the 33// evidence base for that decision. 34// 35// Status: SEED. 2026-05-26. Single-stage radix-2; future commits 36// add split-radix (~40% fewer multiplies) + the SIMD-vector emit 37// variant when V-ext silicon candidate lands. 38 39import "nx_syscalls.nx" 40 41const NX_FFT_N: i64 = 256 42const NX_FFT_LOG2_N: i64 = 8 43const NX_FFT_Q15_ONE: i64 = 32767 44const NX_FFT_ITERS: i64 = 1000 45 46// Twiddle step per Q15. cos(2*pi/256) ~ 0.99969 -> 32757; sin -> 804. 47const NX_FFT_COS_STEP: i64 = 32757 48const NX_FFT_SIN_STEP: i64 = 804 49 50// ===== Signed arithmetic right shift (NishiLang >> is logical only) ================================================= 51// 52// FFT Q15 multiplies produce 32-bit signed intermediate; we shift 53// right by 15 to renormalise. Logical shift wrong-signs negatives; 54// fill top bits with the sign on the negative path. 55 56func nx_fft_sar15(v: i64) -> i64 { 57 if v >= 0 { return v >> 15 } 58 // Negative: shift logically, then OR in sign-extension bits. 59 let logical: i64 = v >> 15 60 let sext: i64 = (0 - 1) << 49 // top 15 bits set 61 return logical | sext 62} 63 64// ===== Twiddle table generation ================================================= 65// 66// Builds cos_tab + sin_tab for k = 0..N/2-1 via the angle-addition 67// recurrence. ~256 Q15 multiplies; runs once per benchmark process. 68 69func nx_fft_gen_twiddles(cos_tab: *i64, sin_tab: *i64) -> i64 { 70 cos_tab[0] = NX_FFT_Q15_ONE 71 sin_tab[0] = 0 72 var k: i64 = 1 73 let half: i64 = NX_FFT_N >> 1 74 while k < half { 75 let prev_cos: i64 = cos_tab[k - 1] 76 let prev_sin: i64 = sin_tab[k - 1] 77 cos_tab[k] = nx_fft_sar15((prev_cos * NX_FFT_COS_STEP) - (prev_sin * NX_FFT_SIN_STEP)) 78 sin_tab[k] = nx_fft_sar15((prev_sin * NX_FFT_COS_STEP) + (prev_cos * NX_FFT_SIN_STEP)) 79 k = k + 1 80 } 81 return 0 82} 83 84// ===== 8-bit bit reverse ================================================= 85// 86// FFT decimation-in-time requires bit-reverse permutation. For 87// N=256 the index is 8 bits; reverse via the standard bit-swap chain. 88 89func nx_fft_bit_reverse_8(x: i64) -> i64 { 90 var v: i64 = x & 0xff 91 // Swap adjacent pairs of 4 bits, then 2 bits, then 1 bit. 92 v = ((v & 0xf0) >> 4) | ((v & 0x0f) << 4) 93 v = ((v & 0xcc) >> 2) | ((v & 0x33) << 2) 94 v = ((v & 0xaa) >> 1) | ((v & 0x55) << 1) 95 return v & 0xff 96} 97 98// ===== Bit-reverse permutation ================================================= 99// 100// In-place; only swap when reversed index > current (each pair 101// touched exactly once). 102 103func nx_fft_bit_reverse_permute(re: *i64, im: *i64) -> i64 { 104 var i: i64 = 0 105 while i < NX_FFT_N { 106 let j: i64 = nx_fft_bit_reverse_8(i) 107 if j > i { 108 let tr: i64 = re[i] 109 re[i] = re[j] 110 re[j] = tr 111 let ti: i64 = im[i] 112 im[i] = im[j] 113 im[j] = ti 114 } 115 i = i + 1 116 } 117 return 0 118} 119 120// ===== Cooley-Tukey radix-2 DIT FFT ================================================= 121// 122// Eight butterfly stages for N=256. Each stage doubles the group 123// size m; the twiddle stride is N/m so larger groups use closer-to- 124// k=0 twiddle values. 125 126func nx_fft_radix2(re: *i64, im: *i64, cos_tab: *i64, sin_tab: *i64) -> i64 { 127 nx_fft_bit_reverse_permute(re, im) 128 var stage: i64 = 1 129 while stage <= NX_FFT_LOG2_N { 130 let m: i64 = 1 << stage 131 let half_m: i64 = m >> 1 132 let step: i64 = NX_FFT_N / m 133 var group: i64 = 0 134 while group < NX_FFT_N { 135 var k: i64 = 0 136 while k < half_m { 137 let t_idx: i64 = k * step 138 let wc: i64 = cos_tab[t_idx] 139 let ws: i64 = sin_tab[t_idx] 140 let a_re: i64 = re[group + k] 141 let a_im: i64 = im[group + k] 142 let b_re: i64 = re[group + k + half_m] 143 let b_im: i64 = im[group + k + half_m] 144 // t = w * b (complex mul); Q15 * Q15 -> Q30, shift right 15. 145 let t_re: i64 = nx_fft_sar15((wc * b_re) - (ws * b_im)) 146 let t_im: i64 = nx_fft_sar15((wc * b_im) + (ws * b_re)) 147 // Butterfly 148 re[group + k] = a_re + t_re 149 im[group + k] = a_im + t_im 150 re[group + k + half_m] = a_re - t_re 151 im[group + k + half_m] = a_im - t_im 152 k = k + 1 153 } 154 group = group + m 155 } 156 stage = stage + 1 157 } 158 return 0 159} 160 161// ===== Bench main ================================================= 162 163func main() -> i64 { 164 // Allocate i64-per-bin arrays for re + im (NishiLang has no i32 165 // narrow type; the Q15 values fit comfortably in i64). 166 let re: *i64 = (sys_mmap(8 * NX_FFT_N)) as *i64 167 let im: *i64 = (sys_mmap(8 * NX_FFT_N)) as *i64 168 let cos_tab: *i64 = (sys_mmap(8 * (NX_FFT_N >> 1))) as *i64 169 let sin_tab: *i64 = (sys_mmap(8 * (NX_FFT_N >> 1))) as *i64 170 171 // Generate twiddle tables once. 172 nx_fft_gen_twiddles(cos_tab, sin_tab) 173 174 // ----- main loop ----- 175 var iter: i64 = 0 176 while iter < NX_FFT_ITERS { 177 // Reseed input each iteration so the compiler can't const-fold 178 // the FFT away. Real signal = sin wave at bin 5 + DC offset. 179 // 5/256 of full circle per sample = 7.03 degrees. Hard-code 180 // a tiny sin via 7 LUT entries indexed by i mod 36 (~10 deg 181 // resolution); accuracy doesn't matter for a perf bench. 182 var i: i64 = 0 183 while i < NX_FFT_N { 184 // Pseudo-sinusoid: tile the index 0..7 against a small LUT 185 // baked here so the bench is purely deterministic. 186 let phase: i64 = (i * 5) & 0x7 187 var sample: i64 = 0 188 if phase == 0 { sample = 0 } 189 if phase == 1 { sample = 11585 } 190 if phase == 2 { sample = 23170 } 191 if phase == 3 { sample = 30273 } 192 if phase == 4 { sample = 32767 } 193 if phase == 5 { sample = 30273 } 194 if phase == 6 { sample = 23170 } 195 if phase == 7 { sample = 11585 } 196 // Add tiny DC offset; iteration-varying so compiler 197 // can't lift sample loop out of the iter loop. 198 re[i] = sample + (iter & 0xff) 199 im[i] = 0 200 i = i + 1 201 } 202 203 nx_fft_radix2(re, im, cos_tab, sin_tab) 204 205 iter = iter + 1 206 } 207 208 // Defeat DCE: fold a byte from re[0] into the exit code (masked). 209 return (re[0] as i64) & 0 210}