Simulating articulatory trajectories with phonological feature interpolation

Abstract

As a first step towards a complete computational model of speech learning involving perception-production loops, we investigate the forward mapping between pseudo-motor commands and articulatory trajectories. Two phonological feature sets, based respectively on generative and articulatory phonology, are used to encode a phonetic target sequence. Different interpolation techniques are compared to generate smooth trajectories in these feature spaces, with a potential optimisation of the target value and timing to capture co-articulation effects. We report the Pearson correlation between a linear projection of the generated trajectories and articulatory data derived from a multi-speaker dataset of electromagnetic articulography (EMA) recordings. A correlation of 0.67 is obtained with an extended feature set based on generative phonology and a linear interpolation technique. We discuss the implications of our results for our understanding of the dynamics of biological motion.

1 Introduction

Recent advances in self-supervised learning (SSL) have led to progress in various speech processing tasks [1, 2] and language modelling from speech units [3]. These SSL models require increasingly greater amounts of (unlabelled) data to capture as much acoustic variance as possible. Moreover, their capacity to learn high-level language representations hinges on the quality of the underlying speech units [4], which are not linguistically interpretable [2, 5]. Importantly, these representations remain sensitive to contextual effects such as co-articulation [6] making them sub-optimal to efficiently code context-invariant phonological units.

According to the motor [7] or perceptuo-motor theories [8] of speech perception, humans’ ‘quest for invariance’ [9] is done by recovering motor or articulatory representations from an auditory input. These representations are supposed to be less variable than their acoustic counterparts. Several studies have found neuro-physiological correlates of this mental ‘sensory-to-motor’ inverse mapping in speech perception [10]. Incorporating such motor representation into SSL speech models could potentially improve their performance (e.g. noise robustness, low-resource downstream tasks, etc.) and lead to more plausible computational models of speech and language acquisition.

In a simplified perception-production loop of speech motor control [11] (see the bottom left of Figure 1), the motor commands are derived from the sensory, acoustic signal and used to generate the underlying articulatory trajectories (via a so-called forward model).

Figure 1: Simplified diagram of a speech perception-production loop (to the left). The focus of this work lies in the forward model and the linear probing (to the right).

Several computational models of such speech perception-production loop have been proposed in the literature [12, 13, 14, 15]. However, in most studies, the motor or articulatory representations are derived from a specific speaker or a specific articulatory model. As a result, these models, while enabling to study some of the underlying processes of speech acquisition, perception and motor control by simulation, are not designed to scale up to large numbers of speakers or languages.

As a first step towards a more universal SSL speech model integrating motor or articulatory representations, we focus in the present study on the forward model, i.e. from the motor commands to the generation of articulatory trajectories. First, we investigate different feature sets to encode a given phonetic target sequence, by relying either on the generative phonology (GP, [16]) or on the articulatory phonology (AP, [17]) theories. Importantly, a phonetic target is here encoded in terms of phonologically-motivated and articulatory-related categories (e.g. the place of articulation for a consonant in GP, the location and degree of a constriction in AP).

To generate continuous (and smooth) trajectories in these feature spaces (with a forward model), we test different interpolation methods, which differ from the dynamic properties desired at each phonetic target (e.g. zero and/or continuous velocity at each target). To account for the uncertainty of our timing heuristic and the (potential) target undershoot phenomenon, we also consider variants performing target optimisation both in space (find the offset from the ideal feature value, e.g. lips only partly closed) and time (reach a target sooner or later). Moreover, the proposed approach can deal with unspecified features for which the value depends on the context (for instance, the position of a constriction modulated by the vocalic context). The generated trajectories in the GP or AP feature spaces are evaluated using a linear probing technique. A linear model is learnt between the generated trajectories on the one hand, and the parameters of an articulatory model [18] built from a multi-speaker electromagnetic articulography (EMA) dataset on the other. Such an evaluation was also used in [19] to probe HuBERT representations, among other SSL models.

The main contributions and findings of the paper are the following: (i) we propose a general methodology to probe pseudo-motor commands and forward models in a computational model of a speech perception-production loop; (ii) we show that features derived from generative phonology (GP) correlate better with real articulatory recordings than articulatory phonology (AP) ones; (iii) a bit surprisingly, we show that a linear interpolation between these features better captures the dynamics of real articulatory data compared to a more complex one (spline based), with constraints on the velocity and/or acceleration at each phonetic target; (iv) we show that the use of unspecified (context-dependent) phonological features improves performance, probably by allowing the forward model to better account for natural co-articulation patterns.

The code with the data processing and interpolation methods can be found at [to ensure author anonymity, the link to the resource will be added after the review process].

2 Methodology

2.1 Phonological feature set

Two phonological feature sets were used. The first one, based on generative phonology and referred to as the GP feature set, was proposed by [20, Chapter 4]. It describes phonetic targets in terms of $26$ manner, laryngeal and place features. The GP feature set is considered under two variants: with unknown support and binary. Following [20, p. 91], some phonemes have zero-value features, notably because their values depend on their local context within an utterance or simply because they are irrelevant to the underlying phoneme. Hence, the ternary-valued (including the zero values) GP feature set is considered as is (to be used by interpolating methods handling unknown values) but also in a binary form, in which a zero value is considered as being the ‘absence’ of the given feature (thus, negatively valued).

The second feature set is based on articulatory phonology (AP), and the location and degrees of constriction of 5 major articulators in the vocal tract. AP-based features have been successfully used in automatic speech recognition, first within a Bayesian framework [21] to deal with pronunciation variability in spontaneous speech, and then in a DNN-based system [22] to increase the robustness to noise. The AP feature values come in the form of categorical distributions over totally ordered categories [21, p. 126]. Feature values are typically Dirac distributions, except for some phoneme features that depend on the phonemic context. Similarly to the GP feature set, we have an AP unknown variant by considering the non-Dirac distributed feature values as unknown values to be found contextually. This feature set is then used in a scalar version, in which each feature-value category is mapped to a real value (scalar AP: 8 features); and a one-hot version, in which the feature-value categories are one-hot encoded (one-hot AP: 32 features).

Table 1: Articulatory score: average Pearson correlation coefficients of the $6$ articulatory parameters and $6$ speakers. The scores correspond to each interpolation method’s best configuration. Default: binary features without optimisation. Variants: unknown features^μ, timing optimisation^†, timing and position optimisation^‡. (Standard error across the 6 different speakers was found to be $0.01$ on average.)

Feature set	# Features	Method	Score $\uparrow$
GP + one-hot phoneme	$73$	piecewise-cst	$\cellcolor{high!58!low!70}0.595$
		linearμ	$\mathbf{\cellcolor{high!100!low!70}0.679}$
		cubic Hermiteμ	$\cellcolor{high!94!low!70}0.668$
		natural cubic‡	$\cellcolor{high!92!low!70}0.663$
one-hot AP + one-hot phoneme	94	linearμ	$\cellcolor{high!92!low!70}0.663$
		cubic Hermiteμ	$\cellcolor{high!84!low!70}0.648$
		natural cubicμ	$\cellcolor{high!74!low!70}0.628$
scalar AP + one-hot phoneme	70	linearμ	$\cellcolor{high!88!low!70}0.656$
		cubic Hermiteμ	$\cellcolor{high!81!low!70}0.642$
		natural cubicμ	$\cellcolor{high!72!low!70}0.624$
one-hot phoneme	$47$	piecewise-cst	$\cellcolor{high!55!low!70}0.589$
		linear	$\cellcolor{high!83!low!70}0.645$
		cubic Hermite‡	$\cellcolor{high!81!low!70}0.642$
		natural cubic‡	$\cellcolor{high!75!low!70}0.630$
GP	$26$	piecewise-cst	$\cellcolor{high!40!low!70}0.559$
		linearμ	$\cellcolor{high!75!low!70}0.630$
		cubic Hermiteμ(†)	$\cellcolor{high!71!low!70}0.622$
		natural cubic‡	$\cellcolor{high!75!low!70}0.629$
one-hot AP	$32$	linearμ	$\cellcolor{high!64!low!70}0.608$
		cubic Hermiteμ‡	$\cellcolor{high!58!low!70}0.596$
		natural cubicμ	$\cellcolor{high!29!low!70}0.538$
scalar AP	$8$	linearμ	$\cellcolor{high!16!low!70}0.511$
		cubic Hermiteμ‡	$\cellcolor{high!14!low!70}0.506$
		natural cubicμ	$\cellcolor{high!0!low!70}0.479$

To ensure that the feature-level information for phonemes is relevant, we also use a feature set with one-hot phoneme encodings. In total, we evaluate seven feature sets: one-hot phonemes, scalar AP (also enriched with one-hot phonemes), one-hot AP (also enriched with one-hot phonemes), binary GP and unknown-supporting GP (also enriched with one-hot phonemes).

2.2 Dataset

The articulatory data comes from the publicly available MOCHA-TIMIT dataset. It provides electromagnetic articulography (EMA) recordings for 460 short sentences read by 8 British English speakers along 12 dimensions (2D midsagittal coordinates for 6 articulators: tongue tip, tongue body, tongue dorsum, lower incisor, upper lip and lower lip.) In this study, we consider 6 speakers, namely fsew0, msak0, ffes0, mjjn0, faet0 and maps0, because a sequence of waveform files did not match the given transcriptions for the other two speakers. For each speaker, the 460 utterances are split into 410 for training (out of which 20 are randomly drawn for development) and 50 for testing.

The EMA data is low-pass filtered at $50\text{\,}\mathrm{Hz}$ and down-sampled from $500\text{\,}\mathrm{Hz}$ to $100\text{\,}\mathrm{Hz}$. As with [23], raw EMA data is then converted into an easier-to-interpret and lower-dimensional set of 6 ‘articulatory parameters’ (jaw height, tongue body, tongue back, tongue tip, lip protrusion and lip height) using a linear decomposition technique, which is often referred to as ‘guided PCA’. It aims to decouple the jaw from tongue and lip movements and extract independent degrees of freedom from the vocal tract.

Each audio recording was segmented at the phonetic level using the Montreal Forced Aligner. Based on the resulting phonetic segmentation, the initial and final utterance silences were removed and the utterances with non-silence boundary phones were discarded. The filtered phonetic segmentations were later mapped into featural segmentations by replacing the phones with the phonological features of their underlying phonemes (lookup table mapping). Finally, we inferred timings for the phonological targets from the time midpoint of each phoneme in the featural segmentation.

2.3 Forward model

Different forward mapping techniques are tested to generate continuous trajectories from the discrete pseudo-motor commands provided by GP and AP (Section 2.1). As a first baseline, we use the piecewise-constant interpolation, which keeps all the phonological features constant for the duration given by the phonetic segmentation.

To test a smoothness degree that better fits the articulatory space, we test linear and cubic interpolation methods. Specifically, we consider two cubic methods: the cubic Hermite spline and the natural cubic spline. The former enforces zero velocity at all targets, and continuity of both position and velocity (so the acceleration is possibly discontinuous at the targets); whereas the latter enforces continuity of all position, velocity and acceleration, with zero acceleration at the initial and final targets.

Formally, let $d\in\mathbf{N}_{>0}$ be the number of phonological features. For a given utterance, let $K\in\mathbf{N}_{>0}$ be its number of non-boundary (or intermediate) targets. Then, its featural segmentation is denoted by $(\mathbf{X},\mathbf{Y})\in\mathbf{R}^{(K+2)\times d}\times\mathbf{R}^{(K+2)\times 2}$, where $\mathbf{X}$ contains the $K+2$ target positions, and $\mathbf{Y}$ the targets’ time intervals. In this work, we remove the boundary silences, so $y_{1,\ast}=\mathbf{0}_{2}$, $y_{K+2,\ast}=y_{K+1,2}\mathbf{1}_{2}$, and $x_{1,\ast}=x_{K+2,\ast}=\mathbf{0}_{d}$. From this, we deduce a vector of midpoint target timings $\mathbf{t}\triangleq\frac{1}{2}Y\mathbf{1}_{2}$. The (base) interpolating function for the utterance $(\mathbf{X},\mathbf{Y})$, expressed as $f(\tau;\mathbf{X},\mathbf{t})$, $0\leq\tau\leq t_{K+2}$, thus satisfies

$$f(t_{k};\mathbf{X},\mathbf{t})=x_{k,\ast},$$ (1)

for all $1\leq k\leq K+2$.

To tackle the uncertainty of the midpoint-timing heuristic assumed for the base interpolating functions and to allow for target undershoot, we include two additive optimisations over time and space. The optimised interpolating function $f$ learns the target positions and timings $(\mathbf{X}^{\prime},\mathbf{t}^{\prime})$ such that

$$f(\tau;\mathbf{X},\mathbf{t})=g(\tau,\mathbf{X}^{\prime},\mathbf{t}^{\prime}),0\leq\tau\leq t_{K+2},$$ (2)

where the base interpolating function $g$ satisfies Equation 1, by minimising the objective function

$$L_{\lambda}(\mathbf{X}^{\prime},\mathbf{t}^{\prime})\triangleq \int_{0}^{t_{K+2}}\lVert g^{\prime\prime}(\tau;\mathbf{X}^{\prime},\mathbf{t}^{\prime})\rVert_{2}^{2}\,d\tau$$

$$+\lambda\sum_{k=2}^{K+1}\lVert g(t_{k}^{\prime};\mathbf{X}^{\prime},\mathbf{t}^{\prime})-x_{k,\ast}\rVert_{2}^{2}.$$ (3)

We optimise the objective function $L_{\lambda}$ by on-utterance gradient descent with the initialisation $(\mathbf{X}^{\prime},\mathbf{t}^{\prime})=(\mathbf{X},\mathbf{t})$.

2.4 Linear probing

To evaluate the different interpolation methods described above, we adopted the same metric as used in [19], referred to as the articulatory score.

Let $S_{i}=\{(\mathbf{X}_{j},\mathbf{t}_{j},\mathbf{Z}_{j})_{j}\}$ be the set of utterances for the $i$th speaker, with $\mathbf{Z}_{j}\in\mathbf{R}^{n_{\mathbf{Z}_{j}}\times 6}$ being the $n_{\mathbf{Z}_{j}}$ articulatory parameters of the $j$th utterance. Then, for each interpolation method $f$ and speaker $S_{i}$, we learn a linear transformation $h_{i}$ that minimises the reconstruction loss from the interpolated articulatory trajectories and the expected articulatory parameters as follows

$$\mathcal{L}_{i}=\frac{1}{\lvert S_{i}\rvert}\sum_{(\mathbf{X},\mathbf{t},\mathbf{Z})\in S_{i}}\frac{1}{n_{\mathbf{Z}}}\sum_{k=1}^{n{{}_{\mathbf{Z}}}}\left\lVert h_{i}\left(f\left(\frac{k}{100};\mathbf{X},\mathbf{t}\right)\right)-z_{k,\ast}\right\rVert_{2}^{2}.$$ (4)

This is done via gradient descent with a learning rate of $0.001$ via the Adam optimiser, with $\beta=(0.9,0.999)$. We run the learning procedure for $100$ epochs unless it is early stopped when the validation loss stops decreasing (patience fixed at $5$).

For the optimised cubic interpolations, we first do a grid search of the hyper-parameters used in the on-utterance target optimisation: (i) timing learning rate within $\{$1\text{\times}{10}^{-6}$,$5\text{\times}{10}^{-6}$,$1\text{\times}{10}^{-5}$,$5\text{\times}{10}^{-5}$,$1\text{\times}{10}^{-4}$\}$, (ii) position learning rate within $\{$1\text{\times}{10}^{-3}$,$1\text{\times}{10}^{-2}$,$1\text{\times}{10}^{-1}$\}$, and (iii) loss weight parameter $\lambda\in\{0,$1\text{\times}{10}^{3}$,$1\text{\times}{10}^{4}$,$1\text{\times}{10}^{5}$,$1\text{\times}{10}^{6}$,$1\text{\times}{10}^{7}$\}$. For each interpolation method, the hyper-parameter configuration with the highest articulatory score in the development set was selected. Finally, we compute the Pearson correlation coefficient (PCC) between the learnt linear projections and the articulatory parameters. The final articulatory score of an interpolation method is then the average PCC of all articulatory parameters and speakers.

3 Results

We run the linear, cubic Hermite spline and natural cubic spline interpolation methods (Section 2.3) on the seven feature sets derived from one-hot phoneme encoding, AP and GP theories (Section 2.1). The piecewise-constant interpolation can only be run on fully specified feature sets, here the one-hot phoneme and binary GP feature sets.

Table 2: Articulatory score for each speaker on the GP feature set enriched with one-hot phonemes.

Method	msak0	fsew0	ffes0	mjjn0	maps0	faet0	Average
piecewise-cst	$\cellcolor{high!47!low!70}0.634$	$\cellcolor{high!39!low!70}0.616$	$\cellcolor{high!27!low!70}0.590$	$\cellcolor{high!28!low!70}0.591$	$\cellcolor{high!3!low!70}0.537$	$\cellcolor{high!34!low!70}0.604$	$\cellcolor{high!30!low!70}0.595$
linear	$\mathbf{\cellcolor{high!90!low!70}0.729}$	$\mathbf{\cellcolor{high!79!low!70}0.704}$	$\mathbf{\cellcolor{high!62!low!70}0.666}$	$\mathbf{\cellcolor{high!58!low!70}0.658}$	$\mathbf{\cellcolor{high!42!low!70}0.623}$	$\mathbf{\cellcolor{high!74!low!70}0.693}$	$\mathbf{\cellcolor{high!68!low!70}0.679}$
cubic Hermite	$\cellcolor{high!85!low!70}0.718$	$\cellcolor{high!75!low!70}0.695$	$\cellcolor{high!57!low!70}0.655$	$\cellcolor{high!54!low!70}0.649$	$\cellcolor{high!37!low!70}0.611$	$\cellcolor{high!69!low!70}0.681$	$\cellcolor{high!63!low!70}0.668$
natural cubic	$\cellcolor{high!82!low!70}0.711$	$\cellcolor{high!71!low!70}0.686$	$\cellcolor{high!55!low!70}0.652$	$\cellcolor{high!46!low!70}0.632$	$\cellcolor{high!38!low!70}0.613$	$\cellcolor{high!70!low!70}0.685$	$\cellcolor{high!60!low!70}0.663$

Table 3: Articulatory score for each articulatory parameter on the GP feature set enriched with one-hot phonemes.

Method	Jaw height	Tongue body	Tongue dorsum	Tongue tip	Lip protrusion	Lip height	Average
piecewise-cst	$\cellcolor{high!53!low!70}0.646$	$\cellcolor{high!56!low!70}0.653$	$\cellcolor{high!6!low!70}0.543$	$\cellcolor{high!4!low!70}0.539$	$\cellcolor{high!1!low!70}0.532$	$\cellcolor{high!58!low!70}0.658$	$\cellcolor{high!30!low!70}0.595$
linear	$\mathbf{\cellcolor{high!84!low!70}0.715}$	$\mathbf{\cellcolor{high!100!low!70}0.750}$	$\mathbf{\cellcolor{high!43!low!70}0.625}$	$\mathbf{\cellcolor{high!44!low!70}0.627}$	$\mathbf{\cellcolor{high!41!low!70}0.621}$	$\mathbf{\cellcolor{high!94!low!70}0.736}$	$\mathbf{\cellcolor{high!68!low!70}0.679}$
cubic Hermite	$\cellcolor{high!79!low!70}0.703$	$\cellcolor{high!96!low!70}0.742$	$\cellcolor{high!40!low!70}0.618$	$\cellcolor{high!39!low!70}0.616$	$\cellcolor{high!36!low!70}0.610$	$\cellcolor{high!87!low!70}0.722$	$\cellcolor{high!63!low!70}0.668$
natural cubic	$\cellcolor{high!81!low!70}0.708$	$\cellcolor{high!92!low!70}0.733$	$\cellcolor{high!34!low!70}0.604$	$\cellcolor{high!35!low!70}0.606$	$\cellcolor{high!39!low!70}0.615$	$\cellcolor{high!83!low!70}0.713$	$\cellcolor{high!60!low!70}0.663$

Table 1 reports the articulatory score by feature space and interpolation method on a held-out test. The scores correspond to each interpolation method’s best configuration on each feature set. We observe that the scalar AP proves to be (very) difficult to interpolate on, but replacing the fixed values (probably needing learning) with equidistant values in the form of one-hot encodings helps close the gap to the GP features. Surprisingly, the one-hot phoneme encodings are the best single feature set.

Given the potential complementarity of information between the GP/AP feature sets and the one-hot phoneme encodings, we probe the GP and AP features enriched with the latter. The mix turns out to be beneficial for all the interpolation methods, although we lose the reduced number of interpretable features sought for the inverse models in perspective.

Tables 3 and 3 show the articulatory scores per speaker and articulatory parameter, respectively, on the best feature space, namely GP features enriched with one-hot phonemes. In both cases, the ranking induced by the average score (linear $\succ$ cubic Hermite $\succ$ natural cubic $\succ$ piecewise constant) is met throughout the conditions, bar the two speakers maps0 and faet0, the jaw height and the lip protrusion (natural cubic $\succ$ cubic Hermite). Interestingly, from Tables 1, 3 and 3, it is clear that the linear interpolation method better exploits the given phonological spaces, regardless of the feature nature, speaker or articulatory parameter.

In Table 1, we see that most of the best scores reported were obtained on features with unknown support. The results in Table 4 support the hypothesis that, in general, keeping and interpolating unknown feature values is better than associating them with a fixed value. On the other hand, the effect of the target timing and/or position optimisation depends on the interpolation method. For instance, when we optimise the timings on the cubic Hermite spline, the articulatory score does not improve, and the spatial optimisation has the same (negative) impact. This is why we see few cubic Hermite interpolations with target optimisations in Table 1.

Table 4: Comparison of GP + one-hot phoneme feature set variants and the effect of target optimisation. The two leftmost scores correspond to the binary and unknown-supporting feature-set variants without optimisation, and the right scores to the timing-only and the timing-and-position optimisations on the underlined feature sets.

Method	Non-optimised		Optimised
	Binary	Unknown	Time	Time & space
linear	$0.659$	$\mathbf{0.679}$
cubic Hermite	$0.645$	$\underline{\mathbf{0.668}}$	$0.660$	$0.649$
natural cubic	$\underline{0.623}$	$0.638$	$0.624$	$\mathbf{0.663}$

4 Conclusion

In this study, we have analysed phonological features as potential pseudo-motor commands in a computational model of a speech perception-production loop. We found that: (i) smooth trajectories on generative phonology features correlate better with articulatory parameters than those on articulatory phonology ones, with a correlation coefficient of $0.67$ when GP features are enriched with one-hot phoneme encodings, (ii) a linear forward model better captures the dynamics of real articulatory data, but target optimisation (in terms of timing and/or position) helps a cubic model to reduce the gap, (iii) with the AP features, a better correlation coefficient is obtained with a one-hot encoding, in which all the values for a given feature are equidistant to one another, rather than with a fixed, scalar continuous one, (iv) interpolating unknown (or context-dependent) features is better than associating a fixed value with them.

Since interpolating under-specified dimensions of articulatory targets appears to lead to a better fit, future work could try to push this strategy further by incorporating more under-specified dimensions in featural segmentations, thus enabling the smoothness of the forward model’s trajectories to better model co-articulation.

Further work should also investigate the reason why linear interpolation of articulatory targets outperforms smoother cubic spline interpolation. This is surprising since articulatory trajectories are not linear, cubic splines are excellent interpolators and the biological motion literature suggests that smooth trajectories should fit articulatory data well [24]. A possible interpretation is that unwarranted assumptions in our analyses cause the observed advantage of linear interpolation methods. Based on our results, the advantage of linear interpolation does not appear sensitive to assumptions regarding target definition (generative vs articulatory, specification, timing, position). It may be the case, however, that the assumption of a fixed inventory of targets at the phonemic level is too optimistic, even for a single speaker, at least without controlling further variables that may modulate target parameters, such as prosodic effects. To test this hypothesis, simulated trajectories could be used to determine if ignoring prosodic effects would predict an advantage for linear interpolation even when using a cubic spline model for the dynamics. Alternatively, our results may indicate that classical results on biological motion, obtained in highly controlled settings, do not accurately characterize the dynamics of biological motion in less restricted environments. To test this hypothesis, our methodology could be applied to more controlled trajectories (e.g. isolated syllables), where it should find that smoother trajectories provide a better fit than linear interpolation. This would validate our methodology, which could then be leveraged to better understand the dynamics of biological motion in the wild.