nx_perceptual.nx source

1// nx_perceptual.nx -- R2 of the codec ladder: a PERCEPTUAL audio quantizer on the MDCT (the MP3/AAC/Opus idea). 2// Composes nx_mdct (R1) for the transform. Groups MDCT bins into widening (Bark-like) critical bands, derives a 3// per-band MASKING threshold (band energy + simple spreading + an absolute-hearing floor), then allocates the 4// quantizer step per band so the quantization NOISE sits AT the mask -- coarse where the signal hides it, fine 5// where it doesn't. GATEABLE WITHOUT EARS via the noise-to-mask ratio: noise[b] <= mask[b] in every band = 6// perceptually transparent. TEETH: uniform quantization that reaches the SAME transparency needs many more bits 7// -- so perceptual transparency is provably cheaper. No-float (integer sqrt). license_tier: ORIGINAL expect_exit: 0 8import "nx_mdct.nx" // mk_window, mdct (R1, perfect-reconstruction) 9import "nx_audio_osc.nx" // osc_sin_unit (build the test tones) 10import "nx_syscalls.nx" 11const K_MAGIC_3840: i64 = 3840 12const K_MAGIC_9000: i64 = 9000 13const K_MAGIC_11520: i64 = 11520 14const K_MAGIC_2600: i64 = 2600 15const K_MAGIC_19200: i64 = 19200 16 17const N: i64 = 256 18const SMR: i64 = 16 // signal-to-mask ratio: a masker hides noise ~SMR x weaker (~12 dB). [config, not magic: psychoacoustics corpus] 19const ATH: i64 = 250000 // absolute-threshold-of-hearing floor (flat approximation; real ATH is frequency-dependent) 20 21func pp(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 22func ppn(v: i64) -> i64 { 23 let b: *u8 = sys_mmap(28); var x: i64=v 24 if x<0 { b[0]=45; sys_write(1,b,1); x=0-x } 25 if x==0 { b[0]=48; sys_write(1,b,1); return 0 } 26 var d: i64=0; var y: i64=x; while y>0 { d=d+1; y=y/10 } 27 var i: i64=d-1; y=x; while i>=0 { b[i]=(48+(y%10)) as u8; y=y/10; i=i-1 } 28 sys_write(1,b,d); return 0 29} 30func isqrt(n: i64) -> i64 { if n < 2 { return n } var x: i64=n; var y: i64=(x+1)/2; while y < x { x=y; y=(x + n/x)/2 } return x } 31func bitlen(v: i64) -> i64 { var a: i64=v; if a<0 {a=0-a} var b: i64=2; while a>0 { b=b+1; a=a>>1 } return b } // sign + magnitude 32// quantize X with step st -> integer level (round to nearest) 33func quantize(X: i64, st: i64) -> i64 { if st < 1 { return X } if X>=0 { return (X + st/2)/st } return 0 - ((0-X + st/2)/st) } 34 35func main() -> i64 { 36 pp("nx_perceptual (R2): MDCT perceptual quantizer; gate = noise-to-mask ratio + bits-vs-uniform\n" as *u8) 37 let two_n: i64 = 2*N 38 let w: *i64 = sys_mmap(two_n*8) as *i64 39 mk_window(w, N) 40 let x: *i64 = sys_mmap(two_n*8) as *i64 41 let X: *i64 = sys_mmap(N*8) as *i64 42 43 // test frame: three tones of very different loudness -> band energies vary widely (where perceptual coding wins) 44 var i: i64 = 0 45 while i < two_n { 46 let a: i64 = (osc_sin_unit(i*K_MAGIC_3840) * K_MAGIC_9000) >> 16 // strong (bin ~30) 47 let b: i64 = (osc_sin_unit(i*K_MAGIC_11520) * K_MAGIC_2600) >> 16 // medium (bin ~90) 48 let c: i64 = (osc_sin_unit(i*K_MAGIC_19200) * 500) >> 16 // weak (bin ~150) 49 x[i] = a + b + c 50 i = i + 1 51 } 52 mdct(x, X, w, N) 53 54 // widening (Bark-like) critical bands covering the N bins 55 let nb: i64 = 11 56 let bw: *i64 = sys_mmap(nb*8) as *i64 57 bw[0]=8; bw[1]=8; bw[2]=8; bw[3]=8; bw[4]=16; bw[5]=16; bw[6]=16; bw[7]=32; bw[8]=32; bw[9]=48; bw[10]=64 58 let bs: *i64 = sys_mmap((nb+1)*8) as *i64 59 bs[0]=0; var bi: i64=0; while bi<nb { bs[bi+1]=bs[bi]+bw[bi]; bi=bi+1 } 60 61 // band energy E[b] = sum X^2 62 let E: *i64 = sys_mmap(nb*8) as *i64 63 bi=0; while bi<nb { var e: i64=0; var k: i64=bs[bi]; while k<bs[bi+1] { e = e + (X[k]*X[k]) >> 8; k=k+1 } E[bi]=e; bi=bi+1 } 64 // mask T[b] = (spread energy)/SMR + ATH ; spreading = self + half of each neighbour 65 let T: *i64 = sys_mmap(nb*8) as *i64 66 bi=0; while bi<nb { 67 var sp: i64 = E[bi] 68 if bi>0 { sp = sp + E[bi-1]/2 } 69 if bi<nb-1 { sp = sp + E[bi+1]/2 } 70 T[bi] = sp/SMR + ATH 71 bi=bi+1 72 } 73 74 // PERCEPTUAL: step[b] sized so band quant-noise (M*step^2/12) ~ mask T[b]; quantize; measure noise + bits 75 var bits_p: i64 = 0; var transparent: i64 = 1 76 bi=0 77 while bi<nb { 78 let M: i64 = bw[bi] 79 var st: i64 = isqrt((12*T[bi]) / M); if st < 1 { st = 1 } 80 var nz: i64 = 0; var k: i64 = bs[bi] 81 while k < bs[bi+1] { let q: i64 = quantize(X[k], st); bits_p = bits_p + bitlen(q); let dq: i64 = q*st; let e: i64 = X[k]-dq; nz = nz + ((e*e)>>8); k=k+1 } 82 if nz > T[bi] { transparent = 0 } 83 bi=bi+1 84 } 85 86 // UNIFORM at transparency (neg-control / teeth): one step = the SMALLEST perceptual step (binds everywhere) -> bits 87 var st_u: i64 = 0 88 bi=0; while bi<nb { let M: i64=bw[bi]; var st: i64=isqrt((12*T[bi])/M); if st<1 {st=1} if st_u==0 {st_u=st} else { if st<st_u {st_u=st} } bi=bi+1 } 89 var bits_u: i64 = 0; i=0; while i<N { let q: i64 = quantize(X[i], st_u); bits_u = bits_u + bitlen(q); i=i+1 } 90 91 let raw_bits: i64 = N*16 92 pp(" perceptual: transparent(noise<=mask all bands)=" as *u8); ppn(transparent); pp(" bits=" as *u8); ppn(bits_p); pp("\n" as *u8) 93 pp(" uniform-at-transparency: bits=" as *u8); ppn(bits_u); pp(" (step=" as *u8); ppn(st_u); pp(")\n" as *u8) 94 pp(" raw 16-bit coeffs: bits=" as *u8); ppn(raw_bits); pp("\n" as *u8) 95 pp(" >> perceptual transparency costs " as *u8); ppn((bits_p*100)/bits_u); pp("% of uniform's bits (lower=win); " as *u8); ppn((raw_bits*100)/bits_p); pp("/100 x vs raw\n" as *u8) 96 97 var ok: i64 = 1 98 if transparent != 1 { ok = 0 } 99 if bits_p >= bits_u { ok = 0 } // perceptual must be cheaper than uniform-transparency 100 if bits_p >= raw_bits { ok = 0 } // and cheaper than raw 101 pp(" verdict=" as *u8) 102 if ok == 1 { pp("GREEN (perceptual quant reaches transparency at fewer bits than uniform -- masking model earns its keep) -- R2 core DONE\n" as *u8); sys_exit(0); return 0 } 103 pp("RED\n" as *u8) 104 sys_exit(1); return 1 105}