`nx_squeaker.nx` source

↩ module page · 290 lines · 14542 B
1// nx_squeaker.nx -- THE TOY.  A dog-toy ball that makes its noise because a dog's jaw compresses it.
2//
3// Operator 2026-07-25: "a dog toy ball that is designed to make the noise when compressed by a dog's
4// mouth".  This is the physical model of that object, from the bite down:
5//
6//     jaw force -> shell compression -> air displaced -> chamber pressure -> reed valve flutters
7//     -> flow through the neck -> radiated sound
8//
9// Every stage is a real physical quantity in manufacturable units (newtons, millimetres, cubic
10// millimetres), because the whole point of rung R5 is to INVERT this model and hand back a geometry
11// somebody can actually mould.  A model in arbitrary units cannot produce a spec.
12//
13// WHY THE PITCH GLIDES, AND WHY THAT IS THE PRODUCT.  The chamber is a Helmholtz resonator:
14//     w = c * sqrt(A_neck / (V_chamber * L_neck))
15// A dog's bite REMOVES volume from the chamber, and w rises as V falls.  So the sound a dog produces
16// by biting is not a fixed beep -- it is a rising squeal whose contour is drawn by the dog's own jaw.
17// Prey distress calls rise in pitch too.  That correspondence is the reason this object can sound like
18// something worth catching instead of like a rubber toy, and it is measurable: the gate checks that a
19// harder bite raises the pitch and that a bigger chamber lowers it, both in the direction the physics
20// demands rather than the direction I would like.
21//
22// THE REED IS THE SAME PHYSICS AS A VOCAL FOLD.  A squeaker tongue and a larynx are both elastic
23// valves driven unstably by the flow past them, so this organ does not implement a second oscillator:
24// it drives nx_vox_source with the chamber pressure and takes its flow back.  That reuse is not just
25// DRY, it is why the toy inherits nonlinear regimes -- the roughness, the pitch breaks and the
26// cycle-to-cycle variation that stop a dog habituating -- from a model already gated at 13/13.
27//
28// Integer fixed-point, deterministic.  license_tier: ORIGINAL  expect_exit: 0
29
30import "nx_syscalls.nx"
31import "nx_isqrt.nx"
32import "nx_vox_source.nx"
33
34const SQ_Q:   i64 = 65536
35const SQ_C:   i64 = 343000        // speed of sound, mm/s
36
37// ---------------------------------------------------------------- parameters, in buildable units
38const SQ_P_V0:        i64 = 0    // mm^3   chamber volume at rest
39const SQ_P_AN_X100:   i64 = 1    // mm^2   neck cross-section, x100
40const SQ_P_LN_X100:   i64 = 2    // mm     neck length, x100
41const SQ_P_DVDF_X100: i64 = 3    // mm^3/N volume displaced per newton of bite, x100 (shell durometer)
42const SQ_P_BITE_N:    i64 = 4    // N      peak bite force
43const SQ_P_RISE_MS:   i64 = 5    // ms     jaw closing time
44const SQ_P_HOLD_MS:   i64 = 6    // ms
45const SQ_P_REL_MS:    i64 = 7    // ms
46const SQ_P_ZETA:      i64 = 8    // Q16    acoustic damping of the resonator
47const SQ_P_REED_F0:   i64 = 9    // Hz     reed natural frequency
48const SQ_P_REED_DET:  i64 = 10   // Hz     second-reed detune -> nonlinear regimes
49const SQ_P_REED_ZETA: i64 = 11   // Q16
50const SQ_P_FLOWGAIN:  i64 = 12   // Q16    reed opening -> neck conductance
51const SQ_P_PGAIN:     i64 = 13   // Q16    chamber pressure -> reed drive
52const SQ_P_OVS:       i64 = 14
53const SQ_P_SEED:      i64 = 15
54const SQ_P_JITTER:    i64 = 16   // Q16
55const SQ_P_N:         i64 = 17
56
57func sq_defaults(p: *i64) -> i64 {
58    // A 65 mm chew ball with a small squeaker chamber, medium-soft shell.
59    p[SQ_P_V0]        = 12000     // mm^3, the squeaker chamber (not the whole ball)
60    p[SQ_P_AN_X100]   = 1200      // 12.00 mm^2 neck
61    p[SQ_P_LN_X100]   = 600       // 6.00 mm neck
62    p[SQ_P_DVDF_X100] = 1500      // 15.00 mm^3 per newton
63    p[SQ_P_BITE_N]    = 200       // N -- a medium gundog's carnassial bite
64    p[SQ_P_RISE_MS]   = 60
65    p[SQ_P_HOLD_MS]   = 120
66    p[SQ_P_REL_MS]    = 90
67    p[SQ_P_ZETA]      = 3276      // 0.05
68    p[SQ_P_REED_F0]   = 1800      // Hz
69    p[SQ_P_REED_DET]  = 0
70    p[SQ_P_REED_ZETA] = 4587      // 0.07
71    p[SQ_P_FLOWGAIN]  = 65536
72    // 128x, MEASURED not chosen: the chamber runs at ~0.018 in normalised units and the reed needs
73    // ~1.0 to leave its rest position, so this is a unit conversion between two subsystems.  Swept
74    // 1x..128x with the gate: reed closures per bite rise 1 -> 3 and then fall again, and the
75    // nonlinear regime is NOT recovered at any value -- see the note on Q7 in nx_squeaker_gate.
76    p[SQ_P_PGAIN]     = 3604480
77    p[SQ_P_OVS]       = 4
78    p[SQ_P_SEED]      = 20260725
79    p[SQ_P_JITTER]    = 1966      // 0.03
80    return 0
81}
82
83// ---------------------------------------------------------------- state
84const SQ_S_P:   i64 = 0     // chamber overpressure (Q16, normalised)
85const SQ_S_U:   i64 = 1     // neck volume flow (Q16)
86const SQ_S_VOL: i64 = 2     // current chamber volume, mm^3
87const SQ_S_N:   i64 = 3
88
89// Helmholtz angular frequency for the CURRENT volume, returned in Q16 radians per sub-sample.
90//   w = c * sqrt(A / (V * L))
91// Computed as sqrt(A*1e12/(V*L)) so the square root runs on an integer large enough to keep precision;
92// the result is then in units of 1e-6 per mm and scales straight into rad/s by the speed of sound.
93func sq_omega(p: *i64, vol: i64, sr: i64) -> i64 {
94    if vol <= 0 { return 0 }
95    let denom: i64 = vol * p[SQ_P_LN_X100]
96    if denom <= 0 { return 0 }
97    let z: i64 = (p[SQ_P_AN_X100] * 1000000000000) / denom
98    let rt: i64 = nx_isqrt(z)                       // sqrt(A/(V*L)) in 1e-6 per mm
99    // w_per_sample_q16 = c * rt * 1e-6 * 65536 / sr
100    return (SQ_C * rt * SQ_Q) / (1000000 * sr)
101}
102
103// Bite force at sub-sample i, in Q16 NEWTONS.  Trapezoid: the jaw closes, holds, releases.
104//
105// The Q16 matters.  Returning whole newtons makes the force a staircase that steps by 1 N, so the air
106// being pumped -- which is the DERIVATIVE of that force -- is zero on almost every sample and a spike
107// on the rest.  The toy then emits a train of clicks instead of a squeal, and the reed, seeing no
108// sustained pressure, flutters twice in 300 ms.  A 200 N bite over 60 ms advances by 0.035 N per
109// sample; anything coarser than that is not a bite, it is a hammer.
110func sq_bite(p: *i64, i: i64, sr: i64) -> i64 {
111    let peak: i64 = p[SQ_P_BITE_N] * SQ_Q
112    let rise: i64 = (p[SQ_P_RISE_MS] * sr) / 1000
113    let hold: i64 = (p[SQ_P_HOLD_MS] * sr) / 1000
114    let rel:  i64 = (p[SQ_P_REL_MS] * sr) / 1000
115    if i < rise {
116        if rise <= 0 { return peak }
117        return (peak / rise) * i + ((peak % rise) * i) / rise
118    }
119    if i < rise + hold { return peak }
120    if i < rise + hold + rel {
121        if rel <= 0 { return 0 }
122        let left: i64 = rise + hold + rel - i
123        return (peak / rel) * left + ((peak % rel) * left) / rel
124    }
125    return 0
126}
127
128// Render `n` samples of RADIATED sound.  Radiation from a small aperture goes as the time derivative
129// of the volume flow, so the output is dU/dt rather than U itself -- using U directly would publish a
130// sound with far too much low end and would flatter every spectral axis in the ruler.
131//
132// diag[0] = peak chamber pressure, diag[1] = final volume, diag[2] = reed closures.
133const SQ_D_PPEAK: i64 = 0
134const SQ_D_VEND:  i64 = 1
135const SQ_D_CYC:   i64 = 2
136const SQ_D_N:     i64 = 3
137
138func sq_render(out: *i64, n: i64, rate: i64, p: *i64, diag: *i64) -> i64 {
139    let s: *i64 = sys_mmap(SQ_S_N * 8) as *i64
140    let vp: *i64 = sys_mmap(VF_P_N * 8) as *i64
141    let vs: *i64 = sys_mmap(VF_S_N * 8) as *i64
142    let vk: *i64 = sys_mmap(VF_K_N * 8) as *i64
143
144    // the reed, expressed as a vocal-fold pair: same equations, different tissue
145    vfold_defaults(vp)
146    vp[VF_P_F0_L]   = p[SQ_P_REED_F0]
147    vp[VF_P_F0_R]   = p[SQ_P_REED_F0] + p[SQ_P_REED_DET]
148    vp[VF_P_ZETA]   = p[SQ_P_REED_ZETA]
149    vp[VF_P_OVS]    = 1                       // this organ owns the oversampling loop
150    vp[VF_P_SEED]   = p[SQ_P_SEED]
151    vp[VF_P_JITTER] = p[SQ_P_JITTER]
152    vp[VF_P_PS]     = 0
153    vp[VF_P_G0_2]   = (vp[VF_P_G0] * 3) / 5      // convergent rest shape: the reed self-starts
154    vfold_init(vp, vk, vs, rate * p[SQ_P_OVS])
155
156    s[SQ_S_P] = 0
157    s[SQ_S_U] = 0
158    s[SQ_S_VOL] = p[SQ_P_V0]
159    diag[SQ_D_PPEAK] = 0
160    diag[SQ_D_CYC] = 0
161
162    let ovs: i64 = p[SQ_P_OVS]
163    let sr: i64 = rate * ovs
164    var prevU: i64 = 0
165    var prevDv: i64 = 0
166    var wasopen: i64 = 0
167    var vmin: i64 = p[SQ_P_V0]
168    var sub: i64 = 0
169    var i: i64 = 0
170    while i < n {
171        var acc: i64 = 0
172        var o: i64 = 0
173        while o < ovs {
174            // ---- the jaw.  Volume displaced tracks force; its RATE of change is the air being pumped.
175            let f: i64 = sq_bite(p, sub, sr)
176            let dv: i64 = (p[SQ_P_DVDF_X100] * f) / 100          // Q16 mm^3 displaced
177            var vol: i64 = p[SQ_P_V0] - (dv / SQ_Q)
178            if vol < p[SQ_P_V0] / 20 { vol = p[SQ_P_V0] / 20 }   // the shell cannot be crushed to nothing
179            s[SQ_S_VOL] = vol
180            if vol < vmin { vmin = vol }
181            let qin: i64 = dv - prevDv                            // Q16 mm^3 per sub-sample
182
183            // ---- the reed, driven by the chamber it sits in.  Clamped to a sane operating range:
184            //      the fold model expects a pressure near 1.0 and rails above ~16, and an unclamped
185            //      chamber drove it to 19 -- the reed stuck wide open and fluttered twice in 300 ms.
186            // s[SQ_S_P] is Q32; VF_P_PS is Q16.  Shifting down is NOT optional: when the resonator was
187            // moved to Q32 this line was left as-is, so the reed was handed a pressure 65536x too large
188            // and sat pinned at the clamp for every sample of every run.  That single missing shift is
189            // what made SQ_P_PGAIN a no-op across a 64x sweep, and it put the whole toy into a
190            // chamber-refill relaxation oscillation near 110 Hz whose frequency was independent of the
191            // reed's tuning -- the reed could not flutter because it was never allowed to move freely.
192            var pd: i64 = ((s[SQ_S_P] >> 16) * p[SQ_P_PGAIN]) / SQ_Q
193            if pd < 0 { pd = 0 }
194            // Ceiling sits below the fold model's own divergence clamp (16.0) so the reed cannot rail,
195            // but NOT so low that it swallows the parameter: at 4.0 the chamber saturated it on every
196            // sample and a 64x sweep of SQ_P_PGAIN produced byte-identical audio -- a live-looking knob
197            // wired to nothing.  A clamp that eats its input is a bug, not a safety feature.
198            if pd > 12 * SQ_Q { pd = 12 * SQ_Q }
199            vp[VF_P_PS] = pd
200            let rflow: i64 = vfold_substep(vs, vp, vk)
201            // reed opening gates the neck: no opening, no path out
202            let cond: i64 = (rflow * p[SQ_P_FLOWGAIN]) / SQ_Q
203
204            // ---- Helmholtz resonator at the CURRENT volume.  This is where the glide comes from.
205            //      P and U are carried at Q32, not Q16.  At Q16 the pressure increment per sub-step
206            //      truncates to the integer 1 for a chamber this size, so the resonator has barely one
207            //      significant figure and the toy renders as a handful of quantisation steps -- the
208            //      same precision failure the vocal-fold velocity had, in a different organ.
209            let w: i64 = sq_omega(p, vol, sr)
210            let r: i64 = (2 * p[SQ_P_ZETA] * w) / SQ_Q
211            // P' = w^2 * (qin - U) ; U' = P - r*U      (=> U'' + rU' + w^2 U = w^2 qin)
212            let net: i64 = (qin * SQ_Q) - s[SQ_S_U]
213            let dp: i64 = ((w * w) * (net >> 16)) >> 16
214            s[SQ_S_P] = s[SQ_S_P] + dp
215            var du: i64 = s[SQ_S_P] - ((r * (s[SQ_S_U] >> 16)))
216            // the reed valve scales the flow that can actually leave
217            du = (du / (SQ_Q + cond)) * cond
218            s[SQ_S_U] = s[SQ_S_U] + du
219
220            // NEVER-BRICK ANALOG: bound the acoustic state so a bad geometry produces a quiet wrong
221            // answer rather than an unbounded one.
222            let acl: i64 = 4096 * SQ_Q * SQ_Q
223            if s[SQ_S_P] > acl { s[SQ_S_P] = acl }
224            if s[SQ_S_P] < 0 - acl { s[SQ_S_P] = 0 - acl }
225            if s[SQ_S_U] > acl { s[SQ_S_U] = acl }
226            if s[SQ_S_U] < 0 - acl { s[SQ_S_U] = 0 - acl }
227
228            var pa: i64 = s[SQ_S_P] >> 16
229            if pa < 0 { pa = 0 - pa }
230            if pa > diag[SQ_D_PPEAK] { diag[SQ_D_PPEAK] = pa }
231            var isopen: i64 = 0
232            if rflow > 0 { isopen = 1 }
233            if wasopen == 1 { if isopen == 0 { diag[SQ_D_CYC] = diag[SQ_D_CYC] + 1 } }
234            wasopen = isopen
235
236            acc = acc + ((s[SQ_S_U] - prevU) >> 16)   // dU/dt = radiated pressure
237            prevU = s[SQ_S_U]
238            prevDv = dv
239            sub = sub + 1
240            o = o + 1
241        }
242        out[i] = acc / ovs
243        i = i + 1
244    }
245    diag[SQ_D_VEND] = vmin
246    return 0
247}
248
249// Render and normalise to a 16-bit audio range.  Returns the raw peak (0 = the toy made no sound).
250func sq_render_audio(out: *i64, n: i64, rate: i64, p: *i64, amp: i64, diag: *i64) -> i64 {
251    sq_render(out, n, rate, p, diag)
252    var sum: i64 = 0
253    var i: i64 = 0
254    while i < n { sum = sum + out[i]; i = i + 1 }
255    let dc: i64 = sum / n
256    var mx: i64 = 0
257    i = 0
258    while i < n {
259        out[i] = out[i] - dc
260        var a: i64 = out[i]
261        if a < 0 { a = 0 - a }
262        if a > mx { mx = a }
263        i = i + 1
264    }
265    if mx <= 0 { return 0 }
266    i = 0
267    while i < n { out[i] = (out[i] * amp) / mx; i = i + 1 }
268    return mx
269}
270
271// The Helmholtz frequency this geometry gives at rest and at full bite, in Hz.  Published alongside
272// any design so the glide is a stated number rather than a claim.
273// Computed straight in Hz from the geometry: f = c*sqrt(A/(V*L))/(2*pi).  Deriving it from the
274// per-sample Q16 omega instead loses the whole value to integer division and silently reports 0 Hz.
275func sq_f_at(p: *i64, vol: i64) -> i64 {
276    if vol <= 0 { return 0 }
277    let denom: i64 = vol * p[SQ_P_LN_X100]
278    if denom <= 0 { return 0 }
279    let z: i64 = (p[SQ_P_AN_X100] * 1000000000000) / denom
280    let rt: i64 = nx_isqrt(z)                       // sqrt(A/(V*L)) in 1e-6 per mm
281    return (SQ_C * rt) / 6283185                    // c*rt/1e6/(2*pi), folded into one divisor
282}
283func sq_vol_bitten(p: *i64) -> i64 {
284    let dv: i64 = (p[SQ_P_DVDF_X100] * p[SQ_P_BITE_N]) / 100
285    var vol: i64 = p[SQ_P_V0] - dv
286    if vol < p[SQ_P_V0] / 20 { vol = p[SQ_P_V0] / 20 }
287    return vol
288}
289func sq_f_rest(p: *i64) -> i64 { return sq_f_at(p, p[SQ_P_V0]) }
290func sq_f_bitten(p: *i64) -> i64 { return sq_f_at(p, sq_vol_bitten(p)) }
nx_squeaker.nx source

`nx_squeaker.nx` source
