GameCult/AquaSynth
Watch
9
0
Fork
You've already forked AquaSynth
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Articulatory tract: Story vowels, nasal branch, stress cue, plosive burst, independent Burg check (unheard) #19

Merged
aqua merged 20 commits from aqua/tract-articulatory into main 2026-11-10 18:13:33 +00:00
Conversation 1 Commits 20 Files changed 25 +1964 -156
25 changed files with 1964 additions and 156 deletions
Download patch file Download diff file Expand all files Collapse all files

135
docs/research/articulatory-vocal-tract-sources.md Normal file
View file

@ -0,0 +1,135 @@
# Articulatory vocal tract: sources and what AquaSynth takes from them
Read 2026-11-08 for follow-up `aqua-body:follow_up:chat-s4456-4456`. The full fact sheet with quoted values and the list of
what could not be retrieved is the Eyes pass `tract-literature-20261108` (research notes in the aqua mind;
working copy `work/tract-literature-facts.md` in the aqua workspace). This page keeps only what the tract design rests on.
Nothing built from it has been heard by anyone; every number below is a measurement or a published value, not a listening result.
## Sources and what each decides
| Source | What we take | What we do not take |
|---|---|---|
| Kelly & Lochbaum 1962; Smith, *Physical Audio Signal Processing* (KL junctions); Liljencrants 1985 / Story 1995 (section length `c / 2Fs`, i.e. one sample per section); Mullen, Howard & Murphy 2006 (boundary reflections, wall loss 0.90-0.98 sets bandwidth) | The ladder: one cell per sample, junction `r = (A[i-1]-A[i])/(A[i-1]+A[i])`, glottal reflection about +0.75, lip reflection about -0.85, a per-pass loss factor for bandwidth | Fractional-delay cells (Välimäki & Karjalainen, Mathur & Story): needed for continuous length change, not built; the tract length is fixed per patch |
| Pink Trombone (Neil Thapen 2017, MIT; read in chdh/pink-trombone-mod `Tract.ts`) | Constants (glottal 0.75, lip -0.85, 0.999 per pass), the 3-port velum junction `rL=(2A_k-S)/S`, `rR`, `rN`, with `S = A_k + A_(k+1) + A_velum` | Its 44-cell tract at twice the sample rate and its Gaussian tongue model as the vowel source: the vowel shapes now come from measured area functions |
| Story, Titze & Hoffman 1996/1998 (MRI area functions; numeric tables paywalled) and **Story 2008, JASA 123(1):327-335, PMC2377017, Table I** (44-section area functions of 11 vowels incl. /ɪ/ /ʊ/, one adult male, 2002 scan) | The 11 vowel area functions, as data: `docs/research/data/story2008-area-functions.csv` | Story's 1994 set (plot only). The sizes are one speaker's, not a population mean (see the measured gap below) |
| Story 2005, JASA 117(5):3231-3254 (area function model: vowel substrate `(pi/4)[Omega + q1 phi1 + q2 phi2]^2`, consonant constriction `C(i) = 1 - exp(-ln16((i-lc)/rc)^2)`, nasal port at about section 22) | The principle: a vowel is a point in a low-dimensional space of whole area functions, blended in diameter space, with a consonant constriction laid on top of it. We blend the 11 measured shapes by weights instead of PCA modes (Omega, phi1, phi2 are in the paper's Table III; not transcribed) | Its 14-parameter control set |
| VocalTractLab (Birkholz; GPL; 2.4 manual, 2005 dissertation) | Nasal cavity area function (Dang & Honda subject, 19 sections x 0.6 cm = 11.4 cm, 0.76-4.24 cm^2); the observation that nasal antiresonances come from the oral cavity as a side branch of the open velum; c = 350 m/s; the wall-loss and radiation (R||L) constants as the reference for the independent check | Its 3D articulators and glottis models; the four paranasal sinus Helmholtz resonators (zeros at 552-1682 Hz) and piriform fossae are not modelled |
| Maeda 1979/1990 (seven jaw/tongue/lip/larynx parameters from cineradiographs; `A = alpha d^beta`) | Structure only: the argument for articulator-level parameters over formant knobs | Per-vowel parameter values (not retrievable); no Maeda parameterisation is built |
| Nasal coupling: Fujimura 1962 (nasal murmur: very low N1 near 250-400 Hz, zeros), Dang & Honda 1994/1995, Story 2005 tier IV | A velum that opens the 3-port junction; the murmur's low first formant and the oral side-branch zero as the measured behaviour to look for | Per-consonant zero ranges for /m n ŋ/: no readable primary source |
| LF glottal flow (Fant, Liljencrants & Lin 1985; Fant 1995 `Rd`; Gobl 2021; Wang & Gobl 2023 regressions `Rkp = 0.118 Rd + 0.224`, `Rap = 0.048 Rd - 0.01`) | Reference for the source spectrum: Rd is a tension axis (low = tense), tilt cutoff falls as Rd rises. The current source is a band-limited saw through a one-pole lowpass whose cutoff is the `tilt` control: it matches the LF pulse's spectral tilt only qualitatively | An LF pulse generator: not built. If the source is moved to LF, `Rd` replaces `tilt` |
| Hillenbrand et al. 1995, JASA 97(5):3099-3111 (male means computed from the public data of github.com/santiagobarreda/hillenbrand_et_al_1995, MIT) and Peterson & Barney 1952 (means from phonTools `pb52`) | The reference formants F1/F2/F3 per vowel for the verification, table below | Table V itself was not read; P&B vowel labels are inferred from the data codes |
| Stress: Sluijter & van Heuven 1996 (ICSLP table + JASA abstract); Fry 1955/58, Lieberman 1960 (secondary reports only) | Stressed syllables are longer (stressed 220 ms against 160 ms unstressed in the initial-stress pairs; a stressed/unstressed ratio of 1.68 in a native-speaker study), 1-3 dB louder overall, f0 higher by about 6 percent when unaccented; duration is the most reliable cue, f0 the weakest in unaccented words | Fry's and Lieberman's own magnitudes (not retrieved) |
| Plosives: Lisker & Abramson 1964 (VOT, secondary table), Chodroff & Wilson 2014 (burst spectra, /b/ COG 2833 Hz, peak 1643 Hz), Stevens & Blumstein (labial onset spectrum diffuse-falling) | A /b/ burst is a short low-frequency-weighted transient with a short voice-onset lag (0-5 ms mean 1 ms for English /b/) | A /b/ voice-bar level in dB and burst duration norms: no source retrieved, so the AquaSynth levels for them are design choices, labelled as such |
## Male reference formants used for verification (Hz)
Hillenbrand 1995 male means (45 men, steady state; computed from the public data): /i/ 343/2323/3001, /ɪ/ 429/2034/2687,
/e/ 476/2090/2692, /ɛ/ 588/1803/2604, /æ/ 591/1930/2595, /ɑ/ 757/1309/2535, /ɔ/ 656/1023/2521, /o/ 498/910/2459,
/ʊ/ 469/1123/2435, /u/ 380/992/2355, /ʌ/ 621/1181/2548. Peterson & Barney 1952 male means (33 men): /i/ 267/2294/2937,
/ɪ/ 392/1993/2569, /ɛ/ 526/1854/2481, /æ/ 664/1727/2420, /ʌ/ 631/1192/2377, /ɑ/ 718/1091/2442, /ɔ/ 568/836/2403,
/ʊ/ 437/1023/2245, /u/ 307/876/2239.
## The gap the area functions must be read against
Story's shapes belong to one speaker. Computed independently (frequency-domain tube chain, `tools/formant-check/tube_formants.py`,
no AquaSynth code), the tabulated area functions at their own lengths give for example /i/ 261/2139, /ʌ/ 561/827, /ɑ/ 666/933,
/u/ 243/684 Hz, against the population means above (/i/ 343/2323, /ʌ/ 621/1181, /ɑ/ 757/1309, /u/ 380/992). So a render that
follows Story's area functions can be correct and still sit 10-30 percent from a population mean. The verification therefore
reports two distances: render against the area function's own predicted formants (does the ladder implement the tube), and
render against the population means (is the speaker typical). Only the first is a defect test.
## Ladder against the tube chain (measured 2026-11-08, unheard)
The ladder's lip termination is the open-end reflectance R(z) = -g (1-b)/(1-b z^-1) (Smith PASP; Valimaki and Karjalainen 1994; Levine and
Schwinger): g = 1 is the radiation load's `reflection` magnitude (sign stays negative), and b follows the lip cell's area, with a = sqrt(A/pi) cm,
q = (a fs / c)^2, c = 35000 cm/s, b/(1-b)^2 = q, b = ((2q+1) - sqrt(4q+1)) / (2q), clamped to [0, .9]. The low-frequency loss 1-|R|^2 is then
(ka)^2 and the DC group delay b/(1-b) samples is the end correction. A constant -.82 reflection, the earlier termination, had no F2 peak for
the rounded and back vowels and no end correction.
`tools/formant-check/ladder_formants.py --filtered` runs the same ladder in numpy (24 cells, 17.52 cm, 350 m/s, loss .998, glottal reflection .75).
Against the independent tube chain (`tube_formants.py`, 17.52 cm; baffled-piston R||L radiation) the numpy ladder reads at g = 1 (F1/F2 error %):
i +1.4/+0.4, ih +1.7/+1.5, e +2.3/+4.2, eh +4.0/+3.1, ae +2.9/+3.5, uh +5.6/+4.3, ah +5.1/+6.6, aw +8.6/+3.7, o +8.2/+5.0, oo_short +5.5/+8.2,
u +6.0/+8.2. Lowering g does not help (g = .95: aw F1 +9.6; g = .9: the aw F2 peak is lost). The numpy ladder therefore reads systematically high
against the tube chain: F1 +4 to +9%, F2 +3 to +8% for the back vowels. That bias is diagnosed, not isolated: raising the glottal reflection to .9 moves only /aw/ F1
(8.6 to 6.2%); a 1.4 times larger end correction leaves a worst case of 6.9%; the remaining suspects (cell-length quantisation of a 24.03-cell
tract, and that the reference itself uses a different radiation model) are untested. It is the open item.
Held tolerances (`VocalTractFormantTests`, harmonic-line reading of the 10 Hz render): the front vowels i ih e eh ae within 6% on F1 and F2;
the back and rounded vowels within 10%, which is the measured disagreement between two lip radiation models (the ladder's one-pole
Levine-Schwinger-form reflectance, the tube chain's baffled piston R||L) plus the ladder's cell quantisation. All eleven meet their bands. For
a while uh, ah, aw, o and oo_short read F1 +14.8, +10.8, +12.4, +25.8 and +13.0% against the tube chain and the tests were skipped with those
numbers; the numpy ladder read those F1 20-85 Hz lower than the render. The cause was the velopharyngeal junction in `FaustExport.VocalLadder.cs`:
it summed the areas of oral cells 11 and 12 and used their reflection coefficients for the junction that sits between cells 10 and 11, so the
pharynx (cells 10-12 of /o/: .52 .33 .52 cm^2) carried the next step's reflection with the wrong sign. Isolated by running the lip wave of the rendered
ladder under an impulse against the numpy copy (the first echo already differed: -0.30 against -0.23), with the velum shut. Fixed, the render agrees
with the numpy ladder to 10 Hz (o 520/830, uh 580/840, ah 690/980, aw 550/720, oo_short 500/920, u 290/820 at the harmonic lines) and the five tests
run. Read by Burg LPC on the renders (Soul pass of 2026-11-09, vowel window about 100-200 ms), F1 error against the tube chain fell from 5c32bf6 to the fix 05efdc9 as uh +5.3 to -0.5, ah +2.5 to -4.3,
aw +4.2 to +0.7, o +16.7 to +5.7, oo_short +12.7 to +4.9 percent (a reconcile rerun on the same renders, 2026-11-10, reproduces these to within 1 point in the 100-200 ms window, aw least: +3.6 to -0.1).
The figures depend on the window: in the 70-150 ms window of the table below the same renders read uh +4.4 to -2.1, ah -0.6 to -6.6, aw +1.3 to -3.4, o +16.3 to +2.3, oo_short +8.8 to +1.7, and the two vowels
the fix did not repair, u and i, read +10.7 to +13.8 and +15.6 to +17.1 (u +4.6 to +6.4 and i +10.8 to +12.2 at 100-200 ms). So the junction fix removed the large positive bias of o and oo_short in either window (ah turns from near zero or slightly high to 4-7 percent low),
and the earlier "-4.3..+5.9 after" range is the five repaired vowels plus u read at one window: it does not cover i and u, which stay 5-17 percent high by Burg and are the LPC's resolution limit at 140 Hz
(see the table note); ah undershoots by 4 to 7 percent after the fix. The earlier "12-26 percent high becomes 1.5-9 percent high" came from the harmonic-line picker, not from LPC. What remains against the tube chain
(F1 +1.5 to +9%, F2 +2 to +8% at the harmonic lines) is the two radiation models' disagreement above.
Burg LPC (`burg_formants.py`) on the WordRender renders at 140 Hz (held vowels): /i/ in nibu 279/2105/2924, /u/ in nibu 283/836/2637, /ih/ in bin
359/2012/2536, /oo_short/ (bʊn) in bun 543/899/2781, against the tube chain i 251/2063/2877, u 271/761/2574, ih 346/1987/2491, oo_short 478/850/2760 and
Hillenbrand male means i 343/2323/3001, u 380/992/2355, ih 429/2034/2687, oo_short 469/1123/2435.
## Not retrieved (from the fact sheet)
Story-Titze-Hoffman 1996/1998 numeric tables; Fant 1960, Stevens-House, Wakita, Ladefoged tables; Maeda per-vowel values;
Dang & Honda 1994 text; per-consonant nasal zero ranges; Fry and Lieberman magnitudes; Klatt 1975; burst duration norms and the
/b/ voice-bar level; Liljencrants 1985; the JOS book tract chapters. The Story 2008 table came through a page summariser, not a PDF
text extraction; every column was checked by computing its formants (all eleven give vowel-typical F1/F2 orderings), which
catches a gross transcription error, not a small one.
## Render verification, stress cue and plosive burst (2026-11-08, unheard)
Renders: `WordRender` at f0 140 Hz, one 180 ms vowel per file (`v_<id>`), plus /ɑmɑ/, /ɑnɑ/, /ˈnibu/ with and without primary stress, /bɪn/, /bʊn/.
Estimator: `tools/formant-check/verify_renders.py` over `burg_formants.py` (Burg LPC order 12 at 10 kHz, 30 ms frames, vowel window 70-150 ms), a separate
Python process that shares no code with the C# tests or the ladder; references are Hillenbrand 1995 and Peterson & Barney 1952 male means (above) and the
tube chain at 17.52 cm (`tube_formants.py`). Files: `/opt/commons/handoff/aqua-tract-articulatory-20261108` (wav/, scripts/, manifest.json, verify.txt, verify.json).
| vowel | Burg F1/F2 | tube chain F1/F2 | render vs tube F1/F2 % | Hillenbrand F1/F2 | render vs Hillenbrand % |
|---|---|---|---|---|---|
| i | 294/2080 | 251/2063 | +17.0/+0.8 | 343/2323 | -14.3/-10.4 |
| ɪ | 375/2024 | 346/1987 | +8.5/+1.8 | 429/2034 | -12.5/-0.5 |
| e | 448/1925 | 432/1827 | +3.7/+5.4 | 476/2090 | -5.9/-7.9 |
| ɛ | 551/1603 | 527/1571 | +4.5/+2.0 | 588/1803 | -6.3/-11.1 |
| æ | 739/1550 | 729/1523 | +1.4/+1.8 | 591/1930 | +25.0/-19.7 |
| ʌ | 538/822 | 549/809 | -2.1/+1.6 | 621/1181 | -13.4/-30.4 |
| ɑ | 607/950 | 650/913 | -6.6/+4.0 | 757/1309 | -19.8/-27.4 |
| ɔ | 499/715 | 516/692 | -3.4/+3.4 | 656/1023 | -24.0/-30.1 |
| o | 488/806 | 477/787 | +2.3/+2.4 | 498/910 | -2.0/-11.4 |
| ʊ | 486/894 | 478/850 | +1.7/+5.2 | 469/1123 | +3.7/-20.4 |
| u | 308/811 | 271/761 | +13.8/+6.6 | 380/992 | -18.9/-18.2 |
(Table re-rendered 2026-11-09 after the junction fix, the nasal and stress changes below; the first table, 2026-11-08, had o, ʊ and u 9-16% high.)
Against its own area function the render sits within 6.6 percent on F1 for every vowel but /i/ (+17%) and /u/ (+13.8%), whose F1 near 270-300 Hz is under the LPC's resolution at a
140 Hz line spacing (the 10 Hz harmonic-line reading in `VocalTractFormantTests` puts them at 260 and 290 Hz, +3.6 and +7.0%); F2 of all eleven is within 6.6 percent. Against the population
means the render is 8-30 percent off for the back vowels' F2, which is the one-speaker gap described above, not a defect test.
Nasal (steady parts: the consonant's middle 60 percent against the mean power of the flanking vowels, 70 ms before and 90 ms from 30 ms after): /ɑmɑ/ -11.8 dB, /ɑnɑ/ -9.6 dB, Burg N1 295 and 332 Hz,
inside Fujimura's low-N1 range. The target band is 8 to 16 dB under the vowels, a DESIGN value: no primary level was retrieved (secondary accounts of Fujimura 1962 and later put nasal murmurs
well under their vowels). Before: -13.2 and -3.8 dB by the same measure (the old whole-segment rms read -2.8 and +1.8 dB). Three causes, in order: the /n/ oral closure began to open at half the gesture,
a third before the velum and the contact (held to 0.8 now); then, with the junction on its own interface the nose took its share and both murmurs rose to -4.2 and -2.2 dB; the nasal path now carries
the damping of the nasal cavity (soft walls, turbinates, sinuses): loss .975 per cell, a bandwidth of about 380 Hz at this sample rate (`IpaGestureExperiment.NasalPathLoss`; the bandwidth is a design value,
none retrieved). /n/ stays 2 dB louder than /m/: the cause is not isolated (a shorter oral side branch is the physical candidate).
Stress (cue values in `WeksaWordPlan.CueFor`): primary stress renders longer by 75 ms for the vowel (design 68 ms), f0 148.1 against 139.9 Hz (ratio 1.059; design 1.06). Level: +2.0 dB mean over the five vowels
i a u e o (`PrimaryStressRendersAboutTwoDecibelsLouderOnAverageOverVowels`, 2.05 dB in the last run) with the amplitude factor at 1.30 (was 1.26: +1.7 dB). The f0 rise alone moves the harmonics against F1 and
changes the level by -1.3 to +1.8 dB by vowel, so a single vowel is a noisy measure of the amplitude cue: /i/ in /nibu/ (steady window 150-250 ms) reads +0.21 dB in the last render (+0.99 before), isolated /i/
+0.79 dB, /e/ +4.1 dB. The mean is what is held; the scatter between vowels is the harmonic/formant interaction of a 140 Hz voice, not the cue.
Plosive: word-initial /b/ closure is the tract wall's voice bar, -29.2 dB re the following vowel, /p/ -40.1 dB (`PlosiveClosureCarriesTheVoiceBarOnlyWhenVoiced`); /bɪn/ closure peak 0.0036 (was 0.0000, below 16-bit
resolution), release 60-95 ms peak 0.0695, vowel 0.129. The wall radiation reads the first-cell pressure of the ladder (R0 + L0), band-limits it to 90-300 Hz and adds it at gain 0.7 (a DESIGN value for a voice-bar
level of -30 dB; the primary level was not retrieved, and the 2.4 manual's skin radiation, Birkholz & Zhang 2025, is the reference to read). In a vowel the term moves the f0 line by under 1 dB; the 200 and 300 Hz harmonics move up to +1.5 dB on one open vowel (wall-only measurement, Soul pass). The band is two second-order sections, not a brick wall: with the mouth and the velum shut the bar sits about 25 dB down above 300 Hz. The burst
itself (difference of renders with and without `release_flow`) peaks at about 0.006-0.008, -24 dB re the vowel peak; the release peak in the table above is the voicing breaking through as the contact opens. Burst gain 0.06 is a DESIGN VALUE.
Weak spots found in the renders: (1) the vowel level is not steady, it rises for about 100 ms after onset (ɑ: rms 0.028 at 80 ms to 0.079 at 160 ms) so any window early in a vowel
reads low; the cause is not isolated (surface smoothing and the pressure ramp are the suspects). (2) /n/ murmur 2 dB above /m/ (above). (3) The voice bar's level is a design value and its band
(90-300 Hz) assumes f0 of 90 Hz and up. (4) Only f0 and level of the first vowel are checked per cue; the cues' perceptual strength is unknown. Everything here is a measurement, nothing was heard.
Mutations (hand-applied, restored, each run against the named tests): pitch ratio ignored, burst removed, burst not shaped by opening, burst gain x20 (0.06 to 1.2), stress duration ignored,
stress pitch ignored, stress amplitude ignored, secondary equal to primary, probe timeline sample rate re-hardcoded to 48000: nine of nine killed. Script: `work/mut-b3.sh` in the aqua workspace.

49
docs/research/data/story2008-area-functions.csv Normal file
View file

@ -0,0 +1,49 @@
# Story 2008, JASA 123(1):327-335, doi 10.1121/1.2805683, PMC2377017, Table I. Areas in cm^2, section 1 = glottis, section 44 = lips.
# Transcribed by an Eyes pass from the PMC page through a page summariser (not pdftotext). tools/formant-check/tube_formants.py checks every column for plausible formants.
section,i,ɪ,e,ɛ,æ,ʌ,ɑ,ɔ,o,ʊ,u
1,0.51,0.28,0.29,0.37,0.31,0.23,0.56,0.27,0.38,0.37,0.54
2,0.59,0.21,0.26,0.30,0.21,0.34,0.62,0.43,0.45,0.38,0.61
3,0.62,0.21,0.30,0.24,0.18,0.47,0.66,0.54,0.57,0.49,0.66
4,0.72,0.30,0.40,0.23,0.23,0.60,0.78,0.67,0.77,0.62,0.75
5,1.24,0.47,0.55,0.29,0.33,0.77,0.97,0.83,1.31,0.85,1.13
6,2.30,0.71,0.74,0.41,0.50,1.06,1.16,0.92,1.92,1.28,1.99
7,3.30,1.12,0.99,0.58,0.78,1.26,1.12,0.89,1.74,1.62,2.83
8,3.59,1.48,1.09,0.82,0.96,1.09,0.82,0.73,1.11,1.47,2.90
9,3.22,1.35,0.90,0.97,0.85,0.80,0.55,0.55,0.75,1.04,2.52
10,2.86,1.05,0.69,0.82,0.63,0.65,0.45,0.44,0.59,0.81,2.40
11,3.00,0.92,0.77,0.62,0.46,0.56,0.37,0.37,0.57,1.03,2.83
12,3.61,0.92,1.31,0.54,0.36,0.47,0.29,0.28,0.68,1.44,3.56
13,4.39,1.19,2.13,0.54,0.33,0.37,0.21,0.18,0.73,1.49,3.99
14,4.95,1.94,2.74,0.68,0.46,0.24,0.15,0.14,0.67,1.28,3.89
15,5.17,2.83,3.03,1.09,0.73,0.17,0.16,0.15,0.58,1.06,3.50
16,5.16,3.31,3.23,1.62,1.00,0.18,0.25,0.14,0.49,0.85,3.04
17,5.18,3.48,3.33,2.03,1.30,0.23,0.34,0.13,0.44,0.69,2.64
18,5.26,3.60,3.27,2.35,1.66,0.27,0.43,0.14,0.42,0.56,2.44
19,5.20,3.64,3.09,2.51,1.97,0.28,0.54,0.18,0.49,0.42,2.31
20,5.02,3.49,2.84,2.39,2.06,0.29,0.61,0.20,0.53,0.27,2.07
21,4.71,3.20,2.66,2.22,2.03,0.33,0.67,0.23,0.38,0.27,1.80
22,4.13,2.90,2.46,2.13,2.01,0.52,0.98,0.49,0.30,0.41,1.52
23,3.43,2.59,2.14,2.00,1.89,0.97,1.76,1.03,0.45,0.51,1.14
24,2.83,2.21,1.79,1.78,1.66,1.50,2.75,1.58,0.61,0.49,0.74
25,2.32,1.87,1.44,1.58,1.49,1.91,3.52,2.06,0.71,0.47,0.42
26,1.83,1.54,1.17,1.43,1.42,2.23,4.08,2.61,0.79,0.50,0.22
27,1.46,1.20,1.00,1.31,1.37,2.65,4.74,3.35,0.86,0.58,0.14
28,1.23,0.92,0.88,1.23,1.34,3.29,5.61,4.34,1.01,0.80,0.20
29,1.08,0.74,0.80,1.24,1.41,4.13,6.60,5.51,1.41,1.19,0.47
30,0.94,0.59,0.81,1.38,1.58,5.00,7.61,6.70,2.09,1.62,0.89
31,0.80,0.52,0.85,1.61,1.82,5.77,8.48,7.75,3.00,2.27,1.15
32,0.67,0.54,0.84,1.82,2.19,6.33,9.06,8.63,4.10,3.24,1.42
33,0.55,0.59,0.86,1.96,2.63,6.61,9.29,9.29,5.16,4.16,2.17
34,0.46,0.65,0.96,2.01,2.97,6.63,9.26,9.59,6.22,5.00,3.04
35,0.40,0.71,1.18,2.00,3.17,6.45,9.06,9.42,7.34,5.70,3.69
36,0.36,0.67,1.35,1.95,3.40,6.04,8.64,8.78,8.15,6.11,4.70
37,0.35,0.61,1.48,1.77,3.56,5.39,7.91,7.82,8.61,6.21,5.74
38,0.35,0.57,1.62,1.48,3.57,4.42,6.98,6.50,8.37,6.29,5.41
39,0.38,0.50,1.49,1.30,3.58,3.29,6.02,4.95,6.76,6.24,3.82
40,0.51,0.48,1.29,1.21,3.44,2.37,5.13,3.47,4.37,4.91,2.34
41,0.74,0.54,1.24,1.10,3.15,1.74,4.55,2.15,2.30,2.61,1.35
42,0.92,0.73,1.17,0.99,3.38,1.36,4.52,1.38,1.06,1.09,0.65
43,0.96,0.93,1.04,0.91,3.99,1.17,4.71,1.11,0.58,0.63,0.29
44,0.91,0.82,0.95,0.78,4.17,0.99,4.72,0.90,0.47,0.59,0.16
section_length_cm,0.384,0.376,0.386,0.393,0.366,0.390,0.388,0.395,0.417,0.440,0.445
vtl_cm,16.90,16.55,16.98,17.30,16.11,17.14,17.09,17.40,18.33,19.34,19.59
1 # Story 2008, JASA 123(1):327-335, doi 10.1121/1.2805683, PMC2377017, Table I. Areas in cm^2, section 1 = glottis, section 44 = lips.
2 # Transcribed by an Eyes pass from the PMC page through a page summariser (not pdftotext). tools/formant-check/tube_formants.py checks every column for plausible formants.
3 section,i,ɪ,e,ɛ,æ,ʌ,ɑ,ɔ,o,ʊ,u
4 1,0.51,0.28,0.29,0.37,0.31,0.23,0.56,0.27,0.38,0.37,0.54
5 2,0.59,0.21,0.26,0.30,0.21,0.34,0.62,0.43,0.45,0.38,0.61
6 3,0.62,0.21,0.30,0.24,0.18,0.47,0.66,0.54,0.57,0.49,0.66
7 4,0.72,0.30,0.40,0.23,0.23,0.60,0.78,0.67,0.77,0.62,0.75
8 5,1.24,0.47,0.55,0.29,0.33,0.77,0.97,0.83,1.31,0.85,1.13
9 6,2.30,0.71,0.74,0.41,0.50,1.06,1.16,0.92,1.92,1.28,1.99
10 7,3.30,1.12,0.99,0.58,0.78,1.26,1.12,0.89,1.74,1.62,2.83
11 8,3.59,1.48,1.09,0.82,0.96,1.09,0.82,0.73,1.11,1.47,2.90
12 9,3.22,1.35,0.90,0.97,0.85,0.80,0.55,0.55,0.75,1.04,2.52
13 10,2.86,1.05,0.69,0.82,0.63,0.65,0.45,0.44,0.59,0.81,2.40
14 11,3.00,0.92,0.77,0.62,0.46,0.56,0.37,0.37,0.57,1.03,2.83
15 12,3.61,0.92,1.31,0.54,0.36,0.47,0.29,0.28,0.68,1.44,3.56
16 13,4.39,1.19,2.13,0.54,0.33,0.37,0.21,0.18,0.73,1.49,3.99
17 14,4.95,1.94,2.74,0.68,0.46,0.24,0.15,0.14,0.67,1.28,3.89
18 15,5.17,2.83,3.03,1.09,0.73,0.17,0.16,0.15,0.58,1.06,3.50
19 16,5.16,3.31,3.23,1.62,1.00,0.18,0.25,0.14,0.49,0.85,3.04
20 17,5.18,3.48,3.33,2.03,1.30,0.23,0.34,0.13,0.44,0.69,2.64
21 18,5.26,3.60,3.27,2.35,1.66,0.27,0.43,0.14,0.42,0.56,2.44
22 19,5.20,3.64,3.09,2.51,1.97,0.28,0.54,0.18,0.49,0.42,2.31
23 20,5.02,3.49,2.84,2.39,2.06,0.29,0.61,0.20,0.53,0.27,2.07
24 21,4.71,3.20,2.66,2.22,2.03,0.33,0.67,0.23,0.38,0.27,1.80
25 22,4.13,2.90,2.46,2.13,2.01,0.52,0.98,0.49,0.30,0.41,1.52
26 23,3.43,2.59,2.14,2.00,1.89,0.97,1.76,1.03,0.45,0.51,1.14
27 24,2.83,2.21,1.79,1.78,1.66,1.50,2.75,1.58,0.61,0.49,0.74
28 25,2.32,1.87,1.44,1.58,1.49,1.91,3.52,2.06,0.71,0.47,0.42
29 26,1.83,1.54,1.17,1.43,1.42,2.23,4.08,2.61,0.79,0.50,0.22
30 27,1.46,1.20,1.00,1.31,1.37,2.65,4.74,3.35,0.86,0.58,0.14
31 28,1.23,0.92,0.88,1.23,1.34,3.29,5.61,4.34,1.01,0.80,0.20
32 29,1.08,0.74,0.80,1.24,1.41,4.13,6.60,5.51,1.41,1.19,0.47
33 30,0.94,0.59,0.81,1.38,1.58,5.00,7.61,6.70,2.09,1.62,0.89
34 31,0.80,0.52,0.85,1.61,1.82,5.77,8.48,7.75,3.00,2.27,1.15
35 32,0.67,0.54,0.84,1.82,2.19,6.33,9.06,8.63,4.10,3.24,1.42
36 33,0.55,0.59,0.86,1.96,2.63,6.61,9.29,9.29,5.16,4.16,2.17
37 34,0.46,0.65,0.96,2.01,2.97,6.63,9.26,9.59,6.22,5.00,3.04
38 35,0.40,0.71,1.18,2.00,3.17,6.45,9.06,9.42,7.34,5.70,3.69
39 36,0.36,0.67,1.35,1.95,3.40,6.04,8.64,8.78,8.15,6.11,4.70
40 37,0.35,0.61,1.48,1.77,3.56,5.39,7.91,7.82,8.61,6.21,5.74
41 38,0.35,0.57,1.62,1.48,3.57,4.42,6.98,6.50,8.37,6.29,5.41
42 39,0.38,0.50,1.49,1.30,3.58,3.29,6.02,4.95,6.76,6.24,3.82
43 40,0.51,0.48,1.29,1.21,3.44,2.37,5.13,3.47,4.37,4.91,2.34
44 41,0.74,0.54,1.24,1.10,3.15,1.74,4.55,2.15,2.30,2.61,1.35
45 42,0.92,0.73,1.17,0.99,3.38,1.36,4.52,1.38,1.06,1.09,0.65
46 43,0.96,0.93,1.04,0.91,3.99,1.17,4.71,1.11,0.58,0.63,0.29
47 44,0.91,0.82,0.95,0.78,4.17,0.99,4.72,0.90,0.47,0.59,0.16
48 section_length_cm,0.384,0.376,0.386,0.393,0.366,0.390,0.388,0.395,0.417,0.440,0.445
49 vtl_cm,16.90,16.55,16.98,17.30,16.11,17.14,17.09,17.40,18.33,19.34,19.59

56
src/AquaSynth.Core/FaustExport.VocalLadder.cs
View file

@ -4,15 +4,15 @@ namespace AquaSynth.Dsl;
// The vocal_network resonator: a Kelly-Lochbaum ladder per waveguide path, one cell per sample (cell length c / sample
// rate), written as a Faust `letrec` so every wave variable is a one-sample state. Cell areas come from the morphology
// rest shape deformed by the path's area control (the same Gaussian model the acoustic graph uses). A branch port whose
// rest shape deformed by the path's area control, or from its measured vowel shapes (DeformedAreaExpression owns both). A branch port whose
// two paths are both in the network becomes a three-port parallel junction (Pink Trombone's velum: ports are the oral
// cells either side of the join and the branch path's first cell, whose area is the branch opening).
public static partial class FaustEmitter
{
// Reflection at the closed glottal end of a path that carries a source.
private const float LadderGlottalReflection = 0.75f;
// Full-open velum diameter in cm, as in Pink Trombone's velum target.
private const float VelumOpenDiameterCentimeters = 0.4f;
// Speed of sound in cm/s used only to size the lip end correction from the lip radius (the ladder's own cell length uses the path's speed).
private const float LipEndCorrectionSoundSpeedCentimeters = 35000f;
private sealed record LadderPath(
WaveguidePath Path,
@ -37,9 +37,10 @@ public static partial class FaustEmitter
/// <summary>
/// Emits the area, reflection and junction definitions and the recursive ladder for <paramref name="members"/>, and returns,
/// per path, the identifier of its lip output (the right-going wave leaving its last cell).
/// per path, the identifier of its lip output (the right-going wave leaving its last cell) and, for each path that ends at a glottis,
/// the identifier of the pressure in its first cell (the tract wall radiation reads it).
/// </summary>
private static IReadOnlyDictionary<string, string> EmitVocalLadder(
private static (IReadOnlyDictionary<string, string> Lips, IReadOnlyDictionary<string, string> GlottalPressures) EmitVocalLadder(
StringBuilder source,
SynthPatch patch,
VocalNetwork network,
@ -68,12 +69,14 @@ public static partial class FaustEmitter
continue;
}
// The velopharyngeal port area, cm^2: opening * coupling * the design maximum (NasalTractAreas); it is the nasal ladder's first cell.
var velum = $"{name}_primitive_branch_{SafeIdentifier(branch.Name)}_velum_area";
source.AppendLine($" {velum} = max(0.000001, pow({admittances[branch.Name]} * {F(VelumOpenDiameterCentimeters)}, 2.0));");
var sum = $"({oral.A(cell)} + {oral.A(cell + 1)} + {velum})";
source.AppendLine($" {velum} = max(0.000001, {admittances[branch.Name]} * {F(NasalTractAreas.VelopharyngealPortMaximumAreaSquareCentimeters)});");
// The junction sits on the interface between oral cells cell - 1 and cell: R(cell) is entered from R(cell - 1), L(cell - 1) from L(cell).
var sum = $"({oral.A(cell - 1)} + {oral.A(cell)} + {velum})";
var key = $"{name}_primitive_branch_{SafeIdentifier(branch.Name)}";
source.AppendLine($" {key}_refl_left = (2.0 * {oral.A(cell)} - {sum}) / {sum};");
source.AppendLine($" {key}_refl_right = (2.0 * {oral.A(cell + 1)} - {sum}) / {sum};");
source.AppendLine($" {key}_refl_left = (2.0 * {oral.A(cell - 1)} - {sum}) / {sum};");
source.AppendLine($" {key}_refl_right = (2.0 * {oral.A(cell)} - {sum}) / {sum};");
source.AppendLine($" {key}_refl_nose = (2.0 * {velum} - {sum}) / {sum};");
junctions.Add(new LadderJunction(oral, nasal, cell, velum, $"{key}_refl_left", $"{key}_refl_right", $"{key}_refl_nose"));
}
@ -93,7 +96,7 @@ public static partial class FaustEmitter
continue;
}
var area = DeformedAreaExpression(parameters, member.Area.Shape, member.Area.Deformation, $"/vocal/areas/{member.AreaIndex}/area", (i + 0.5f) / member.Cells);
var area = DeformedAreaExpression(parameters, member.Area.Shape, member.Area.Deformation, $"/vocal/areas/{member.AreaIndex}/area", (i + 0.5f) / member.Cells, member.Area.Vowels);
source.AppendLine($" {member.A(i)} = max(0.000001, {area});");
}
@ -125,25 +128,40 @@ public static partial class FaustEmitter
: $" '{member.L(i)} = ({member.L(i + 1)} + {member.Refl(i + 1)} * ({member.R(i)} + {member.L(i + 1)})) * {loss};");
}
equations.AppendLine($" '{member.L(last)} = {member.R(last)} * ({lipReflection}) * {loss};");
// Open-end reflectance R(z) = -g (1-b)/(1-b z^-1): g is the radiation load's reflection (DC gain, negative), b follows the
// lip cell's area so that the low-frequency loss 1-|R|^2 ~ (ka)^2 and the DC group delay b/(1-b) is the end correction
// (Levine & Schwinger; Smith PASP; Valimaki & Karjalainen 1994). With a = sqrt(A/pi) in cm, q = (a fs / c)^2, c = 35000 cm/s,
// b/(1-b)^2 = q gives b = ((2q+1) - sqrt(4q+1)) / (2q). The pole is one-sample state in the same slot the constant used.
var q = $"{member.Prefix}_lip_q";
var b = $"{member.Prefix}_lip_b";
source.AppendLine($" {q} = max(0.000001, pow(sqrt({member.A(last)} / ma.PI) * ma.SR / {F(LipEndCorrectionSoundSpeedCentimeters)}, 2.0));");
source.AppendLine($" {b} = min(0.9, max(0.0, ((2.0 * {q} + 1.0) - sqrt(4.0 * {q} + 1.0)) / (2.0 * {q})));");
equations.AppendLine($" '{member.L(last)} = ({b} * {member.L(last)} + (1.0 - {b}) * ({lipReflection}) * {member.R(last)}) * {loss};");
}
var outputs = string.Join(", ", members.Select(member => member.R(member.Cells - 1)));
// Outputs: every path's lip wave, then the first-cell pressure R0 + L0 of each path that ends at a glottis (a branch path's first
// cell is the junction, not a glottis).
var glottal = members.Where(member => !targets.ContainsKey(member)).ToArray();
var outputs = string.Join(", ", members.Select(member => member.R(member.Cells - 1)).Concat(glottal.Select(member => $"({member.R(0)} + {member.L(0)})")));
var outputCount = members.Count + glottal.Length;
var ladder = $"{name}_primitive_ladder";
source.AppendLine($" {ladder} = ({outputs}) letrec {{");
source.Append(equations);
source.AppendLine(" };");
var lips = new Dictionary<string, string>(StringComparer.OrdinalIgnoreCase);
for (var m = 0; m < members.Count; m++)
var glottalPressures = new Dictionary<string, string>(StringComparer.OrdinalIgnoreCase);
for (var m = 0; m < outputCount; m++)
{
var lip = $"{members[m].Prefix}_lip";
source.AppendLine(members.Count == 1
? $" {lip} = {ladder};"
: $" {lip} = {ladder} : ba.selector({m}, {members.Count});");
lips[members[m].Path.Name] = lip;
var isLip = m < members.Count;
var owner = isLip ? members[m] : glottal[m - members.Count];
var tap = $"{owner.Prefix}_{(isLip ? "lip" : "glottal_pressure")}";
source.AppendLine(outputCount == 1
? $" {tap} = {ladder};"
: $" {tap} = {ladder} : ba.selector({m}, {outputCount});");
(isLip ? lips : glottalPressures)[owner.Path.Name] = tap;
}
return lips;
return (lips, glottalPressures);
}
}

84
src/AquaSynth.Core/FaustExport.cs
View file

@ -647,8 +647,8 @@ public static partial class FaustEmitter
}
}
var lips = members.Count == 0
? new Dictionary<string, string>()
var (lips, glottalPressures) = members.Count == 0
? (new Dictionary<string, string>(), new Dictionary<string, string>())
: EmitVocalLadder(source, patch, network, name, members, inputs, admittances, parameters, warnings);
foreach (var member in members)
@ -684,6 +684,17 @@ public static partial class FaustEmitter
radiated.Add(PrimitiveRadiationExpression(source, patch, name, load, wave, parameters, options, keepAlive));
}
}
if (glottalPressures.TryGetValue(path.Name, out var glottalPressure))
{
// Tract wall radiation: the pressure inside the tract, low-passed, leaves through the cheeks and neck whether or not the mouth
// is open. It is the only sound a closed mouth makes (the voice bar of a voiced stop) and sits far under the mouth's output
// in a vowel.
var wall = $"{member.Prefix}_wall_radiation";
source.AppendLine($" {wall} = {ProbeSignal(options, $"/debug/{name}/path/{path.Name}/wall_radiation", $"(({glottalPressure}) : fi.highpass(2, {WallRadiationHighpassHz}) : fi.lowpass(2, {WallRadiationLowpassHz})) * {WallRadiationGain}", -4, 4)};");
keepAlive.Add(wall);
radiated.Add(wall);
}
}
var output = radiated.Count == 0 ? "0.0" : string.Join(" + ", radiated);
@ -714,6 +725,8 @@ public static partial class FaustEmitter
var opening = parameters.Expression(OwnerField(owner, "opening"), port.Opening);
var impedance = parameters.Expression(OwnerField(owner, "impedance"), port.Impedance);
var flowScale = parameters.Expression(OwnerField(owner, "flow_scale"), port.FlowScale);
// f0 is the voice frequency times the source's ratio (the stress cue), smoothed over 10 ms so a step in the ratio is a glide.
var pitchRatio = $"({parameters.Expression(OwnerField(owner, "pitch_ratio"), port.PitchRatio)} : si.smooth(ba.tau2pole(0.01)))";
var load = $"{voiceName}_primitive_source_{safe}_load_pressure";
var voiced = $"{voiceName}_primitive_source_{safe}_voiced";
var flow = $"{voiceName}_primitive_source_{safe}_flow";
@ -721,7 +734,7 @@ public static partial class FaustEmitter
// Harmonic glottal source: a band-limited saw through a one-pole lowpass. Tilt 0 is the brightest source and 1 the
// steepest rolloff; the cutoff sweeps five octaves, 6400 Hz down to 200 Hz, and sits below Nyquist.
var cutoff = $"min(0.45 * ma.SR, 200.0 * pow(2.0, (1.0 - clip01({parameters.Expression(OwnerField(owner, "tilt"), port.Tilt)})) * 5.0))";
source.AppendLine($" {voiced} = (os.sawtooth({frequency}) : fi.lowpass(1, {cutoff})) * clip01({opening}) * clip01(0.25 + {tension});");
source.AppendLine($" {voiced} = (os.sawtooth(({frequency}) * {pitchRatio}) : fi.lowpass(1, {cutoff})) * clip01({opening}) * clip01(0.25 + {tension});");
source.AppendLine($" {flow} = {ProbeSignal(options, $"/debug/{voiceName}/source/{port.Name}/flow", $"({flowScale}) * (({voiced}) * ma.tanh((({pressure}) * (0.5 + {tension})) / max(0.05, {impedance})) + 0.08 * no.noise * {parameters.Expression(OwnerField(owner, "noise"), port.Noise)})", -2, 2)};");
keepAlive.Add(load);
keepAlive.Add(voiced);
@ -744,13 +757,19 @@ public static partial class FaustEmitter
var opening = $"clip01({parameters.Expression(OwnerField(owner, "opening"), contact.Opening)})";
var resistance = parameters.Expression(OwnerField(owner, "resistance"), contact.Resistance);
var stored = parameters.Expression(OwnerField(owner, "stored_pressure"), contact.StoredPressure);
var releaseGain = parameters.Expression(OwnerField(owner, "release_flow"), contact.ReleaseFlow);
var released = $"{voiceName}_primitive_contact_{safe}_released_flow";
var reservoir = $"{voiceName}_primitive_contact_{safe}_reservoir";
var burst = $"{voiceName}_primitive_contact_{safe}_burst";
source.AppendLine($" {reservoir} = {ProbeSignal(options, $"/debug/{voiceName}/contact/{contact.Name}/reservoir", $"(({stored}) + (1.0 - {opening}) * {resistance})", 0, 2)};");
source.AppendLine($" {released} = {ProbeSignal(options, $"/debug/{voiceName}/contact/{contact.Name}/released_flow", $"({reservoir}) * {opening} * (1.0 - min(0.95, {resistance} * 0.5))", -2, 2)};");
// Release burst: noise whose amplitude is the stored pressure through the opening area as it passes between shut and open
// (4 * o * (1 - o) peaks at half open and vanishes at both ends), so a shut or fully open contact is silent.
source.AppendLine($" {burst} = {ProbeSignal(options, $"/debug/{voiceName}/contact/{contact.Name}/burst", $"no.noise * ({releaseGain}) * ({stored}) * 4.0 * {opening} * (1.0 - {opening})", -2, 2)};");
keepAlive.Add(reservoir);
keepAlive.Add(burst);
keepAlive.Add(released);
return released;
return $"({released} + {burst})";
}
// Open fraction of a branch port: the ladder reads it as the velum area when the branch joins two paths of the network,
@ -792,6 +811,16 @@ public static partial class FaustEmitter
// /u/'s F1 and put a 450 Hz hump in its envelope. The aperture stays the gain consonant closures use.
private const string RadiationHighpassHz = "200.0";
// Tract wall radiation (the voice bar). Corners: skin and soft tissue pass only the low end of the voice (below about 200-300 Hz), and
// a second-order highpass at 90 Hz keeps the slow pressure drift of a shutting tract out (a 40 Hz first-order corner left it in). The gain
// is a DESIGN VALUE set so a word-initial voiced closure reads 30 dB under the following vowel (/p/ 41 dB); the wall term moves a
// vowel's f0 line by under 1 dB and its 200 and 300 Hz harmonics by up to 1.5 dB on one open vowel (no primary level retrieved;
// docs/research/articulatory-vocal-tract-sources.md). Both corners are second-order sections: the closed-tract bar is about 25 dB down
// above 300 Hz, not cut off there.
private const string WallRadiationHighpassHz = "90.0";
private const string WallRadiationLowpassHz = "300.0";
private const string WallRadiationGain = "0.7";
private static string PrimitiveRadiationExpression(
StringBuilder source,
SynthPatch patch,
@ -1671,16 +1700,20 @@ public static partial class FaustEmitter
(segment.StartPosition + segment.EndPosition) * 0.5f);
/// <summary>
/// The one owner of the Gaussian tongue/constriction/lip area deformation: cross-section area (cm^2) of a tract at a
/// normalized position, as a Faust expression over the control surfaces under <paramref name="areaPath"/>. Acoustic
/// graph terminals and segments and the vocal ladder's cells all read their area here.
/// The one owner of the tract area model: cross-section area (cm^2) of a tract at a normalized position, as a Faust
/// expression over the control surfaces under <paramref name="areaPath"/>. Acoustic graph terminals and segments and the
/// vocal ladder's cells all read their area here. The base diameter is either the rest shape deformed by the Gaussian
/// tongue and lip terms, or, when <paramref name="vowels"/> are given, the log-area blend of the measured vowel shapes
/// weighted by the vowel_* surfaces (no tongue term; lip_opening is then a cap on the lip cells). The constriction term
/// applies to either base.
/// </summary>
private static string DeformedAreaExpression(
ParameterMap parameters,
TractAreaFunction shape,
AcousticAreaControl? area,
string areaPath,
float position)
float position,
IReadOnlyList<VowelAreaShape>? vowels = null)
{
var restDiameter = shape.DiameterAt(position);
if (area is null)
@ -1690,17 +1723,29 @@ public static partial class FaustEmitter
var sections = Math.Max(2, shape.Sections);
var index = F(position * (sections - 1));
var tongueIndex = ScaledAreaIndex(parameters.Expression(OwnerField(areaPath, "tongue_index"), area.TongueIndex), area.IndexScale);
var tongueDiameter = parameters.Expression(OwnerField(areaPath, "tongue_diameter"), area.TongueDiameter);
var tongueWidth = $"max(0.001, {F(area.TongueWidth)} * {F(sections)})";
var constrictionIndex = ScaledAreaIndex(parameters.Expression(OwnerField(areaPath, "constriction_index"), area.ConstrictionIndex), area.IndexScale);
var constrictionDiameter = parameters.Expression(OwnerField(areaPath, "constriction_diameter"), area.ConstrictionDiameter);
var constrictionWidth = $"max(0.001, {F(area.ConstrictionWidth)} * {F(sections)})";
var lipOpening = parameters.Expression(OwnerField(areaPath, "lip_opening"), area.LipOpening);
var lipWeight = MathF.Exp(0 - MathF.Pow((1 - position) / Math.Max(0.001f, area.LipWidth), 2));
var tongueWeight = $"exp(0.0 - pow(({index} - ({tongueIndex})) / ({tongueWidth}), 2.0))";
var constrictionWeight = $"exp(0.0 - pow(({index} - ({constrictionIndex})) / ({constrictionWidth}), 2.0))";
var baseDiameter = $"({F(restDiameter)} + (({tongueDiameter}) - {F(restDiameter)}) * ({tongueWeight}) + (({lipOpening}) - {F(restDiameter)}) * {F(lipWeight)})";
string baseDiameter;
if (vowels is { Count: > 0 })
{
// The cap applies to the cells inside the lip Gaussian's half-power region; elsewhere it is unbounded.
var blended = VowelBlendDiameter(parameters, vowels, areaPath, position);
baseDiameter = lipWeight >= 0.5f ? $"min({blended}, ({lipOpening}))" : blended;
baseDiameter = $"({baseDiameter})";
}
else
{
var tongueIndex = ScaledAreaIndex(parameters.Expression(OwnerField(areaPath, "tongue_index"), area.TongueIndex), area.IndexScale);
var tongueDiameter = parameters.Expression(OwnerField(areaPath, "tongue_diameter"), area.TongueDiameter);
var tongueWidth = $"max(0.001, {F(area.TongueWidth)} * {F(sections)})";
var tongueWeight = $"exp(0.0 - pow(({index} - ({tongueIndex})) / ({tongueWidth}), 2.0))";
baseDiameter = $"({F(restDiameter)} + (({tongueDiameter}) - {F(restDiameter)}) * ({tongueWeight}) + (({lipOpening}) - {F(restDiameter)}) * {F(lipWeight)})";
}
var diameter = $"max(0.0, min({baseDiameter}, ({constrictionDiameter}) + max(0.0, ({baseDiameter}) - ({constrictionDiameter})) * (1.0 - ({constrictionWeight}))))";
return $"(({diameter}) * ({diameter}))";
@ -1708,6 +1753,19 @@ public static partial class FaustEmitter
scale == 1 ? expression : $"(({expression}) * {F(scale)})";
}
// Diameter of the blend exp(sum_k w_k ln A_k / max(1e-3, sum_k w_k)) of the vowel shapes at one position, A_k read at that position
// by linear interpolation in area. Weights below zero count as zero; when they sum to under 1e-3 the blend is the geometric mean of
// the shapes (select2 picks it), so zero weights are a defined shape rather than NaN.
private static string VowelBlendDiameter(ParameterMap parameters, IReadOnlyList<VowelAreaShape> vowels, string areaPath, float position)
{
var logs = vowels.Select(vowel => MathF.Log(Math.Max(0.000001f, StoryVowelAreas.AreaAt(vowel.Areas, position)))).ToArray();
var weights = vowels.Select(vowel => $"max(0.0, {parameters.Expression(OwnerField(areaPath, StoryVowelAreas.SurfaceField(vowel.Name)), 0)})").ToArray();
var total = $"({string.Join(" + ", weights)})";
var numerator = $"({string.Join(" + ", weights.Select((weight, k) => $"{weight} * ({F(logs[k])})"))})";
var geometricMean = F(MathF.Exp(logs.Average()));
return $"sqrt(select2({total} < 0.001, exp({numerator} / max(0.001, {total})), {geometricMean}))";
}
private static WaveClockPolicy ResolveWaveClock(
AcousticPortNetwork network,
IReadOnlyDictionary<string, WaveClockPolicy> waveClocks,

55
src/AquaSynth.Core/IpaGestureExperiment.cs
View file

@ -8,7 +8,8 @@ public sealed record IpaGestureExperimentTarget(
string Ipa,
string Descriptor,
float DurationSeconds = 0.18f,
string Notes = "");
string Notes = "",
float PitchRatio = 1f);
public sealed record IpaGestureExperimentVariant(
string Id,
@ -91,7 +92,8 @@ public static class IpaGestureExperiment
string roundId,
IEnumerable<IpaGestureExperimentTarget> targets,
IEnumerable<IpaGestureExperimentVariant> variants,
int timelineBlocks = 12)
int timelineBlocks = 12,
float sampleRate = 48000)
{
ArgumentException.ThrowIfNullOrWhiteSpace(rootDirectory);
ArgumentException.ThrowIfNullOrWhiteSpace(roundId);
@ -128,7 +130,7 @@ public static class IpaGestureExperiment
File.WriteAllText(scriptPath, script, Encoding.UTF8);
var patch = PatchScript.Parse(script);
var timeline = ProbeTimelineReport.Build(patch, NetworkName, Math.Max(1, timelineBlocks));
var timeline = ProbeTimelineReport.Build(patch, NetworkName, Math.Max(1, timelineBlocks), sampleRate);
var timelinePath = Path.Combine(timelineDirectory, $"{candidateId}.csv");
File.WriteAllText(timelinePath, ProbeTimelineReport.ToCsv(timeline), Encoding.UTF8);
@ -224,19 +226,29 @@ public static class IpaGestureExperiment
}
private static string GestureLine(IpaGestureExperimentTarget target, float start, float duration, float intensity) =>
$"phoneme_gesture name={SafeDslName(target.Id)} ipa={target.Ipa} {(target.Descriptor.Length == 0 ? "" : $"descriptor={target.Descriptor.Replace(' ', '_')} ")}start={Format(start)} dur={Format(duration)} intensity={Format(intensity)}";
$"phoneme_gesture name={SafeDslName(target.Id)} ipa={target.Ipa} {(target.Descriptor.Length == 0 ? "" : $"descriptor={target.Descriptor.Replace(' ', '_')} ")}start={Format(start)} dur={Format(duration)} intensity={Format(intensity)}{(target.PitchRatio == 1f ? "" : $" pitch_ratio={Format(target.PitchRatio)}")}";
// One cell per sample at 48 kHz: 17.52 cm * 48000 / 350 m/s = 24.03 cells, a ladder 17.5 cm long (the mean of Story's vocal tract lengths).
private const int TractSpeedMetersPerSecond = 350;
// Per-cell pass factor of the nasal path. The nasal cavity is the most heavily damped part of the tract (soft walls, turbinates, the paranasal
// sinuses), so a nasal murmur has broad formants and sits well under its vowels. A factor of .975 per 0.73 cm cell is a bandwidth of about
// 380 Hz, a DESIGN VALUE set from the murmur level: with the junction on its own interface and .997 (the oral tract's .998 less a little) the
// murmur read -4.2 dB (/m/) and -2.2 dB (/n/) re the flanking vowels; the target is 8 to 16 dB under (no primary level or bandwidth retrieved).
private const string NasalPathLoss = ".975";
private static string Scaffold(IReadOnlyList<string> gestureLines, float sustain, float frequency) => $$"""
patch gain=.45
morphology name=oral length_cm=17 diameters=.55,.75,1.05,1.45,1.25,.8 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=1 lip_opening=1.2
morphology name=nasal length_cm=12 diameters=.04,.28,.52,.72
waveguide_path name=oral_path morphology=oral strategy=thiran order=1 max_delay=4096 loss=.998
waveguide_path name=nasal_path morphology=nasal strategy=thiran order=1 max_delay=4096 loss=.997
morphology name=oral length_cm={{Format(StoryVowelAreas.MeanLengthCentimeters)}} diameters=.55,.75,1.05,1.45,1.25,.8 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=1 constriction_width=.03 lip_opening=4
vowel_set area_function=oral source={{StoryVowelAreas.Source}}
morphology name=nasal length_cm={{Format(NasalTractAreas.LengthCentimeters)}} areas={{NasalTractAreas.AreasList()}}
waveguide_path name=oral_path morphology=oral speed={{TractSpeedMetersPerSecond}} strategy=thiran order=1 max_delay=4096 loss=.998
waveguide_path name=nasal_path morphology=nasal speed={{TractSpeedMetersPerSecond}} strategy=thiran order=1 max_delay=4096 loss={{NasalPathLoss}}
source_port name=folds path=oral_path kind=glottal position=0 pressure=.66 tension=.54 opening=.42 noise=.04 impedance=.32 flow_scale=.12
branch_port name=velopharynx from=oral_path from_position=.45 to=nasal_path opening=.01 coupling=1
branch_port name=velopharynx from=oral_path from_position=.45 to=nasal_path opening=0 coupling=1
constriction_contact name=contact path=oral_path position=.92 opening=1 resistance=.45 stored_pressure=.12
radiation_load name=mouth path=oral_path kind=lip position=1 aperture=.8 reflection=-.82 impedance=.28
radiation_load name=nostril path=nasal_path kind=lip position=1 aperture=.8 reflection=-.82 impedance=.28
radiation_load name=mouth path=oral_path kind=lip position=1 aperture=.8 reflection=-1 impedance=.28
radiation_load name=nostril path=nasal_path kind=lip position=1 aperture=.8 reflection=-1 impedance=.28
probe_timeline name=flow network=voice blocks=12 block_size=64
vocal_network name=voice paths=oral_path,nasal_path sources=folds contacts=contact branches=velopharynx radiation=mouth,nostril probe=flow
{{string.Join("\n", gestureLines)}}
@ -286,7 +298,7 @@ public static class IpaGestureExperiment
IReadOnlyList<ProbeTimelineSample> timeline)
{
var tokens = DescriptorTokens(target.Descriptor);
var expected = ExpectedSurfaces(tokens);
var expected = ExpectedSurfaces(tokens, target.Descriptor);
var touched = patch.ControlSplines
.Where(spline => spline.Enabled)
.Select(spline => spline.SurfacePath)
@ -295,7 +307,7 @@ public static class IpaGestureExperiment
var coverage = expected.Length == 0
? 1
: expected.Count(surface => touched.Contains(surface)) / (float)expected.Length;
var direction = DirectionScore(tokens, patch.ControlSplines);
var direction = DirectionScore(tokens, target.Descriptor, patch.ControlSplines);
var contour = ContourTimingScore(patch.ControlSplines);
var primitive = PrimitiveTimelineScore(tokens, timeline);
var score = 0.25f * coverage + 0.25f * direction + 0.20f * contour + 0.20f * primitive;
@ -313,7 +325,7 @@ public static class IpaGestureExperiment
new(candidate.Id, target.Id, "gesture", name, Math.Clamp(value, 0, 1));
}
private static string[] ExpectedSurfaces(HashSet<string> tokens)
private static string[] ExpectedSurfaces(HashSet<string> tokens, string descriptor)
{
var surfaces = new List<string>
{
@ -326,11 +338,13 @@ public static class IpaGestureExperiment
{
surfaces.AddRange(
[
"/vocal/areas/0/area/tongue_index",
"/vocal/areas/0/area/tongue_diameter",
"/vocal/areas/0/area/lip_opening",
"/vocal/radiation/0/aperture"
]);
if (StoryVowelAreas.IdForDescriptor(descriptor) is { } vowelId)
{
surfaces.Add($"/vocal/areas/0/area/{StoryVowelAreas.SurfaceField(vowelId)}");
}
}
else
{
@ -367,7 +381,7 @@ public static class IpaGestureExperiment
return surfaces.Distinct(StringComparer.OrdinalIgnoreCase).ToArray();
}
private static float DirectionScore(HashSet<string> tokens, IReadOnlyList<ControlSpline> splines)
private static float DirectionScore(HashSet<string> tokens, string descriptor, IReadOnlyList<ControlSpline> splines)
{
var checks = new List<float>();
if (tokens.Contains("nasal"))
@ -390,8 +404,11 @@ public static class IpaGestureExperiment
if (tokens.Contains("vowel"))
{
checks.Add(Exists(splines, "/vocal/areas/0/area/tongue_index"));
checks.Add(Exists(splines, "/vocal/areas/0/area/tongue_diameter"));
if (StoryVowelAreas.IdForDescriptor(descriptor) is { } vowelId)
{
checks.Add(Exists(splines, $"/vocal/areas/0/area/{StoryVowelAreas.SurfaceField(vowelId)}"));
}
checks.Add(Exists(splines, "/vocal/radiation/0/aperture"));
}

9
src/AquaSynth.Core/Model.cs
View file

@ -443,11 +443,15 @@ public sealed record AcousticPortNetwork(
public IReadOnlyList<string> Connections { get; init; } = Connections ?? Array.Empty<string>();
}
/// <summary>A named measured tract shape: areas in cm^2 from the glottis to the lips, any number of sections.</summary>
public sealed record VowelAreaShape(string Name, IReadOnlyList<float> Areas);
public sealed record AreaFunction(
string Name,
TractAreaFunction Shape,
AcousticAreaControl? Deformation = null,
int EmitSections = 8);
int EmitSections = 8,
IReadOnlyList<VowelAreaShape>? Vowels = null);
public sealed record WaveguidePath(
string Name,
@ -469,7 +473,8 @@ public sealed record SourcePort(
float Noise = 0.08f,
float Impedance = 0.35f,
float FlowScale = 0.02f,
float Tilt = 0.5f);
float Tilt = 0.5f,
float PitchRatio = 1f);
public sealed record ConstrictionContact(
string Name,

33
src/AquaSynth.Core/NasalTractAreas.cs Normal file
View file

@ -0,0 +1,33 @@
namespace AquaSynth.Dsl;
/// <summary>
/// The owner of the nasal tract area function and of the velopharyngeal port's maximum area. Data: Birkholz P, 2005, "3D-Artikulatorische
/// Sprachsynthese", dissertation Universitaet Rostock, Table 2.3 (the nasal cavity of one subject, after Dang J, Honda K, 1994, "Morphological
/// and acoustical analysis of the nasal and the paranasal cavities", JASA 96(4):2088-2100): 19 sections of 0.6 cm (11.4 cm), areas in cm^2,
/// velopharyngeal end first, nostrils last. A test compares <see cref="Areas"/> to the cited values.
/// Not modelled: the paranasal sinus Helmholtz resonators (zeros at 552-1682 Hz in the same source) and the piriform fossae. The table is a
/// one-subject measurement and is unheard here: nothing in this file says how any nasal sounds.
/// </summary>
public static class NasalTractAreas
{
public const string Source = "birkholz2005-table2.3";
public const float SectionLengthCentimeters = 0.6f;
public const float LengthCentimeters = 11.4f;
public static IReadOnlyList<float> Areas { get; } =
[
1.63f, 2.07f, 2.72f, 3.59f, 4.24f, 3.26f, 3.04f, 3.04f, 2.72f, 2.5f,
2.39f, 2.39f, 1.85f, 0.76f, 1.41f, 1.74f, 1.30f, 1.74f, 0.76f,
];
/// <summary>
/// Maximum cross-section area of the velopharyngeal port in cm^2: the ladder's port area is opening * coupling * this. A DESIGN VALUE:
/// no source for the range of the port area was retrieved (Birkholz's velum port is a controllable area; its range is not in the
/// material read). 0.8 cm^2 is a wide-open port, below the first nasal section (1.63 cm^2) so the port, not the table, is the narrowest
/// first cell. Pink Trombone's earlier constant, a 0.4 cm diameter (0.16 cm^2 area at full opening and coupling), was likewise unsourced.
/// </summary>
public const float VelopharyngealPortMaximumAreaSquareCentimeters = 0.8f;
/// <summary>The scaffold's morphology line value for areas=, invariant-culture comma separated.</summary>
public static string AreasList() => string.Join(",", Areas.Select(area => area.ToString("0.##", System.Globalization.CultureInfo.InvariantCulture)));
}

171
src/AquaSynth.Core/PatchScript.cs
View file

@ -11,9 +11,16 @@ public static class PatchScript
/// <summary>
/// The DSL's own IPA table: the phones a phoneme_gesture knows without an explicit descriptor. A hand-written script
/// that names another phone falls back to an alveolar approximant; a caller that must not be approximated asks here first.
/// Vowels come from <see cref="StoryVowelAreas"/>, the one table that also says which measured tract shape each selects.
/// </summary>
public static bool TryIpaDescriptor(string ipa, out string descriptor)
{
if (StoryVowelAreas.TryVowelForIpa(ipa, out var vowelDescriptor, out _))
{
descriptor = vowelDescriptor;
return true;
}
var known = ipa switch
{
"p" => "voiceless_bilabial_plosive",
@ -26,11 +33,6 @@ public static class PatchScript
"ʒ" => "voiced_postalveolar_fricative",
"ʍ" => "voiceless_labial_velar_fricative",
"w" => "voiced_labial_velar_approximant",
"i" => "voiced_close_front_unrounded_vowel",
"u" => "voiced_close_back_rounded_vowel",
"a" or "ɑ" => "voiced_open_back_unrounded_vowel",
"æ" => "voiced_open_front_unrounded_vowel",
"ə" => "voiced_mid_central_unrounded_vowel",
_ => null
};
descriptor = known ?? "";
@ -245,6 +247,10 @@ public static class PatchScript
FlushPendingOperatorGraph();
AddAreaFunction(fields, line);
break;
case "vowel_set":
FlushPendingOperatorGraph();
AddVowelSet(fields, line);
break;
case "waveguide_path":
FlushPendingOperatorGraph();
AddWaveguidePath(fields, line);
@ -741,6 +747,43 @@ public static class PatchScript
RegisterAreaFunctionSurfaces(_areaFunctions.Count - 1, area);
}
// vowel_set area_function=<morphology name> source=story2008: attaches the eleven measured vowel shapes to a morphology and
// registers a weight surface /vocal/areas/<i>/area/vowel_<id> in [0,1] for each (default 0; all zero is the geometric mean shape).
// The morphology's tongue surfaces stay registered but no longer shape the tract. Must precede the phoneme_gestures that drive it.
private void AddVowelSet(IReadOnlyDictionary<string, string> fields, int line)
{
var name = RequiredAny(fields, ["area_function", "morphology", "area", "shape"], line);
var source = Required(fields, "source", line);
if (!source.Equals(StoryVowelAreas.Source, StringComparison.OrdinalIgnoreCase))
{
throw new PatchScriptException(line, $"unknown vowel source `{source}` (only `{StoryVowelAreas.Source}`)");
}
if (!_areaFunctionsByName.TryGetValue(name, out var area))
{
throw new PatchScriptException(line, $"unknown area function `{name}`");
}
if (area.Deformation is null)
{
throw new PatchScriptException(line, $"area function `{name}` has no area control (constriction and lip surfaces) for vowels to work with");
}
if (area.Vowels is not null)
{
throw new PatchScriptException(line, $"area function `{name}` already has vowel shapes");
}
var index = _areaFunctions.IndexOf(area);
var withVowels = area with { Vowels = StoryVowelAreas.All.Select(vowel => new VowelAreaShape(vowel.Id, vowel.Areas)).ToArray() };
_areaFunctions[index] = withVowels;
_areaFunctionsByName[name] = withVowels;
foreach (var vowel in StoryVowelAreas.All)
{
RegisterNormalizedPrimitiveSurface(OwnerField($"/vocal/areas/{index}/area", StoryVowelAreas.SurfaceField(vowel.Id)), name, StoryVowelAreas.SurfaceField(vowel.Id), 0);
}
}
private void AddWaveguidePath(IReadOnlyDictionary<string, string> fields, int line)
{
var name = Required(fields, "name", line);
@ -895,7 +938,8 @@ public static class PatchScript
GetBoundFloat(fields, line, 0.08f, $"/vocal/sources/{_sourcePorts.Count}/noise", "noise", "aspiration"),
GetBoundFloat(fields, line, 0.35f, $"/vocal/sources/{_sourcePorts.Count}/impedance", "impedance", "load", "source_impedance"),
GetBoundFloat(fields, line, 0.02f, $"/vocal/sources/{_sourcePorts.Count}/flow_scale", "flow_scale", "scale"),
GetBoundFloat(fields, line, 0.5f, $"/vocal/sources/{_sourcePorts.Count}/tilt", "tilt", "spectral_tilt"));
GetBoundFloat(fields, line, 0.5f, $"/vocal/sources/{_sourcePorts.Count}/tilt", "tilt", "spectral_tilt"),
GetBoundFloat(fields, line, 1f, $"/vocal/sources/{_sourcePorts.Count}/pitch_ratio", "pitch_ratio", "f0_ratio"));
_sourcePorts.Add(port);
_sourcePortsByName[name] = port;
RegisterSourcePortSurfaces(_sourcePorts.Count - 1, port);
@ -1342,6 +1386,7 @@ public static class PatchScript
RegisterNormalizedPrimitiveSurface(OwnerField(owner, "impedance"), port.Name, "impedance", port.Impedance, 0.05f, 2);
RegisterNormalizedPrimitiveSurface(OwnerField(owner, "flow_scale"), port.Name, "flow_scale", port.FlowScale, 0, 0.2f);
RegisterNormalizedPrimitiveSurface(OwnerField(owner, "tilt"), port.Name, "tilt", port.Tilt);
RegisterNormalizedPrimitiveSurface(OwnerField(owner, "pitch_ratio"), port.Name, "pitch_ratio", port.PitchRatio, 0.5f, 2f);
}
private void RegisterAreaFunctionSurfaces(int index, AreaFunction area)
@ -1459,7 +1504,7 @@ public static class PatchScript
}
}
private void EmitPhonemeGestureSplines(string name, string descriptor, float start, float duration, float intensity)
private void EmitPhonemeGestureSplines(string name, string descriptor, float start, float duration, float intensity, float pitchRatio = 1f)
{
var tokens = DescriptorTokens(descriptor);
var voiced = tokens.Contains("voiced") || tokens.Contains("vowel") || tokens.Contains("nasal") || tokens.Contains("approximant");
@ -1474,16 +1519,36 @@ public static class PatchScript
AddAuthoredSpline($"{name}_source_pressure", "/vocal/sources/0/pressure", start, duration, (0, 0.4f), (0.2f, pressure), (1, pressure * 0.85f));
AddAuthoredSpline($"{name}_source_tension", "/vocal/sources/0/tension", start, duration, (0, voiced ? 0.55f : 0.18f), (1, voiced ? 0.62f : 0.20f));
AddAuthoredSpline($"{name}_source_opening", "/vocal/sources/0/opening", start, duration, (0, voiced ? 0.45f : 0.25f), (1, voiced ? 0.48f : 0.28f));
// Intensity is an amplitude factor on the source's flow scale (the level lever: pressure saturates in the tanh), and the f0 ratio
// is held per gesture; both default to the source's own value.
HoldFractionOfDefault($"{name}_source_level", "/vocal/sources/0/flow_scale", start, duration, intensity);
HoldOnSurface($"{name}_source_pitch", "/vocal/sources/0/pitch_ratio", start, duration, pitchRatio);
AddAuthoredSpline($"{name}_source_noise", "/vocal/sources/0/noise", start, duration, (0, tokens.Contains("fricative") ? 0.45f : 0.05f), (1, tokens.Contains("fricative") ? 0.55f : 0.05f));
if (tokens.Contains("vowel"))
{
EmitVowelGesture(name, tokens, start, duration);
EmitVowelGesture(name, descriptor, start, duration);
return;
}
AddAuthoredSpline($"{name}_constriction_index", "/vocal/areas/0/area/constriction_index", start, duration, (0, place), (1, place));
AddAuthoredSpline($"{name}_constriction_diameter", "/vocal/areas/0/area/constriction_diameter", start, duration, (0, constriction), (0.5f, constriction), (1, Math.Min(0.65f, constriction + 0.2f)));
if (_areaFunctions.Count > 0 && _areaFunctions[0].Vowels is not null && _areaFunctions[0].Deformation is { } control)
{
// Under vowel_set the constriction is placed by distance along the tract: the surface is a section index, so a normalized
// place fraction would land at fraction * sections / (sections - 1) of the tract, past the lips for a labial.
var sections = Math.Max(2, _areaFunctions[0].Shape.Sections);
var position = ConsonantPositionFromGlottis(tokens, place, _areaFunctions[0].Shape.LengthCentimeters);
HoldOnSurface($"{name}_constriction_index", "/vocal/areas/0/area/constriction_index", start, duration, position * (sections - 1) / control.IndexScale);
}
else
{
AddAuthoredSpline($"{name}_constriction_index", "/vocal/areas/0/area/constriction_index", start, duration, (0, place), (1, place));
}
// A nasal stop holds its oral closure for as long as the velum and the contact stay put (to 0.8 of the gesture): the murmur is the
// closed-mouth sound, and an alveolar closure that began to open at 0.5 let the oral tract radiate again over the last third (/n/ read
// -3.8 dB under its vowels against /m/ -13.2 dB, whose lips hold to 0.7).
var closureHold = tokens.Contains("nasal") ? 0.8f : 0.5f;
AddAuthoredSpline($"{name}_constriction_diameter", "/vocal/areas/0/area/constriction_diameter", start, duration, (0, constriction), (closureHold, constriction), (1, Math.Min(0.65f, constriction + 0.2f)));
if (tokens.Contains("bilabial") || tokens.Contains("labiodental") || tokens.Contains("labial") || tokens.Contains("rounded"))
{
@ -1501,13 +1566,16 @@ public static class PatchScript
{
AddAuthoredSpline($"{name}_contact", "/vocal/contacts/0/opening", start, duration, (0, 0.02f), (0.65f, 0.02f), (1, 0.9f));
AddAuthoredSpline($"{name}_stored_pressure", "/vocal/contacts/0/stored_pressure", start, duration, (0, 0.4f), (0.65f, 0.75f), (1, 0.05f));
AddAuthoredSpline($"{name}_velum", "/vocal/branches/0/opening", start, duration, (0, 0.01f), (1, 0.01f));
// The release burst is turbulence at the opening constriction: its gain is held for the plosive and the exporter shapes
// it by the stored pressure and the opening (see PrimitiveContactExpression).
HoldOnSurface($"{name}_release_flow", "/vocal/contacts/0/release_flow", start, duration, PlosiveBurstGain);
AddAuthoredSpline($"{name}_velum", "/vocal/branches/0/opening", start, duration, (0, VelumClosed), (1, VelumClosed));
}
else if (tokens.Contains("fricative"))
{
AddAuthoredSpline($"{name}_contact", "/vocal/contacts/0/opening", start, duration, (0, 0.18f), (1, 0.22f));
AddAuthoredSpline($"{name}_contact_resistance", "/vocal/contacts/0/resistance", start, duration, (0, 0.75f), (1, 0.85f));
AddAuthoredSpline($"{name}_velum", "/vocal/branches/0/opening", start, duration, (0, 0.01f), (1, 0.01f));
AddAuthoredSpline($"{name}_velum", "/vocal/branches/0/opening", start, duration, (0, VelumClosed), (1, VelumClosed));
}
else
{
@ -1515,40 +1583,36 @@ public static class PatchScript
}
}
// Vowel tract shapes for the Kelly-Lochbaum ladder on the experiment scaffold (six rest diameters, 17 cm, 48 kHz): tongue
// section index (0..5), tongue and lip diameters in cm. Each was found by grid search over the exported ladder's impulse
// response (peak picking) toward textbook F1/F2, not by hand: measured F1/F2 in Hz in the trailing comment.
// Features are (front 0..1, openness 0..1, rounded 0/1); a descriptor selects the nearest anchor, there is no interpolation
// because the shapes are not monotone in formant space.
private static readonly (float Front, float Open, float Round, float TongueIndex, float TongueCm, float LipCm)[] VowelShapes =
[
(1.0f, 0.0f, 0, 3.00f, 0.15f, 2.20f), // i 271/2313
(0.0f, 0.0f, 1, 2.25f, 0.45f, 0.05f), // u 292/876
(0.0f, 1.0f, 0, 1.50f, 0.35f, 1.00f), // a 703/1084
(1.0f, 1.0f, 0, 0.50f, 0.30f, 1.00f), // ae 665/1755
(1.0f, 0.5f, 0, 2.50f, 1.80f, 2.20f), // e 535/1799
(0.0f, 0.5f, 1, 1.50f, 0.10f, 1.60f), // o 454/962
(0.5f, 0.5f, 0, 2.00f, 1.20f, 0.60f), // schwa 497/1503
];
private void EmitVowelGesture(string name, HashSet<string> tokens, float start, float duration)
// A vowel gesture raises the weight of the descriptor's Story vowel to 1 and holds every other vowel weight at 0 (the schwa and
// any vowel descriptor without a Story shape hold all of them at 0: the geometric-mean shape). The constriction and the lip cap
// stay open, the velum closed. On a morphology without vowel_set the vowel surfaces do not exist and nothing is held.
private void EmitVowelGesture(string name, string descriptor, float start, float duration)
{
var front = tokens.Contains("front") ? 1f : tokens.Contains("back") ? 0f : 0.5f;
var open = tokens.Contains("close") ? 0f : tokens.Contains("mid") ? 0.5f : tokens.Contains("open") ? 1f : 0.5f;
var round = tokens.Contains("rounded") && !tokens.Contains("unrounded") ? 1f : 0f;
var shape = VowelShapes
.OrderBy(candidate => MathF.Abs(candidate.Front - front) + MathF.Abs(candidate.Open - open) + MathF.Abs(candidate.Round - round))
.First();
HoldOnSurface($"{name}_tongue_index", "/vocal/areas/0/area/tongue_index", start, duration, shape.TongueIndex);
HoldOnSurface($"{name}_tongue_diameter", "/vocal/areas/0/area/tongue_diameter", start, duration, shape.TongueCm);
HoldOnSurface($"{name}_lip_area", "/vocal/areas/0/area/lip_opening", start, duration, shape.LipCm);
// A vowel leaves the tract unconstricted: hold the constriction at the surface maximum.
var target = StoryVowelAreas.IdForDescriptor(descriptor);
foreach (var vowel in StoryVowelAreas.All)
{
var weight = vowel.Id == target ? 1f : 0f;
HoldOnSurface($"{name}_{StoryVowelAreas.SurfaceField(vowel.Id)}", $"/vocal/areas/0/area/{StoryVowelAreas.SurfaceField(vowel.Id)}", start, duration, weight);
}
// A vowel leaves the tract unconstricted: hold the constriction and the lip cap at the surface maximum.
AddAuthoredSpline($"{name}_constriction_diameter", "/vocal/areas/0/area/constriction_diameter", start, duration, (0, 1), (1, 1));
AddAuthoredSpline($"{name}_lip_area", "/vocal/areas/0/area/lip_opening", start, duration, (0, 1), (1, 1));
// The ladder's lip cell already carries the lip area into the resonance and the level. Aperture is the radiating
// gain the consonant closures use (the lip cell cannot shut), so a vowel keeps it fully open; scaling it by the lip
// area as well counted the rounding twice and left /u/ 30 dB under /i/.
AddAuthoredSpline($"{name}_lip_radiation", "/vocal/radiation/0/aperture", start, duration, (0, 1), (1, 1));
AddAuthoredSpline($"{name}_velum", "/vocal/branches/0/opening", start, duration, (0, 0.01f), (1, 0.01f));
AddAuthoredSpline($"{name}_velum", "/vocal/branches/0/opening", start, duration, (0, VelumClosed), (1, VelumClosed));
// The contact is fully open for a vowel; a contact left half open by the previous stop would keep making release noise.
AddAuthoredSpline($"{name}_contact", "/vocal/contacts/0/opening", start, duration, (0, 1), (1, 1));
}
// Holds a surface at its default value times a factor (clamped to the surface range).
private void HoldFractionOfDefault(string name, string surfacePath, float start, float duration, float factor)
{
if (!_controlSurfacesByPath.TryGetValue(surfacePath, out var surface)) return;
var normalized = Math.Clamp(surface.DefaultNormalized * factor, 0, 1);
AddAuthoredSpline(name, surfacePath, start, duration, (0, normalized), (1, normalized));
}
// Holds a surface at a value in its own unit (cm, section index); splines carry the surface's normalized range.
@ -1589,6 +1653,29 @@ public static class PatchScript
.Split(['_', ',', '+', ' '], StringSplitOptions.RemoveEmptyEntries | StringSplitOptions.TrimEntries)
.ToHashSet(StringComparer.OrdinalIgnoreCase);
// A shut velum is opening 0: the port area is linear in the opening (opening * coupling * NasalTractAreas maximum), so 0.01 would
// leave a 0.008 cm^2 leak that shows in the oral formants.
private const float VelumClosed = 0f;
// Gain of the plosive release burst (contact release_flow, surface range 0..2). A DESIGN VALUE set by measurement against the
// following vowel; no published burst level was retrieved.
private const float PlosiveBurstGain = 0.06f;
// Distance from the lip plane to the alveolar ridge in cm. A DESIGN VALUE; no source for the distance was retrieved. 2.555 cm puts the
// closure on the centre of ladder cell 20 of 24 of the scaffold tract (17.52 cm at 48 kHz and 350 m/s: (20.5 / 24) * 17.52 = 14.965 cm
// from the glottis), the cell grid being what a narrow constriction (scaffold constriction_width .03, e-fold 0.63 cm) can close; it is
// the ladder's counterpart of the tube check's closed end after section 37 of 44 (14.73 cm).
private const float AlveolarRidgeCentimetersFromLips = 2.555f;
// Normalized position (0 glottis, 1 lips) of a consonant's constriction in the vowel_set path: labials at the lip plane, alveolars
// AlveolarRidgeCentimetersFromLips behind it, every other place at its PlaceTarget fraction of the tract length.
private static float ConsonantPositionFromGlottis(HashSet<string> tokens, float placeFraction, float lengthCentimeters)
{
if (tokens.Contains("bilabial") || tokens.Contains("labiodental") || tokens.Contains("labial")) return 1;
if (tokens.Contains("alveolar")) return 1 - AlveolarRidgeCentimetersFromLips / lengthCentimeters;
return placeFraction;
}
private static float PlaceTarget(HashSet<string> tokens)
{
if (tokens.Contains("bilabial") || tokens.Contains("labiodental") || tokens.Contains("labial")) return 0.96f;
@ -2595,13 +2682,14 @@ public static class PatchScript
var start = GetFloat(fields, line, 0, "start", "time", "at");
var duration = Math.Max(0.001f, GetFloat(fields, line, 0.12f, "dur", "duration", "length"));
var intensity = Math.Clamp(GetFloat(fields, line, 1, "intensity", "gain", "strength"), 0, 2);
var pitchRatio = Math.Clamp(GetFloat(fields, line, 1, "pitch_ratio", "f0_ratio"), 0.5f, 2f);
var descriptor = GestureDescriptor(ipa, fields);
var name = TryGetAny(fields, ["name", "n"], out var nameText)
? nameText
: $"ipa_{SafeName(ipa)}_{_phonemeGestures.Count}";
_phonemeGestures.Add(new PhonemeGesture(name, ipa, descriptor, start, duration, intensity));
EmitPhonemeGestureSplines(name, descriptor, start, duration, intensity);
EmitPhonemeGestureSplines(name, descriptor, start, duration, intensity, pitchRatio);
}
private void AddControlCurve(IReadOnlyDictionary<string, string> fields, int line)
@ -4014,6 +4102,7 @@ public static class PatchScript
"nasal_branch" or "nose_branch" or "nasal" => "nasal_branch",
"tract_motion" or "shape_motion" or "slew" => "tract_motion",
"path" or "acoustic_path" or "morphology_path" => "acoustic_path",
"vowel_set" or "vowelset" => "vowel_set",
"waveguide_path" or "waveguidepath" or "wg_path" => "waveguide_path",
"source_port" or "acoustic_source" or "port_source" => "source_port",
"constriction_contact" or "contact_port" or "contact" => "constriction_contact",

9
src/AquaSynth.Core/ProbeTimelineReport.cs
View file

@ -7,7 +7,7 @@ public sealed record ProbeTimelineSample(int Block, string Primitive, string Sig
public static class ProbeTimelineReport
{
public static IReadOnlyList<ProbeTimelineSample> Build(SynthPatch patch, string networkName, int blocks = 4)
public static IReadOnlyList<ProbeTimelineSample> Build(SynthPatch patch, string networkName, int blocks = 4, float sampleRate = 48000)
{
var network = patch.VocalNetworks.First(network => network.Name.Equals(networkName, StringComparison.OrdinalIgnoreCase));
var areas = patch.AreaFunctions.ToDictionary(area => area.Name, StringComparer.OrdinalIgnoreCase);
@ -27,7 +27,7 @@ public static class ProbeTimelineReport
for (var block = 0; block < Math.Max(1, blocks); block++)
{
var timeSeconds = block * blockSize / 48000f;
var timeSeconds = block * blockSize / sampleRate;
var outgoingByPath = new Dictionary<string, float>(StringComparer.OrdinalIgnoreCase);
foreach (var pathName in network.Paths)
{
@ -37,7 +37,7 @@ public static class ProbeTimelineReport
var pathIndex = PathIndex(patch.WaveguidePaths, path);
var speed = SurfaceValue(surfacesByField, splinesBySurface, $"/vocal/paths/{pathIndex}/speed", path.PropagationSpeedMetersPerSecond, timeSeconds);
var loss = SurfaceValue(surfacesByField, splinesBySurface, $"/vocal/paths/{pathIndex}/loss", path.Loss, timeSeconds);
var delay = area.Shape.LengthMeters / Math.Max(1, speed) * 48000;
var delay = area.Shape.LengthMeters / Math.Max(1, speed) * sampleRate;
var incoming = PathIncomingFlow(patch, network, sources, contacts, branches, path.Name, surfacesByField, splinesBySurface, timeSeconds);
var outgoing = incoming * Math.Clamp(loss, 0, 1);
outgoingByPath[path.Name] = outgoing;
@ -80,6 +80,9 @@ public static class ProbeTimelineReport
Add(samples, block, $"contact:{contact.Name}", "resistance", resistance);
Add(samples, block, $"contact:{contact.Name}", "reservoir", reservoir);
Add(samples, block, $"contact:{contact.Name}", "released_flow", reservoir * opening * (1 - Math.Min(0.95f, resistance * 0.5f)));
var releaseGain = SurfaceValue(surfacesByField, splinesBySurface, $"{owner}/release_flow", contact.ReleaseFlow, timeSeconds);
// The exporter's burst is white noise of this amplitude (noise itself is not mirrored).
Add(samples, block, $"contact:{contact.Name}", "burst_amplitude", releaseGain * stored * 4 * opening * (1 - opening));
}
foreach (var branchName in network.Branches)

93
src/AquaSynth.Core/StoryVowelAreas.cs Normal file
View file

@ -0,0 +1,93 @@
namespace AquaSynth.Dsl;
/// <summary>
/// The owner of the measured vowel tract shapes and of the IPA vowel -> Story vowel mapping. Data: Story BH, 2008, "Comparison of
/// magnetic resonance imaging-based vocal tract area functions obtained from the same speaker in 1994 and 2002", JASA 123(1):327-335,
/// doi 10.1121/1.2805683, PMC2377017, Table I: eleven vowels, 44 sections glottis -> lips, areas in cm^2. The cited copy is
/// docs/research/data/story2008-area-functions.csv; a test compares every cell of this table to that file (transcription guard).
/// The table is a one-speaker measurement and is unheard here: nothing in this file says how any vowel sounds.
///
/// Story's vocal tract lengths run 16.1-19.6 cm (the CSV carries each). The ladder has one length per patch, so a patch that uses
/// these shapes sets its oral morphology to <see cref="MeanLengthCentimeters"/> and every vowel is read at that length; the sections are
/// resampled onto the ladder's cells by <see cref="AreaAt"/>.
///
/// Control surface ids (ascii, used in /vocal/areas/&lt;i&gt;/area/vowel_&lt;id&gt;): i ih e eh ae uh ah aw o oo_short u for
/// i ɪ e ɛ æ ʌ ɑ ɔ o ʊ u. The schwa has no Story shape: it is the all-weights-zero state, the geometric mean of the eleven shapes.
/// </summary>
public static class StoryVowelAreas
{
public const string Source = "story2008";
public const float MeanLengthCentimeters = 17.52f;
public const int Sections = 44;
public sealed record Vowel(string Ipa, string Id, string Descriptor, IReadOnlyList<float> Areas);
public static string SurfaceField(string id) => $"vowel_{id}";
/// <summary>
/// Vowels other than the eleven shapes that the DSL's IPA table accepts, and the shape each selects. ɑ and a share a shape; the schwa
/// maps to null (no weight raised: the geometric-mean shape).
/// </summary>
public static bool TryVowelForIpa(string ipa, out string descriptor, out string? id)
{
if (ipa == "a") ipa = "ɑ";
if (ipa == "ə")
{
descriptor = SchwaDescriptor;
id = null;
return true;
}
var vowel = All.FirstOrDefault(candidate => candidate.Ipa == ipa);
descriptor = vowel?.Descriptor ?? "";
id = vowel?.Id;
return vowel is not null;
}
public const string SchwaDescriptor = "voiced_mid_central_unrounded_vowel";
/// <summary>
/// The explicit descriptor -> Story vowel id lookup (there is no feature-distance guessing). A vowel descriptor with no entry,
/// and the schwa's, return null: the geometric-mean shape.
/// </summary>
public static string? IdForDescriptor(string descriptor) =>
All.FirstOrDefault(vowel => vowel.Descriptor.Equals(descriptor.Replace('-', '_').ToLowerInvariant(), StringComparison.Ordinal))?.Id;
/// <summary>
/// Area in cm^2 of <paramref name="areas"/> at a normalized glottis -> lips position: linear interpolation between section centres
/// (section j is centred at (j + 0.5) / n), held at the first and last section beyond them.
/// </summary>
public static float AreaAt(IReadOnlyList<float> areas, float position)
{
var x = Math.Clamp(position * areas.Count - 0.5f, 0, areas.Count - 1);
var left = (int)MathF.Floor(x);
var right = Math.Min(areas.Count - 1, left + 1);
return areas[left] + (areas[right] - areas[left]) * (x - left);
}
public static readonly IReadOnlyList<Vowel> All =
[
new("i", "i", "voiced_close_front_unrounded_vowel", [
0.51f, 0.59f, 0.62f, 0.72f, 1.24f, 2.30f, 3.30f, 3.59f, 3.22f, 2.86f, 3.00f, 3.61f, 4.39f, 4.95f, 5.17f, 5.16f, 5.18f, 5.26f, 5.20f, 5.02f, 4.71f, 4.13f, 3.43f, 2.83f, 2.32f, 1.83f, 1.46f, 1.23f, 1.08f, 0.94f, 0.80f, 0.67f, 0.55f, 0.46f, 0.40f, 0.36f, 0.35f, 0.35f, 0.38f, 0.51f, 0.74f, 0.92f, 0.96f, 0.91f]),
new("ɪ", "ih", "voiced_near_close_near_front_unrounded_vowel", [
0.28f, 0.21f, 0.21f, 0.30f, 0.47f, 0.71f, 1.12f, 1.48f, 1.35f, 1.05f, 0.92f, 0.92f, 1.19f, 1.94f, 2.83f, 3.31f, 3.48f, 3.60f, 3.64f, 3.49f, 3.20f, 2.90f, 2.59f, 2.21f, 1.87f, 1.54f, 1.20f, 0.92f, 0.74f, 0.59f, 0.52f, 0.54f, 0.59f, 0.65f, 0.71f, 0.67f, 0.61f, 0.57f, 0.50f, 0.48f, 0.54f, 0.73f, 0.93f, 0.82f]),
new("e", "e", "voiced_close_mid_front_unrounded_vowel", [
0.29f, 0.26f, 0.30f, 0.40f, 0.55f, 0.74f, 0.99f, 1.09f, 0.90f, 0.69f, 0.77f, 1.31f, 2.13f, 2.74f, 3.03f, 3.23f, 3.33f, 3.27f, 3.09f, 2.84f, 2.66f, 2.46f, 2.14f, 1.79f, 1.44f, 1.17f, 1.00f, 0.88f, 0.80f, 0.81f, 0.85f, 0.84f, 0.86f, 0.96f, 1.18f, 1.35f, 1.48f, 1.62f, 1.49f, 1.29f, 1.24f, 1.17f, 1.04f, 0.95f]),
new("ɛ", "eh", "voiced_open_mid_front_unrounded_vowel", [
0.37f, 0.30f, 0.24f, 0.23f, 0.29f, 0.41f, 0.58f, 0.82f, 0.97f, 0.82f, 0.62f, 0.54f, 0.54f, 0.68f, 1.09f, 1.62f, 2.03f, 2.35f, 2.51f, 2.39f, 2.22f, 2.13f, 2.00f, 1.78f, 1.58f, 1.43f, 1.31f, 1.23f, 1.24f, 1.38f, 1.61f, 1.82f, 1.96f, 2.01f, 2.00f, 1.95f, 1.77f, 1.48f, 1.30f, 1.21f, 1.10f, 0.99f, 0.91f, 0.78f]),
new("æ", "ae", "voiced_near_open_front_unrounded_vowel", [
0.31f, 0.21f, 0.18f, 0.23f, 0.33f, 0.50f, 0.78f, 0.96f, 0.85f, 0.63f, 0.46f, 0.36f, 0.33f, 0.46f, 0.73f, 1.00f, 1.30f, 1.66f, 1.97f, 2.06f, 2.03f, 2.01f, 1.89f, 1.66f, 1.49f, 1.42f, 1.37f, 1.34f, 1.41f, 1.58f, 1.82f, 2.19f, 2.63f, 2.97f, 3.17f, 3.40f, 3.56f, 3.57f, 3.58f, 3.44f, 3.15f, 3.38f, 3.99f, 4.17f]),
new("ʌ", "uh", "voiced_open_mid_back_unrounded_vowel", [
0.23f, 0.34f, 0.47f, 0.60f, 0.77f, 1.06f, 1.26f, 1.09f, 0.80f, 0.65f, 0.56f, 0.47f, 0.37f, 0.24f, 0.17f, 0.18f, 0.23f, 0.27f, 0.28f, 0.29f, 0.33f, 0.52f, 0.97f, 1.50f, 1.91f, 2.23f, 2.65f, 3.29f, 4.13f, 5.00f, 5.77f, 6.33f, 6.61f, 6.63f, 6.45f, 6.04f, 5.39f, 4.42f, 3.29f, 2.37f, 1.74f, 1.36f, 1.17f, 0.99f]),
new("ɑ", "ah", "voiced_open_back_unrounded_vowel", [
0.56f, 0.62f, 0.66f, 0.78f, 0.97f, 1.16f, 1.12f, 0.82f, 0.55f, 0.45f, 0.37f, 0.29f, 0.21f, 0.15f, 0.16f, 0.25f, 0.34f, 0.43f, 0.54f, 0.61f, 0.67f, 0.98f, 1.76f, 2.75f, 3.52f, 4.08f, 4.74f, 5.61f, 6.60f, 7.61f, 8.48f, 9.06f, 9.29f, 9.26f, 9.06f, 8.64f, 7.91f, 6.98f, 6.02f, 5.13f, 4.55f, 4.52f, 4.71f, 4.72f]),
new("ɔ", "aw", "voiced_open_mid_back_rounded_vowel", [
0.27f, 0.43f, 0.54f, 0.67f, 0.83f, 0.92f, 0.89f, 0.73f, 0.55f, 0.44f, 0.37f, 0.28f, 0.18f, 0.14f, 0.15f, 0.14f, 0.13f, 0.14f, 0.18f, 0.20f, 0.23f, 0.49f, 1.03f, 1.58f, 2.06f, 2.61f, 3.35f, 4.34f, 5.51f, 6.70f, 7.75f, 8.63f, 9.29f, 9.59f, 9.42f, 8.78f, 7.82f, 6.50f, 4.95f, 3.47f, 2.15f, 1.38f, 1.11f, 0.90f]),
new("o", "o", "voiced_close_mid_back_rounded_vowel", [
0.38f, 0.45f, 0.57f, 0.77f, 1.31f, 1.92f, 1.74f, 1.11f, 0.75f, 0.59f, 0.57f, 0.68f, 0.73f, 0.67f, 0.58f, 0.49f, 0.44f, 0.42f, 0.49f, 0.53f, 0.38f, 0.30f, 0.45f, 0.61f, 0.71f, 0.79f, 0.86f, 1.01f, 1.41f, 2.09f, 3.00f, 4.10f, 5.16f, 6.22f, 7.34f, 8.15f, 8.61f, 8.37f, 6.76f, 4.37f, 2.30f, 1.06f, 0.58f, 0.47f]),
new("ʊ", "oo_short", "voiced_near_close_near_back_rounded_vowel", [
0.37f, 0.38f, 0.49f, 0.62f, 0.85f, 1.28f, 1.62f, 1.47f, 1.04f, 0.81f, 1.03f, 1.44f, 1.49f, 1.28f, 1.06f, 0.85f, 0.69f, 0.56f, 0.42f, 0.27f, 0.27f, 0.41f, 0.51f, 0.49f, 0.47f, 0.50f, 0.58f, 0.80f, 1.19f, 1.62f, 2.27f, 3.24f, 4.16f, 5.00f, 5.70f, 6.11f, 6.21f, 6.29f, 6.24f, 4.91f, 2.61f, 1.09f, 0.63f, 0.59f]),
new("u", "u", "voiced_close_back_rounded_vowel", [
0.54f, 0.61f, 0.66f, 0.75f, 1.13f, 1.99f, 2.83f, 2.90f, 2.52f, 2.40f, 2.83f, 3.56f, 3.99f, 3.89f, 3.50f, 3.04f, 2.64f, 2.44f, 2.31f, 2.07f, 1.80f, 1.52f, 1.14f, 0.74f, 0.42f, 0.22f, 0.14f, 0.20f, 0.47f, 0.89f, 1.15f, 1.42f, 2.17f, 3.04f, 3.69f, 4.70f, 5.74f, 5.41f, 3.82f, 2.34f, 1.35f, 0.65f, 0.29f, 0.16f]),
];
}

28
src/AquaSynth.Core/WeksaPhones.cs
View file

@ -71,19 +71,30 @@ public static class WeksaPhonesReader
/// <summary>
/// A Weksa phones line as a word the tract can render. Every phone must be in the DSL's IPA table
/// (<see cref="PatchScript.TryIpaDescriptor"/>); a phone outside it is listed in <see cref="Unsupported"/> and the plan
/// carries no targets, so nothing is approximated. Duration is the vowel/consonant length AquaSynth owns; stress moves
/// only the gesture intensity, the one thing the word script uses stress for.
/// carries no targets, so nothing is approximated. Base durations are the vowel/consonant lengths AquaSynth owns; stress
/// carries its three acoustic cues (<see cref="StressCue"/>): the marked phone is longer, louder and higher in f0.
/// </summary>
public sealed record WeksaWordPlan(IReadOnlyList<(IpaGestureExperimentTarget Target, float Intensity)> Phones, IReadOnlyList<WeksaPhone> Unsupported)
{
public const float VowelSeconds = 0.18f, ConsonantSeconds = 0.095f;
public bool Renderable => Unsupported.Count == 0;
public static float StressIntensity(WeksaStress stress) => stress switch
/// <summary>Duration factor, amplitude factor on the source flow and f0 ratio of a stressed phone against an unstressed one.</summary>
public readonly record struct StressCue(float Duration, float Amplitude, float PitchRatio);
/// <summary>
/// Primary stress: duration 220/160 = 1.375 (stressed against unstressed syllable in initial-stress pairs, Sluijter and van
/// Heuven 1996), amplitude 1.30 = 2.3 dB on the source flow, f0 +6 percent (their unaccented f0 difference). The amplitude is set
/// by the rendered result, not the nominal 2 dB (their 1-3 dB overall-intensity difference): the f0 rise alone moves the harmonics against
/// F1 and changes the rendered level by -1.3 to +1.8 dB by vowel (i -0.7, a -0.3, u -1.3, e -0.8, o +1.8), so 1.26 rendered +1.7 dB on
/// average over those five and 1.30 renders +2.0. Secondary stress is placed halfway, in dB and in percent: a DESIGN VALUE, no
/// secondary-stress measurement was retrieved. Unheard.
/// </summary>
public static StressCue CueFor(WeksaStress stress) => stress switch
{
WeksaStress.Primary => 1.2f,
WeksaStress.Secondary => 1.1f,
_ => 1f
WeksaStress.Primary => new(220f / 160f, 1.30f, 1.06f),
WeksaStress.Secondary => new(1.19f, 1.14f, 1.03f),
_ => new(1f, 1f, 1f)
};
public static WeksaWordPlan From(WeksaPhones word, string label)
@ -93,8 +104,9 @@ public sealed record WeksaWordPlan(IReadOnlyList<(IpaGestureExperimentTarget Tar
var targets = word.Phones.Select((p, i) =>
{
PatchScript.TryIpaDescriptor(p.Ipa, out var descriptor);
var seconds = descriptor.EndsWith("_vowel", StringComparison.Ordinal) ? VowelSeconds : ConsonantSeconds;
return (new IpaGestureExperimentTarget($"{label}_{i}_{p.Ipa}", p.Ipa, "", seconds), StressIntensity(p.Stress));
var cue = CueFor(p.Stress);
var seconds = (descriptor.EndsWith("_vowel", StringComparison.Ordinal) ? VowelSeconds : ConsonantSeconds) * cue.Duration;
return (new IpaGestureExperimentTarget($"{label}_{i}_{p.Ipa}", p.Ipa, "", seconds, PitchRatio: cue.PitchRatio), cue.Amplitude);
}).ToArray();
return new WeksaWordPlan(targets, []);
}

13
src/AquaSynth.Faust/IpaTrialOrchestrator.cs
View file

@ -139,7 +139,8 @@ public static class IpaTrialOrchestrator
targetSet.Id,
targetSet.Targets.Select(target => target.Target),
targetSet.Variants,
options.TimelineBlocks);
options.TimelineBlocks,
options.SampleRate);
var analysis = IpaGestureExperiment.AnalyzeRound(round);
var targetById = targetSet.Targets.ToDictionary(target => target.Target.Id, StringComparer.OrdinalIgnoreCase);
@ -459,7 +460,7 @@ public static class IpaTrialOrchestrator
IReadOnlyList<ProbeTimelineSample> timeline)
{
var tokens = DescriptorTokens(target.Descriptor);
var expected = ExpectedSurfaces(tokens);
var expected = ExpectedSurfaces(tokens, target.Descriptor);
var touched = patch.ControlSplines
.Where(spline => spline.Enabled)
.Select(spline => spline.SurfacePath)
@ -478,7 +479,7 @@ public static class IpaTrialOrchestrator
return Math.Clamp(.5f * coverage + .25f * motion + .25f * primitive, 0, 1);
}
private static string[] ExpectedSurfaces(HashSet<string> tokens)
private static string[] ExpectedSurfaces(HashSet<string> tokens, string descriptor)
{
var surfaces = new List<string>
{
@ -491,11 +492,13 @@ public static class IpaTrialOrchestrator
{
surfaces.AddRange(
[
"/vocal/areas/0/area/tongue_index",
"/vocal/areas/0/area/tongue_diameter",
"/vocal/areas/0/area/lip_opening",
"/vocal/radiation/0/aperture"
]);
if (StoryVowelAreas.IdForDescriptor(descriptor) is { } vowelId)
{
surfaces.Add($"/vocal/areas/0/area/{StoryVowelAreas.SurfaceField(vowelId)}");
}
}
else
{

88
tests/AquaSynth.Dsl.Tests/LipReflectanceTests.cs Normal file
View file

@ -0,0 +1,88 @@
using System.Text.RegularExpressions;
using AquaSynth.Dsl;
namespace AquaSynth.Dsl.Tests;
// The ladder's lip termination is the open-end reflectance R(z) = -g (1-b)/(1-b z^-1) with b/(1-b)^2 = (a fs / c)^2, a = sqrt(A_lip / pi) in cm,
// c = 35000 cm/s (FaustExport.VocalLadder.cs). These tests evaluate the emitted Faust expression itself and check it against the derivation.
public sealed class LipReflectanceTests
{
private const int Rate = 48000;
private const string Script = """
morphology name=oral length_cm=17.52 diameters=.55,.75,1.05,1.45,1.25,.8 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=4 lip_opening=4
waveguide_path name=oral_path morphology=oral speed=350 loss=.998
source_port name=folds path=oral_path pressure=.7 tension=.55 opening=.5 noise=0 impedance=.3 flow_scale=.12 tilt=0
radiation_load name=mouth path=oral_path aperture=1 reflection=-1 impedance=.28
vocal_network name=voice paths=oral_path sources=folds radiation=mouth
vocal network=voice freq=100 gain=.5 sustain=1
""";
private static string Emit() => FaustEmitter.Emit(PatchScript.Parse(Script), new FaustExportOptions("lip_pole", SampleRate: Rate)).Source;
// b by bisection on the derivation b/(1-b)^2 = q, independent of the closed form.
private static double SolvePole(double areaCm2)
{
var a = Math.Sqrt(areaCm2 / Math.PI);
var q = Math.Pow(a * Rate / 35000.0, 2);
double lo = 0, hi = 1 - 1e-12;
for (var i = 0; i < 200; i++)
{
var mid = (lo + hi) / 2;
if (mid / Math.Pow(1 - mid, 2) < q) lo = mid; else hi = mid;
}
return Math.Clamp((lo + hi) / 2, 0, .9);
}
[Fact]
public void LipReflectanceKeepsReflectionAsDcGainAndFeedsBackThePoleInTheLastCellEquation()
{
var source = Emit();
var last = Regex.Match(source, @"'(\w+_L\d+) = \((\w+_lip_b) \* (\w+_L\d+) \+ \(1\.0 - (\w+_lip_b)\) \* \((.+?)\) \* (\w+_R\d+)\) \* clip01");
Assert.True(last.Success, "last-cell equation is not the one-pole reflectance");
Assert.Equal(last.Groups[1].Value, last.Groups[3].Value);
Assert.Equal(last.Groups[2].Value, last.Groups[4].Value);
Assert.False(string.IsNullOrWhiteSpace(last.Groups[5].Value));
}
[Theory]
[InlineData(.3)]
[InlineData(1.0)]
[InlineData(3.14159265)]
[InlineData(8.0)]
public void EmittedPoleSolvesTheLowFrequencyLossDerivationForTheLipCellArea(double areaCm2)
{
var source = Emit();
var q = Regex.Match(source, @"^\s*(\w+)_lip_q = (.+);$", RegexOptions.Multiline);
var b = Regex.Match(source, @"^\s*(\w+)_lip_b = (.+);$", RegexOptions.Multiline);
Assert.True(q.Success && b.Success);
var prefix = q.Groups[1].Value;
// The lip pole reads the last cell's area and nothing else of the tract.
var lastCell = Regex.Matches(source, $@"^\s*{prefix}_A(\d+) = ", RegexOptions.Multiline).Count - 1;
Assert.Contains($"{prefix}_A{lastCell}", q.Groups[2].Value);
var program = $"import(\"stdfaust.lib\");\n{prefix}_A{lastCell} = {areaCm2.ToString(System.Globalization.CultureInfo.InvariantCulture)};\n" +
$"{prefix}_lip_q = {q.Groups[2].Value};\n{prefix}_lip_b = {b.Groups[2].Value};\nprocess = {prefix}_lip_b;\n";
using var compiler = new AquaSynthPatchCompiler(new AquaSynthNativeOptions(SampleRate: Rate));
if (!compiler.TryCompileSource(new AquaSynthCompileIdentity("lip_pole", "lip_pole", ""), program, .05f, out var patch, out var error))
{
if (error?.Contains("not found", StringComparison.OrdinalIgnoreCase) == true) return;
Assert.Fail($"lip pole program failed: {error}");
}
using (patch)
{
var samples = patch!.Render();
Assert.Equal(SolvePole(areaCm2), samples[^1], 3);
}
}
[Fact]
public void PoleGrowsWithTheLipAreaAndEndCorrectionOfAOneCentimetreRadiusIsAboutHalfARadius()
{
// b/(1-b) samples is the DC group delay: for a = 1 cm at 48 kHz it is 0.96 samples = 0.7 cm, which is the 0.6 a end correction scale.
Assert.True(SolvePole(.5) < SolvePole(2) && SolvePole(2) < SolvePole(6));
var b = SolvePole(Math.PI);
Assert.InRange(b / (1 - b) * 35000.0 / Rate, .5, .8);
}
}

113
tests/AquaSynth.Dsl.Tests/NasalCouplingTests.cs Normal file
View file

@ -0,0 +1,113 @@
using System.Globalization;
using System.Text.RegularExpressions;
using AquaSynth.Dsl;
namespace AquaSynth.Dsl.Tests;
// The nasal side branch: the cited nasal tract, the velopharyngeal port area, and where the scaffold's nasal and bilabial/alveolar gestures
// close the mouth. The murmur itself is measured in VocalTractFormantTests (renders).
public sealed class NasalCouplingTests
{
// Birkholz 2005 Table 2.3 after Dang & Honda 1994, copied here independently of NasalTractAreas: 19 sections of 0.6 cm, cm^2,
// velopharyngeal end first.
private static readonly float[] CitedAreas =
[
1.63f, 2.07f, 2.72f, 3.59f, 4.24f, 3.26f, 3.04f, 3.04f, 2.72f, 2.5f, 2.39f, 2.39f, 1.85f, 0.76f, 1.41f, 1.74f, 1.30f, 1.74f, 0.76f,
];
[Fact]
public void TheNasalTractTableIsTheCitedNineteenSections()
{
Assert.Equal(CitedAreas, NasalTractAreas.Areas.ToArray());
Assert.Equal(11.4f, NasalTractAreas.LengthCentimeters, 4);
Assert.Equal(NasalTractAreas.LengthCentimeters, CitedAreas.Length * NasalTractAreas.SectionLengthCentimeters, 3);
}
[Fact]
public void TheScaffoldNasalMorphologyIsTheCitedTableAtItsLength()
{
var patch = PatchScript.Parse(IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget("v_m", "m", "", .2f), 1f)]));
var nasal = patch.AreaFunctions.Single(area => area.Name == "nasal").Shape;
Assert.Equal(11.4f, nasal.LengthCentimeters, 3);
Assert.Equal(CitedAreas.Length, nasal.Diameters.Count);
for (var i = 0; i < CitedAreas.Length; i++)
{
Assert.Equal(CitedAreas[i], nasal.Diameters[i] * nasal.Diameters[i], 2);
}
}
// The emitted velum area expression, evaluated with the admittance (opening * coupling) it reads set to a constant: the port area is
// linear in the opening with the cited maximum, 0.8 cm^2.
[Theory]
[InlineData(1.0, 0.8)]
[InlineData(0.7, 0.56)]
[InlineData(0.25, 0.2)]
public void TheVelopharyngealPortAreaIsTheOpeningTimesTheMaximumPortArea(double opening, double expectedAreaSquareCentimeters)
{
var source = FaustEmitter.Emit(
PatchScript.Parse(IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget("v_m", "m", "", .2f), 1f)])),
new FaustExportOptions("velum", SampleRate: 48000)).Source;
var area = Regex.Match(source, @"^\s*(\w+)_velum_area = (.+);$", RegexOptions.Multiline);
Assert.True(area.Success);
var admittance = Regex.Match(area.Groups[2].Value, @"\b(\w+_admittance)\b");
Assert.True(admittance.Success, area.Groups[2].Value);
var program = $"import(\"stdfaust.lib\");\n{admittance.Groups[1].Value} = {opening.ToString(CultureInfo.InvariantCulture)};\n" +
$"area = {area.Groups[2].Value};\nprocess = area;\n";
using var compiler = new AquaSynthPatchCompiler(new AquaSynthNativeOptions(SampleRate: 48000));
if (!compiler.TryCompileSource(new AquaSynthCompileIdentity("velum", "velum", ""), program, .05f, out var compiled, out var error))
{
if (error?.Contains("not found", StringComparison.OrdinalIgnoreCase) == true) return;
Assert.Fail($"velum program failed: {error}");
}
using (compiled)
{
Assert.Equal(expectedAreaSquareCentimeters, compiled!.Render()[^1], 3);
}
}
private static float Spline(SynthPatch patch, string namePrefix, string surfacePath) =>
patch.ControlSplines.Single(spline => spline.Name.StartsWith(namePrefix, StringComparison.Ordinal) && spline.SurfacePath == surfacePath).Points[0].Value;
// Distance of a vowel_set consonant constriction from the lip plane in cm, from the constriction_index spline: the surface is a section
// index in [0, sections] and the area model reads position = index / (sections - 1).
private static double CentimetresFromLips(string ipa)
{
var patch = PatchScript.Parse(IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget("v_c", ipa, "", .2f), 1f)]));
var oral = patch.AreaFunctions.Single(area => area.Name == "oral");
var index = Spline(patch, "v_c_constriction_index", "/vocal/areas/0/area/constriction_index") * oral.Shape.Sections;
return (1 - index / (oral.Shape.Sections - 1)) * oral.Shape.LengthCentimeters;
}
[Fact]
public void UnderVowelSetMClosesAtTheLipsAndNAtTheAlveolarRidgeInCentimetres()
{
Assert.InRange(CentimetresFromLips("m"), -.05, .05);
Assert.InRange(CentimetresFromLips("n"), 2.4, 2.7);
Assert.True(CentimetresFromLips("n") > CentimetresFromLips("m") + 2);
}
[Fact]
public void WithoutVowelSetTheConstrictionStaysOnTheNormalizedPlaceFraction()
{
// The earlier surface semantics are kept for hand-written scripts: the place fraction is the normalized spline value itself.
var patch = PatchScript.Parse("""
morphology name=oral length_cm=17 diameters=.6,.8,1.2,1.6,1.3,.9 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=1 lip_opening=1.5
morphology name=nasal length_cm=12 diameters=.05,.35,.6,.8
waveguide_path name=oral_path morphology=oral loss=.998
waveguide_path name=nasal_path morphology=nasal loss=.997
source_port name=folds path=oral_path pressure=.7 tension=.55 opening=.45 noise=.05 impedance=.3
branch_port name=velopharynx from=oral_path from_position=.45 to=nasal_path opening=0 coupling=1
constriction_contact name=contact path=oral_path position=.92 opening=.5 resistance=.4 stored_pressure=.1
radiation_load name=mouth path=oral_path aperture=.8 reflection=-.82 impedance=.28
vocal_network name=voice paths=oral_path,nasal_path sources=folds contacts=contact branches=velopharynx radiation=mouth
phoneme_gesture name=en ipa=n start=0 dur=.16 intensity=.8
phoneme_gesture name=em ipa=m start=.2 dur=.16 intensity=.8
vocal network=voice freq=150 gain=.2 sustain=.6
""");
Assert.Equal(.78f, Spline(patch, "en_constriction_index", "/vocal/areas/0/area/constriction_index"));
Assert.Equal(.96f, Spline(patch, "em_constriction_index", "/vocal/areas/0/area/constriction_index"));
}
}

5
tests/AquaSynth.Dsl.Tests/PatchScriptTests.cs
View file

@ -582,6 +582,7 @@ public sealed class PatchScriptTests
{
var patch = PatchScript.Parse("""
morphology name=oral length_cm=17 diameters=.6,.8,1.2,1.6,1.3,.9 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=1 lip_opening=1.5
vowel_set area_function=oral source=story2008
morphology name=nasal length_cm=12 diameters=.05,.35,.6,.8
waveguide_path name=oral_path morphology=oral loss=.998
waveguide_path name=nasal_path morphology=nasal loss=.997
@ -602,7 +603,9 @@ public sealed class PatchScriptTests
Assert.Contains(patch.ControlSplines, spline => spline.Name.StartsWith("ma_contact", StringComparison.Ordinal) && spline.SurfacePath == "/vocal/contacts/0/opening");
Assert.Contains(patch.ControlSplines, spline => spline.Name.StartsWith("hw_source_noise", StringComparison.Ordinal) && spline.SurfacePath == "/vocal/sources/0/noise");
Assert.Contains(patch.ControlSplines, spline => spline.Name.StartsWith("hw_constriction_index", StringComparison.Ordinal) && spline.SurfacePath == "/vocal/areas/0/area/constriction_index");
Assert.Contains(patch.ControlSplines, spline => spline.Name.StartsWith("vowel_tongue_index", StringComparison.Ordinal) && spline.SurfacePath == "/vocal/areas/0/area/tongue_index");
Assert.Contains(patch.ControlSplines, spline => spline.Name.StartsWith("vowel_vowel_ah", StringComparison.Ordinal) && spline.SurfacePath == "/vocal/areas/0/area/vowel_ah" && spline.Points.All(point => point.Value == 1));
Assert.Contains(patch.ControlSplines, spline => spline.Name.StartsWith("vowel_vowel_i", StringComparison.Ordinal) && spline.SurfacePath == "/vocal/areas/0/area/vowel_i" && spline.Points.All(point => point.Value == 0));
Assert.DoesNotContain(patch.ControlSplines, spline => spline.SurfacePath == "/vocal/areas/0/area/tongue_index" && spline.Name.StartsWith("vowel_", StringComparison.Ordinal));
var catalog = ControlSurfaceCatalog.FromPatch(patch);
var timeline = catalog.CreateTimeline();

31
tests/AquaSynth.Dsl.Tests/ProbeTimelineSampleRateTests.cs Normal file
View file

@ -0,0 +1,31 @@
using AquaSynth.Dsl;
namespace AquaSynth.Dsl.Tests;
public sealed class ProbeTimelineSampleRateTests
{
private const string TimelineScript = """
morphology name=oral length_cm=17.52 diameters=.55,.75,1.05,1.45,1.25,.8 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=4 lip_opening=4
waveguide_path name=oral_path morphology=oral speed=350 loss=.998
source_port name=folds path=oral_path pressure=.7 tension=.55 opening=.5 noise=0 impedance=.3 flow_scale=.12
radiation_load name=mouth path=oral_path aperture=.8 reflection=-.82 impedance=.28
probe_timeline name=flow blocks=4 block_size=64
vocal_network name=voice paths=oral_path sources=folds radiation=mouth probe=flow
vocal network=voice freq=140 gain=.2 sustain=.3
""";
[Fact]
public void ProbeTimelineReadsTheSampleRateItIsGiven()
{
var patch = PatchScript.Parse(TimelineScript + "\ncontrol_spline name=speed_step surface=/vocal/paths/0/speed interp=hold points=0:.25,0.0028:.5");
float Delay(int block, float? rate) => (rate is { } r ? ProbeTimelineReport.Build(patch, "voice", blocks: 3, sampleRate: r) : ProbeTimelineReport.Build(patch, "voice", blocks: 3))
.First(s => s.Signal == "delay_samples" && s.Block == block).Value;
// Block 0 (t = 0) reads speed 325 m/s: the delay is length / speed * rate, so 44.1 kHz and 48 kHz stand in the ratio of the rates.
Assert.Equal(0.1752f / 325 * 48000, Delay(0, 48000), 3);
Assert.Equal(0.1752f / 325 * 44100, Delay(0, 44100), 3);
Assert.Equal(Delay(0, 48000), Delay(0, null), 3);
// Block 2 sits at 64 * 2 / rate = 2.67 ms at 48 kHz (before the 2.8 ms speed step: 325 m/s) and 2.90 ms at 44.1 kHz (after it: 550 m/s).
Assert.Equal(0.1752f / 325 * 48000, Delay(2, 48000), 3);
Assert.Equal(0.1752f / 550 * 44100, Delay(2, 44100), 3);
}
}

74
tests/AquaSynth.Dsl.Tests/StoryVowelAreasTests.cs Normal file
View file

@ -0,0 +1,74 @@
using System.Globalization;
using AquaSynth.Dsl;
namespace AquaSynth.Dsl.Tests;
// StoryVowelAreas.cs is the code owner of Story (2008) Table I; docs/research/data/story2008-area-functions.csv is the cited copy.
// Nothing here says how any vowel sounds: these tests guard the transcription, the descriptor lookup and the sample-rate threading.
public sealed class StoryVowelAreasTests
{
private static string RepoFile(params string[] relative)
{
for (var dir = new DirectoryInfo(AppContext.BaseDirectory); dir is not null; dir = dir.Parent)
{
if (File.Exists(Path.Combine(dir.FullName, "AquaSynth.sln")) || Directory.Exists(Path.Combine(dir.FullName, "docs", "research")))
{
return Path.Combine([dir.FullName, .. relative]);
}
}
throw new InvalidOperationException("repository root not found above " + AppContext.BaseDirectory);
}
[Fact]
public void EveryCellOfTheCodeTableEqualsTheCitedCsv()
{
var lines = File.ReadAllLines(RepoFile("docs", "research", "data", "story2008-area-functions.csv"))
.Where(line => !line.StartsWith('#') && line.Length > 0).ToArray();
var header = lines[0].Split(',');
Assert.Equal("section", header[0]);
Assert.Equal(StoryVowelAreas.All.Select(v => v.Ipa), header.Skip(1));
var sections = lines.Skip(1).Where(line => char.IsDigit(line[0])).Select(line => line.Split(',')).ToArray();
Assert.Equal(StoryVowelAreas.Sections, sections.Length);
for (var v = 0; v < StoryVowelAreas.All.Count; v++)
{
Assert.Equal(StoryVowelAreas.Sections, StoryVowelAreas.All[v].Areas.Count);
for (var s = 0; s < sections.Length; s++)
{
Assert.Equal(float.Parse(sections[s][v + 1], CultureInfo.InvariantCulture), StoryVowelAreas.All[v].Areas[s]);
}
}
}
[Fact]
public void AsciiIdsAreUniqueAndDescriptorsSelectTheirVowelExactly()
{
Assert.Equal(11, StoryVowelAreas.All.Select(v => v.Id).Distinct().Count());
Assert.Equal(11, StoryVowelAreas.All.Select(v => v.Descriptor).Distinct().Count());
foreach (var vowel in StoryVowelAreas.All)
{
Assert.True(vowel.Descriptor.EndsWith("_vowel", StringComparison.Ordinal));
Assert.Equal(vowel.Id, StoryVowelAreas.IdForDescriptor(vowel.Descriptor));
Assert.True(PatchScript.TryIpaDescriptor(vowel.Ipa, out var descriptor));
Assert.Equal(vowel.Descriptor, descriptor);
}
Assert.Equal("ah", StoryVowelAreas.IdForDescriptor("voiced_open_back_unrounded_vowel"));
PatchScript.TryIpaDescriptor("a", out var a);
Assert.Equal("ah", StoryVowelAreas.IdForDescriptor(a));
// The schwa and a vowel descriptor with no Story shape both select no vowel: the geometric-mean shape.
PatchScript.TryIpaDescriptor("ə", out var schwa);
Assert.Null(StoryVowelAreas.IdForDescriptor(schwa));
Assert.Null(StoryVowelAreas.IdForDescriptor("voiced_close_central_unrounded_vowel"));
}
[Fact]
public void ResamplingInterpolatesBetweenSectionCentresAndHoldsTheEnds()
{
float[] areas = [1, 3, 5, 7];
Assert.Equal(1, StoryVowelAreas.AreaAt(areas, 0));
Assert.Equal(7, StoryVowelAreas.AreaAt(areas, 1));
Assert.Equal(3, StoryVowelAreas.AreaAt(areas, 0.375f), 4);
Assert.Equal(4, StoryVowelAreas.AreaAt(areas, 0.5f), 4);
}
}

206
tests/AquaSynth.Dsl.Tests/StressAndReleaseRenderTests.cs Normal file
View file

@ -0,0 +1,206 @@
using AquaSynth.Dsl;
using AquaSynth.Faust;
using Xunit.Abstractions;
namespace AquaSynth.Dsl.Tests;
// Renders the two cues the tract gained in this pass through the native Faust path: the source's f0 ratio (the stress pitch cue) and
// the plosive release burst. Each asserts a measured signal, not the emitted text.
[Collection("NativeFaust")]
public sealed class StressAndReleaseRenderTests(ITestOutputHelper output)
{
private const int Rate = 48000;
private static float[]? Render(string script, float seconds, Func<SynthPatch, SynthPatch>? edit = null)
{
var parsed = PatchScript.Parse(script);
var source = FaustEmitter.Emit(edit is null ? parsed : edit(parsed), new FaustExportOptions("stress_release", SampleRate: Rate)).Source;
using var compiler = new AquaSynthPatchCompiler(new AquaSynthNativeOptions(SampleRate: Rate));
if (!compiler.TryCompileSource(new AquaSynthCompileIdentity("stress_release", "stress_release", ""), source, seconds, out var patch, out var error))
{
if (error?.Contains("Faust toolchain not found", StringComparison.OrdinalIgnoreCase) == true ||
error?.Contains("Faust DLL not found", StringComparison.OrdinalIgnoreCase) == true)
{
return null;
}
Assert.Fail($"compile failed: {error}");
}
using (patch) return patch!.Render(seconds);
}
private static double Rms(float[] x, double from, double to)
{
var a = (int)(from * Rate);
var b = Math.Min(x.Length, (int)(to * Rate));
double sum = 0;
for (var i = a; i < b; i++) sum += (double)x[i] * x[i];
return Math.Sqrt(sum / Math.Max(1, b - a));
}
private static double Peak(float[] x, double from, double to) =>
x.Skip((int)(from * Rate)).Take((int)((to - from) * Rate)).Max(Math.Abs);
// f0 by the highest autocorrelation peak between 60 and 400 Hz.
private static double F0(float[] x, double from, double to)
{
var seg = x.Skip((int)(from * Rate)).Take((int)((to - from) * Rate)).Select(v => (double)v).ToArray();
double best = 0; var bestLag = 0;
for (var lag = Rate / 400; lag <= Rate / 60; lag++)
{
double c = 0;
for (var i = 0; i + lag < seg.Length; i++) c += seg[i] * seg[i + lag];
if (c > best) { best = c; bestLag = lag; }
}
return (double)Rate / bestLag;
}
[Fact]
public void SourcePitchRatioScalesTheRenderedF0()
{
var f0 = (float ratio) =>
{
var script = IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget("w_0_a", "a", "", 0.5f, PitchRatio: ratio), 1f)], 110);
var x = Render(script, 0.6f);
return x is null ? double.NaN : F0(x, 0.2, 0.5);
};
var unity = f0(1f);
var up = f0(1.5f);
output.WriteLine($"f0 at ratio 1: {unity:F1}, ratio 1.5: {up:F1}");
if (double.IsNaN(unity)) return;
Assert.InRange(unity, 108, 112);
Assert.InRange(up, 163, 167);
}
[Fact]
public void IntensityIsAnAmplitudeFactorOnTheRenderedLevel()
{
double Level(float intensity)
{
var script = IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget("w_0_a", "a", "", 0.5f), intensity)], 110);
var x = Render(script, 0.6f);
return x is null ? double.NaN : Rms(x, 0.3, 0.5);
}
var unity = Level(1f);
var stressed = Level(1.26f);
if (double.IsNaN(unity)) return;
var db = 20 * Math.Log10(stressed / unity);
output.WriteLine($"level at intensity 1.26 against 1: {db:F2} dB (factor wanted {20 * Math.Log10(1.26):F2})");
Assert.InRange(db, 1.7, 2.3);
}
[Fact]
public void PlosiveReleaseBurstIsAboveTheClosureAndBelowTheFollowingVowel()
{
var script = IpaGestureExperiment.BuildWordScript(
[(new IpaGestureExperimentTarget("w_0_b", "b", "", 0.095f), 1f), (new IpaGestureExperimentTarget("w_1_a", "a", "", 0.18f), 1f)], 140);
var withBurst = Render(script, 0.4f);
var withoutBurst = Render(script, 0.4f, p => p with { ControlSplines = p.ControlSplines.Where(s => !s.SurfacePath.EndsWith("/release_flow", StringComparison.Ordinal)).ToArray() });
if (withBurst is null || withoutBurst is null) return;
// Everything but the burst's noise is deterministic and shared, so the difference of the two renders is the burst alone.
var burst = withBurst.Zip(withoutBurst, (x, y) => x - y).ToArray();
const double releaseFrom = 0.095 * 0.65, vowelFrom = 0.095;
var vowelPeak = Peak(withBurst, vowelFrom + 0.06, vowelFrom + 0.18);
var burstPeak = Peak(burst, releaseFrom, vowelFrom + 0.01);
var burstBeforeRelease = Peak(burst, 0, releaseFrom - 0.005);
var burstAfter = Peak(burst, vowelFrom + 0.04, 0.4);
output.WriteLine($"burst peak {burstPeak:F4} in the release, {burstBeforeRelease:F5} in the closure, {burstAfter:F5} in the vowel; vowel peak {vowelPeak:F4}");
Assert.True(burstPeak > 0.004, "the burst is a transient of audible size in the release window");
Assert.True(burstBeforeRelease < burstPeak / 20, "the closure is silent: no burst before the release");
Assert.True(burstAfter < burstPeak / 5, "the burst is over once the vowel is open");
Assert.True(burstPeak < vowelPeak / 2, "the burst stays well under the following vowel's peak (the b level fix)");
}
// /ɑ C ɑ/ at the scaffold's own durations (vowel 0.18 s, consonant 0.095 s). Returns the render and the consonant's [start, end).
private (float[]? Samples, double From, double To) VowelConsonantVowel(string ipa, float frequency = 140)
{
var script = IpaGestureExperiment.BuildWordScript(
[(new IpaGestureExperimentTarget("w_0_a", "a", "", 0.18f), 1f), (new IpaGestureExperimentTarget("w_1_c", ipa, "", 0.095f), 1f), (new IpaGestureExperimentTarget("w_2_a", "a", "", 0.18f), 1f)], frequency);
return (Render(script, 0.5f), 0.18, 0.275);
}
// Level of the consonant's middle 60 percent against the mean power of the steady parts of the two flanking vowels (the 70 ms before the
// consonant and 90 ms from 30 ms after it), in dB.
private static double ConsonantDbReVowels(float[] x, double from, double to)
{
var v1 = Rms(x, from - 0.08, from - 0.01);
var v2 = Rms(x, to + 0.03, to + 0.12);
var vowels = Math.Sqrt((v1 * v1 + v2 * v2) / 2);
var span = to - from;
return 20 * Math.Log10(Math.Max(1e-9, Rms(x, from + 0.2 * span, to - 0.2 * span)) / vowels);
}
// Fujimura 1962 and later accounts put the nasal murmur well under the flanking vowels; the target band is 8 to 16 dB under (a DESIGN
// value: no primary level was retrieved, docs/research/articulatory-vocal-tract-sources.md). Two defects held it out of the band: /n/
// read -3.8 dB while its alveolar closure began to open at half the gesture, a third before the velum and the contact did (now held to
// 0.8: /n/ -8.8, /m/ -12.8), and once the nasal junction sat on its own interface the nose took its share and both read -2.2 and -4.2 dB
// until the nasal path carried its real damping (loss .975 per cell; /m/ -11.8, /n/ -9.6).
[Theory]
[InlineData("m")]
[InlineData("n")]
public void NasalMurmurIsWellBelowTheFlankingVowels(string ipa)
{
var (x, from, to) = VowelConsonantVowel(ipa);
if (x is null) return;
var db = ConsonantDbReVowels(x, from, to);
output.WriteLine($"/{ipa}/ murmur {db:F1} dB re the flanking vowels");
Assert.InRange(db, -16, -8);
}
// The primary-stress cue as the plan applies it (amplitude and f0 ratio, same vowel length) against the unstressed vowel, over five vowels:
// the f0 rise alone moves the harmonics against F1 and the level by -1.3 to +1.8 dB by vowel, so the amplitude factor is held to the mean
// rendered result, 2 dB (Sluijter and van Heuven 1996: 1 to 3 dB overall). Measured 1.7 dB with the amplitude at 1.26 (outside the band), 2.0 dB at 1.30.
[Fact]
public void PrimaryStressRendersAboutTwoDecibelsLouderOnAverageOverVowels()
{
var cue = WeksaWordPlan.CueFor(WeksaStress.Primary);
var diffs = new List<double>();
foreach (var vowel in new[] { "i", "a", "u", "e", "o" })
{
double Level(float ratio, float intensity)
{
var script = IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget($"w_0_{vowel}", vowel, "", 0.25f, PitchRatio: ratio), intensity)], 140, 0.2f);
var x = Render(script, 0.45f);
return x is null ? double.NaN : Rms(x, 0.12, 0.22);
}
var flat = Level(1f, 1f);
if (double.IsNaN(flat)) return;
var db = 20 * Math.Log10(Level(cue.PitchRatio, cue.Amplitude) / flat);
output.WriteLine($"/{vowel}/ stressed against unstressed: {db:F2} dB");
diffs.Add(db);
}
output.WriteLine($"mean {diffs.Average():F2} dB");
Assert.InRange(diffs.Average(), 1.8, 2.3);
}
// Standard deviation, so the slow offset a shutting tract leaves behind does not count as sound.
private static double Deviation(float[] x, double from, double to)
{
var seg = x.Skip((int)(from * Rate)).Take((int)((to - from) * Rate)).Select(v => (double)v).ToArray();
var mean = seg.Average();
return Math.Sqrt(seg.Sum(v => (v - mean) * (v - mean)) / seg.Length);
}
// A word-initial plosive's closure (0.03 to 0.058 s of a 0.095 s consonant shut to 0.062 s) against the vowel after it (0.15 to 0.25 s). The
// mouth is shut, so all that gets out is the tract wall radiation of the voiced source: the voice bar, about 30 dB under the vowel for /b/
// (a DESIGN value, no primary level retrieved; Stevens 1998 per secondary sources has voicing weaker as transglottal pressure falls), and
// near-silent for /p/, whose source is weak. The closure was 0.0000 (below 16 bits) before the wall radiation.
[Theory]
[InlineData("b", -34, -26)]
[InlineData("p", -60, -38)]
public void PlosiveClosureCarriesTheVoiceBarOnlyWhenVoiced(string ipa, double lowDb, double highDb)
{
var script = IpaGestureExperiment.BuildWordScript(
[(new IpaGestureExperimentTarget("w_0_c", ipa, "", 0.095f), 1f), (new IpaGestureExperimentTarget("w_1_a", "a", "", 0.18f), 1f)], 140);
var x = Render(script, 0.4f);
if (x is null) return;
var db = 20 * Math.Log10(Math.Max(1e-9, Deviation(x, 0.03, 0.058)) / Deviation(x, 0.15, 0.25));
output.WriteLine($"/{ipa}/ closure {db:F1} dB re the following vowel");
Assert.InRange(db, lowDb, highDb);
}
}

310
tests/AquaSynth.Dsl.Tests/VocalTractFormantTests.cs
View file

@ -14,24 +14,25 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
private const int Rate = 48000;
private const float Pitch = 10f;
private static string Script(float tongueIndex, float tongueDiameter, float lipOpening, float velum = .01f) => $$"""
private static string Script(float tongueIndex, float tongueDiameter, float lipOpening, float velum = 0f) => $$"""
morphology name=oral length_cm=17 diameters=.55,.75,1.05,1.45,1.25,.8 tongue_index={{tongueIndex}} tongue_diameter={{tongueDiameter}} constriction_index=4 constriction_diameter=4 lip_opening={{lipOpening}}
morphology name=nasal length_cm=12 diameters=.04,.28,.52,.72
waveguide_path name=oral_path morphology=oral loss=.998
waveguide_path name=nasal_path morphology=nasal loss=.997
source_port name=folds path=oral_path pressure=.7 tension=.55 opening=.5 noise=0 impedance=.3 flow_scale=.12 tilt=0
branch_port name=velopharynx from=oral_path from_position=.45 to=nasal_path opening={{velum}} coupling=1
radiation_load name=mouth path=oral_path aperture=.8 reflection=-.82 impedance=.28
radiation_load name=nostril path=nasal_path aperture=.8 reflection=-.82 impedance=.28
radiation_load name=mouth path=oral_path aperture=.8 reflection=-1 impedance=.28
radiation_load name=nostril path=nasal_path aperture=.8 reflection=-1 impedance=.28
vocal_network name=voice paths=oral_path,nasal_path sources=folds branches=velopharynx radiation=mouth,nostril
vocal network=voice freq={{Pitch}} min_freq=5 gain=.5 sustain=1.5
""";
// Three seconds of the 10 Hz source would be wasteful: one second after the attack is ten whole periods, so every
// harmonic of Pitch sits exactly on a DFT bin.
private float[]? Render(string script)
private float[]? Render(string script, Func<SynthPatch, SynthPatch>? edit = null)
{
var source = FaustEmitter.Emit(PatchScript.Parse(script), new FaustExportOptions("vocal_formants", SampleRate: Rate)).Source;
var parsed = PatchScript.Parse(script);
var source = FaustEmitter.Emit(edit is null ? parsed : edit(parsed), new FaustExportOptions("vocal_formants", SampleRate: Rate)).Source;
using var compiler = new AquaSynthPatchCompiler(new AquaSynthNativeOptions(SampleRate: Rate));
if (!compiler.TryCompileSource(new AquaSynthCompileIdentity("vocal_formants", "vocal_formants", ""), source, 1.4f, out var patch, out var error))
{
@ -163,7 +164,7 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
[Fact]
public void LadderStaysBoundedForASecondWithTheVelumOpen()
{
foreach (var (ti, td, lo, velum) in new[] { (3f, .15f, 2.2f, .01f), (1.5f, .5f, 2.2f, 1f), (2f, 1.2f, .6f, 1f) })
foreach (var (ti, td, lo, velum) in new[] { (3f, .15f, 2.2f, 0f), (1.5f, .5f, 2.2f, 1f), (2f, 1.2f, .6f, 1f) })
{
var window = Render(Script(ti, td, lo, velum));
if (window is null) return;
@ -178,22 +179,79 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
private static double Rms(float[] window) => Math.Sqrt(window.Sum(x => (double)x * x) / window.Length);
// Rms of the part of the window above a frequency (Parseval over the bins of a zero-padded transform). The tract wall radiates a voice bar
// below about 300 Hz with the mouth shut and the velum shut, so "the nose is the only door" is read where the wall barely reaches.
private static double RmsAbove(float[] window, double hz) => RmsBetween(window, hz, Rate / 2.0);
private static double RmsBetween(float[] window, double lowHz, double highHz)
{
var n = 1 << 17;
var re = new double[n];
var im = new double[n];
for (var i = 0; i < window.Length; i++) re[i] = window[i] * (0.5 - 0.5 * Math.Cos(2 * Math.PI * i / window.Length));
Fft(re, im);
double power = 0;
for (var k = (int)(lowHz * n / Rate); k < Math.Min(n / 2, (int)(highHz * n / Rate)); k++) power += re[k] * re[k] + im[k] * im[k];
var windowPower = Enumerable.Range(0, window.Length).Sum(i => Math.Pow(0.5 - 0.5 * Math.Cos(2 * Math.PI * i / window.Length), 2));
return Math.Sqrt(2 * power / n / windowPower);
}
// The wall term's band, read from the render of a shut mouth and a shut velum (the wall is then the only path out). Behaviour, not
// brick walls: the 300 Hz low-pass is a second-order roll-off and the 90 Hz high-pass a second-order shelf, so these compare band levels.
private float[]? ShutMouthShutVelum() => Render(MouthShutScript(0f));
private static string MouthShutScript(float velum) => Script(1.5f, .5f, 2.2f, velum)
.Replace("constriction_diameter=4", "constriction_diameter=.001", StringComparison.Ordinal)
.Replace("diameters=.04,.28,.52,.72", "diameters=.9,1.3,1.5,1.2", StringComparison.Ordinal);
// The velum-open junction's nasal coefficient, (2 Yv - sum Y) / sum Y: with the mouth shut the first resonance the nose lets out sits at
// 520 Hz (400 Hz when the coefficient is (Yv - sum Y) / sum Y, and output above 600 Hz falls from 39 to 24 times the shut-velum level).
[Fact]
public void TheNasalReflectionCoefficientPlacesTheFirstResonanceOfAShutMouthWithTheVelumOpen()
{
var open = Render(MouthShutScript(1f));
if (open is null) return;
var first = Peaks(open)[0];
output.WriteLine($"mouth shut, velum open: first peak {first:F0} Hz");
Assert.InRange(first, 460, 580);
}
[Fact]
public void TheWallBarRollsOffBetweenItsTwoCornersAndOpensAboveTheSlowDrift()
{
var shut = ShutMouthShutVelum();
if (shut is null) return;
var drift = RmsBetween(shut, 5, 40);
var bar = RmsBetween(shut, 90, 200);
var upper = RmsBetween(shut, 230, 300);
var above = RmsAbove(shut, 600);
output.WriteLine($"wall bar band levels: 5-40 Hz {drift:E2}, 90-200 Hz {bar:E2}, 230-300 Hz {upper:E2}, above 600 Hz {above:E2}");
// The high-pass keeps the slow pressure drift out: measured 5-40 Hz at 0.4 of the bar, and 3.2 of it with the corner at 20 Hz.
Assert.True(drift < bar, $"sub-90 Hz drift {drift:E2} vs the 90-200 Hz bar {bar:E2}");
// The low-pass at 300 Hz still passes the 230-300 Hz harmonics: measured 0.14 of the bar, and 0.07 of it with the corner at 150 Hz.
Assert.True(upper > bar / 10, $"230-300 Hz {upper:E2} vs the 90-200 Hz bar {bar:E2}");
// The bar is not a brick wall but it is a bar: the low end stands well over what leaks above 600 Hz.
Assert.True(bar > 10 * above, $"bar {bar:E2} vs above 600 Hz {above:E2}");
}
[Fact]
public void WithTheMouthClosedOnlyAnOpenVelumLetsSoundOutThroughTheNose()
{
// constriction_diameter .001 at tract index 4 shuts the oral ladder; the only way out is the junction into the nasal ladder.
string Closed(float velum) => Script(1.5f, .5f, 2.2f, velum)
.Replace("constriction_diameter=4", "constriction_diameter=.001", StringComparison.Ordinal)
.Replace("diameters=.04,.28,.52,.72", "diameters=.9,1.3,1.5,1.2", StringComparison.Ordinal);
var shut = Render(Closed(.01f));
var open = Render(Closed(1f));
// constriction_diameter .001 at tract index 4 shuts the oral ladder; the only way out is the junction into the nasal ladder. A shut
// velum is opening 0: the port area is linear in the opening, and a flow source into a nearly lossless closed tract leaves through even
// a 0.008 cm^2 port (opening .01 measured rms 1.0e-2 against 2.1e-2 open).
var shut = Render(MouthShutScript(0f));
var open = Render(MouthShutScript(1f));
if (shut is null || open is null) return;
output.WriteLine($"mouth shut: velum closed rms {Rms(shut):E2}, velum open rms {Rms(open):E2}; open peaks {string.Join(" ", Peaks(open).Select(p => p.ToString("F0")))}");
Assert.True(Rms(open) > 5 * Rms(shut), $"nasal output rms {Rms(open):E2} vs velum-closed {Rms(shut):E2}");
output.WriteLine($"mouth shut, above 600 Hz: velum closed rms {RmsAbove(shut, 600):E2}, velum open rms {RmsAbove(open, 600):E2}; open peaks {string.Join(" ", Peaks(open).Select(p => p.ToString("F0")))}");
Assert.True(RmsAbove(open, 600) > 5 * RmsAbove(shut, 600), $"nasal output rms {RmsAbove(open, 600):E2} vs velum-closed {RmsAbove(shut, 600):E2}");
}
// The authored vowel gestures (PatchScript's VowelShapes, reached through the IPA word script) rendered end to end through the
// The authored vowel gestures (the Story vowel weights, reached through the IPA word script) rendered end to end through the
// exporter, merged splines and the ladder; the source is a bright 10 Hz saw so the tract's response is sampled densely.
private float[]? GestureWindow(string ipa)
{
@ -270,17 +328,17 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
return peaks.ToArray();
}
// Measured through the exported ladder on 2026-11-08 (impulse peaks i 271/2313, u 292/876, a 703/1084; the textbook values are
// i 270/2290, u 300/870, a 730/1090). The envelope of the 10 Hz render reads i 270/2380 and a 690/1060.
// Tube-chain F1/F2 of the Story 2008 area functions at 17.52 cm (tools/formant-check/tube_formants.py, an independent frequency-domain
// transfer-matrix estimate: i 251/2063, ɑ 650/913, u 271/761). The envelope of the 10 Hz render is read against those with 10% bands.
public static TheoryData<string, double, double> GestureVowels => new()
{
{ "i", 271, 2313 },
{ "a", 703, 1084 },
{ "i", 251, 2063 },
{ "a", 650, 913 },
};
[Theory]
[MemberData(nameof(GestureVowels))]
public void AuthoredVowelGesturesPlaceF1AndF2NearTheTextbook(string ipa, double f1, double f2)
public void AuthoredVowelGesturesPlaceF1AndF2NearTheStoryTubeChain(string ipa, double f1, double f2)
{
var peaks = GesturePeaks(ipa);
if (peaks is null) return;
@ -290,14 +348,14 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
}
[Fact]
public void AuthoredRoundedBackVowelPlacesF1AndF2NearTheTextbook()
public void AuthoredRoundedBackVowelPlacesF1AndF2NearTheStoryTubeChain()
{
// With a radiating highpass that rose with aperture (540 Hz at aperture 1) the envelope read u's F1 at 450 Hz.
var peaks = GesturePeaks("u");
if (peaks is null) return;
Assert.InRange(peaks[0], 292 * .85, 350);
Assert.InRange(peaks[1], 876 * .9, 876 * 1.1);
Assert.InRange(peaks[0], 271 * .85, 350);
Assert.InRange(peaks[1], 761 * .9, 761 * 1.1);
}
[Fact]
@ -309,7 +367,7 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
if (i is null || u is null || a is null) return;
Assert.True(a[0] > i[0] + 250, $"F1 i {i[0]} a {a[0]}");
// F2: i 2380 > a 1060 > u 880.
// F2 (tube chain): i 2063 > a 913 > u 761.
Assert.True(i[1] > a[1] + 800 && a[1] > u[1] + 100, $"F2 i {i[1]} a {a[1]} u {u[1]}");
}
@ -317,7 +375,7 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
public void RoundedBackVowelIsAudibleNextToTheCloseFrontVowel()
{
// The rendered /ˈnibu/ had /u/ 29 dB under /i/: the gesture scaled the radiating aperture (an output gain) by the lip area
// the ladder already carries. Measured here: -4.3 dB; with the old aperture scale -22.7 dB.
// the ladder already carries. Measured with the earlier grid-searched shapes: -4.3 dB; with the old aperture scale -22.7 dB.
var i = GestureWindow("i");
var u = GestureWindow("u");
if (i is null || u is null) return;
@ -330,4 +388,206 @@ public sealed class VocalTractFormantTests(ITestOutputHelper output)
Assert.InRange(gap, -10, 3);
Assert.True(uPeaks[1] < 1000, $"u F2 {uPeaks[1]}");
}
// ---- Story 2008 measured vowel shapes, set through the vowel_* weight surfaces ----
private static string StoryScript(IReadOnlyDictionary<string, float> weights, float lipOpening = 4f) => $$"""
morphology name=oral length_cm=17.52 diameters=.55,.75,1.05,1.45,1.25,.8 tongue_index=3 tongue_diameter=1.4 constriction_index=4 constriction_diameter=4 lip_opening={{lipOpening}}
vowel_set area_function=oral source=story2008
morphology name=nasal length_cm=12 diameters=.04,.28,.52,.72
waveguide_path name=oral_path morphology=oral speed=350 loss=.998
waveguide_path name=nasal_path morphology=nasal speed=350 loss=.997
source_port name=folds path=oral_path pressure=.7 tension=.55 opening=.5 noise=0 impedance=.3 flow_scale=.12 tilt=0
branch_port name=velopharynx from=oral_path from_position=.45 to=nasal_path opening=0 coupling=1
radiation_load name=mouth path=oral_path aperture=1 reflection=-1 impedance=.28
radiation_load name=nostril path=nasal_path aperture=.8 reflection=-1 impedance=.28
vocal_network name=voice paths=oral_path,nasal_path sources=folds branches=velopharynx radiation=mouth,nostril
{{string.Join("\n", weights.Select(w => $"control_spline name=w_{w.Key} surface=/vocal/areas/0/area/vowel_{w.Key} interp=linear points=0:{w.Value.ToString(System.Globalization.CultureInfo.InvariantCulture)},4:{w.Value.ToString(System.Globalization.CultureInfo.InvariantCulture)}"))}}
vocal network=voice freq={{Pitch}} min_freq=5 gain=.5 sustain=1.5
""";
private double[]? StoryPeaks(string id, float lipOpening = 4f)
{
var window = Render(StoryScript(new Dictionary<string, float> { [id] = 1f }, lipOpening));
if (window is null) return null;
var peaks = Peaks(window);
output.WriteLine($"story {id}: peaks " + string.Join(" ", peaks.Select(p => p.ToString("F0"))) + " lpc " + string.Join(" ", LpcPeaks(window).Select(p => p.ToString("F0"))));
Assert.True(peaks.Length >= 2, $"{id}: fewer than two spectral peaks");
return peaks;
}
// F1/F2 (F3) of the Story 2008 area functions at 17.52 cm from tools/formant-check/tube_formants.py (an independent
// transfer-matrix estimate in the frequency domain): python3 tools/formant-check/tube_formants.py docs/research/data/story2008-area-functions.csv --length-cm 17.52
//
// Bands: the front vowels (i ih e eh ae) are held to 6% on F1 and F2. The back and rounded vowels (uh ah aw o oo_short u) are held to 10%:
// the measured disagreement between two lip radiation models (the ladder's one-pole Levine-Schwinger-form reflectance with g = 1, the tube
// chain's baffled piston R||L) plus the ladder's cell quantisation. The nasal junction used to take the reflection coefficients of the
// interface one cell further on (sum over cells 11 and 12 instead of 10 and 11), which put the wrong step in the pharynx and read the
// back vowels' F1 12 to 26% high (uh +14.8%, ah +10.8%, aw +12.4%, o +25.8%, oo_short +13.0% against this table, five tests skipped on those
// numbers); the ladder now agrees with the numpy copy (tools/formant-check/ladder_formants.py --filtered --lip-reflection -1) to 10 Hz.
public static TheoryData<string, double, double, double> StoryTubeChain => new()
{
{ "i", 251, 2063, .06 }, { "ih", 346, 1987, .06 }, { "e", 432, 1827, .06 }, { "eh", 527, 1571, .06 }, { "ae", 729, 1523, .06 },
{ "uh", 549, 809, .10 }, { "ah", 650, 913, .10 }, { "aw", 516, 692, .10 }, { "o", 477, 787, .10 }, { "oo_short", 478, 850, .10 },
{ "u", 271, 761, .10 },
};
[Theory]
[MemberData(nameof(StoryTubeChain))]
public void LadderReproducesTheTubeChainFormantsOfEachStoryVowel(string id, double f1, double f2, double band)
{
var peaks = StoryPeaks(id);
if (peaks is null) return;
Assert.InRange(peaks[0], f1 * (1 - band), f1 * (1 + band));
Assert.InRange(peaks[1], f2 * (1 - band), f2 * (1 + band));
}
[Fact]
public void NearCloseVowelsSitApartFromTheCloseVowelsInF1()
{
var i = StoryPeaks("i");
var ih = StoryPeaks("ih");
var u = StoryPeaks("u");
var oo = StoryPeaks("oo_short");
if (i is null || ih is null || u is null || oo is null) return;
Assert.True(ih[0] >= i[0] + 60, $"F1 i {i[0]} ih {ih[0]}");
Assert.True(oo[0] >= u[0] + 60, $"F1 u {u[0]} oo_short {oo[0]}");
}
[Fact]
public void BlendingTwoVowelsHalfAndHalfIsTheCellwiseGeometricMeanOfTheirAreas()
{
var blend = Render(StoryScript(new Dictionary<string, float> { ["i"] = .5f, ["u"] = .5f }));
if (blend is null) return;
// Reference: one custom shape on the 24-cell grid (round(17.52e-2 * 48000 / 350)), the geometric mean of i and u at each cell
// centre, selected by its default weight 0 (a lone shape at zero weight is itself). The blend must render the same lines.
const int Cells = 24;
var i = StoryVowelAreas.All.Single(v => v.Id == "i").Areas;
var u = StoryVowelAreas.All.Single(v => v.Id == "u").Areas;
var mean = Enumerable.Range(0, Cells).Select(c => (float)Math.Sqrt(StoryVowelAreas.AreaAt(i, (c + .5f) / Cells) * (double)StoryVowelAreas.AreaAt(u, (c + .5f) / Cells))).ToArray();
var reference = Render(StoryScript(new Dictionary<string, float>()), patch => patch with
{
AreaFunctions = patch.AreaFunctions.Select(a => a.Name == "oral" ? a with { Vowels = [new VowelAreaShape("mean", mean)] } : a).ToArray()
});
if (reference is null) return;
var blendPeaks = Peaks(blend);
var referencePeaks = Peaks(reference);
output.WriteLine($"blend peaks {string.Join(" ", blendPeaks.Select(p => p.ToString("F0")))}; reference {string.Join(" ", referencePeaks.Select(p => p.ToString("F0")))}");
Assert.Equal(referencePeaks.Take(3), blendPeaks.Take(3));
var iPeaks = StoryPeaks("i")!;
Assert.True(blendPeaks[1] < iPeaks[1] - 100, $"blend F2 {blendPeaks[1]} against i {iPeaks[1]}");
}
[Fact]
public void WeightsThatSumToZeroGiveTheGeometricMeanShapeNotNaN()
{
var zero = Render(StoryScript(new Dictionary<string, float>()));
if (zero is null) return;
Assert.All(zero, sample => Assert.True(float.IsFinite(sample)));
Assert.True(Rms(zero) > 1e-4, $"silent: {Rms(zero)}");
// Independent route to the same shape: all eleven weights at 1 take the normalised-blend branch, whose result is the geometric mean.
var reference = Render(StoryScript(StoryVowelAreas.All.ToDictionary(v => v.Id, _ => 1f)));
if (reference is null) return;
Assert.Equal(Peaks(reference).Take(3), Peaks(zero).Take(3));
}
[Fact]
public void LipOpeningCapsTheLipCellOfAStoryVowel()
{
var open = StoryPeaks("ah", 4f);
var capped = StoryPeaks("ah", .4f);
if (open is null || capped is null) return;
output.WriteLine($"ah lip cap 4 cm F1 {open[0]}; 0.4 cm F1 {capped[0]}");
Assert.True(capped[0] < open[0] - 80, $"lip cap .4 cm F1 {capped[0]} vs open {open[0]}");
}
// ---- Nasal coupling: the murmur of /m/ and /n/ against the nostril transfer of tools/formant-check/tube_formants.py --nasal ----
// The independent tube chain (Story 2008 geometric-mean oral shape, port after section 22 of 44, A_vp 0.56 cm^2 = the gesture's 0.7 opening,
// the 19-section nasal table, rigid closed oral end after section 44 for /m/ and section 37 for /n/, 17.52 cm): N1 325 / 352 Hz, transfer
// zeros (short circuit of the oral side branch) 882 / 1136 Hz; the quarter-wave c / 4L of the uniform branch would say 999 / 1465 Hz.
private const double TubeM1 = 325, TubeN1 = 352, TubeMZero = 882, TubeNZero = 1136;
// Smoothed line spectrum in dB (flattened by the saw's 1/f) of a steady 10 Hz render, indexed by line number k (k * Pitch Hz).
private static double[] SmoothedDb(float[] window)
{
var lines = Lines(window);
var db = lines.Select(line => 20 * Math.Log10(Math.Max(1e-12, line))).ToArray();
return db.Select((_, k) => k < 2 || k > db.Length - 3 ? db[k] : db.Skip(k - 2).Take(5).Average()).ToArray();
}
// First dip above the first peak: the first local minimum of the smoothed spectrum, 200-2000 Hz, at least 3 dB under the lower of the
// maxima on either side of it.
private static double FirstDip(double[] db)
{
for (var k = (int)(200 / Pitch); k < (int)(2000 / Pitch); k++)
{
if (!(db[k] < db[k - 1] && db[k] <= db[k + 1])) continue;
var left = db[k]; for (var j = k; j > 0 && db[j] >= left; j--) left = db[j];
var right = db[k]; for (var j = k; j < db.Length - 1 && db[j] >= right; j++) right = db[j];
if (Math.Min(left, right) - db[k] >= 3) return k * Pitch;
}
return double.NaN;
}
private (double N1, double Dip)? Murmur(string ipa)
{
var window = GestureWindow(ipa);
if (window is null) return null;
var n1 = Peaks(window)[0];
var dip = FirstDip(SmoothedDb(window));
output.WriteLine($"/{ipa}/ murmur: N1 {n1:F0} Hz, first dip {dip:F0} Hz");
return (n1, dip);
}
[Fact]
public void TheMurmurOfMAndNHasItsFirstPeakBelow400HzNearTheTubeChain()
{
var m = Murmur("m");
var n = Murmur("n");
if (m is null || n is null) return;
Assert.InRange(m.Value.N1, 250, 400);
Assert.InRange(n.Value.N1, 250, 400);
Assert.InRange(m.Value.N1, TubeM1 * .85, TubeM1 * 1.15);
Assert.InRange(n.Value.N1, TubeN1 * .85, TubeN1 * 1.15);
}
[Fact]
public void TheFirstZeroOfTheNasalMurmurSitsWhereTheOralSideBranchPutsItAndRisesWhenTheBranchIsShorter()
{
var m = Murmur("m");
var n = Murmur("n");
if (m is null || n is null) return;
Assert.InRange(m.Value.Dip, TubeMZero * .85, TubeMZero * 1.15);
Assert.InRange(n.Value.Dip, TubeNZero * .85, TubeNZero * 1.15);
// /n/ closes the mouth 2.6 cm from the lips, /m/ at the lips: a shorter branch beyond the port, a higher zero.
Assert.True(n.Value.Dip >= 1.15 * m.Value.Dip, $"/n/ dip {n.Value.Dip} Hz against /m/ dip {m.Value.Dip} Hz");
}
[Fact]
public void AnOralClosureWithTheVelumShutLeavesTheMurmurSilentAndOpeningTheVelumLetsItOut()
{
var open = GestureWindow("m");
if (open is null) return;
// The same /m/ with the velum spline held shut: the gesture's own velum control is the only door, so remove it.
var shut = Render(IpaGestureExperiment.BuildWordScript([(new IpaGestureExperimentTarget("v_m", "m", "", 1.5f), 1f)], Pitch, 0.1f)
.Replace("noise=.04", "noise=0 tilt=1", StringComparison.Ordinal)
.Replace($"freq={Pitch}", $"freq={Pitch} min_freq=5", StringComparison.Ordinal),
patch => patch with { ControlSplines = patch.ControlSplines.Where(spline => spline.SurfacePath != "/vocal/branches/0/opening").ToArray() });
if (shut is null) return;
output.WriteLine($"/m/ above 600 Hz: velum open rms {RmsAbove(open, 600):E2}, velum shut rms {RmsAbove(shut, 600):E2}");
Assert.True(RmsAbove(open, 600) > 10 * RmsAbove(shut, 600), $"open {RmsAbove(open, 600):E2} shut {RmsAbove(shut, 600):E2}");
}
}

81
tests/AquaSynth.Dsl.Tests/WeksaPhonesTests.cs
View file

@ -61,36 +61,81 @@ public sealed class WeksaPhonesTests
[Fact]
public void ReportsPhonesTheTractCannotMakeAndRendersNothing()
{
var word = WeksaPhonesReader.Parse(Line($"[{Ph("n", "none")},{Ph("ɪ", "primary")},{Ph("b", "none")},{Ph("ʊ", "none")}]"));
var word = WeksaPhonesReader.Parse(Line($"[{Ph("n", "none")},{Ph("ŋ", "primary")},{Ph("b", "none")},{Ph("θ", "none")}]"));
var plan = WeksaWordPlan.From(word, "x");
Assert.False(plan.Renderable);
Assert.Equal(["ɪ", "ʊ"], plan.Unsupported.Select(p => p.Ipa));
Assert.Equal(["ŋ", "θ"], plan.Unsupported.Select(p => p.Ipa));
Assert.Empty(plan.Phones);
}
[Fact]
public void StressMovesOnlyTheIntensityOfTheMarkedPhone()
public void NearCloseVowelsRenderAsVowelsWithTheirOwnStoryShape()
{
string Script(WeksaStress onI) => IpaGestureExperiment.BuildWordScript(
WeksaWordPlan.From(new WeksaPhones("en-US", "ˈnibu",
[new("n", WeksaStress.None), new("i", onI), new("b", WeksaStress.None), new("u", WeksaStress.None)]), "w").Phones.ToArray());
var word = WeksaPhonesReader.Parse(Line($"[{Ph("b", "none")},{Ph("ɪ", "primary")},{Ph("n", "none")},{Ph("b", "none")},{Ph("ʊ", "none")}]"));
var plan = WeksaWordPlan.From(word, "x");
Assert.True(plan.Renderable);
Assert.Equal([WeksaWordPlan.ConsonantSeconds, WeksaWordPlan.VowelSeconds * WeksaWordPlan.CueFor(WeksaStress.Primary).Duration, WeksaWordPlan.ConsonantSeconds, WeksaWordPlan.ConsonantSeconds, WeksaWordPlan.VowelSeconds],
plan.Phones.Select(p => p.Target.DurationSeconds));
var patch = PatchScript.Parse(IpaGestureExperiment.BuildWordScript(plan.Phones.ToArray(), 140));
var ih = patch.ControlSplines.Where(s => s.Name.StartsWith("x_1_ɪ_vowel_", StringComparison.Ordinal)).ToDictionary(s => s.SurfacePath, s => s.Points[0].Value);
Assert.Equal(1f, ih["/vocal/areas/0/area/vowel_ih"]);
Assert.Equal(0f, ih["/vocal/areas/0/area/vowel_i"]);
var oo = patch.ControlSplines.Where(s => s.Name.StartsWith("x_4_ʊ_vowel_", StringComparison.Ordinal)).ToDictionary(s => s.SurfacePath, s => s.Points[0].Value);
Assert.Equal(1f, oo["/vocal/areas/0/area/vowel_oo_short"]);
Assert.Equal(0f, oo["/vocal/areas/0/area/vowel_u"]);
}
var none = Script(WeksaStress.None).Split('\n');
var primary = Script(WeksaStress.Primary).Split('\n');
Assert.Equal(none.Length, primary.Length);
var differing = Enumerable.Range(0, none.Length).Where(i => none[i] != primary[i]).ToArray();
Assert.Single(differing);
Assert.Contains("ipa=i ", primary[differing[0]]);
var strip = (string l) => Regex.Replace(l, @"intensity=\S+", "");
Assert.Equal(strip(none[differing[0]]), strip(primary[differing[0]]));
Assert.Contains("intensity=1.2", primary[differing[0]]);
private static WeksaWordPlan Plan(WeksaStress onI) => WeksaWordPlan.From(new WeksaPhones("en-US", "ˈnibu",
[new("n", WeksaStress.None), new("i", onI), new("b", WeksaStress.None), new("u", WeksaStress.None)]), "w");
[Fact]
public void StressIsLongerLouderAndHigherOnTheMarkedPhoneOnly()
{
var none = Plan(WeksaStress.None).Phones;
var primary = Plan(WeksaStress.Primary).Phones;
var secondary = Plan(WeksaStress.Secondary).Phones;
foreach (var k in new[] { 0, 2, 3 })
{
Assert.Equal(none[k].Target, primary[k].Target);
Assert.Equal(none[k].Intensity, primary[k].Intensity);
}
Assert.Equal(WeksaWordPlan.VowelSeconds, none[1].Target.DurationSeconds);
Assert.Equal(1f, none[1].Target.PitchRatio);
Assert.Equal(1f, none[1].Intensity);
Assert.Equal(WeksaWordPlan.VowelSeconds * 220f / 160f, primary[1].Target.DurationSeconds, 5);
Assert.Equal(1.06f, primary[1].Target.PitchRatio);
Assert.Equal(1.30f, primary[1].Intensity);
// Secondary sits strictly between unstressed and primary on every cue.
Assert.InRange(secondary[1].Target.DurationSeconds, none[1].Target.DurationSeconds + 1e-4f, primary[1].Target.DurationSeconds - 1e-4f);
Assert.InRange(secondary[1].Target.PitchRatio, 1.001f, 1.059f);
Assert.InRange(secondary[1].Intensity, 1.001f, 1.299f);
}
[Fact]
public void StressReachesTheScriptAsDurationF0RatioAndSourceLevelAndShiftsLaterPhones()
{
var patch = (WeksaStress s) => PatchScript.Parse(IpaGestureExperiment.BuildWordScript(Plan(s).Phones.ToArray(), 140));
var none = patch(WeksaStress.None);
var primary = patch(WeksaStress.Primary);
float First(SynthPatch p, string prefix, string surface) =>
p.ControlSplines.Single(sp => sp.Name.StartsWith(prefix, StringComparison.Ordinal) && sp.SurfacePath == surface).Points[0].Value;
// pitch_ratio surface range 0.5..2: normalized (1.06 - .5) / 1.5; unstressed holds 1 -> (1 - .5) / 1.5.
Assert.Equal((1.06f - .5f) / 1.5f, First(primary, "w_1_i_source_pitch", "/vocal/sources/0/pitch_ratio"), 4);
Assert.Equal((1f - .5f) / 1.5f, First(none, "w_1_i_source_pitch", "/vocal/sources/0/pitch_ratio"), 4);
Assert.Equal(First(none, "w_1_i_source_level", "/vocal/sources/0/flow_scale") * 1.30f, First(primary, "w_1_i_source_level", "/vocal/sources/0/flow_scale"), 4);
// The b after the stressed vowel starts later by the vowel's extra length.
float Start(SynthPatch p) => p.PhonemeGestures.Single(g => g.Ipa == "b").StartSeconds;
Assert.Equal(WeksaWordPlan.VowelSeconds * (220f / 160f - 1f), Start(primary) - Start(none), 4);
}
[Fact]
public void PhonesTheWordPlanAcceptsAreExactlyTheDslTable()
{
foreach (var ipa in new[] { "n", "i", "b", "u", "a", "ə" })
Assert.True(PatchScript.TryIpaDescriptor(ipa, out _));
Assert.False(PatchScript.TryIpaDescriptor("ɪ", out _));
foreach (var ipa in new[] { "n", "i", "b", "u", "a", "ə", "ɪ", "ʊ", "e", "ɛ", "æ", "ʌ", "ɑ", "ɔ", "o" })
Assert.True(PatchScript.TryIpaDescriptor(ipa, out var descriptor) && descriptor.EndsWith("_vowel", StringComparison.Ordinal) == "iuaəɪʊeɛæʌɑɔo".Contains(ipa));
foreach (var ipa in new[] { "ŋ", "θ", "ɹ", "l" })
Assert.False(PatchScript.TryIpaDescriptor(ipa, out _));
}
}

73
tools/formant-check/burg_formants.py Normal file
View file

@ -0,0 +1,73 @@
#!/usr/bin/env python3
"""Independent formant estimator for AquaSynth renders. Shares no code with the C# tests or with the ladder.
Burg linear prediction (not autocorrelation LPC, not the tract's own impulse response): stable on short windows and
unbiased by the window edge. Usage:
burg_formants.py FILE.wav T0 T1 [--order N] [--frame-ms 30] [--json]
Resamples to 10 kHz (FFT-domain), pre-emphasises 0.97, Hamming 30 ms frames hop 10 ms in [T0,T1], Burg order 12 default,
poles with bandwidth < 400 Hz and 150 Hz < f < 4.5 kHz are formants. Prints the per-formant median over frames and the
interquartile spread. It measures; it does not listen.
"""
import sys, json, wave, argparse
import numpy as np
def read_wav(path):
with wave.open(path) as w:
sr, n, ch, sw = w.getframerate(), w.getnframes(), w.getnchannels(), w.getsampwidth()
raw = w.readframes(n)
assert sw == 2, "16-bit PCM only"
x = np.frombuffer(raw, dtype=np.int16).astype(np.float64) / 32768.0
return (x.reshape(-1, ch).mean(axis=1) if ch > 1 else x), sr
def resample(x, sr, target):
if sr == target:
return x
n = int(round(len(x) * target / sr))
X = np.fft.rfft(x)
m = n // 2 + 1
Y = np.zeros(m, dtype=complex)
k = min(m, len(X))
Y[:k] = X[:k]
return np.fft.irfft(Y, n) * (n / len(x))
def burg(x, order):
f = x.copy(); b = x.copy(); a = np.array([1.0])
for _ in range(order):
fp, bp = f[1:], b[:-1]
k = -2.0 * np.dot(fp, bp) / (np.dot(fp, fp) + np.dot(bp, bp) + 1e-30)
a = np.concatenate([a, [0.0]]); a = a + k * a[::-1]
f, b = fp + k * bp, bp + k * fp
return a
def formants(frame, fs, order):
a = burg(frame * np.hamming(len(frame)), order)
r = np.roots(a); r = r[np.imag(r) > 1e-4]
fq = np.angle(r) * fs / (2 * np.pi); bw = -np.log(np.abs(r)) * fs / np.pi
return sorted(float(f) for f, w in zip(fq, bw) if 150 < f < 4500 and w < 400)
def main():
p = argparse.ArgumentParser()
p.add_argument("wav"); p.add_argument("t0", type=float); p.add_argument("t1", type=float)
p.add_argument("--order", type=int, default=12); p.add_argument("--frame-ms", type=float, default=30)
p.add_argument("--json", action="store_true")
a = p.parse_args()
x, sr = read_wav(a.wav)
fs = 10000
y = resample(x, sr, fs)
y = np.append(y[0], y[1:] - 0.97 * y[:-1])
n = int(a.frame_ms * fs / 1000); hop = int(0.010 * fs)
rows = []
for s in range(int(a.t0 * fs), int(a.t1 * fs) - n, hop):
fr = y[s:s + n]
if fr.std() < 1e-6: continue
rows.append(formants(fr, fs, a.order)[:4])
out = {"file": a.wav, "t0": a.t0, "t1": a.t1, "order": a.order, "frames": len(rows)}
for i in range(4):
v = [r[i] for r in rows if len(r) > i]
out[f"F{i+1}"] = None if len(v) < max(3, len(rows) // 2) else {"median": round(float(np.median(v)), 1), "q1": round(float(np.percentile(v, 25)), 1), "q3": round(float(np.percentile(v, 75)), 1), "n": len(v)}
if a.json: print(json.dumps(out))
else:
print(f"{a.wav} [{a.t0},{a.t1}] order {a.order} frames {len(rows)}: " + " ".join(f"F{i+1}={out[f'F{i+1}']['median']} ({out[f'F{i+1}']['q1']}-{out[f'F{i+1}']['q3']})" if out[f"F{i+1}"] else f"F{i+1}=-" for i in range(4)))
if __name__ == "__main__":
main()

56
tools/formant-check/ladder_formants.py Normal file
View file

@ -0,0 +1,56 @@
#!/usr/bin/env python3
"""Impulse-response formants of the AquaSynth Kelly-Lochbaum ladder, in numpy, for the Story 2008 area functions. It copies the update
equations of FaustExport.VocalLadder.cs (glottal reflection .75, per-cell loss, one cell per sample, the lip reflection a constant) so
that the ladder's formants can be compared with tube_formants.py (frequency domain, frequency-dependent radiation) without Faust.
ladder_formants.py CSV [--cells 24] [--lip-reflection -0.82] [--loss 0.998] [--vowel V] [--filtered]
Prints the first three spectral peaks above 150 Hz at 48 kHz. Unheard: this says where resonances sit, not how a vowel sounds.
"""
import argparse, csv
import numpy as np
def load(path):
rows = [l for l in open(path, encoding='utf-8') if not l.startswith('#')]
r = list(csv.reader(rows))
return {h: np.array([float(x[j + 1]) for x in r[1:] if x[0].isdigit()]) for j, h in enumerate(r[0][1:])}
def at(a, p): # linear in area between section centres, as StoryVowelAreas.AreaAt
x = np.clip(p * len(a) - 0.5, 0, len(a) - 1); l = int(np.floor(x)); r = min(len(a) - 1, l + 1)
return a[l] + (a[r] - a[l]) * (x - l)
def lip_pole(area_cm2, fs=48000.0, c=35000.0):
"""One-pole coefficient b of the open-end reflectance -g(1-b)/(1-b z^-1): b/(1-b)^2 = (a fs / c)^2 for the lip radius a
(Levine-Schwinger low-ka loss (ka)^2 = 1-|R|^2), clamped to [0, .9] as the Faust emission does."""
a = np.sqrt(max(area_cm2, 1e-6) / np.pi); q = (a * fs / c) ** 2
if q < 1e-6: return 0.0
return float(np.clip(((2 * q + 1) - np.sqrt(4 * q + 1)) / (2 * q), 0.0, 0.9))
def impulse(areas, lip, loss, glot=0.75, n=1 << 14, filtered=False):
N = len(areas)
k = [0] + [(areas[i - 1] - areas[i]) / (areas[i - 1] + areas[i]) for i in range(1, N)]
R = np.zeros(N); L = np.zeros(N); out = np.zeros(n)
b = lip_pole(areas[-1]) if filtered else 0.0
for t in range(n):
Rn = np.zeros(N); Ln = np.zeros(N)
Rn[0] = (L[0] * glot + (1.0 if t == 0 else 0.0)) * loss
for i in range(1, N): Rn[i] = (R[i - 1] - k[i] * (R[i - 1] + L[i])) * loss
for i in range(N - 1): Ln[i] = (L[i + 1] + k[i + 1] * (R[i] + L[i + 1])) * loss
Ln[N - 1] = (b * L[N - 1] + (1 - b) * lip * R[N - 1]) * loss
R, L = Rn, Ln; out[t] = R[N - 1]
return out
def peaks(out, fs=48000, n=3):
S = np.abs(np.fft.rfft(out)); f = np.fft.rfftfreq(len(out), 1 / fs)
idx = [i for i in range(1, len(S) - 1) if S[i] > S[i - 1] and S[i] >= S[i + 1] and f[i] > 150 and S[i] > 0.05 * S.max()]
return [round(f[i]) for i in idx][:n]
if __name__ == '__main__':
p = argparse.ArgumentParser()
p.add_argument('csv'); p.add_argument('--cells', type=int, default=24)
p.add_argument('--lip-reflection', type=float, default=-0.82); p.add_argument('--loss', type=float, default=0.998)
p.add_argument('--vowel'); p.add_argument('--filtered', action='store_true', help='lip reflectance is the one-pole open-end filter with DC gain |--lip-reflection|, pole from the lip area')
a = p.parse_args()
for v, areas in load(a.csv).items():
if a.vowel and v != a.vowel: continue
cells = np.array([at(areas, (i + .5) / a.cells) for i in range(a.cells)])
print(v, peaks(impulse(cells, a.lip_reflection, a.loss, filtered=a.filtered)))

66
tools/formant-check/spectrum_dips.py Normal file
View file

@ -0,0 +1,66 @@
#!/usr/bin/env python3
"""Averaged-spectrum peaks and dips of a stretch of a WAV, independent of the AquaSynth code (numpy only).
spectrum_dips.py WAV T0 T1 [--f0 HZ] [--fmax 4000] [--nperseg 1024]
Welch estimate (Hann, 75 % overlap, zero-padded to 4x) of the mono signal between T0 and T1 seconds. Prints local maxima and minima in dB (relative
to the largest) below --fmax. With --f0 the spectrum of a periodic source is a line spectrum, so the extrema are read from the harmonic envelope
instead: the largest bin within a quarter of f0 of each harmonic k*f0; frequencies are then multiples of f0 and the resolution is f0.
Used on the nasal murmur of a WordRender render against the nostril transfer of tube_formants.py --nasal.
"""
import sys, argparse, wave
import numpy as np
def read_wav(path):
with wave.open(path) as w:
n, ch, width, rate = w.getnframes(), w.getnchannels(), w.getsampwidth(), w.getframerate()
raw = w.readframes(n)
if width != 2: raise SystemExit(f"{path}: only 16-bit PCM")
x = np.frombuffer(raw, dtype='<i2').astype(float) / 32768.0
return x.reshape(-1, ch)[:, 0], rate
def welch(x, rate, nperseg):
nperseg = min(nperseg, len(x))
win = np.hanning(nperseg)
step = max(1, nperseg // 4)
nfft = 4 * nperseg
acc = np.zeros(nfft // 2 + 1); count = 0
for s in range(0, len(x) - nperseg + 1, step):
seg = x[s:s + nperseg] * win
acc += np.abs(np.fft.rfft(seg, nfft)) ** 2; count += 1
return np.fft.rfftfreq(nfft, 1 / rate), acc / max(count, 1), count
def harmonic_envelope(freqs, power, f0, fmax):
out = []
k = 1
while k * f0 < fmax:
sel = (freqs >= k * f0 - f0 / 4) & (freqs <= k * f0 + f0 / 4)
out.append((k * f0, power[sel].max()))
k += 1
return np.array([o[0] for o in out]), np.array([o[1] for o in out])
def extrema(freqs, db):
res = []
for i in range(1, len(db) - 1):
if db[i] > db[i - 1] and db[i] >= db[i + 1]: res.append((freqs[i], db[i], 'peak'))
elif db[i] < db[i - 1] and db[i] <= db[i + 1]: res.append((freqs[i], db[i], 'dip'))
return res
def main():
p = argparse.ArgumentParser()
p.add_argument('wav'); p.add_argument('t0', type=float); p.add_argument('t1', type=float)
p.add_argument('--f0', type=float); p.add_argument('--fmax', type=float, default=4000); p.add_argument('--nperseg', type=int, default=1024)
a = p.parse_args()
x, rate = read_wav(a.wav)
seg = x[int(a.t0 * rate):int(a.t1 * rate)]
freqs, power, count = welch(seg, rate, a.nperseg)
keep = (freqs > 30) & (freqs < a.fmax)
freqs, power = freqs[keep], power[keep]
if a.f0: freqs, power = harmonic_envelope(freqs, power, a.f0, a.fmax)
db = 10 * np.log10(np.maximum(power, 1e-30) / power.max())
ex = extrema(freqs, db)
print(f"{a.wav} {a.t0}-{a.t1}s rate {rate} Welch segments {count}" + (f" harmonic envelope f0 {a.f0}" if a.f0 else ""))
print(" peaks: " + " ".join(f"{f:.0f}({d:.0f}dB)" for f, d, k in ex if k == 'peak'))
print(" dips: " + " ".join(f"{f:.0f}({d:.0f}dB)" for f, d, k in ex if k == 'dip'))
if __name__ == '__main__':
main()

186
tools/formant-check/tube_formants.py Normal file
View file

@ -0,0 +1,186 @@
#!/usr/bin/env python3
"""Frequency-domain formants of a published area function. Independent of the AquaSynth ladder (no time-domain scattering,
no Faust): a chain of lossy uniform tubes (ABCD matrices) between a high-impedance glottis and a piston-in-baffle radiation load.
tube_formants.py CSV [--vowel V] [--length-cm L] [--c 350] [--fmax 5000]
CSV is docs/research/data/story2008-area-functions.csv. --length-cm rescales every vowel to one tract length (the ladder
runs one length for all vowels); default is each vowel's own VTL. Prints F1..F4 as peaks of |U_lips / U_glottis|.
tube_formants.py CSV --nasal [--vowel V|mean] [--closure-section N|none] [--a-vp 0.56] [--port-section 22] [--ladder-closure-cm CM [--lip-cap-cm D]] [--length-cm 17.52]
Nasal mode: |U_nostril / U_glottis| of the oral tract (Story area function, `mean` = the cellwise geometric mean of the eleven) with the
velopharyngeal port at section --port-section of 44 (Story 2005 tier IV), port area --a-vp (cm^2) as the first nasal section, the 19-section
nasal table after it (Birkholz 2005 Table 2.3 after Dang & Honda 1994, 0.6 cm sections), the nostril as a piston-in-baffle load, and the oral
tract beyond the port either closed rigid at the end of section --closure-section (the oral stop of /m n/) or open at the lips (nasalised vowel).
Prints N1 (lowest peak), the local minima with their depth, and the oral side-branch quarter-wave prediction c / (4 L_branch).
"""
import sys, argparse
import numpy as np
RHO = 1.14e-3 # g/cm^3 as in VocalTractLab's dissertation constants (SI 1.14 kg/m^3)
def load(path):
rows = [l.strip() for l in open(path) if l.strip() and not l.startswith('#')]
head = rows[0].split(',')[1:]
areas = {v: [] for v in head}
extra = {}
for r in rows[1:]:
c = r.split(',')
if c[0].isdigit():
for v, x in zip(head, c[1:]): areas[v].append(float(x))
else:
extra[c[0]] = dict(zip(head, map(float, c[1:])))
return head, {v: np.array(a) for v, a in areas.items()}, extra
def alpha_of(f, wall=0.0):
# attenuation per cm: viscous/thermal ~ f^0.5 plus wall term, giving F1 bandwidths of order 50-100 Hz
return 4e-5 * np.sqrt(f) + wall + 1.5e-3
def abcd(areas, section_cm, c, f, wall=0.0):
"""Chain matrix of lossy uniform tubes, (P, U)_in = M (P, U)_out, cgs."""
w = 2 * np.pi * f
alpha = alpha_of(f, wall)
k = w / c * (1 - 1j * alpha * c / w) if w > 0 else 0
M = np.eye(2, dtype=complex)
for A in areas:
Z = RHO * c / A
kl = k * section_cm
M = M @ np.array([[np.cos(kl), 1j * Z * np.sin(kl)], [1j * np.sin(kl) / Z, np.cos(kl)]])
return M
def radiation(A_l, c, f):
"""Low-ka piston-in-baffle radiation impedance (Flanagan R||L), matched impedance above ka 1.2."""
w = 2 * np.pi * f
ka = w / c * np.sqrt(A_l / np.pi)
return (RHO * c / A_l) * (ka ** 2 / 2 + 1j * 8 * ka / (3 * np.pi)) if ka < 1.2 else (RHO * c / A_l) * (1 - 0j)
def transfer(areas, section_cm, c, freqs, wall=0.0):
"""|U_lip/U_glottis| for tubes of the given areas (cm^2) and one section length (cm). cgs units."""
out = np.zeros(len(freqs))
for n, f in enumerate(freqs):
M = abcd(areas, section_cm, c, f, wall)
Zr = radiation(areas[-1], c, f)
# U_glottis = 1 (open-circuit glottis); lip P_l = Zr*U_l; U_g = M10*P_l + M11*U_l = U_l*(M10*Zr + M11)
out[n] = 1 / abs(M[1, 0] * Zr + M[1, 1])
return out
# Nasal tract area function, cm^2, velopharyngeal end first (Birkholz 2005 Table 2.3 after Dang & Honda 1994), 0.6 cm sections.
NASAL_TABLE = [1.63, 2.07, 2.72, 3.59, 4.24, 3.26, 3.04, 3.04, 2.72, 2.5, 2.39, 2.39, 1.85, 0.76, 1.41, 1.74, 1.30, 1.74, 0.76]
NASAL_SECTION_CM = 0.6
def nasal_transfer(oral_areas, oral_section_cm, port_section, closure_section, a_vp, c, freqs):
"""(|U_nostril / U_glottis|, |oral branch admittance|). The transfer's zeros are the maxima of the second:
P_j = 0 where the mouth side of the junction is a short circuit. The oral tract is split after `port_section` sections: sections before it are the pharynx side, sections after
it the mouth side (closed rigid end after `closure_section` sections, or open at the lips when None). The nasal branch is the table with its
first section replaced by the port area a_vp (the ladder's first nasal cell is the port)."""
pharynx = oral_areas[:port_section]
end = len(oral_areas) if closure_section is None else closure_section
mouth = oral_areas[port_section:end]
nasal = [a_vp] + NASAL_TABLE[1:]
out = np.zeros(len(freqs))
yo = np.zeros(len(freqs))
for n, f in enumerate(freqs):
Mp = abcd(pharynx, oral_section_cm, c, f)
Mn = abcd(nasal, NASAL_SECTION_CM, c, f)
Zn = radiation(NASAL_TABLE[-1], c, f)
yn = (Mn[1, 0] * Zn + Mn[1, 1]) / (Mn[0, 0] * Zn + Mn[0, 1]) # nasal input admittance U/P
Mm = abcd(mouth, oral_section_cm, c, f)
if closure_section is None:
Zm = radiation(oral_areas[-1], c, f)
ym = (Mm[1, 0] * Zm + Mm[1, 1]) / (Mm[0, 0] * Zm + Mm[0, 1])
else:
ym = Mm[1, 0] / Mm[0, 0] # U_out = 0
yo[n] = abs(ym)
# glottis: U_g = M10 P_j + M11 U_1 = 1 with U_1 = (yn + ym) P_j
Pj = 1 / (Mp[1, 0] + Mp[1, 1] * (yn + ym))
out[n] = abs(Pj / (Mn[0, 0] * Zn + Mn[0, 1])) # U_nostril = P_j / (M00 Zn + M01)
return out, yo
def extrema(freqs, mag):
"""Local maxima and minima of mag with parabolic refinement in log magnitude: [(f, dB, 'peak'|'dip')]."""
res = []
for i in range(1, len(mag) - 1):
kind = 'peak' if mag[i] > mag[i - 1] and mag[i] >= mag[i + 1] else 'dip' if mag[i] < mag[i - 1] and mag[i] <= mag[i + 1] else None
if kind is None: continue
y0, y1, y2 = np.log(mag[i - 1]), np.log(mag[i]), np.log(mag[i + 1])
d = 0.5 * (y0 - y2) / (y0 - 2 * y1 + y2)
res.append((freqs[i] + d * (freqs[1] - freqs[0]), 20 * np.log10(mag[i]), kind))
return res
def ladder_closure(A, centre_fraction, constriction_cm=0.12, width=0.54, lip_cap_cm=None):
"""The AquaSynth ladder's own consonant shape on a 44-section area function, so the tube chain can be run on the constricted tract the ladder
is given (a Gaussian narrowing, not a rigid closed end): d' = min(d, c + max(0, d - c)(1 - w)), w = exp(-((5 p - 5 centre) / width)^2) with p the
section centre as a fraction of the tract (5 = sections - 1 of the 6-diameter morphology, width in those index units); lip_cap_cm caps the last section's diameter."""
n = len(A)
out = np.array(A, dtype=float)
for i in range(n):
p = (i + 0.5) / n
w = np.exp(-(((p - centre_fraction) * 5) / width) ** 2)
d = np.sqrt(out[i])
out[i] = min(d, constriction_cm + max(0.0, d - constriction_cm) * (1 - w)) ** 2
if lip_cap_cm is not None:
out[-1] = min(np.sqrt(out[-1]), lip_cap_cm) ** 2
return out
def nasal_report(args, head, areas, extra):
if args.vowel == 'mean':
A = np.exp(np.mean([np.log(areas[v]) for v in head], axis=0))
else:
A = areas[args.vowel]
L = args.length_cm or 17.52
sec = L / len(A)
if args.ladder_closure_cm is not None: # centre of the Gaussian narrowing, cm from the glottis
A = ladder_closure(A, args.ladder_closure_cm / L, constriction_cm=args.closure_diameter_cm, lip_cap_cm=args.lip_cap_cm)
closure = None if args.closure_section in (None, 'none') else int(args.closure_section)
freqs = np.arange(20, args.fmax, 2.0)
m, yo = nasal_transfer(A, sec, args.port_section, closure, args.a_vp, args.c * 100, freqs)
ex = extrema(freqs, m)
peaks_ = [e for e in ex if e[2] == 'peak']
# depth of a dip: below the smaller of its neighbouring peaks
lines = []
for j, e in enumerate(ex):
if e[2] != 'dip': continue
nb = [x[1] for x in (ex[j - 1] if j > 0 else None, ex[j + 1] if j + 1 < len(ex) else None) if x is not None and x[2] == 'peak']
lines.append(f"dip {e[0]:.0f} Hz ({e[1] - min(nb):+.1f} dB vs lower neighbour peak)" if nb else f"dip {e[0]:.0f} Hz")
zeros = [freqs[i] for i in range(1, len(yo) - 1) if yo[i] > yo[i - 1] and yo[i] >= yo[i + 1]]
where = 'open mouth' if closure is None else f"closed after section {closure}"
print(f"{args.vowel} port@{args.port_section}/{len(A)} A_vp={args.a_vp} L={L}cm {where}")
print(" peaks: " + " ".join(f"{e[0]:.0f}({e[1]:.0f}dB)" for e in peaks_[:5]))
print(" local minima (the first is the first null): " + "; ".join(lines[:5]))
print(" transfer zeros (oral branch short-circuit, maxima of |Y_oral|): " + " ".join(f"{z:.0f}" for z in zeros[:4]))
if closure is not None:
Lb = (closure - args.port_section) * sec
print(f" oral side branch {Lb:.2f} cm: quarter-wave {args.c * 100 / (4 * Lb):.0f} Hz, three-quarter {3 * args.c * 100 / (4 * Lb):.0f} Hz")
def peaks(freqs, mag, n=4):
idx = [i for i in range(1, len(mag) - 1) if mag[i] > mag[i - 1] and mag[i] >= mag[i + 1]]
res = []
for i in idx:
y0, y1, y2 = np.log(mag[i - 1]), np.log(mag[i]), np.log(mag[i + 1])
d = 0.5 * (y0 - y2) / (y0 - 2 * y1 + y2) # parabolic refinement
res.append(freqs[i] + d * (freqs[1] - freqs[0]))
return res[:n]
def main():
p = argparse.ArgumentParser()
p.add_argument('csv'); p.add_argument('--vowel'); p.add_argument('--length-cm', type=float)
p.add_argument('--c', type=float, default=350.0); p.add_argument('--fmax', type=float, default=5000)
p.add_argument('--nasal', action='store_true'); p.add_argument('--closure-section', default=None)
p.add_argument('--ladder-closure-cm', type=float); p.add_argument('--lip-cap-cm', type=float); p.add_argument('--closure-diameter-cm', type=float, default=0.12)
p.add_argument('--a-vp', type=float, default=0.56); p.add_argument('--port-section', type=int, default=22)
a = p.parse_args()
head, areas, extra = load(a.csv)
if a.nasal:
nasal_report(a, head, areas, extra)
return
freqs = np.arange(50, a.fmax, 5.0)
for v in head:
if a.vowel and v != a.vowel: continue
A = areas[v]
L = a.length_cm or extra['vtl_cm'][v]
sec = L * 100 / 100 / len(A) # cm per section
sec = L / len(A)
m = transfer(A, sec, a.c * 100, freqs)
print(f"{v} L={L:.2f}cm " + " ".join(f"F{i+1}={f:.0f}" for i, f in enumerate(peaks(freqs, m))))
if __name__ == '__main__':
main()

92
tools/formant-check/verify_renders.py Normal file
View file

@ -0,0 +1,92 @@
#!/usr/bin/env python3
"""Independent check of the WordRender outputs against published formants and the stress/burst cues. Shares no code with the C#
tests or the ladder; formants come from burg_formants.py (Burg LPC), f0 from autocorrelation, levels from the WAV.
verify_renders.py RENDER_DIR [--json OUT.json]
RENDER_DIR holds one subdirectory per word (v_<id>, ama, ana, nibu, nibu_flat, bin, bun) as written by the WordRender runs of
work/b3-render-all.sh. Prints a table: per vowel Burg F1/F2/F3 at the vowel's steady part against Hillenbrand 1995 male means and
Peterson & Barney 1952 male means (docs/research/articulatory-vocal-tract-sources.md) and the Story tube-chain prediction
(tube_formants.py, 17.52 cm), then the nasal murmur, the stress cue and the b burst. It measures; it does not listen.
"""
import sys, os, json, wave, subprocess, argparse
import numpy as np
HERE = os.path.dirname(os.path.abspath(__file__))
CSV = os.path.join(HERE, "..", "..", "docs", "research", "data", "story2008-area-functions.csv")
# Male means, F1/F2/F3 Hz (docs/research/articulatory-vocal-tract-sources.md).
HILL = {"i": (343, 2323, 3001), "ih": (429, 2034, 2687), "e": (476, 2090, 2692), "eh": (588, 1803, 2604), "ae": (591, 1930, 2595),
"uh": (621, 1181, 2548), "ah": (757, 1309, 2535), "aw": (656, 1023, 2521), "o": (498, 910, 2459), "oo_short": (469, 1123, 2435), "u": (380, 992, 2355)}
PB = {"i": (267, 2294, 2937), "ih": (392, 1993, 2569), "eh": (526, 1854, 2481), "ae": (664, 1727, 2420), "uh": (631, 1192, 2377),
"ah": (718, 1091, 2442), "aw": (568, 836, 2403), "oo_short": (437, 1023, 2245), "u": (307, 876, 2239)}
STORY_IDS = {"i": "i", "ih": "ɪ", "e": "e", "eh": "ɛ", "ae": "æ", "uh": "ʌ", "ah": "ɑ", "aw": "ɔ", "o": "o", "oo_short": "ʊ", "u": "u"}
def wav(path):
with wave.open(path) as w:
x = np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(float) / 32768.0
return x, w.getframerate()
def burg(path, t0, t1):
out = subprocess.run([sys.executable, os.path.join(HERE, "burg_formants.py"), path, str(t0), str(t1), "--json"], capture_output=True, text=True, check=True).stdout
j = json.loads(out)
return [j[f"F{i}"]["median"] if j[f"F{i}"] else None for i in (1, 2, 3)]
def tube(vowel):
out = subprocess.run([sys.executable, os.path.join(HERE, "tube_formants.py"), CSV, "--vowel", STORY_IDS[vowel], "--length-cm", "17.52"], capture_output=True, text=True).stdout
nums = [float(tok) for line in out.splitlines() for tok in line.replace("=", " ").replace("Hz", " ").split() if tok.replace(".", "", 1).isdigit()]
return out.strip().splitlines()[-1] if out.strip() else "?"
def rms(x, sr, a, b):
s = x[int(a * sr):int(b * sr)]
return float(np.sqrt((s ** 2).mean())) if len(s) else 0.0
def f0(x, sr, a, b):
s = x[int(a * sr):int(b * sr)]; s = s - s.mean()
lo, hi = sr // 400, sr // 60
c = [np.dot(s[:-l], s[l:]) for l in range(lo, hi)]
return sr / (lo + int(np.argmax(c)))
def pct(v, ref): return None if v is None else round(100.0 * (v - ref) / ref, 1)
def main():
ap = argparse.ArgumentParser(); ap.add_argument("dir"); ap.add_argument("--json"); a = ap.parse_args()
res = {"vowels": {}}
print("vowel Burg F1/F2/F3 (window) | Hillenbrand F1/F2 err% | P&B F1/F2 err% | tube chain")
for v in HILL:
p = os.path.join(a.dir, f"v_{v}", f"v_{v}.wav")
if not os.path.exists(p): print(v, "MISSING"); continue
f = burg(p, 0.07, 0.15) # 180 ms vowel; the onset ramp is excluded
h = HILL[v]; pb = PB.get(v)
res["vowels"][v] = {"burg": f, "hillenbrand": h, "pb": pb, "err_hill_pct": [pct(f[0], h[0]), pct(f[1], h[1])], "err_pb_pct": [pct(f[0], pb[0]), pct(f[1], pb[1])] if pb else None}
print(f"{v:9s}{f} | {h[:2]} {res['vowels'][v]['err_hill_pct']} | {pb[:2] if pb else '-'} {res['vowels'][v]['err_pb_pct'] if pb else ''} | {tube(v)}")
# nasal murmur: low first formant and weaker level than the flanking vowels
for w, c in (("ama", "m"), ("ana", "n")):
p = os.path.join(a.dir, w, f"{w}.wav")
if not os.path.exists(p): continue
x, sr = wav(p)
# steady parts only: the consonant's middle 60 percent against the mean power of the 70 ms before it and the 90 ms from 30 ms after it
c0, c1 = .18, .275
v1, v2 = rms(x, sr, c0 - .08, c0 - .01), rms(x, sr, c1 + .03, c1 + .12)
m = rms(x, sr, c0 + .2 * (c1 - c0), c1 - .2 * (c1 - c0))
db = round(20 * np.log10(max(m, 1e-9) / np.sqrt((v1 ** 2 + v2 ** 2) / 2)), 1)
res[w] = {"rms": [v1, m, v2], "murmur_db_re_flanking_vowels": db, "burg_murmur": burg(p, .20, .26)}
print(f"{w}: steady rms V/{c}/V {['%.4f' % q for q in (v1, m, v2)]} murmur {db} dB re the flanking vowels; Burg of the murmur {res[w]['burg_murmur']}")
# stress
p0, p1 = os.path.join(a.dir, "nibu_flat", "nibu_flat.wav"), os.path.join(a.dir, "nibu", "nibu.wav")
if os.path.exists(p0) and os.path.exists(p1):
x0, sr = wav(p0); x1, _ = wav(p1)
d = 0.18 * 220 / 160 - 0.18
# vowel i spans 0.095-0.275 (flat) / 0.095-0.3225 (primary); compare steady parts
rf, rp = rms(x0, sr, .15, .25), rms(x1, sr, .15, .25)
ff, fp = f0(x0, sr, .15, .25), f0(x1, sr, .15, .25)
res["stress"] = {"rms_flat": rf, "rms_primary": rp, "level_db": round(20 * np.log10(rp / rf), 2), "f0_flat": ff, "f0_primary": fp, "f0_ratio": round(fp / ff, 3),
"length_flat_s": len(x0) / sr, "length_primary_s": len(x1) / sr, "length_diff_s": round(len(x1) / sr - len(x0) / sr, 3), "expected_diff_s": round(d, 3)}
print("stress:", res["stress"])
# b release in bin: closure/ burst / vowel
p = os.path.join(a.dir, "bin", "bin.wav")
if os.path.exists(p):
x, sr = wav(p)
step = int(sr * 0.002); env = [(i / sr, float(np.abs(x[i:i + step]).max())) for i in range(0, len(x) - step, step)]
res["bin_env_peak_per_2ms"] = [(round(t * 1000), round(v, 4)) for t, v in env[:70]]
pk = max(v for _, v in env)
print("bin: peak", round(pk, 4), "b closure (0-60 ms) max", round(max(v for t, v in env if t < .06), 4), "release 60-95 ms max", round(max(v for t, v in env if .06 <= t < .1), 4), "vowel 110-260 ms max", round(max(v for t, v in env if .11 <= t < .26), 4))
if a.json: json.dump(res, open(a.json, "w"), indent=1, default=float)
if __name__ == "__main__": main()