diff --git a/docs/research/articulatory-vocal-tract-sources.md b/docs/research/articulatory-vocal-tract-sources.md new file mode 100644 index 0000000..61a0f69 --- /dev/null +++ b/docs/research/articulatory-vocal-tract-sources.md @@ -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_`), 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. diff --git a/docs/research/data/story2008-area-functions.csv b/docs/research/data/story2008-area-functions.csv new file mode 100644 index 0000000..0a6e572 --- /dev/null +++ b/docs/research/data/story2008-area-functions.csv @@ -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 diff --git a/src/AquaSynth.Core/FaustExport.VocalLadder.cs b/src/AquaSynth.Core/FaustExport.VocalLadder.cs index ae06b70..c3a9ca0 100644 --- a/src/AquaSynth.Core/FaustExport.VocalLadder.cs +++ b/src/AquaSynth.Core/FaustExport.VocalLadder.cs @@ -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 /// /// Emits the area, reflection and junction definitions and the recursive ladder for , 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). /// - private static IReadOnlyDictionary EmitVocalLadder( + private static (IReadOnlyDictionary Lips, IReadOnlyDictionary 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(StringComparer.OrdinalIgnoreCase); - for (var m = 0; m < members.Count; m++) + var glottalPressures = new Dictionary(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); } } diff --git a/src/AquaSynth.Core/FaustExport.cs b/src/AquaSynth.Core/FaustExport.cs index d18bf3c..ca98116 100644 --- a/src/AquaSynth.Core/FaustExport.cs +++ b/src/AquaSynth.Core/FaustExport.cs @@ -647,8 +647,8 @@ public static partial class FaustEmitter } } - var lips = members.Count == 0 - ? new Dictionary() + var (lips, glottalPressures) = members.Count == 0 + ? (new Dictionary(), new Dictionary()) : 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); /// - /// 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 . 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 . 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 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. /// private static string DeformedAreaExpression( ParameterMap parameters, TractAreaFunction shape, AcousticAreaControl? area, string areaPath, - float position) + float position, + IReadOnlyList? 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 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 waveClocks, diff --git a/src/AquaSynth.Core/IpaGestureExperiment.cs b/src/AquaSynth.Core/IpaGestureExperiment.cs index 6916988..b343df3 100644 --- a/src/AquaSynth.Core/IpaGestureExperiment.cs +++ b/src/AquaSynth.Core/IpaGestureExperiment.cs @@ -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 targets, IEnumerable 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 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 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 tokens) + private static string[] ExpectedSurfaces(HashSet tokens, string descriptor) { var surfaces = new List { @@ -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 tokens, IReadOnlyList splines) + private static float DirectionScore(HashSet tokens, string descriptor, IReadOnlyList splines) { var checks = new List(); 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")); } diff --git a/src/AquaSynth.Core/Model.cs b/src/AquaSynth.Core/Model.cs index 5052134..175bdfa 100644 --- a/src/AquaSynth.Core/Model.cs +++ b/src/AquaSynth.Core/Model.cs @@ -443,11 +443,15 @@ public sealed record AcousticPortNetwork( public IReadOnlyList Connections { get; init; } = Connections ?? Array.Empty(); } +/// A named measured tract shape: areas in cm^2 from the glottis to the lips, any number of sections. +public sealed record VowelAreaShape(string Name, IReadOnlyList Areas); + public sealed record AreaFunction( string Name, TractAreaFunction Shape, AcousticAreaControl? Deformation = null, - int EmitSections = 8); + int EmitSections = 8, + IReadOnlyList? 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, diff --git a/src/AquaSynth.Core/NasalTractAreas.cs b/src/AquaSynth.Core/NasalTractAreas.cs new file mode 100644 index 0000000..572fb93 --- /dev/null +++ b/src/AquaSynth.Core/NasalTractAreas.cs @@ -0,0 +1,33 @@ +namespace AquaSynth.Dsl; + +/// +/// 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 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. +/// +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 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, + ]; + + /// + /// 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. + /// + public const float VelopharyngealPortMaximumAreaSquareCentimeters = 0.8f; + + /// The scaffold's morphology line value for areas=, invariant-culture comma separated. + public static string AreasList() => string.Join(",", Areas.Select(area => area.ToString("0.##", System.Globalization.CultureInfo.InvariantCulture))); +} diff --git a/src/AquaSynth.Core/PatchScript.cs b/src/AquaSynth.Core/PatchScript.cs index 622bd7c..ff2f3db 100644 --- a/src/AquaSynth.Core/PatchScript.cs +++ b/src/AquaSynth.Core/PatchScript.cs @@ -11,9 +11,16 @@ public static class PatchScript /// /// 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 , the one table that also says which measured tract shape each selects. /// 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= source=story2008: attaches the eleven measured vowel shapes to a morphology and + // registers a weight surface /vocal/areas//area/vowel_ 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 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 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 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 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 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 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", diff --git a/src/AquaSynth.Core/ProbeTimelineReport.cs b/src/AquaSynth.Core/ProbeTimelineReport.cs index d5ab2fb..42bf89f 100644 --- a/src/AquaSynth.Core/ProbeTimelineReport.cs +++ b/src/AquaSynth.Core/ProbeTimelineReport.cs @@ -7,7 +7,7 @@ public sealed record ProbeTimelineSample(int Block, string Primitive, string Sig public static class ProbeTimelineReport { - public static IReadOnlyList Build(SynthPatch patch, string networkName, int blocks = 4) + public static IReadOnlyList 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(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) diff --git a/src/AquaSynth.Core/StoryVowelAreas.cs b/src/AquaSynth.Core/StoryVowelAreas.cs new file mode 100644 index 0000000..51ef70b --- /dev/null +++ b/src/AquaSynth.Core/StoryVowelAreas.cs @@ -0,0 +1,93 @@ +namespace AquaSynth.Dsl; + +/// +/// 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 and every vowel is read at that length; the sections are +/// resampled onto the ladder's cells by . +/// +/// Control surface ids (ascii, used in /vocal/areas/<i>/area/vowel_<id>): 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. +/// +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 Areas); + + public static string SurfaceField(string id) => $"vowel_{id}"; + + /// + /// 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). + /// + 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"; + + /// + /// 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. + /// + public static string? IdForDescriptor(string descriptor) => + All.FirstOrDefault(vowel => vowel.Descriptor.Equals(descriptor.Replace('-', '_').ToLowerInvariant(), StringComparison.Ordinal))?.Id; + + /// + /// Area in cm^2 of 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. + /// + public static float AreaAt(IReadOnlyList 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 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]), + ]; +} diff --git a/src/AquaSynth.Core/WeksaPhones.cs b/src/AquaSynth.Core/WeksaPhones.cs index 9939c79..4e21898 100644 --- a/src/AquaSynth.Core/WeksaPhones.cs +++ b/src/AquaSynth.Core/WeksaPhones.cs @@ -71,19 +71,30 @@ public static class WeksaPhonesReader /// /// A Weksa phones line as a word the tract can render. Every phone must be in the DSL's IPA table /// (); a phone outside it is listed in 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 (): the marked phone is longer, louder and higher in f0. /// public sealed record WeksaWordPlan(IReadOnlyList<(IpaGestureExperimentTarget Target, float Intensity)> Phones, IReadOnlyList Unsupported) { public const float VowelSeconds = 0.18f, ConsonantSeconds = 0.095f; public bool Renderable => Unsupported.Count == 0; - public static float StressIntensity(WeksaStress stress) => stress switch + /// Duration factor, amplitude factor on the source flow and f0 ratio of a stressed phone against an unstressed one. + public readonly record struct StressCue(float Duration, float Amplitude, float PitchRatio); + + /// + /// 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. + /// + 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, []); } diff --git a/src/AquaSynth.Faust/IpaTrialOrchestrator.cs b/src/AquaSynth.Faust/IpaTrialOrchestrator.cs index 5f6e2fe..4454451 100644 --- a/src/AquaSynth.Faust/IpaTrialOrchestrator.cs +++ b/src/AquaSynth.Faust/IpaTrialOrchestrator.cs @@ -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 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 tokens) + private static string[] ExpectedSurfaces(HashSet tokens, string descriptor) { var surfaces = new List { @@ -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 { diff --git a/tests/AquaSynth.Dsl.Tests/LipReflectanceTests.cs b/tests/AquaSynth.Dsl.Tests/LipReflectanceTests.cs new file mode 100644 index 0000000..f2cb069 --- /dev/null +++ b/tests/AquaSynth.Dsl.Tests/LipReflectanceTests.cs @@ -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); + } +} diff --git a/tests/AquaSynth.Dsl.Tests/NasalCouplingTests.cs b/tests/AquaSynth.Dsl.Tests/NasalCouplingTests.cs new file mode 100644 index 0000000..cac20de --- /dev/null +++ b/tests/AquaSynth.Dsl.Tests/NasalCouplingTests.cs @@ -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")); + } +} diff --git a/tests/AquaSynth.Dsl.Tests/PatchScriptTests.cs b/tests/AquaSynth.Dsl.Tests/PatchScriptTests.cs index 5e4eec0..0a0e6df 100644 --- a/tests/AquaSynth.Dsl.Tests/PatchScriptTests.cs +++ b/tests/AquaSynth.Dsl.Tests/PatchScriptTests.cs @@ -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(); diff --git a/tests/AquaSynth.Dsl.Tests/ProbeTimelineSampleRateTests.cs b/tests/AquaSynth.Dsl.Tests/ProbeTimelineSampleRateTests.cs new file mode 100644 index 0000000..9d381d9 --- /dev/null +++ b/tests/AquaSynth.Dsl.Tests/ProbeTimelineSampleRateTests.cs @@ -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); + } +} diff --git a/tests/AquaSynth.Dsl.Tests/StoryVowelAreasTests.cs b/tests/AquaSynth.Dsl.Tests/StoryVowelAreasTests.cs new file mode 100644 index 0000000..e38739b --- /dev/null +++ b/tests/AquaSynth.Dsl.Tests/StoryVowelAreasTests.cs @@ -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); + } +} diff --git a/tests/AquaSynth.Dsl.Tests/StressAndReleaseRenderTests.cs b/tests/AquaSynth.Dsl.Tests/StressAndReleaseRenderTests.cs new file mode 100644 index 0000000..bebfa1c --- /dev/null +++ b/tests/AquaSynth.Dsl.Tests/StressAndReleaseRenderTests.cs @@ -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? 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(); + 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); + } +} diff --git a/tests/AquaSynth.Dsl.Tests/VocalTractFormantTests.cs b/tests/AquaSynth.Dsl.Tests/VocalTractFormantTests.cs index ca0f694..775c07d 100644 --- a/tests/AquaSynth.Dsl.Tests/VocalTractFormantTests.cs +++ b/tests/AquaSynth.Dsl.Tests/VocalTractFormantTests.cs @@ -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? 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 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 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 { [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 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 { ["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()), 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())); + 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}"); + } } diff --git a/tests/AquaSynth.Dsl.Tests/WeksaPhonesTests.cs b/tests/AquaSynth.Dsl.Tests/WeksaPhonesTests.cs index 27ef6d6..34349af 100644 --- a/tests/AquaSynth.Dsl.Tests/WeksaPhonesTests.cs +++ b/tests/AquaSynth.Dsl.Tests/WeksaPhonesTests.cs @@ -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 _)); } } diff --git a/tools/formant-check/burg_formants.py b/tools/formant-check/burg_formants.py new file mode 100644 index 0000000..9f72aed --- /dev/null +++ b/tools/formant-check/burg_formants.py @@ -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() diff --git a/tools/formant-check/ladder_formants.py b/tools/formant-check/ladder_formants.py new file mode 100644 index 0000000..eae87f8 --- /dev/null +++ b/tools/formant-check/ladder_formants.py @@ -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))) diff --git a/tools/formant-check/spectrum_dips.py b/tools/formant-check/spectrum_dips.py new file mode 100644 index 0000000..b3083e0 --- /dev/null +++ b/tools/formant-check/spectrum_dips.py @@ -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='= 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() diff --git a/tools/formant-check/tube_formants.py b/tools/formant-check/tube_formants.py new file mode 100644 index 0000000..9697744 --- /dev/null +++ b/tools/formant-check/tube_formants.py @@ -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() diff --git a/tools/formant-check/verify_renders.py b/tools/formant-check/verify_renders.py new file mode 100644 index 0000000..80c009b --- /dev/null +++ b/tools/formant-check/verify_renders.py @@ -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_, 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="