nx_vox_source.nx source

1// nx_vox_source.nx -- a SELF-OSCILLATING physical voice source. Rung R2 of the bioacoustics lane. 2// 3// WHAT MAKES THIS DIFFERENT FROM EVERY OTHER AUDIO ORGAN IN THE TREE. nx_audio_bark, nx_tts_formant, 4// nx_phoneme_synth and friends are all SOURCE-FILTER synthesisers: something generates a pulse train at 5// a frequency you tell it, and a filter bank shapes it. The frequency is an input. Nothing oscillates 6// on its own, so nothing can do what real animal larynges do -- jump registers, split into two pitches, 7// break into chaos -- unless a human scripts each event by hand. That is why nx_audio_bark's header 8// says "tune by ear": there is nothing to measure because there is no physics. 9// 10// Here the frequency is an OUTPUT. This is a mass-spring-damper model of two vocal folds, driven by 11// nothing but a constant subglottal pressure, and it oscillates because of the aerodynamics: 12// 13// * Each fold is two coupled masses (lower m1, upper m2) -- the standard two-mass idealisation of the 14// body-cover structure of a real fold. 15// * The glottal areas are a1 = g0 + xL1 + xR1 and a2 = g0 + xL2 + xR2. 16// * Bernoulli gives the driving pressure P1 = Ps * (1 - (amin/a1)^2) where amin is the narrowest 17// point. When the glottis is CONVERGENT (a2 < a1) the narrowest point is downstream, pressure 18// inside the glottis stays high, and the folds are pushed APART. When it is DIVERGENT (a1 < a2) 19// the flow separates at the entry, P1 collapses to zero, and the springs pull the folds shut. 20// * Because m2 lags m1, the glottis is convergent while opening and divergent while closing. Energy 21// goes in during opening and not during closing. That asymmetry IS the oscillation -- it is why a 22// larynx works, and it is the whole reason this file exists. 23// 24// EVERYTHING THE RULER MEASURES FALLS OUT OF THE PHYSICS, unscripted: 25// * f0 emerges from mass and stiffness; it is never written down anywhere. 26// * SUBHARMONICS appear when the two folds entrain at 2:1 instead of 1:1. 27// * BIPHONATION appears when left and right folds are detuned enough to stop entraining and each 28// runs at its own rate. This is why the model is left/right asymmetric rather than a 29// single fold mirrored -- a symmetric model can never produce it. 30// * CHAOS appears at high subglottal pressure, where the collision nonlinearity dominates. 31// * NON-REPETITION comes from turbulent pressure jitter, so no two cycles are ever identical. A dog 32// habituates to a sample loop precisely because a loop cannot do this. 33// 34// The claim "these phenomena are emergent" is not an assertion in a comment -- nx_vox_source_gate 35// SWEEPS asymmetry and pressure and shows the regimes appearing in the order the physics predicts, as 36// measured by nx_bioacoustic_bench, which was built and frozen BEFORE this file existed. 37// 38// Integer fixed-point throughout (Q16), explicit Euler at an oversampled step. Deterministic: the same 39// parameters give bit-identical audio, which is what makes a regression baseline meaningful. 40// license_tier: ORIGINAL expect_exit: 0 41 42import "nx_syscalls.nx" 43import "nx_isqrt.nx" 44 45const VF_Q: i64 = 65536 46const VF_PI2: i64 = 411774 // 2*pi in Q16 -- same constant the rest of the audio family uses 47 48// ---------------------------------------------------------------- parameters (data-driven, rule 11) 49const VF_P_F0_L: i64 = 0 // Hz natural frequency of the LEFT fold 50const VF_P_F0_R: i64 = 1 // Hz natural frequency of the RIGHT fold (detune -> biphonation) 51const VF_P_ZETA: i64 = 2 // Q16 damping ratio 52const VF_P_KCOUP: i64 = 3 // Q16 spring coupling between the lower and upper mass of one fold 53const VF_P_G0: i64 = 4 // Q16 rest half-gap; negative means the folds rest pressed together 54const VF_P_PS: i64 = 5 // Q16 subglottal pressure -- the ONLY energy input 55const VF_P_KCOLL: i64 = 6 // Q16 extra stiffness when the folds collide (the closure nonlinearity) 56const VF_P_F2RATIO: i64 = 7 // Q16 upper-mass natural frequency as a fraction of the lower 57const VF_P_DRIVE: i64 = 8 // Q16 aerodynamic force gain 58const VF_P_OVS: i64 = 9 // integration steps per audio sample 59const VF_P_SEED: i64 = 10 // jitter RNG seed 60const VF_P_JITTER: i64 = 11 // Q16 turbulent pressure fluctuation -- why no two cycles repeat 61// Rest gap of the UPPER aperture. A valve whose two apertures are identical at rest is perfectly 62// symmetric, so a1 == a2, the Bernoulli term is exactly zero, and nothing can start it: pressure alone 63// produces silence forever. The standalone model got away with this because it is handed a small 64// asymmetric nudge at t=0 and its pressure is on from the first sample, so the nudge grows before it 65// damps. Drive the same model from a chamber that takes 60 ms to pressurise and the nudge is long dead 66// by the time there is anything to amplify -- which is exactly how the squeaker came out silent. 67// Real folds and real reed tongues both rest in a CONVERGENT shape; this parameter is that shape, and 68// it makes the model self-starting from pressure alone. Zero means "same as VF_P_G0". 69const VF_P_G0_2: i64 = 12 // Q16 70const VF_P_N: i64 = 13 71 72// ---------------------------------------------------------------- state 73const VF_S_XL1: i64 = 0 74const VF_S_VL1: i64 = 1 75const VF_S_XL2: i64 = 2 76const VF_S_VL2: i64 = 3 77const VF_S_XR1: i64 = 4 78const VF_S_VR1: i64 = 5 79const VF_S_XR2: i64 = 6 80const VF_S_VR2: i64 = 7 81const VF_S_RNG: i64 = 8 82const VF_S_N: i64 = 9 83 84// ---------------------------------------------------------------- precomputed coefficients 85// UNITS. The first version of this model stored displacement and velocity both in Q16 and added the 86// aerodynamic drive directly as a pressure. Those are not the same units: stiffness enters the update 87// as (w*dt)^2, which at 420 Hz and a 384 kHz sub-rate is about 5e-5, while the drive was order 1. 88// The spring was therefore ~20000x too weak to ever pull the folds back, they flew to the safety clamp 89// and stayed there, and the model emitted a DC rail with jitter on it. It still scored 460 on the 90// ruler and still reported ax_vocal=1000, because a railed signal is locally smooth -- which is exactly 91// why the ground-truth probe (count glottal closures directly) had to exist before any of this was 92// believed. openduty_permille=999 and closures=1 is what a "working" voice source looked like. 93// 94// Now every force is expressed as an ACCELERATION in the same units as the stiffness term, and the 95// drive is stated as "the displacement this pressure would hold the fold at", which is unit-correct by 96// construction. Velocity is carried in Q32 because the per-sample acceleration of a 420 Hz oscillator 97// at Q16 rounds to the integer 3 and quantisation eats the dynamics. 98const VF_K_WL1: i64 = 0 // Q16 w*dt for the lower-left mass 99const VF_K_WL2: i64 = 1 100const VF_K_WR1: i64 = 2 101const VF_K_WR2: i64 = 3 102const VF_K_RL1: i64 = 4 // Q16 2*zeta*w*dt 103const VF_K_RL2: i64 = 5 104const VF_K_RR1: i64 = 6 105const VF_K_RR2: i64 = 7 106const VF_K_N: i64 = 8 107const VF_XCLAMP: i64 = 16 * 65536 // +-16.0 displacement; a divergent solve is bounded, not loud 108 109func vfold_defaults(p: *i64) -> i64 { 110 p[VF_P_F0_L] = 420 111 p[VF_P_F0_R] = 420 // symmetric by default -> ordinary 1:1 phonation 112 p[VF_P_ZETA] = 6553 // 0.10 113 p[VF_P_KCOUP] = 32768 // 0.50 114 p[VF_P_G0] = 6553 // 0.10 -- a small resting gap 115 p[VF_P_PS] = 65536 // 1.00 116 p[VF_P_KCOLL] = 262144 // 4.00 -- collision is much stiffer than the free spring 117 p[VF_P_F2RATIO] = 45875 // 0.70 -- the upper mass is lighter and slower 118 p[VF_P_DRIVE] = 65536 // 1.00 119 p[VF_P_OVS] = 4 120 p[VF_P_SEED] = 20260725 121 p[VF_P_JITTER] = 1310 // 0.02 -- 2% turbulent pressure fluctuation 122 p[VF_P_G0_2] = 0 // 0 = mirror VF_P_G0 (symmetric rest shape) 123 return 0 124} 125 126// Natural frequency -> per-sample stiffness and damping. A mass-spring with stiffness K (Q16) advanced 127// by x += v, v += -K*x - R*v oscillates at sqrt(K/Q) radians PER SAMPLE, so K = Q * (2*pi*f/SR)^2. 128func vfold_init(p: *i64, k: *i64, s: *i64, rate: i64) -> i64 { 129 let sr: i64 = rate * p[VF_P_OVS] 130 let wl: i64 = (VF_PI2 * p[VF_P_F0_L]) / sr 131 let wr: i64 = (VF_PI2 * p[VF_P_F0_R]) / sr 132 let wl2: i64 = (wl * p[VF_P_F2RATIO]) / VF_Q 133 let wr2: i64 = (wr * p[VF_P_F2RATIO]) / VF_Q 134 k[VF_K_WL1] = wl 135 k[VF_K_WR1] = wr 136 k[VF_K_WL2] = wl2 137 k[VF_K_WR2] = wr2 138 k[VF_K_RL1] = (2 * p[VF_P_ZETA] * wl) / VF_Q 139 k[VF_K_RR1] = (2 * p[VF_P_ZETA] * wr) / VF_Q 140 k[VF_K_RL2] = (2 * p[VF_P_ZETA] * wl2) / VF_Q 141 k[VF_K_RR2] = (2 * p[VF_P_ZETA] * wr2) / VF_Q 142 var i: i64 = 0 143 while i < VF_S_N { s[i] = 0; i = i + 1 } 144 // A tiny asymmetric nudge, not a driving term: the model must AMPLIFY this on its own or it is not 145 // a self-oscillator. vs gate assertion S2 sets Ps = 0 and requires the nudge to die away. 146 s[VF_S_XL1] = 256 147 s[VF_S_XR1] = 192 148 s[VF_S_RNG] = p[VF_P_SEED] 149 return 0 150} 151 152func vfold_rng(s: *i64) -> i64 { 153 var r: i64 = s[VF_S_RNG] 154 r = r ^ (r << 13) 155 r = r ^ (r >> 7) 156 r = r ^ (r << 17) 157 s[VF_S_RNG] = r 158 if r < 0 { return 0 - r } 159 return r 160} 161 162// One integration sub-step. Returns the glottal flow for this sub-step (Q16). 163func vfold_substep(s: *i64, p: *i64, k: *i64) -> i64 { 164 // --- glottal areas. a1 is the lower (inlet) aperture, a2 the upper (outlet). 165 var g02: i64 = p[VF_P_G0_2] 166 if g02 == 0 { g02 = p[VF_P_G0] } 167 let a1: i64 = p[VF_P_G0] + s[VF_S_XL1] + s[VF_S_XR1] 168 let a2: i64 = g02 + s[VF_S_XL2] + s[VF_S_XR2] 169 var amin: i64 = a1 170 if a2 < amin { amin = a2 } 171 172 // --- turbulent pressure jitter. Physically this is broadband noise from the jet; practically it 173 // is the reason a real animal never emits the same cycle twice. 174 var ps: i64 = p[VF_P_PS] 175 if p[VF_P_JITTER] > 0 { 176 let j: i64 = vfold_rng(s) % (2 * p[VF_P_JITTER] + 1) 177 ps = ps + (((j - p[VF_P_JITTER]) * ps) / VF_Q) 178 } 179 180 // --- Bernoulli drive. Zero while the glottis is shut: no flow, no driving pressure. 181 var p1: i64 = 0 182 var flow: i64 = 0 183 if amin > 0 { 184 if a1 > 0 { 185 let ratio: i64 = (amin * VF_Q) / a1 186 let sq: i64 = (ratio * ratio) / VF_Q 187 p1 = (ps * (VF_Q - sq)) / VF_Q 188 } 189 // U = amin * sqrt(2*Ps/rho); constants fold into DRIVE, so U ~ amin * sqrt(ps) 190 var pp: i64 = ps 191 if pp < 0 { pp = 0 } 192 flow = (amin * nx_isqrt(pp * VF_Q)) / VF_Q 193 } 194 // The drive is expressed as the DISPLACEMENT this pressure would hold a fold at, so that turning it 195 // into an acceleration below is just the stiffness term evaluated at that displacement -- units 196 // cannot drift apart again. 197 let xdrive: i64 = (p[VF_P_DRIVE] * p1) / VF_Q 198 199 // --- collision. When a fold crosses the midline the tissue is compressed, far stiffer than the 200 // free spring. Stated as a MULTIPLE of the fold's own stiffness, so it scales with pitch. 201 // This is the strongest nonlinearity here and the route to chaos at high pressure. 202 var xcoll1: i64 = 0 203 var xcoll2: i64 = 0 204 if a1 < 0 { xcoll1 = (p[VF_P_KCOLL] * (0 - a1)) / VF_Q } 205 if a2 < 0 { xcoll2 = (p[VF_P_KCOLL] * (0 - a2)) / VF_Q } 206 207 // Coupling between the lower and upper mass, also as a fraction of stiffness. 208 let cplL: i64 = (p[VF_P_KCOUP] * (s[VF_S_XL1] - s[VF_S_XL2])) / VF_Q 209 let cplR: i64 = (p[VF_P_KCOUP] * (s[VF_S_XR1] - s[VF_S_XR2])) / VF_Q 210 211 // --- accelerations, all in Q32 (velocity units). acc = w^2 * (xdrive - x) etc, computed as 212 // (W*W*X)>>16 which keeps full precision instead of rounding w^2 to the integer 3. 213 let wL1: i64 = k[VF_K_WL1] 214 let wL2: i64 = k[VF_K_WL2] 215 let wR1: i64 = k[VF_K_WR1] 216 let wR2: i64 = k[VF_K_WR2] 217 218 let eL1: i64 = xdrive + xcoll1 - s[VF_S_XL1] - cplL 219 let eL2: i64 = xcoll2 - s[VF_S_XL2] + cplL 220 let eR1: i64 = xdrive + xcoll1 - s[VF_S_XR1] - cplR 221 let eR2: i64 = xcoll2 - s[VF_S_XR2] + cplR 222 223 var aL1: i64 = ((wL1 * wL1) * eL1) >> 16 224 var aL2: i64 = ((wL2 * wL2) * eL2) >> 16 225 var aR1: i64 = ((wR1 * wR1) * eR1) >> 16 226 var aR2: i64 = ((wR2 * wR2) * eR2) >> 16 227 aL1 = aL1 - ((k[VF_K_RL1] * s[VF_S_VL1]) >> 16) 228 aL2 = aL2 - ((k[VF_K_RL2] * s[VF_S_VL2]) >> 16) 229 aR1 = aR1 - ((k[VF_K_RR1] * s[VF_S_VR1]) >> 16) 230 aR2 = aR2 - ((k[VF_K_RR2] * s[VF_S_VR2]) >> 16) 231 232 // Semi-implicit (symplectic) Euler: advance velocity first, then position with the NEW velocity. 233 // Plain explicit Euler pumps energy into an oscillator every cycle and is unconditionally unstable 234 // here; this ordering is stable for w*dt < 2, which the oversampling guarantees. 235 s[VF_S_VL1] = s[VF_S_VL1] + aL1 236 s[VF_S_VL2] = s[VF_S_VL2] + aL2 237 s[VF_S_VR1] = s[VF_S_VR1] + aR1 238 s[VF_S_VR2] = s[VF_S_VR2] + aR2 239 s[VF_S_XL1] = s[VF_S_XL1] + (s[VF_S_VL1] >> 16) 240 s[VF_S_XL2] = s[VF_S_XL2] + (s[VF_S_VL2] >> 16) 241 s[VF_S_XR1] = s[VF_S_XR1] + (s[VF_S_VR1] >> 16) 242 s[VF_S_XR2] = s[VF_S_XR2] + (s[VF_S_VR2] >> 16) 243 244 // NEVER-BRICK ANALOG (rule 26 applied to numerics): an integrator can diverge, and a divergent 245 // state would emit full-scale noise into whatever is downstream. Clamp so a bad parameter set 246 // produces a bounded, obviously-wrong sound instead of an unbounded one. The gate checks that a 247 // healthy configuration never touches this clamp -- a model that lives on the rail is not a model. 248 // Displacements are Q16 and velocities Q32, so they need DIFFERENT bounds -- clamping the Q32 249 // velocities to a Q16 bound would silently throttle every fold to a crawl. 250 let vclamp: i64 = VF_XCLAMP * VF_Q 251 var i: i64 = 0 252 while i < 8 { 253 var lim: i64 = vclamp 254 if (i % 2) == 0 { lim = VF_XCLAMP } 255 if s[i] > lim { s[i] = lim } 256 if s[i] < 0 - lim { s[i] = 0 - lim } 257 i = i + 1 258 } 259 return flow 260} 261 262// Render `n` samples of glottal flow into out[] (Q16 flow units). Returns the peak absolute value, 263// which is 0 exactly when the model did not oscillate. 264func vfold_render(out: *i64, n: i64, rate: i64, p: *i64, s: *i64, k: *i64) -> i64 { 265 var peak: i64 = 0 266 let ovs: i64 = p[VF_P_OVS] 267 var i: i64 = 0 268 while i < n { 269 var acc: i64 = 0 270 var o: i64 = 0 271 while o < ovs { acc = acc + vfold_substep(s, p, k); o = o + 1 } 272 let v: i64 = acc / ovs 273 out[i] = v 274 var a: i64 = v 275 if a < 0 { a = 0 - a } 276 if a > peak { peak = a } 277 i = i + 1 278 } 279 return peak 280} 281 282// ---------------------------------------------------------------- self-reported ground truth 283// The model counts its OWN glottal closures, which gives an f0 that owes nothing to any pitch tracker. 284// This exists because the first working-looking version of this organ was not oscillating at all -- 285// it had railed to the safety clamp and was emitting DC with jitter on it, and the ruler still scored 286// it 460 with ax_vocal=1000. A physical model must be able to state what it did independently of the 287// instrument measuring it, or there is no way to tell which of the two is wrong. 288// 289// diag[0] = closures counted, diag[1] = samples with the glottis open, diag[2] = clamp hits. 290// Returns closures per second. 291const VF_D_CLOSURES: i64 = 0 292const VF_D_OPEN: i64 = 1 293const VF_D_CLAMP: i64 = 2 294const VF_D_N: i64 = 3 295 296func vfold_render_diag(out: *i64, n: i64, rate: i64, p: *i64, s: *i64, k: *i64, diag: *i64) -> i64 { 297 diag[VF_D_CLOSURES] = 0 298 diag[VF_D_OPEN] = 0 299 diag[VF_D_CLAMP] = 0 300 let ovs: i64 = p[VF_P_OVS] 301 var wasopen: i64 = 0 302 var i: i64 = 0 303 while i < n { 304 var acc: i64 = 0 305 var o: i64 = 0 306 while o < ovs { acc = acc + vfold_substep(s, p, k); o = o + 1 } 307 out[i] = acc / ovs 308 let a1: i64 = p[VF_P_G0] + s[VF_S_XL1] + s[VF_S_XR1] 309 let a2: i64 = p[VF_P_G0] + s[VF_S_XL2] + s[VF_S_XR2] 310 var amin: i64 = a1 311 if a2 < amin { amin = a2 } 312 var isopen: i64 = 0 313 if amin > 0 { isopen = 1; diag[VF_D_OPEN] = diag[VF_D_OPEN] + 1 } 314 if wasopen == 1 { if isopen == 0 { diag[VF_D_CLOSURES] = diag[VF_D_CLOSURES] + 1 } } 315 wasopen = isopen 316 // A healthy configuration must never reach the divergence clamp. A model that lives on the 317 // rail is not a model, and the gate refuses one. 318 var j: i64 = 0 319 while j < 8 { 320 var lim: i64 = VF_XCLAMP * VF_Q 321 if (j % 2) == 0 { lim = VF_XCLAMP } 322 if s[j] >= lim { diag[VF_D_CLAMP] = diag[VF_D_CLAMP] + 1 } 323 if s[j] <= 0 - lim { diag[VF_D_CLAMP] = diag[VF_D_CLAMP] + 1 } 324 j = j + 1 325 } 326 i = i + 1 327 } 328 if n <= 0 { return 0 } 329 return (diag[VF_D_CLOSURES] * rate) / n 330} 331 332// Convenience: initialise and render in one call, then scale to a 16-bit audio range. 333// The glottal flow is a positive-only pulse train, so its DC component is removed first -- radiation 334// from a real mouth differentiates the flow, and leaving the DC in would waste headroom on an offset 335// no dog can hear. 336func vfold_render_audio(out: *i64, n: i64, rate: i64, p: *i64, amp: i64) -> i64 { 337 let s: *i64 = sys_mmap(VF_S_N * 8) as *i64 338 let k: *i64 = sys_mmap(VF_K_N * 8) as *i64 339 vfold_init(p, k, s, rate) 340 let peak: i64 = vfold_render(out, n, rate, p, s, k) 341 if peak <= 0 { return 0 } 342 var sum: i64 = 0 343 var i: i64 = 0 344 while i < n { sum = sum + out[i]; i = i + 1 } 345 let dc: i64 = sum / n 346 var mx: i64 = 0 347 i = 0 348 while i < n { 349 out[i] = out[i] - dc 350 var a: i64 = out[i] 351 if a < 0 { a = 0 - a } 352 if a > mx { mx = a } 353 i = i + 1 354 } 355 if mx <= 0 { return 0 } 356 i = 0 357 while i < n { out[i] = (out[i] * amp) / mx; i = i + 1 } 358 return mx 359}