nx_vox_tract.nx source

1// nx_vox_tract.nx -- the TIMBRE half of animal-sound synthesis: a vocal tract as a physical tube. 2// 3// nx_vox_source produces a glottal flow. That flow sounds like a buzz until it passes through a 4// throat, and the shape of the throat is most of what makes a rabbit sound like a rabbit and a dog 5// sound like a dog. This organ is that throat, modelled as a chain of cylindrical sections. 6// 7// WHAT IS GENUINELY NEW HERE, AND WHAT IS DELIBERATELY NOT. A tree-wide audit ran before a line was 8// written (operator: "make sure we arent duplicating stuff"), and it found that most of this problem is 9// already solved: 10// * nx_lpc_levinson ALREADY emits PARCOR reflection coefficients in Q30. Reflection coefficients are 11// exactly what a tube chain is made of, so this organ adopts that representation verbatim rather 12// than inventing a second one -- a tract built here is interchangeable with the voice stack's. 13// * nx_lpc_synth ALREADY implements the all-pole filter in direct form. It is imported and USED. 14// What did not exist anywhere in 16000+ files, and is therefore all this organ adds: 15// * vt_area_to_k -- tube GEOMETRY (cross-sectional areas) -> reflection coefficients. The whole 16// point: the existing coefficients are FITTED to a recording, these are DERIVED 17// from a shape you can state, so a species can be described instead of sampled. 18// * vt_k_to_a -- the step-up recursion, so a geometry can drive the existing direct-form filter. 19// * vt_ladder -- the Kelly-Lochbaum scattering ladder. 20// 21// WHY A LADDER AT ALL, IF THE DIRECT FORM ALREADY EXISTS. They compute the same transfer function, so 22// the ladder has to earn its place: it is the PHYSICAL structure. Each coefficient is one tube 23// junction, forward and backward travelling waves are separately available, |k|<1 guarantees stability 24// by construction, and a section can be moved while the filter runs (which is what articulation is). 25// The direct form has none of that. To prove this is a second REALISATION and not a second 26// IMPLEMENTATION, nx_vox_tract_gate asserts the ladder and nx_lpc_synth agree sample-for-sample on the 27// same geometry. If they ever diverge, one of them is wrong and the gate says so. 28// 29// TRACT LENGTH SETS THE FORMANTS, and that is the species signature. For a uniform tube the 30// resonances sit at (2n-1)*c/(4L), so a 150 mm dog throat puts its first formant near 570 Hz and a 31// 15 mm mouse throat puts it near 5700 Hz. Section count is NOT free: a waveguide advances one section 32// per sample, so N = 2*L*rate/c. The organ derives N rather than letting a caller pick a number that 33// silently means a different tube than the one they asked for. 34// 35// license_tier: ORIGINAL expect_exit: 0 36 37import "nx_syscalls.nx" 38import "nx_lpc_synth.nx" // REUSE: the existing direct-form all-pole filter 39 40const VT_Q30: i64 = 1073741824 // 2^30 -- the PARCOR scale nx_lpc_levinson already uses 41const VT_C: i64 = 343000 // speed of sound, mm/s 42const VT_MAXSEC: i64 = 96 // sections; 96 covers a 170 mm tract at 96 kHz 43const VT_SHAPE_N: i64 = 8 // control points in a species shape profile 44 45// ---------------------------------------------------------------- geometry -> reflection coefficients 46// Section count a waveguide must use for a tube of length L at this sample rate: one section per 47// sample of travel time. Returned clamped to VT_MAXSEC. 48func vt_nsec(len_mm: i64, rate: i64) -> i64 { 49 if len_mm <= 0 { return 0 } 50 var n: i64 = (2 * len_mm * rate) / VT_C 51 if n < 2 { n = 2 } 52 if n > VT_MAXSEC { n = VT_MAXSEC } 53 return n 54} 55 56// Reflection coefficient at the junction between section i and i+1: 57// k = (A_i - A_{i+1}) / (A_i + A_{i+1}) 58// Areas may be in any consistent unit -- only their ratio matters, which is why a species can be 59// described by a SHAPE and scaled freely. Writes nsec-1 coefficients in Q30. Returns the count. 60func vt_area_to_k(areas: *i64, nsec: i64, k_out: *i64) -> i64 { 61 if nsec < 2 { return 0 } 62 var i: i64 = 0 63 while i < nsec - 1 { 64 let a0: i64 = areas[i] 65 let a1: i64 = areas[i + 1] 66 let sum: i64 = a0 + a1 67 var k: i64 = 0 68 if sum > 0 { k = ((a0 - a1) * VT_Q30) / sum } 69 // |k| < 1 is what makes a tube chain stable; a degenerate area could push it to the boundary. 70 if k > VT_Q30 - 1 { k = VT_Q30 - 1 } 71 if k < 1 - VT_Q30 { k = 1 - VT_Q30 } 72 k_out[i] = k 73 i = i + 1 74 } 75 return nsec - 1 76} 77 78// Step-up recursion: reflection coefficients -> direct-form predictor coefficients, both Q30. 79// 80// SIGN CONVENTION IS NOT A DETAIL. nx_lpc_synth computes s[n] = e[n] - sum a[j] s[n-j], i.e. 81// A(z) = 1 + sum a_j z^-j, and that convention requires PLUS in the step-up recursion. The other 82// common textbook convention (A(z) = 1 - sum ...) uses minus, and picking it here produced a filter 83// that was internally self-consistent and completely disagreed with the organ it was supposed to 84// reuse -- max sample difference 40768 out of a 32767 range. Both halves were "right"; they were 85// right about different conventions. This one follows nx_lpc_synth because reuse is the point. 86func vt_k_to_a(k: *i64, order: i64, a_out: *i64) -> i64 { 87 if order < 1 { return 0 } 88 let tmp: *i64 = sys_mmap((order + 1) * 8) as *i64 89 var i: i64 = 0 90 while i <= order { a_out[i] = 0; tmp[i] = 0; i = i + 1 } 91 var m: i64 = 1 92 while m <= order { 93 // a^(m)_m = k_m 94 tmp[m - 1] = k[m - 1] 95 var j: i64 = 1 96 while j < m { 97 tmp[j - 1] = a_out[j - 1] + ((k[m - 1] * a_out[m - j - 1]) / VT_Q30) 98 j = j + 1 99 } 100 j = 1 101 while j <= m { a_out[j - 1] = tmp[j - 1]; j = j + 1 } 102 m = m + 1 103 } 104 return order 105} 106 107// ---------------------------------------------------------------- the Kelly-Lochbaum ladder 108// All-pole lattice synthesis. exc[] is the excitation (glottal flow), out[] the radiated signal. 109// b[] carries the backward-travelling wave between samples -- that state IS the air in the tube. 110func vt_ladder(exc: *i64, n: i64, k: *i64, order: i64, out: *i64) -> i64 { 111 if order < 1 { return 0 - 1 } 112 let b: *i64 = sys_mmap((order + 2) * 8) as *i64 113 var i: i64 = 0 114 while i <= order + 1 { b[i] = 0; i = i + 1 } 115 var t: i64 = 0 116 while t < n { 117 var f: i64 = exc[t] 118 // walk the junctions from the lips back to the glottis 119 // Matched to the step-up above and hence to nx_lpc_synth: subtract on the forward path, 120 // add on the backward path. T5 in the gate is what holds these three in agreement. 121 var m: i64 = order 122 while m >= 1 { 123 f = f - ((k[m - 1] * b[m - 1]) / VT_Q30) 124 b[m] = b[m - 1] + ((k[m - 1] * f) / VT_Q30) 125 m = m - 1 126 } 127 b[0] = f 128 out[t] = f 129 t = t + 1 130 } 131 return 0 132} 133 134// ---------------------------------------------------------------- species 135// A species is a tract LENGTH plus a SHAPE (relative areas, lips last). Length sets where the 136// formants sit; shape sets their pattern. Both are data, not code (rule 11) -- the table below is the 137// fallback and data/vox_species.cfg overrides it. 138const VT_SP_LEN: i64 = 0 // mm, glottis to lips 139const VT_SP_F0: i64 = 1 // Hz, typical fundamental for the source driving it 140const VT_SP_SHAPE: i64 = 2 // VT_SHAPE_N relative areas follow 141const VT_SP_N: i64 = 10 // 2 + VT_SHAPE_N 142 143const VT_NSPECIES: i64 = 5 144const VT_SP_RABBIT: i64 = 0 145const VT_SP_RODENT: i64 = 1 146const VT_SP_BIRD: i64 = 2 147const VT_SP_DOG: i64 = 3 148const VT_SP_UNIFORM: i64 = 4 // a plain tube: the control case, no shaping at all 149 150// Shapes are relative areas from glottis (index 0) to lips. A narrow glottal end opening into a wider 151// pharynx and closing again at the lips is the generic mammalian pattern; the differences between these 152// rows are what a listener hears as species. 153// 154// PROVENANCE: these are PLAUSIBLE SHAPES consistent with published tract lengths for each class, NOT 155// measured area functions from imaging or casts. The lane has no licensed morphometric source, so the 156// lengths are defensible and the shapes are illustrative. Stated here so the distinction survives. 157func vt_species_defaults(tbl: *i64) -> i64 { 158 var i: i64 = 0 159 while i < VT_NSPECIES * VT_SP_N { tbl[i] = 0; i = i + 1 } 160 161 // rabbit: ~55 mm tract, distress calls around 1-3 kHz 162 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_LEN] = 55 163 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_F0] = 1400 164 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 0] = 10 165 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 1] = 22 166 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 2] = 34 167 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 3] = 30 168 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 4] = 20 169 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 5] = 14 170 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 6] = 12 171 tbl[VT_SP_RABBIT * VT_SP_N + VT_SP_SHAPE + 7] = 9 172 173 // rodent: ~16 mm tract, squeaks reaching well past 8 kHz -- the terrier target 174 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_LEN] = 16 175 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_F0] = 4200 176 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 0] = 6 177 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 1] = 9 178 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 2] = 14 179 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 3] = 16 180 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 4] = 13 181 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 5] = 10 182 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 6] = 8 183 tbl[VT_SP_RODENT * VT_SP_N + VT_SP_SHAPE + 7] = 7 184 185 // upland bird: ~75 mm trachea, near-uniform bore -- birds shape far less than mammals 186 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_LEN] = 75 187 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_F0] = 2300 188 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 0] = 12 189 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 1] = 14 190 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 2] = 15 191 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 3] = 15 192 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 4] = 15 193 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 5] = 14 194 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 6] = 13 195 tbl[VT_SP_BIRD * VT_SP_N + VT_SP_SHAPE + 7] = 12 196 197 // dog: ~150 mm tract -- present so the toy can be told what NOT to sound like 198 tbl[VT_SP_DOG * VT_SP_N + VT_SP_LEN] = 150 199 tbl[VT_SP_DOG * VT_SP_N + VT_SP_F0] = 450 200 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 0] = 14 201 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 1] = 30 202 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 2] = 46 203 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 3] = 42 204 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 4] = 28 205 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 5] = 22 206 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 6] = 20 207 tbl[VT_SP_DOG * VT_SP_N + VT_SP_SHAPE + 7] = 18 208 209 // uniform tube: every area equal, so every reflection coefficient is zero. The control. 210 tbl[VT_SP_UNIFORM * VT_SP_N + VT_SP_LEN] = 55 211 tbl[VT_SP_UNIFORM * VT_SP_N + VT_SP_F0] = 1400 212 i = 0 213 while i < VT_SHAPE_N { tbl[VT_SP_UNIFORM * VT_SP_N + VT_SP_SHAPE + i] = 20; i = i + 1 } 214 return 0 215} 216 217// Expand a species' 8-point shape into `nsec` section areas by linear interpolation. 218func vt_species_areas(tbl: *i64, sp: i64, nsec: i64, areas: *i64) -> i64 { 219 if nsec < 2 { return 0 } 220 let base: i64 = sp * VT_SP_N + VT_SP_SHAPE 221 var i: i64 = 0 222 while i < nsec { 223 // position along the shape, in units of (VT_SHAPE_N-1) scaled by 1024 for the interpolation 224 let pos: i64 = (i * (VT_SHAPE_N - 1) * 1024) / (nsec - 1) 225 let idx: i64 = pos / 1024 226 let frac: i64 = pos % 1024 227 var a0: i64 = tbl[base + idx] 228 var a1: i64 = a0 229 if idx + 1 < VT_SHAPE_N { a1 = tbl[base + idx + 1] } 230 areas[i] = a0 + ((a1 - a0) * frac) / 1024 231 if areas[i] < 1 { areas[i] = 1 } 232 i = i + 1 233 } 234 return nsec 235} 236 237// First formant of a uniform tube of this length: F1 = c / (4L). Published beside any species so the 238// tract's headline number is checkable by hand rather than taken on trust. 239func vt_f1_uniform(len_mm: i64) -> i64 { 240 if len_mm <= 0 { return 0 } 241 return VT_C / (4 * len_mm) 242} 243 244// Convenience: drive a species tract with an excitation, in one call. 245func vt_speak(tbl: *i64, sp: i64, exc: *i64, n: i64, rate: i64, out: *i64) -> i64 { 246 let len_mm: i64 = tbl[sp * VT_SP_N + VT_SP_LEN] 247 let nsec: i64 = vt_nsec(len_mm, rate) 248 if nsec < 2 { return 0 - 1 } 249 let areas: *i64 = sys_mmap(nsec * 8) as *i64 250 let k: *i64 = sys_mmap(nsec * 8) as *i64 251 vt_species_areas(tbl, sp, nsec, areas) 252 let order: i64 = vt_area_to_k(areas, nsec, k) 253 return vt_ladder(exc, n, k, order, out) 254}