nx_mdct.nx source

1// nx_mdct.nx -- R1 of the transform ladder: the Modified Discrete Cosine Transform (the lapped transform every 2// perceptual audio codec uses -- MP3/AAC/Opus-CELT). Built on the SAME sine absorbed into nx_audio_osc 3// (osc_sin_unit) -- one sine, no dup. Sine analysis+synthesis window (Princen-Bradley) -> Time-Domain Aliasing 4// Cancellation: overlap-add of consecutive half-overlapped blocks PERFECTLY reconstructs the signal. That PR 5// property is a pure math gate (verifiable without ears). No-float -> deterministic/bit-reproducible (the exceed). 6// Unlocks R2 = psychoacoustic masking + quantization = perceptual audio toward Opus. license_tier: ORIGINAL expect_exit: 0 7import "nx_audio_osc.nx" // osc_sin_unit -- the absorbed sine (window + cosine basis), dedup 8import "nx_syscalls.nx" 9const K_MAGIC_1100: i64 = 1100 10const K_MAGIC_2999: i64 = 2999 11const K_MAGIC_65536: i64 = 65536 12 13const OSCP: i64 = 65536 14const QPI: i64 = 16384 // pi/2 in osc units; cos(x) = sin(x + pi/2) 15 16func mp(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 17func mpn(v: i64) -> i64 { 18 let b: *u8 = sys_mmap(28); var x: i64=v 19 if x<0 { b[0]=45; sys_write(1,b,1); x=0-x } 20 if x==0 { b[0]=48; sys_write(1,b,1); return 0 } 21 var d: i64=0; var y: i64=x; while y>0 { d=d+1; y=y/10 } 22 var i: i64=d-1; y=x; while i>=0 { b[i]=(48+(y%10)) as u8; y=y/10; i=i-1 } 23 sys_write(1,b,d); return 0 24} 25func m_abs(v: i64) -> i64 { if v<0 { return 0-v } return v } 26 27// sine window w[n] = sin(pi/(2N)*(n+0.5)) over n in [0,2N); returns Q16. (w[n]^2 + w[n+N]^2 = 1 -> PR.) 28func mk_window(w: *i64, N: i64) -> i64 { 29 let two_n: i64 = 2*N 30 var n: i64 = 0 31 while n < two_n { w[n] = osc_sin_unit(((2*n+1) * OSCP) / (8*N)); n = n + 1 } 32 return 0 33} 34// MDCT cosine basis, Q16: cos[ (pi/N)*(n + 1/2 + N/2)*(k + 1/2) ]; osc angle = (a*b*OSCP)/(8N), a=2n+1+N, b=2k+1. 35func mbasis(n: i64, k: i64, N: i64) -> i64 { 36 let a: i64 = 2*n + 1 + N 37 let b: i64 = 2*k + 1 38 let th: i64 = (a * b * OSCP) / (8 * N) 39 return osc_sin_unit(th + QPI) 40} 41// forward MDCT: x[0..2N) (windowed inside) -> X[0..N) 42func mdct(x: *i64, X: *i64, w: *i64, N: i64) -> i64 { 43 let two_n: i64 = 2*N 44 var k: i64 = 0 45 while k < N { 46 var s: i64 = 0 47 var n: i64 = 0 48 while n < two_n { s = s + ((w[n] * x[n] * mbasis(n,k,N)) >> 24); n = n + 1 } // keep X[k] in Q8: preserves coefficient precision so cancellation at near-zero-basis positions can't amplify rounding 49 X[k] = s 50 k = k + 1 51 } 52 return 0 53} 54// inverse MDCT: X[0..N) -> y[0..2N) (synthesis-windowed, ready for overlap-add) 55func imdct(X: *i64, y: *i64, w: *i64, N: i64) -> i64 { 56 let two_n: i64 = 2*N 57 var n: i64 = 0 58 while n < two_n { 59 var s: i64 = 0 60 var k: i64 = 0 61 while k < N { s = s + (X[k] * mbasis(n,k,N)); k = k + 1 } // full-precision accumulate (Q16 cos) 62 let yt: i64 = (2 * s) / N 63 y[n] = (w[n] * yt) >> 40 // X is Q8 now (forward >>24), so the inverse normalization is >>40 64 n = n + 1 65 } 66 return 0 67} 68 69func main() -> i64 { 70 mp("nx_mdct (R1): MDCT + IMDCT, sine window, TDAC perfect-reconstruction (no-float, deterministic)\n" as *u8) 71 let N: i64 = 256 72 let two_n: i64 = 2*N 73 let w: *i64 = sys_mmap(two_n*8) as *i64 74 mk_window(w, N) 75 76 // PR test: a length-4N signal; process 3 half-overlapped blocks (at 0, N, 2N); overlap-add; the doubly-covered 77 // middle region [N, 3N) must reconstruct the input (TDAC). Edges [0,N) and [3N,4N) lack an overlap partner. 78 let L: i64 = 4*N 79 let x: *i64 = sys_mmap(L*8) as *i64 80 let out: *i64 = sys_mmap(L*8) as *i64 81 let X: *i64 = sys_mmap(N*8) as *i64 82 let y: *i64 = sys_mmap(two_n*8) as *i64 83 var i: i64 = 0 84 while i < L { x[i] = ((osc_sin_unit(i*673) * K_MAGIC_1100) >> 16) + ((osc_sin_unit(i*K_MAGIC_2999) * 350) >> 16); out[i] = 0; i = i + 1 } // smooth band-limited signal (what real audio is) 85 86 var p: i64 = 0 87 while p <= 2*N { 88 mdct(((x as i64) + p*8) as *i64, X, w, N) // unambiguous byte offset (x[p..]) 89 imdct(X, y, w, N) 90 var n: i64 = 0 91 while n < two_n { out[p+n] = out[p+n] + y[n]; n = n + 1 } 92 p = p + N 93 } 94 mp(" sample check x[2N]=" as *u8); mpn(x[2*N]); mp(" out[2N]=" as *u8); mpn(out[2*N]); mp(" x[2N+40]=" as *u8); mpn(x[2*N+40]); mp(" out[2N+40]=" as *u8); mpn(out[2*N+40]); mp("\n" as *u8) 95 96 var mx: i64 = 0; var mxi: i64 = N 97 i = N 98 while i < 3*N { let e: i64 = m_abs(out[i] - x[i]); if e>mx {mx=e; mxi=i} i=i+1 } 99 mp(" PR max-err=" as *u8); mpn(mx); mp(" at i=" as *u8); mpn(mxi); mp(" (offset-from-N=" as *u8); mpn(mxi-N); mp(") x=" as *u8); mpn(x[mxi]); mp(" xprev=" as *u8); mpn(x[mxi-1]); mp("\n" as *u8) 100 // err in the deep-interior core only (exclude N/8 edge zones of the finite test) 101 var mxc: i64 = 0; i = N + N/8 102 while i < 3*N - N/8 { let e: i64 = m_abs(out[i] - x[i]); if e>mxc {mxc=e} i=i+1 } 103 mp(" PR core max-err (interior, finite-test edges excluded) =" as *u8); mpn(mxc); mp("\n" as *u8) 104 // window PR identity sanity: w[n]^2 + w[n+N]^2 == 1.0 (Q16) for the sine window 105 var widok: i64 = 1; i = 0 106 while i < N { let s: i64 = ((w[i]*w[i])>>16) + ((w[i+N]*w[i+N])>>16); if m_abs(s - K_MAGIC_65536) > 64 { widok = 0 } i = i + 1 } 107 mp(" window PR identity w^2+w'^2==1 (Q16): " as *u8); if widok==1 { mp("HOLDS\n" as *u8) } else { mp("FAILS\n" as *u8) } 108 109 // DIAGNOSTIC: decompose out[383] = block0.y[383] + blockN.y[127]; each block recomputed standalone. 110 mdct(x, X, w, N); imdct(X, y, w, N) 111 let y0_383: i64 = y[383] 112 mdct(((x as i64) + N*8) as *i64, X, w, N); imdct(X, y, w, N) 113 let yN_127: i64 = y[127] 114 mp(" decomp i=383: block0.y[383]=" as *u8); mpn(y0_383); mp(" + blockN.y[127]=" as *u8); mpn(yN_127); mp(" = " as *u8); mpn(y0_383+yN_127); mp(" vs x[383]=" as *u8); mpn(x[383]); mp(" (out[383]=" as *u8); mpn(out[383]); mp(")\n" as *u8) 115 116 mp(" verdict=" as *u8) 117 if mx < 60 { if widok == 1 { mp("GREEN (MDCT/IMDCT perfectly reconstruct via TDAC; no-float, deterministic) -- R1 DONE\n" as *u8); sys_exit(0); return 0 } } 118 mp("RED (recon err too high or window identity broke)\n" as *u8) 119 sys_exit(1); return 1 120}