Lu Xie

Pascal Brault

Anne-Lise Thomann

Xiao Yang

Yong Zhang

Guangyi Shang
email: gyshang@buaa.edu.cn

Molecular Dynamics Simulation of Al-Co-Cr-Cu-Fe-Ni high entropy alloy thin film growth

Keywords: MOLECULAR DYNAMICS SIMULATION, COATINGS, HIGH-ENTROPY ALLOYS

Molecular dynamics (MD) simulations are used to study AlCoCrCuFeNi high entropy alloy (HEA) thin film growth on a silicon (100) substrate. Effect of the atomic composition is studied on morphology and atomic scale structure. Input data are chosen to fit with experimental operating conditions of magnetron sputtering deposition process. It is observed that the different structures are determined by the chemical composition and atomic size mismatch. The simulated results are in good agreement with the solid-solution formation rules proposed by Zhang et al [1] for multi-principal component HEAs which based on the two parameters  and , respectively describing describe the comprehensive effect of the atomic-size difference in the n-element alloy and the effects of enthalpy and entropy of mixing on formation of multi-component solid-solutions. When   1.1 and  ≤ 6.6%, the multi-component solid solution phase could form. In contrast, the multi-component alloys forming intermetallic compounds and bulk metallic glasses (BMG) have larger value of  and smaller value of . The value of  for BMG is smaller than that of intermetallic compounds.

Introduction:

High-entropy alloys (HEAs) have attracted increasing attentions because of their unique compositions, microstructures, and adjustable properties. Proposed by Yeh et al [START_REF] Zhang | Microstructures and properties of high-entropy alloys[END_REF][START_REF] Yeh | Nanostructured High-Entropy Alloys with Multiple Principal Elements: Novel Alloy Design Concepts and Outcomes[END_REF], HEAs contain five or more principal elements in equal or near-equal atomic ratios, the corresponding atomic concentrations lie between 5% and 35%. Because of the high entropy of mixing, formation of brittle intermetallic phases can be avoided and simple disordered solid-solution structures like the face-centred cubic (FCC) and bodycentred cubic (BCC) ones can be stabilized [START_REF] Yang | Prediction of high-entropy stabilized solid-solution in multicomponent alloys[END_REF][START_REF] Zhang | Microstructures and properties of high-entropy alloys[END_REF]. The combination of numerous metallic elements with different sizes induces a sluggish atomic diffusion, leading to nanocrystalline or amorphous structures. Bulk HEAs demonstrate superior potential for engineering applications [START_REF] Lu | An assessment on the future development of high-entropy alloys: Summary from a recent workshop[END_REF] due to their high strength [START_REF] Zhou | Solid solution alloys of AlCoCrFeNiTix with excellent room-temperature mechanical properties[END_REF], hardness [START_REF] Zhou | Microstructure and compressive properties of multicomponent Alx(TiVCrMnFeCoNiCu)100-x high-entropy alloys[END_REF], wear resistance [START_REF] Hsu | Effect of iron content on wear behavior of AlCoCrFexMo0.5Ni high-entropy alloys[END_REF], high-temperature softening resistance and oxidation resistance [START_REF] Tsai | Deformation and annealing behaviors of high-entropy alloy Al0.5CoCrCuFeNi[END_REF]. HEAs can be viewed as atomic-scale composites; their final properties come from: 1) the combination of properties of their constituting elements and 2) the stabilized solid solutions and structures. Thus for the design of new alloys dedicated to a given application, the prediction of structural features of the material would be of particular interest.

However, it is quite difficult to synthesize HEAs as bulk materials because very fast cooling rates are required. Compared to many thin film deposition techniques, the sputtering method is an easy way to control the chemical composition, grain size, preferred orientation and hardness by varying the deposition parameters. Few studies on HEA thin films deposited by magnetron sputtering have been published [START_REF] Yeh | Nanostructured High-Entropy Alloys with Multiple Principal Elements: Novel Alloy Design Concepts and Outcomes[END_REF][START_REF] Chang | Influence of substrate bias, deposition temperature and post-deposition annealing on the structure and properties of multi-principal-component (AlCrMoSiTi)N coatings[END_REF][START_REF] Tsai | Diffusion barrier properties of AlMoNbSiTaTiVZr high-entropy alloy layer between copper and silicon[END_REF][START_REF] Chen | Nanostructured nitride films of multielement high-entropy alloys by reactive DC sputtering[END_REF][START_REF] Dolique | Thermal stability of AlCoCrCuFeNi high entropy alloy thin films studied by in-situ XRD analysis[END_REF][START_REF] Dolique | Complex structure/composition relationship in thin films of AlCoCrCuFeNi high entropy alloy[END_REF]. In addition, understanding the growth mode of such complex multi-element films remains a challenging research topic. In this view, carrying out simulation of the process at the atomic level and comparing with the published literatures is of great interest. Molecular Dynamics (MD) simulations are a powerful tool in computational materials science, which can provide a deep insight into the atomic mechanism of ordering of HEA. It can reproduce such processes and provide information not accessible experimentally. In previous works [START_REF] Dolique | Thermal stability of AlCoCrCuFeNi high entropy alloy thin films studied by in-situ XRD analysis[END_REF][START_REF] Dolique | Complex structure/composition relationship in thin films of AlCoCrCuFeNi high entropy alloy[END_REF], the growth of AlCoCrCuFeNi HEA clusters and their annealing properties have been studied. Moreover, phase transition from bcc to fcc structure was evidenced by MD simulation [START_REF] Xie | AlCoCrCuFeNi high entropy alloy cluster growth and annealing on silicon: A classical molecular dynamics simulation study[END_REF][START_REF] Lu | Molecular dynamics simulations of clusters and thin film growth in the context of plasma sputtering deposition[END_REF], in agreement with experimental results obtained on bulk or thin film alloys. In this work, we performed a deeper and more precise work on the microstructure evolution of HEA thin films grown in operating conditions corresponding to magnetron sputter deposition. First, we studied the nucleation and growth processes occurring in HEA thin film formation; Second, the conditions of solid solution formation were explored in terms of the  and , parameters typically used for the characterization of multi-components HEA materials [START_REF] Yang | Prediction of high-entropy stabilized solid-solution in multicomponent alloys[END_REF].

In those simulations, the atomic interactions for the Al-Co-Cr-Cu-Fe-Ni alloy were described by an embedded atom method many-body potential [START_REF] Daw | The embedded-atom method: a review of theory and applications[END_REF][START_REF] Zhou | Misfit-energy-increasing dislocations in vapordeposited CoFe/NiFe multilayers[END_REF]. Tersoff empirical potential was used for describing the substrate Si-Si interactions [START_REF] Graves | Molecular dynamics for low temperature plasma-surface interaction studies[END_REF][START_REF] Tersoff | Modeling solid-state chemistry: Interatomic potentials for multicomponent systems[END_REF]. The interactions of metal atoms with silicon were described using a Lennard-Jones potential and Lorenz-Berthelot mixing rules [START_REF] Xie | AlCoCrCuFeNi high entropy alloy cluster growth and annealing on silicon: A classical molecular dynamics simulation study[END_REF]. The local atomic structures were identified and characterized by calculating the radial distribution functions (RDF) and simulating the corresponding X-ray diffraction (XRD) patterns. The time evolution of the crystallization process was also addressed. Comparison with available experimental data confirmed the relevance of MD simulations to predict HEA thin film features.

Experiments:

Thin HEA thin films have been synthesized at the GREMI laboratory (Orléans, FRANCE) using DC magnetron sputtering technique with mosaic targets (patent no.WO/2008/028981). The experimental set up is widely described in reference 12.

Three targets are placed in the chamber in front of rotating Si (100) substrates in order to ensure large overlapping of the sputtered atom fluxes of pure elements (99.99%).

The six chemical elements are arranged on the targets (Fig. 1) [START_REF] Dolique | High-entropy alloys deposited by magnetron sputtering[END_REF]: Fe, Co and Ni are on target 1, Cu and Cr on target 2 and Al on target 3.

The relative surface fraction of each element on the target was calculated, taking into account the difference in sputtering yields estimated with TRIM software in order to reach the equimolar composition. However it is possible to adjust the stoichiometry of AlCoCrCuFeNi alloy in a narrow range, by varying the electrical powers applied to the magnetron cathodes independently. Note that in the present target configuration, the concentrations of elements present on the same target are linked.

Thin films of more than 1 μm thickness were synthesized on Si (100) substrates. X-ray diffraction (XRD) analysis (Cu Kα radiation, λ = 0.154 nm, Bragg-Brentano geometry) was performed to determine the crystalline phases and the stoichiometry of the thin films was determined by Energy Dispersive Spectroscopy (EDS). In this paper three samples will be under investigation: Al2Co9Cr32Cu39Fe12Ni6, sample 1, Al3Co26Cr15Cu18Fe20Ni18 sample 2 and Al39Co10Cr14Cu18Fe13Ni6 sample 3.

Simulation details:

MD simulations for the present HEA system were carried out in a three dimensional cell. Periodic boundary conditions were used in two horizontal x and y directions so that the simulated systems could be viewed as large planes in the x-y dimension. A free boundary condition was used in the vertical z direction to allow for surface growth.

The substrate is Silicon (100) with dimensions of 10a × 10a × 5a. a is the Si lattice parameter, i.e. 5.4308 Å. The velocities of the deposited atoms are sampled from a Maxwell-Gaussian distribution with a most probable energy of 1 eV, and with randomly selected incident angles. The mean kinetic energy of 1 eV was chosen in order to match at the best the deposition conditions in magnetron sputter technique at low argon pressure., it simulates sputtering experiments at low argon pressure where a small amount of buffered gas only randomizes sputtered atom velocities without not too much energy loss The deposition model of the six atoms is shown in Fig. 2. To prevent the shift of the computational system due to the impact of the adatoms with the top surface, the positions of four bottom atomic plane of atoms (depth a) were fixed during simulations. Besides, the system temperature will rise due to kinetic energy released by impinging atoms. To mimic isothermal growth conditions, an intermediate region above the fix region (depth 2a) is controlled by a Berendsen thermostat and maintained at a desired substrate temperature (here 300K). Atoms located in the surface region (thickness 2a) are left free to move and thus are able to respond to incoming atom impacts.

During MD simulation, 15000 atoms are released to the Si surface for realizing the three deposition sets. Table 1 gives the number of atoms really living in the simulated film and the corresponding composition. The sticking coefficient is found close to 0.9 and the compositions of the simulated samples are closed to the targeted ones.

After simulation, radial distribution functions (RDF) were used to employ the correlation between the atoms and the simulated powder X-ray diffraction patterns (XRDs) to analyze the phase of HEA thin film. The simulated XRD patterns were determined by using Debye functional analysis as described in detail by Kazakov and co-workers [START_REF] Bauchy | Structural, vibrational, and elastic properties of a calcium aluminosilicate glass from molecular dynamics simulations: The role of the potential[END_REF]. I is the intensity of the diffracted coherent radiation and is determined by the Debye formula.

𝐼 𝑘 (𝑏) = ∑ ∑ 𝑓 𝑛 𝑛≠𝑚 𝑛 (𝑏)𝑓 𝑚 (𝑏) 𝑠𝑖𝑛(2𝜋𝑏𝑟 𝑛𝑚 ) 2𝜋𝑏𝑟 𝑛𝑚
where b=2sin(θ)/ λ and λ is the wavelength of the incident radiation, 2θ is the scattering angle, and rnm is the distance between atoms n and m. The functions fn(b) and fm(b) are the scattering factors for atoms n and m, respectively.

Results and discussions:

Simulation snapshots for sample1-3 thin films are displayed in Fig. 3. It is clear that the atoms, in sample 1 and sample 2, are arranged in a crystalline structure. For sample 3, however, an amorphous structure is formed on the overall thickness.

For analyzing deposit properties, first simulated and experimental XRD patterns are compared (figure 4); second, we discuss the process of HEA thin film growth (Fig. 567); and third, we use solid-solution formation rules to characterize the collective behavior of these three samples (Fig. 8,9). It is important to note that it is not possible to directly compare intensity and width of experimental XRD peaks to that of simulated ones. The FWHM of XRD peaks depends on the size of crystalline domains, which of course are larger on 300 nm thick HEA films than on simulated films composed of a limited number of atoms (1500 in our case).

Thus the total number of atoms participating to the diffraction pattern is completely different. The relative intensities of the peaks are also related to the crystalline quality.

It appears from Fig. 4 that experimental S2 and simulated S3 samples exhibit a worth crystallinity than the others. This difference between experiment and simulation remains unexplained but may be attributed to the experimental conditions (energy flux transferred to the surface during the growth for example) that were not exactly during the synthesis of S1, S2 and S3 samples.

The growth process of the HEA thin film

Examination of the snapshots, as displayed in Fig. 5a-① to 5a-⑥, at different times,

shows the growth of Al2Co9Cr32Cu39Fe12Ni6 thin film deposited on Si (100) substrate.

At first step (t = 0.5 ns), it is found that small clusters are formed, which leave the substrate partly uncovered. As the number of deposited atoms increases, a continuous film is formed (from 2.5 ns). It can be clearly seen that the structure is amorphous on the first two snapshots of Fig. 5a, which is also evidenced on the RDF plotted in Fig. 5b. The RDF curve is relatively smooth and does not exhibit any defined peaks at expected bulk crystal neighboring positions shown in table 3. Those positions have been calculated from element crystal values found in reference 14 and are given here for equimolar composition. The RDF shape lies between a crystal and liquid indicating an amorphous structure. Going on with the deposition process, when more atoms are released onto the substrate, the RDF patterns change a little bit: the first peak (at 2.52 Å) becomes narrower, the second peak (above 4 Å) splits and a new peak at r ≈ 2.84 Å gradually arises. This indicates that the thin film changes from amorphous to crystalline structure. This phenomenon can also be observed in Fig. 5a. At ≈ 3.5 ns to 6.0 ns. First, grains with crystallized structure are formed on the top of the film (see Fig. 5a-). Then, new grains nucleate at locations close to the existing ones. Finally, this process propagates across the thin film leading to a single-crystalline structure (see Fig. 5a-).

It is found that the emergence and growth of a critical nucleus plays a key role in the crystallization of the HEA film alloy. The critical nucleus triggers a strong spatial correlation with other crystal nuclei in its surroundings just before the onset of crystallization, which may provide a subsequent rapid crystallization.

The snapshots of sample 2 deposited on Si(100) substrate at different stages are shown in Fig. 6a. It is obvious that sample 1 and sample 2 exhibit the similar final product of magnetron sputter deposition, i.e. a crystallized film (shown in Fig. 5a-⑥ and Fig. 6a-⑥). However, the drastic difference can be observed in the growth process of both alloys. For sample 2, at 2.5 ns, various crystalline clusters with same lattice orientation are formed on the Si(100) substrate (shown in Fig. 6a-), while for sample 1 a continuous amorphous film is observed (shown in Fig. 5a-). These results are consistent with the RDF analysis. It is seen that at 2.5 ns RDF patterns features of sample 2 are typical of a crystallized film: narrower first peak, appearance of a peak at 2.84 Å and splitting of the second peak. This implies that the film in sample 2 has the strongest ability of crystallization, which could only be correlated to the chemical composition. It seems thus, that the atomic size difference has significant effect on the film structure. After several crystalline clusters have been formed, they grow in size along the direction perpendicular and parallel to the Si(100) substrate. Interfaces between adjoining clusters of various orientations are formed. Finally, the different crystalline cluster join together. In order to reduce the energy of the system, the structure of the sample rearranged and the direction of the lattice rotated. This is depicted in the Fig. 6-③ -Fig. 6-⑥. In this simulation the energy of depositing atoms seems to be high enough to overcome the diffusion activation energy and lead to reconstruction the film as a whole with a new lattice orientation. Fig. 7 shows the snapshots and RDF of sample 3 at different times. Comparing with sample 1 and 2, the RDF at all times exhibits the typical features of an amorphous structure, where the first and second peak are broad and higher-order peaks fade away quickly. This is also clear on the corresponding snapshots given in Fig. 7(a). All these points indicate that sample 3 grows as an amorphous film.

RDF features can be explained by the close-packed hard ball model. The structure stabilized in multi-element alloys can be attributed to the effect of mixing numerous elements with different atomic sizes as depicted in Fig. 8. The circles represent the first, second, third, fourth and fifth shells, respectively. It is seen that a crystal structure will have many sharp RDF peaks corresponding to the first, second, third, fourth, fifth shells and so on (Fig. 8(a)).

In an amorphous structure, because of the most compact system rules, the first shell will be stable (only a change due to the presence of neighbor atoms with different atom radius could be detected). While, the second and third shells, and also the fourth and fifth shells, being quite close in space, it is expected that a certain degree of shell diffusion renders them indistinguishable from each other (Fig. 8(b)) [START_REF] Kao | Rapidly solidified structure of alloys with up to eight equal-molar elements-a simulation by molecular dynamics[END_REF]., Thus an amorphous structure exhibits a split second peak reflecting the partial merging of the second and third shells, and the third broad peak reflecting the merging of the third and fourth shells.

Following the hard ball model, atomic radii and deviation from the mean radius for each element in the system are shown in table 4. The deviation of the mean distance of Al is much higher than other element. Considering sample 3 has the largest ratio of Al (39%), it is more easily for the second and third shells to be distinguishable and behaved as a liquid-like structure or amorphous structure (Fig. 8(b)).

Solid-solution formation rules for HEA

Currently, the phase stability of high entropy alloys has been widely investigated [START_REF] Yang | Prediction of high-entropy stabilized solid-solution in multicomponent alloys[END_REF][START_REF] Troparevsky | Criteria for Predicting the Formation of Single-Phase High-Entropy Alloys[END_REF][START_REF] Ye | A geometric model for intrinsic residual strain and phase stability in high entropy alloys[END_REF][START_REF] Guo | Solid solutioning in equiatomic alloys: limit set by topological instability[END_REF] Among them Zhang et al [START_REF] Yang | Prediction of high-entropy stabilized solid-solution in multicomponent alloys[END_REF] proposed that two parameters can be used to characterize the collective behavior of the constituent elements in multi-component The parameters of the three samples grown by MD simulations are calculated, listed in table 5 and mapped in Fig. 9. The values of δ (2.28 %, 2.89 %) and  (2.794, 38.986) of sample 1 and 2 lies in the area where the formation of solid solutions is predicted.

(
Sample 3 with δ = 5.94% and  = 2.595, stays in the overlap area of "Solid Solutions" and "Intermetallics", which indicates both solid solutions and intermetallics are likely to be formed. Thus, solid solution could be predicted to form in all samples, and ordered intermetallics could also appear in sample 3, for which parameters tend to enter the metallic glass zone. However, the simulation results shown in Fig. 3 evidence that, crystalline structures are just observed in sample 1 and 2, and do not appear in sample 3. The sole amorphous structure formed in sample 3could be ascribed to a kinetic factor.

It is well known that deposition by magnetron sputter technique could lead to the formation of metastable phases. With no external heating of the substrate, low ordered phases are usually obtained. The solid solution formation rules as defined from ,  graph, is based on thermodynamics consideration, the kinetic effect is not emphasized, and non-equilibrium phase formation cannot be predicted.

Thus, the rules proposed by ZHANG apply well for conventional bulk alloy equilibrium synthesis techniques.. In the case of thin film growth, expecially when deposition conditions deviate from thermodynamic equilibrium, they could only give some expected trends. For a more precise prediction of the expected phases, simulations of the deposition process, taking into account kinetec effects, as those highlighted by MD simulations, have to be implemented

Conclusions:

MD simulations were carried out for describing the initial growth of AlCoCrCuFeNi HEA thin films. From the simulation results, it has been shown that the number of elements, atomic size difference have significant effects on the atomic configuration.

Crystalline layer grows for low Al content (2 % and 3 %), while, pure amorphous thin film is observed for the highest proportion (39 %) of the atom which radius deviates the most from the mean one calculated for the equimolarity. By calculating the parameters δ and Ω for the reported HEAs, a tendency to evolve from solid solution to bulk metallic glasses was expected, in good agreement with the present MD simulation results.

4. 1

 1 Determination of Simulated Powder X-ray Diffraction Patterns. The experimental XRD θ-2θ peaks are shown in figure 4(a). The following peak positions are observed: 2θ = 43.81 for sample 1 attributed to the fcc (111)  structure, 2θ = 44.41 for sample 2 attributed to the bcc (110)  structure and 2θ = 43.77 for sample 3. However, the broad and low intensity peak in sample 2 could indicate a mixture of both solid solutions and a less ordered phase.The calculated X-ray intensities vs. 2θ of simulated HEA thin films are shown in figure4(b). By comparing both XRD plot sets, it can be observed that the peak positions in the simulated XRD are close to that detected in the experimental XRD patterns (see table2). The simulated peaks at positions 2θ = 43.97•, 44.57• and 42.50• correspond to samples 1~3 respectively, sample 2 being close to the bcc(1 1 0)α structure while sample 3 peak is fairly consistent with fcc (1 1 1) β[START_REF] Xie | AlCoCrCuFeNi high entropy alloy cluster growth and annealing on silicon: A classical molecular dynamics simulation study[END_REF]. In this latter case, because the peak is broader than for Sample 1 and 2, the structure has a much shorter range order, in agreement with the amorphous-like snapshot of this sample While, the peak intensities are very different between simulations and experiments.

  more than 4 elements) alloys: δ (δ = √∑ 𝑐 𝑖 (1 -𝑟 𝑖 /𝑟̅ )2 𝑛 𝑖=1where ci is the atomic percentage of the i th component, and 𝑟̅ is the average atomic radius and ri is the atomic radius), describing the comprehensive effect of the atomic-size difference in the nelement alloy and Ω(Ω = 𝑇 𝑚 Δ𝑆 𝑚𝑖𝑥 /|Δ𝐻 𝑚𝑖𝑥 |, where ΔHmix is enthalpy of mixing, ΔSmix is entropy of mixing and Tm is the average melting temperature for n-elements alloy), describing the effect of ΔHmix and ΔSmix on the formation of multi-component solid-solutions. From the values of parameters δ and Ω for multi-component alloys, a solid-solution formation rule has been proposed by Zhang et al. When   1.1 and  ≤ 6.6%, the multi-component solid solution phase could form. In contrast, for larger value of  and smaller value of , intermetallic compounds and bulk metallic glass(BMG) are formed. The value of  for BMG is smaller than that predicted for intermetallic compounds.

Acknowledgement: This work was supported by the China National Key Basic Research Program 973 (Grant No.2013CB934004), The authors would like to thank FEDER (Fonds Européen de Développement Régional) and the "Région Centre" for financial support, ZY acknowledge the financial support of National Natural Science Foundation of China (Grant No. 51471025 and 51210105006).


Sébastien Demange

Slim Ouni

An episodic memory-based solution for the acoustic-to-articulatory inversion problem

Keywords: numbers: 43.70.Bk, 43.70.Jt, 43.70.Aj

This paper presents an acoustic-to-articulatory inversion method based on an episodic memory. An episodic memory is an interesting model for two reasons. First, it does not rely on any assumptions about the mapping function but rather it relies on real synchronized acoustic and articulatory data streams. Second, the memory inherently represents the real articulatory dynamics as observed. It is argued that the computational models of episodic memory, as they are usually designed, cannot provide a satisfying solution for the acoustic-to-articulatory inversion problem due to the insufficient quantity of training data. Therefore, an episodic memory is proposed, called generative episodic memory (G-Mem), which is able to produce articulatory trajectories that do not belong to the set of episodes the memory is based on. The generative episodic memory is evaluated using two Electromagnetic Articulography (EMA) corpora: one for English and one for French. Comparisons with a codebook-based method and with a classical episodic memory (which is termed concatenative episodic memory) are presented in order to evaluate the proposed generative episodic memory in terms of both its modeling of articulatory dynamics and its generalization capabilities. The results show the effectiveness of the method where an overall root-mean-square error of 1.65 mm and a correlation of .71 are obtained for the G-Mem method. They are comparable to those of methods recently proposed.

I. INTRODUCTION

Recovering the vocal tract shape from speech acoustics could allow a number of breakthroughs in automatic speech processing. For instance, the location of critical articulators ? could be exploited to better characterize a given phoneme, which may improve speech synthesis ? and recognition ? ? ? ? . Indeed, articulatory features vary much more slowly than speech acoustic features, and thus they should be more robust than acoustic parameterizations, especially in noisy environments ? .

Recently, databases of synchronized acoustic and articulatory data streams, using electromagnetic articulography (EMA) for instance, have become available. These corpora enable machine learning algorithms to perform acoustic-to-articulatory regression. The main techniques use support vector machines ? , Gaussian mixture models ? , hidden Markov models (HMM) ? ? , or artificial neural networks ? ? . These methods try to learn the acoustic-toarticulatory mapping function, which is known to be highly non-linear, non-unique (i.e., different vocal tract shapes producing the same acoustics) and thus difficult to model ? ? ? ? ? .

As an alternative, codebook-based approaches ? ? ? ? ? make very few assumptions about the mapping function, but rather rely on a collection of pairs of acoustic-articulatory data.

We believe that the main difficulty of inversion is the lack of a good representation of the dynamics. The non-uniqueness problem ? will very likely vanish if the dynamics are fully integrated within the inversion methods. In fact, Qin & Carreira-Perpiñán ? , and to some extent Neiberg et al. ? , argued that natural human speech is produced with a unique vocal tract configuration and there are few cases of non-uniqueness. Phonetic context naturally imposes constraints on the vocal tract related to coarticulation. Effective modeling of articulatory dynamics seems essential to solving the inversion problem. These dynamics can be modeled with HMMs ? ? and neural networks ? ; they can also be inferred from time derivative features ? . For the codebook-based methods the dynamics are not modeled at all. Instead, continuity constraints are used during inversion. However, despite these constraints, the recovered articulatory trajectories show many discontinuities and need to be smoothed with signal processing techniques.

In this paper, we present a new data-driven approach based on the concept of episodic memory ? . An episodic memory can be considered to be a codebook that includes a temporal dimension. While a codebook models the relationship between two static observations, an episodic memory can model the relationship between two sequences of observations. The main advantage of episodic memory is that it keeps track of the order of the observations and thus preserves the acoustic and articulatory dynamics of each episode.

An episodic memory can deal with the non-linearity of the mapping function by using the one-to-one correspondence between the synchronized acoustic and articulatory observations. It can deal with the non-uniqueness of the mapping by exploiting the articulatory dynamics encoded through the time ordering of both the episodes and the observations they are comprised of. In addition, the need to smooth inferred articulator trajectories is greatly reduced compared to current codebook-based inversion methods.

In the following sections, we first present our motivations for developing a generative as opposed to a concatenative episodic memory. Then, we present our inversion method based on the generative memory model. Finally, inversion results are presented and compared to other methods that have been reported previously.

II. EPISODIC MODELING FOR SPEECH INVERSION

In speech processing, the episodes, which are the units of an episodic memory, are individual acoustic realizations of predefined linguistic units (e.g. phones, diphones, syllables or words). A similar approach has been used in speech recognition ? , referred to as templatebased or exemplar-based speech recognition, and in unit selection speech synthesis ? ? . The memories used for these two problems are concatenative (C-Mem) because the episodes are indivisible and the recognized or synthesized sentences are always comprised of a concatenation of episodes.

We apply an episodic memory model to the acoustic-to-articulatory mapping problem, whereby episodes comprise the synchronized acoustic and articulatory feature sequences for each linguistic unit. Then, for inversion, the memory could exploit precisely the acousticarticulatory relationship, as well as real acoustic and articulatory dynamics.

Speech inversion differs from speech synthesis and recognition because the mapping is between two continuous spaces, while for speech recognition and synthesis, the target or the source space, respectively, is discrete. This is an important issue for an episodic memory model because the memories need to contain many episodes of each linguistic unit in order to achieve adequate coverage of the variability present in speech. Usually, several hours of acoustic speech are required to reach acceptable performance when episodic memory is used for synthesis or recognition. Furthermore, keeping in mind that the articulators can compensate for one another during speech (a phenomenon which introduces non-uniqueness to the solution), an episodic memory intended for speech inversion should be even larger to ensure good coverage of both the acoustic and articulatory variability. This requirement cannot be fulfilled, as almost all articulatory corpora currently available contain less than an hour of speech, usually much less. The biggest corpus, the EMA part of the mngu0 dataset ? , made available very recently and containing 1263 sentences (1 hour 27 minutes of speech) is still probably too small. We think that this lack of data is the main reason why episodic modeling has not yet been used for this problem.

To address this problem, we propose to provide the memory with a mechanism to simulate many more episodes than the ones it contains. We define an episode as synchronized acoustic and articulatory realizations of a linguistic unit (LU ). Let us consider two episodes X and Y of a given linguistic unit. X and Y are almost identical, differing only at the beginning and end, due to coarticulation effects. A C-Mem which does not contain any episode of LU whose left and right contexts are those of X and Y , respectively, will invariably fail to invert acoustic realizations of LU in this context. However, the memory could perform better if it were allowed to go through the first part of X, then to switch to Y at any time during the central, nearly identical part, and finally to go through the final part of Y .

Even though all episodes of a given linguistic unit are not identical, they can exhibit local articulatory similarities. Therefore, we propose to allow the memory to switch between any two episodes X and Y of the same linguistic unit during the inversion at times when observations of X and Y are similar. Care will be taken to produce realistic articulatory dynamics by considering similarities with regard to temporal alignments and with regard to the positions of the articulators. As the proposed memory will be able to generate episodes other than the ones in the database, we will refer to this memory as a generative memory (G-Mem).

III. GENERATIVE MEMORY BASED INVERSION

As before, we define an episode as synchronized acoustic and articulatory realizations of a particular linguistic unit, in this case, phonemes. The phoneme identity will be referred to as the class of the episode. In the following, we consider two given episodes X and Y .

X is a particular realization of a given phoneme expressed as a sequence of K articulatoryacoustic observations X = (x 1 , . . . , x K ). Each observation

x i = (x art i , x ac i ), where i ∈ [1..N ],
is a synchronized pair composed of an articulatory observation x art i and its corresponding acoustic observation x ac i . The scalar articulatory observation x art i can be a given articulator description or even the x-or y-coordinate of a sensor glued onto an articulator as used in our work. Similarly, Y is another realization of the same phoneme expressed as a sequence of N observations: Y = (y 1 , . . . , y N ), where y i = (y art i , y ac i ).

A. Local articulatory similarity

Dynamic Time Warping (DTW) ? is a general algorithm to find the shortest distance D(X, Y ) between two episodes X and Y , which may vary in length. The episodes are warped non-linearly in order to minimize the effects of their temporal variability. Any given mapping leads to a particular alignment path Φ = (Φ 1 , . . . , Φ M ), where M is the number of alignments. Φ i = (Φ x,i , Φ y,i ) is the i th observation pairing along Φ with Φ x,i and Φ y,i as the indices in X and Y of the aligned observations. The distance between X and Y given Φ is the sum of the local distances d(., .) between the aligned observations along Φ. The choice of this local distance d(., .) depends on the nature of the observations used to perform the mapping. Here, the articulatory observations are used because we are focusing on producing realistic articulatory trajectories. In our work, the articulatory observations are the x-and y-coordinates of sensors glued onto articulators in the midsagittal plane. Thus, Euclidean distance is used.

D(X, Y

) is the shortest distance over all Φ:

D(X, Y ) = arg min Φ M i=1 d(x Φ x,i , y Φ y,i ) (1)
Many variations of the algorithm have been proposed ? in order to prevent degenerate paths from occurring. In this work, we applied Itakura's constraints ? . These constraints make the DTW asymmetric, such that each observation in sequence X is aligned with exactly one observation in Y . The mappings of many episodes onto X result in many alignment paths of the same length (equal to the length of X). Therefore, the distances from these episodes to X can be fairly compared and ranked. The other Itakura constraints impose bounds on the temporal deformation, in order to preserve a certain temporal consistency of the aligned observations. Let X art be an articulatory trajectory of an episode X expressed as a sequence of K articulatory positions (x art 1 , x art 2 , . . . , x art K ). We define each articulatory observation x art i+1 as the natural articulatory target (local target) of x art i since it has been observed to follow x art i . In fact, x art i+1 is a specific articulatory position, but we can suppose it could have been slightly different. Indeed, starting from x art i at time i, the articulators could have reached a different position at time i + 1 close to x art i+1 with no significant consequences to the acoustics. Then, for each x i we define an articulatory target interval AT I x i as:

AT I x i = [x art i+1 -δ, x art i+1 + δ] (2)
where δ is a given positive value.

We consider any articulatory position y art j to be similar to any articulatory position x art i , if y j is aligned with x i when mapping Y onto X, and if y art j belongs to the articulatory target interval of x art i-1 .

B. Building the generative episodic memory

We model a G-Mem as an oriented graph G G-M em . The nodes are the pairs of synchronized acoustic and articulatory observations comprising the episodes. The oriented edges indicate the allowed transitions between the articulatory positions the memory can follow during the inversion process. They are created according to Algorithm 1.

Require: Let Γ be the set of all episodes

for all X = (x 1 , . . . , x K ) and Y = (y 1 , . . . , y N ) in Γ, class(Y ) = class(X) do if D(X, Y ) ≤ ∆ * length(X) then for all i = 1 to K -1 do if y art Φ y,i+1 ∈ [x art i+1 -δ, x art i+1 + δ] then Add an edge from x i to y Φ y,i+1 end if end for end if end for Algorithm. 1 Building the graph G G-M em
For any pair of episodes X and Y of the same class we create an oriented edge from a given x i to a given y j if y j is similar to x i+1 from an articulatory point of view (i.e, y art j is similar to x art i+1 ) as defined previously:

Φ y,i+1 = j (3) y art j ∈ AT I x i = [x art i+1 -δ, x art i+1 + δ] (4)
In addition we impose that the asymmetric distance D(X, Y ) falls below a given threshold (proportional to the length of X):

D(X, Y ) ≤ ∆ * length(X) (5)
The goal is to prohibit the memory from switching between episodes, which are globally very different (from an articulatory point of view) because it could lead the memory to produce unrealistic articulatory trajectories. As an example, consider the movements of the tongue tip, which might rise or fall during the production of a given phoneme. Although these trajectories are very different, it is likely that the tongue tip can reach similar positions midway through the fall and rise. Combining the fall and rise could possibly lead to a degenerate trajectory.

At the episode boundaries the memory is only subject to the articulatory continuity requirement expressed by equation ( 4). Let Z = (z 1 , z 2 , . . . , z P ) be the episode, which was observed after X. Then, an edge from x K to the first observation w 1 of any episode W of

any class is created if w art 1 ∈ AT I x K = [z art 1 -δ, z art 1 + δ].
If the episode X is the last of a record, its natural articulatory target is unknown and equation (4) cannot be satisfied;

thus no edge to any other episode is possible. Note that a C-Mem only accounts for these transitions between episode boundaries. So, a C-Mem can be seen as a particular case of G-Mem for which ∆ is set to zero. The bottom graph shows the resulting episodic memory with all possible transitions from X to Y and from Y to X. The resulting transitions within the memory are represented by the plain arrows from any x i to x i+1 and from any y j to y j+1 resulting from the mappings of X and Y onto themselves (not shown in the figure). The dotted arrows from any x i to y j result from the mapping of Y onto X while those from any y j to x i result from the mapping of X onto Y .

C. An example

While a C-Mem would have contained only two articulatory trajectories X art and Y art , the resulting G-Mem can simulate six additional trajectories combining X and Y . Thus, the G-Mem is able to produce the following 8 articulatory trajectories according to a given input acoustic signal to be inverted:

1. x art 1 → x art 2 → x art 3 → x art 4 (X art ) 2. x art 1 → x art 2 → x art 3 → y art 5 3. x art 1 → y art 3 → y art 4 → y art 5 4. x art 1 → y art 3 → y art 4 → x art 4 5. y art 1 → y art 2 → y art 3 → y art 4 → y art 5 (Y art ) 6. y art 1 → y art 2 → y art 3 → y art 4 → x art 4 7. y art 1 → x art 2 → x art 3 → x art 4 8. y art 1 → x art 2 → x art 3 → y art 5
As explained in section III.B, note that when using the G-Mem, all the inverted trajectories start from the beginning of an episode and finish at the end of an episode.

As the nodes of the oriented graph G G-M em are bimodal (composed of an acoustic and an articulatory observation), each path within the graph corresponds to a particular articulatory trajectory and also to the acoustics that would have been produced by the articulatory trajectory. Thus, acoustic-to-articulatory inversion is performed by searching the path within G G-M em that best matches the acoustic speech signal to be inverted. The estimated articulatory trajectory is extracted from the articulatory observations of the visited nodes along this path.

All search paths can start only at nodes that represent the first observation of an episode.

During inversion a breadth-first search is performed, applying the Viterbi algorithm. That is, at each step, the K-best paths obtained at the previous step are propagated through the G G-M em along the oriented edges defining allowed articulatory movements. The K-best propagated paths are kept while the others are discarded and the process is repeated up to the end of the speech signal.

The winning path is the one with the best acoustic score selected from all paths ending at a node that corresponds to the final observation of an episode. The score of each path is expressed as the sum of acoustic distances between the speech frames and the acoustic observations of the visited nodes along the path, computed on a predefined acoustic window W.

IV. DATA AND INVERSION EXPERIMENTS

A. Corpora

All the experiments presented in this work were carried out on the following two corpora of synchronized acoustic speech signal and articulatory trajectories.

a. MOCHA We used the EMA corpus of MOCHA ? . Two speakers of British English, one female (fsew ) and one male (msak ), were recorded reading 460 short phonetically balanced British-TIMIT sentences. The audio is provided as waveforms sampled at 16 kHz, and each EMA receiver position is given as 2D coordinates in the mid-sagittal plane. We used seven sensors located at the lower incisors (li), upper lip (ul), lower lip (ll), tongue tip (tt), tongue body (tb), tongue dorsum (td), and velum (vl). The phonetic segmentation provided with the audio stream was used.

b. mdem We have recorded a French EMA corpus using an articulograph (AG500, Carstens Medizinelektronik). A male French speaker (mdem) was recorded uttering 400 phonetically balanced sentences. The audio is provided as waveforms sampled at 16 kHz, and the EMA data consist of 2D coordinates in the mid-sagittal plane. We used 6 sensors fixed in the midsagittal plane on the lower lip (ll), upper lip (ul), tongue tip (tt), tongue body (tb), tongue dorsum (td), and tongue back dorsum (tbd). The phonetic segmentation was obtained from a word-level transcription of the sentences, a dictionary containing several pronunciation variants for each word and a set of French monophone HMMs trained on several hours of speech and adapted to mdem's voice. The segmentation was obtained by force aligning the phone HMMs onto the acoustics given the sentence word transcription and the pronunciation dictionary.

Each corpus was split into training, development and test sets. For MOCHA, care was taken that the selected utterances for each set corresponded exactly to the ones used by Richmond ? , as this split was also used in Toutios & Margaritis ? and Zhang & Renals ? .

Information about the different sets are given in table I. Note that the durations only account for usable speech (without the start and end silences). Figure 2 shows the distributions of the articulatory samples for each speaker and coil.

B. Feature extraction

The silences occurring at the beginning and end of each recording were first discarded, as the articulators can move unpredictably during such intervals. A Linear Predictive Analy-sis ? was performed on the speech signal using the HTK toolkit ? . 12 cepliftered MF-PLPs ?

and the logarithmic energy of the signal comprised the acoustic feature vector extracted from every 25 ms speech frame shifted by 10 ms.

The articulatory data were first down-sampled from 500 to 100 Hz (for MOCHA) and from 200 Hz to 100 Hz (for mdem) to match the acoustic frame-shift. Then, all trajectories were low-pass filtered, in order to remove the recording noise, using a cut-off frequency of 20

Hz. A final data adjustment was performed in order to take into account the observed drift in the global trajectories throughout the recording sessions (see Richmond ? for details).

According to Richmond ? , these variations can reflect speaker compensation/adaptation to the presence of the coils within the mouth, or may be due to measurement errors. We should note that the true underlying cause of these apparent long-term inconsistencies has not been properly and deeply investigated. In our study, we applied the same procedure proposed by Richmond ? to remove these very long-term variations by subtracting a low-pass filtered version of the trajectory means from the EMA data.

C. Quality measurements

The quality of the recovered articulatory trajectories was evaluated using root-meansquare error (RMSE) that quantifies the difference between the measured x art i and estimated f (x ac i ) articulatory coordinates:

RM SE = 1 N N i=1 (f (x ac i ) -x art i ) 2 (6)
as well as Pearson's correlation between the values in the two trajectories. This is obtained by dividing the covariance of the two trajectories by the product of their standard deviations.

Cor = N i=1 (f (x ac i ) -f (x ac )).(x art i -x art ) N i=1 (f (x ac i ) -f (x ac )) 2 . N i=1 (x art i -x art ) 2 (7)

D. Experiment design

We have implemented three different inversion methods: a codebook-based approach, as described in Suzuki et al. ? , and two memory-based approaches (C-Mem and G-Mem) as described above. We have chosen to compare the memory-based approaches to this codebook-based approach, as they differ only in the manner of inferring the dynamics of the recovered articulatory trajectories. While the memories can model and use observed trajectories, the codebook only relies on continuity constraints.

In fact, this codebook method consists in looking up a set Γ i of the N best entries in the codebook for each given acoustic signal sample Y i . The best entries are the ones, which minimize the acoustic distance to Y i . In fact, this distance is the average acoustic distance computed over a window of a predefined length W. The articulatory trajectory is obtained by looking for the path through Γ i that minimizes the weighted sum of the acoustic distances and squared distances between subsequent articulatory parameters. The reader may refer to Suzuki et al. ? for details of the implementation.

We consider the recovery of the articulatory trajectories for each coil, and along both the

x and y axes, as independent inversion problems. Thus, the experiments presented here consist of fourteen (for fsew and msak ) and twelve (for mdem) distinct inversion problems. For each inversion problem, a dedicated codebook, C-Mem and G-Mem are built and optimized.

The parameters of the codebook are the length of the spectral window, and the weight of the articulatory constraints with regard to the acoustic distances. During the inversion, the 1000 best codebook entries are considered at any given time. For G-Mem, the parameters to be set are δ, the half ATI length, and ∆, the maximum allowed articulatory distance between two episodes X and Y for allowing the memory to switch from one to the other.

Note that the parameter ∆ is equal to zero for a C-Mem, so that each episode can only be combined with itself. An acoustic window was used to compute the acoustic distances. The length of this window is a parameter to be optimized. Euclidean distance is used for both acoustic and articulatory distance calculation.

Two types of experiment have been carried out: with and without phonetic constraints.

Note that the G-Mem used phonetic segmentation only for building the memory, and the segmentation is not needed for inversion. The purpose of adding experiments with phonetic constraints is to ascertain what extra information is added for use of the phonetic knowledge to the different inversion methods. The reference phonetic segmentation is obtained by force aligning HMM acoustic models onto the speech signal according to the transcription of the uttered words. In the phonetically constrained mode, the estimated articulatory position at time t has to come from a codebook entry or an episode of the same phoneme as the one indicated by the segmentation at time t.

For the three systems, the parameters have been jointly optimized through a grid search optimizing RMSE on the development set. Since the articulators move in different ranges, different sets of parameters were obtained for all inversion problems. Globally, the length of the acoustic window W is approximately 150 ms for the codebook and 90 ms for the memories. δ ranges from a few hundredths of a millimeter (for the velum) to at most one millimeter, and ∆ ranges from a tenth of a millimeter to two millimeters for different articulators. Finally, the memory search beam width is set to 10 000 for the memories and 1000 for the codebook.

V. INVERSION RESULTS AND ANALYSIS

A. Results Figure 4 shows the tongue tip movements (thick curves) along the vertical (up/down) axis recovered by each of the three approaches from a two second long speech signal corresponding to the French sentence "juste quelques extrémités de branches gelées" (only a few frozen branch tips). For each of the three graphs, the reference trajectory is provided as the thin curve and the estimation errors are emphasized as filled areas between the estimated and reference trajectories. Though optimized, the dynamic articulatory constraints of the codebook-based approach do not prevent discontinuities, as the recovered articulatory trajectory is very jerky compared to the smooth and continuously varying reference trajectory.

The results obtained by both the C-Mem and G-Mem are visibly better. One might have expected smoother trajectories using the C-Mem, as the result is expressed as a concatenation of natural articulatory episodes. Through deeper analysis of the decoding paths within the C-Mem, we have noticed that most of the time the C-Mem does not contain episodes, which acoustically match the test signal well enough. In order to counteract this lack of good episodes, the C-Mem tends to select many short episodes. Indeed, the sum of the acoustic distances of short episodes, which locally match the test signal well, is usually smaller than the acoustic distance of a longer episode with partial acoustic mismatches. Finally, the G-Mem succeeds in estimating the articulatory movements accurately. The combination of episodes significantly reduces the estimation error. Furthermore, the articulatory dynamics embodied in the G-Mem transition graph contribute to the naturalness of the resulting articulatory movements.

B. Statistical significance

We have conducted a statistical analysis in order to quantify the confidence in our results

and also in the improvements. We have applied a bootstrap method proposed in Bisani & Ney ? that was originally designed for speech recognition performance evaluation, but which is equally applicable here. This method relies on a bootstrap replication by creating N pseudo data sets X i from the original data set X with replacements.

A pseudo data set can thus contain several or no samples of the original data set.

The statistics Stat i are then calculated over each pseudo data set X i . For a large N the distribution of the statistics is approximately Gaussian and thus, the true statistic lies with 99% confidence within the interval Stat i ± 3σ, where Stat i and σ are the mean and standard deviation of the statistics distribution, respectively. See Bisani & Ney ? for more details on this statistical analysis.

We have computed the 99% confidence intervals of the RMSE, generating 100 000 pseudo data sets for each experiment. These intervals are reported in figure 3 as error bars on each bar. These confidence intervals only represent the range of performance we would obtain when applying one of the three methods on a new data set of one of the three speakers.

In addition, table II gives the probabilities of improvement, as well as the 99% confidence intervals of the expected improvement of each method over the others, for the three corpora.

Each cell contains an integer, which is the probability that the method corresponding to the row outperforms the method corresponding to the column, as well as the 99% confidence interval of the expected improvements, expressed in percent. These numbers confirm the performance boost of the memory-based approaches over the codebook-based method and the superiority of the G-Mem over the C-Mem.

C. Smoothness and naturalness of articulatory trajectories

We argued in the introduction that the episodic memories are well-suited models for articulatory inversion, as they are able to preserve the dynamics of the episodes. We offer here an analysis that gives more insight on this point. Many proposed solutions to the inversion problem include low-pass filtering the recovered articulatory trajectories. This filtering aims to remove rapid changes (usually of low amplitude) that come from errors during the acoustic-to-articulatory mapping, and that can be considered noise. Indeed, all articulators have their proper velocity and can move more or less quickly, but they all move continuously. Comparing the recovered trajectories with their smoothed version is therefore one possible way of assessing how well the articulatory dynamics have been approximated.

For each coil, in both directions, we have determined the best cut-off frequency minimizing the RMSE on the test set. Figure 5 summarizes the best case scenario smoothing results.

As in Figure 3, the bars represent the RMSE, and Pearson's correlations are provided below the bars. The numbers in square brackets are the relative percentage improvements over the non-smoothed trajectories.

The most obvious effect is that the codebook-based method significantly benefits from this filtering. Its RMSE improves by approximately 10%, while the improvements for the memory-based approaches do not exceed 3%. The same trend can be observed for the Pearson's correlations. This indicates that the memory-based approaches really do take advantage of the observed dynamics of the episodes. Applying articulatory continuity constraints during the inversion does not yield the same benefits. Note that similar observations have been reported in Toda et al. ? . Indeed, the authors proposed an MLE-based mapping that accounts for correlation between frames. They reported significant improvements over their baseline GMM-based mapping, but also showed that the low-pass filtering effect became negligible. A relative improvement of 9.36% for RMSE was obtained over the baseline GMM-based mapping with a low-pass filter, but only .72% over the MLE-based mapping.

D. Computational resource requirement

Let N s be the number of samples contained in the training set, K s be the search beam width, and B be the average branching factor of the G-Mem, i.e. the average number of transitions allowed from any sample.

At any time t, the codebook method needs to compute the acoustic distance between the current test sample and all the N s training samples. Then, the K s best training samples are selected and the dynamic articulatory constraints are computed between each of them and each of the K s best hypotheses computed at time t -1. The complexity O codebook is then:

O codebook = N s + (K s × K s ) ≈ K 2 s (8)
Using a beam width of 1000 leads to one million dynamic articulatory constraints computed for each test sample.

At any time t, the G-Mem needs to propagate all of the K s best hypotheses computed at time t -1. Setting the average branching factor to B, the complexity O G-M em is:

O G-M em = K s × B (9)
The average branching factor B is about 250 for the G-Mem. Using a beam width of K s = 10 000 leads to 2 500 000 paths to be investigated. Actually, these paths end in a subset of all the N s acoustic/articulatory frames contained in the memory. Thus, the maximum number of acoustic distances to be computed each time is N s . Unlike the codebook approach, propagating the path here does not require any distance computation.

So, the two different beam widths used for the codebook-based and G-Mem-based methods roughly lead to the same computational resource requirements.

Figure 6 shows the evolution of the running time for both the codebook-based and G-Mem-based methods on the msak test set as a function of the search beam width. The plot is logarithmically scaled along both axes. As stated in equations ( 8) and ( 9), we can verify that the complexity is exponential for the codebook and linear for the G-Mem with respect Varying the beam width from 1000 to 8000 does not lead to a significant RMSE improvement.

Therefore, a beam width of 1000 appears to be a good compromise between performance and execution time.

VI. DISCUSSION

We have proposed an episodic memory solution for acoustic-to-articulatory inversion.

This model has some similarities with a codebook; in particular, it relies on a one-to-one correspondence between articulatory and acoustic observations. However, unlike a codebook, an episodic memory accounts for the temporal dimension and is thereby able to preserve the articulatory dynamics of the episodes. We have also proposed an algorithm, which allows the memory to combine different episodes during the inversion to simulate many other episodes than the ones in the database. Through the presented experiments we have shown that the trajectories produced by the memory-based models are better than those produced by the codebook.

In addition, the estimation errors using a G-Mem are very encouraging compared with the state of the art. Hiroya & Honda ? reported RMSEs of 1.50 and 1.73 mm with and without phonetic segmentations, respectively, using a HMM-based production model. However, we cannot directly compare our results with theirs, as they used a Japanese database. On MOCHA, Toda et al. ? used Gaussian mixture models to map the acoustic space onto the articulatory space. They reduced the RMSE from 1.58 to 1.40 mm by applying a maximum likelihood estimation (MLE) of the dynamic features. Zhang & Renals ? obtained an RMSE of 1.71 mm using a trajectory HMM. They included velocity features in their acoustic front end and performed speech recognition prior to inversion to provide their system with a phonetic segmentation. Even without phonetic segmentation, the G-Mem performs slightly better. Moubayed & Ananthakrishnan ? proposed a memory-based method. They used a linear regression on the local neighborhood of the codebook entries to map the acoustic input frames onto the articulatory space. They also used MLE of the dynamic features to improve the results. An RMSE of 1.52 mm was reported using this method. Finally, Richmond ?

reported an RMSE of 1.40 mm using trajectory mixture density neural networks. We share with Zhang & Renals ? , Richmond ? and Moubayed & Ananthakrishnan ? the same train, dev and test sets and all reported RMS errors range from 1.40 to 1.73 mm on this corpus.

We can claim that the proposed G-Mem, with an RMSE of 1.65 mm, performs as well as the machine learning based approaches.

The proposed G-Mem can be improved in the following direction. A local linear regression as proposed by Moubayed & Ananthakrishnan ? could further improve our results, as the G-Mem can produce unseen trajectories but it is unable to precisely map an acoustic frame onto the articulatory space if this acoustic frame does not belong to the memory.

Many studies have shown significant improvements using the reference phonetic segmentation of the acoustic signal to be inverted. This suggests that the inversion takes advantage of phonetic knowledge: if the phonetic content of the speech signal to be inverted is available, the articulatory movements can be recovered more precisely. As for the HMM-based method ? , a first phone recognition pass could be performed with well-trained acoustic HMMs in order to provide the G-Mem with an accurate phonetic segmentation.

As an alternative, a language model could be used to rescore the paths within the memory during the inversion. A dictionary could further constrain the search paths and decrease the RMSE to reach results similar to the ones we have obtained using the reference phonetic segmentation.

Finally, we define the episodes as the acoustic and articulatory realizations of particular phonemes. That is, we implicitly hypothesize that the articulatory trajectories can be segmented into nonoverlapping elementary units that are phones, and that this segmentation is the same for all studied articulators. However, Ananthakrishnan & Engwall ? recently proposed an automatic method for segmenting the articulatory movements into articulatory gestures. This method accounts for the notion of critical articulators. That is, the production of a particular phoneme depends mainly on the movements of few articulators (the critical articulators), which have to reach precise articulatory target positions. The other articulators can move more freely and thus can anticipate the production of the next phoneme to be uttered. Then, the movements of our articulators overlap with each other and contribute to coarticulate all the uttered sounds. We believe that such an articulatory based segmentation might lead to better inversion results than those presented in our study, which used a purely acoustic segmentation.

VII. CONCLUSION

In this paper, we have shown that the concept of episodic memory can help to enhance the articulatory-to-acoustic mapping. It has not been our aim to provide an inversion method that beats the state of the art. Instead, we have focused on the added value of this approach to address the inversion problem. Both C-Mem and G-Mem integrate the articulatory dynamics in the memory structure, which is essential to fully resolve the inversion problem. Indeed the articulatory dynamics reflect the articulatory strategy, which is chosen by a speaker to utter a sequence of phones. Regarding the physical constraints of each articulator as well as the phonetic content of the sentence, the speaker determines the best manner, among many, to articulate and coarticulate each sound. Thus, the articulatory dynamics help to resolve the non-uniqueness of the acoustic-to-articulatory inversion problem. The dynamics are actually those naturally produced by the speaker and not the result of an a posteriori computational smoothing process. Moreover, as the size of the corpus is very limited and does not include all possible articulatory transitions, the G-Mem was provided with a generalization mechanism. The G-Mem is able to combine different observed articulatory strategies according to the speech signal to be inverted. Then, even if a particular phone sequence has never been observed, the memory can predict how the II. Probabilities of improvement and 99% confidence intervals of expected improvements (both expressed in percentage). Each cell contains the probability that the method corresponding to the row outperforms the method corresponding to the column, as well as the 99% confidence interval of the expected improvements.

Figure 1

 1 Figure 1 shows an example of a G-Mem built from a set of two episodes X and Y of the same class. The top two graphs illustrate the DTW mapping of Y onto X, and of X onto Y . Lines are drawn between aligned observations of X and Y . A thick line indicates the

Figure 3

 3 Figure 3 shows the results of the three inversion methods, with (grey) and without (black) phonetic constraints. The bars show the overall RMSE (means over the coils and x-and y-coordinates) in millimeters for each corpus; the RMSE values are indicated above the bars. The respective Pearson's correlations are given below the bars. The figure shows that the memory-based approaches always outperform the codebook-based method. This suggests the articulatory dynamics are modeled more effectively in the memory-based systems than by continuity constraints used by the codebook for speech inversion. It also illustrates how

Figures 7

 7 Figures 7 shows the evolution of the overall RMSE for both the codebook based and G-Mem based methods on the test set of msak as a function of the search beam width.

  FIG. 1 Illustration of a G-Mem built from two episodes X and Y . The top two graphs illustrate the DTW mappings of Y onto X and X onto Y . The thick lines indicate that the local distance between the aligned observations of X and Y is smaller than δ. The bottom graph show the resulting G-Mem. . . . 26

TABLE I .

 I Overview of the corpora.

	Corpora Sets Durations Sentences Phones
		train 16 min 35 sec	368	11179
	fsew	dev 1 min 57 sec	46	1324
		test 2 min 5 sec	46	1457
		train 13 min 59 sec	368	11179
	msak dev 1 min 41 sec	46	1324
		test 1 min 45 sec	46	1457
		train 8 min 24 sec	319	6355
	mdem dev 1 min 2 sec	40	817
		test 1 min 3 sec	40	814
	TABLE

Acknowledgment

The authors are thankful to Yves Laprie for providing the articulograph machine facility.


Sébastien Demange
email: sebastien.demange@loria.fr

Continuous speech recognition using articulatory dynamics

Keywords: speech recognition, articulatory dynamic, episodic memory

In this paper we present a speech recognition system which uses articulatory dynamics. We do not extend the acoustic feature with any explicit articulatory measurements but instead the articulatory dynamics of speech are structurally embodied within episodic memories. The proposed recognizer is made of different memories each specialized for a particular articulator. As all the articulators don't contribute equally to the realization of a particular phoneme, the specialized memories don't perform equally regarding each phoneme. We show, through phone string recognition experiments that combining the recognition hypotheses resulting from the different articulatory specialized memories leads to significant recognition improvements.

Introduction

Over the last three decades the HMM has been established as the reference models for speech recognition. Its mathematical simplicity, the development of automatic and robust training and decoding algorithms have contributed to its success. However, while important progresses have been done in the past, the recognition accuracies seem to converge asymptotically toward a limit which is situated below the human ability.

This ceiling of performance can be partly explained by some well known HMM limitations. First, though necessary for practical reasons, the HMM probabilistic framework relies on unrealistic assumptions such as the independence of the observations and the first order Markov assumption. Second, it provides insufficient capabilities for modeling trajectories and durations which can be related to the dynamic of speech. Finally, an HMM based speech recognizer usually relies only on two sources of knowledge: the acoustic models of the phonemes and a language model providing the probability to observe a particular word sequence.

Speech perception is much more complex and uses other communication modalities which contribute to disambiguate the acoustic signal. In particular, theories in cognitive science underline the role of the articulation in speech perception [START_REF] Fadiga | Speech listening specifically modulates the excitability of tongue muscles: a tms study[END_REF][START_REF] Watkins | Seeing and hearing speech excites the motor system involved in speech production[END_REF]. We propose here a new approach relying on episodic memories which make use of articulatory dynamics knowledge. An episodic memory [START_REF] Tulving | Episodic and semantic memory[END_REF] can be viewed as a collection of past events (also called episodes) which have been experienced and stored by a subject. In response to a given stimulus, episodes (similar to the stimulus) are activated in memory and contribute to recognize (categorize) the stimulus. Such a process is expected to take place during speech perception [START_REF] Goldinger | Echo of echoes? an episodic theory of lexical access[END_REF]. In addition to its biological basement, an episodic memory is able to describe each episode with respect to different heterogeneous modalities.

For example, in this work, we define an episode as a particular realization of a phoneme which has been observed. The realization is described by both acoustic and articulatory observation sequences.

Recently, databases of synchronized acoustic and articulatory data streams using electromagnetic articulography have been possible. However, due to practical difficulties during the data acquisition, most of them is limited to few tens of minutes. Then, only few episodes of each phoneme are available which is quite insufficient for covering the important speech variability. The originality of the proposed memory lies in its ability to combine different episodes during the recognition process based on their articulatory dynamics. The memory is extended with inter-episodes transitions based on articulatory continuity constraints and which respect the original temporal order of the observations. Any path across different episodes of a same phoneme result to a particular acoustic realization which would have been produced by the resulting articulatory observation sequence.

In the next section we present how the memory is built. Section 3 explains how to recover the phone string from an acoustic speech signal using the acoustic-articulatory memory. The corpora, the feature extraction as well as the experiments set-up are presented in section 4 and the recognition results are exposed along the section 5.

Acoustic-articulatory episodic memory

The main advantage of an episodic memory is that it keeps trace of the order of the observations and thus preserves the acoustic and articulatory dynamics of each episode. In order to preserve this property, the inter-episodes transitions have to be defined carefully. Indeed, they're expected to provide the memory with generalization capabilities but must not allow the memory to produce unrealistic pattern from a dynamical point of view. The inter-episode transitions are defined accordingly to the concept of articulatory target interval (ATI).

Articulatory target interval

Let X be a particular articulatory realization of a given phoneme expressed as a sequence of K observations: X = (x1, x2, . . . , xK ). We define each observation xi+1 as the natural articulatory target of xi as it has been observed following xi. In fact, xi+1 is a particular articulatory configuration but we can suppose it could have been slightly different. Indeed, starting from xi at time t, the articulators could have reached a different target at time t + 1 close to xi+1 without affecting the rest of the realization. Then, for each xi we define an articulatory target interval AT Ix i as the interval

[xi+1 -δ, xi+1 + δ],
where δ is a positive value.

Modeling the articulatory dynamics

Inter-episode transitions

Let Y = (y1, y2, . . . , yN ) be a second articulatory realization of the same phoneme as X. We define φ = (Φ1, . . . , ΦM ) as the alignment path corresponding to the shortest distance D(X, Y ) between X and Y obtained by the well-known dynamic time warping algorithm (DTW). Each Φi is a pair of indexes of the elements of X and Y , which are aligned together: Φi = (Φx,i, Φy,i). For example, Φ3 = (4, 5) indicates that the third element of the path is an alignment between x4 the fourth articulatory configuration of X and y5 the fifth articulatory configuration of Y . For our problem, we extended the DTW algorithm with the Itakura constraints [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF] to impose temporal constraints on the alignment paths ensuring that aligned articulatory configurations occur at similar time in their respective episode. Once the DTW distance between X and Y is computed, a transition in the memory from any xi to a yj is created if yj matches the two following conditions:

Φy,i+1 = j (1) yj ∈ AT Ix i (2)
Equation 1 requires yj to be aligned with xi+1 when mapping Y onto X . In other words, it has to be aligned with the following observation of xi. This condition ensures that the transitions are consistent with the temporality of the episodes. Equation 2states that yj has to belong to the AT Ixi. It locally ensures the physical articulatory validity and naturalness of the transition since yj is close to xi+1, which is the natural articulatory target of xi. Note that the articulatory trajectories of two episodes of the same phone can be significantly different due to the coarticulation effects as their phonetic contexts can differ. Combining two episodes, which match only on a very small segment but which drastically differ outside, could result in unrealistic trajectory. To avoid this undesired effect, transitions from X to Y are created only if Y is similar enough to X:

D(X, Y ) ≤ ∆ ( 3
)
where ∆ is a positive value. One can remark that the memory is still conservative as all the original episodes it is made of are preserved since:

D(X, X) = 0 ≤ ∆ (4) Φx,i+1 = i + 1 (5) xi+1 ∈ AT Ix i (6)

Between-episode transitions

The transitions at episode boundaries are only subject to the articulatory continuity requirement expressed by the equation 2. Let Z be the episode observed following X. Then, a transition from xK (the last observation of X) to the first observation w1 of any realization W of any phoneme is created if

w1 ∈ AT Ix K = [z1 -δ, z1 + δ].
If the episode X is the last of a record the natural articulatory target of xK is unknown and equation 2 cannot be verified, thus no transition to any other episode is possible from xK .

Phone string recognition

In practice the memory is modeled as an oriented graph. The nodes are synchronized acoustic and articulatory observations and the edges are the allowed transitions reflecting natural articulatory dynamics. The transitions are created according to the procedure described above in the articulatory space using Euclidean distance. Each path within the graph corresponds to a physically possible articulatory trajectory and its corresponding acoustic counterpart. In addition, the nodes are supplemented with linguistic and temporal attributes which are the phoneme label of the acoustic-articulatory observations and their relative index within their original episode. We define the relative index RI(xi) of any observation xi being part of an episode X made of K observations as i/K. Recognizing speech consists in finding the path within the memory which acoustically best matches the speech signal X to be recognized. All paths can start only at nodes representing the first observation of the original episodes and can end only at nodes representing the last observation of the original episodes. During the recognition a breadth first search is performed applying the Viterbi algorithm. At each step, only the K best paths are propagated through all defined transitions. The score of each path is the sum of the local acoustic distances between the visited nodes and the observations of X. The local acoustic distances are computed over a window using the Euclidean distance.

The recognized phone string is deduced from the winning path. The phone segmentation is deduced from the evolution of the RI of each node along the path. Any part of the path across different episodes of the same phoneme might exhibit a continuous increase of the RIs as the inter-episodes transitions respect the temporality of the original episodes. A significant decrease indicates that the path start a new phoneme as the transitions at the episode boundaries are defined between the last (RI = 1) and the first (small RI) observations of the original episodes. All the observations within a particular segment comes from realizations of the same phoneme, that is they do not necessarily come from a unique episode, but may come from multiple episodes of a same phoneme between which inter-episode transitions have be defined.

Experiment set-up and corpora

Corpora

All the experiments presented in this work have been carried out on two corpora of synchronized acoustic speech signal and articulatory trajectories. Table 1 synthesizes the corpora.

The first corpus is MOCHA [START_REF] Wrench | A multichannel articulatory database and its applicationfor automatic speech recognition[END_REF]. Two speakers, a female (fsew) and a male (msak) British English speakers, were recorded while reading 460 short phonetically balanced British-TIMIT sentences. We use in this work the acoustic and EMA data streams. The acoustic is provided as waveforms sampled at 16 kHz and the EMA data consist in 2D data with coordinates expressed in the mid-sagital plane. Nine sensors are used, located at the bridge of the nose (bn), upper incisors (ui), lower incisors (li), upper lip (ul), lower lip (ll), tongue tip (tt), tongue body (tb), tongue dorsum (td) and velum. The two first are used to normalize the trajectories of the last seven with regard to the head movements. The y axis (upright) passes through the upper incisors (origin of the coordinate system) and bridge nose sensors.

The second one is a corpus we have recorded with an articulograph (AG500, Carstens Medizinelektronik). For now only one speaker, a male French speaker (mdem) has been recorded reading 400 phonetically balanced sentences. The data streams are synchronized acoustic waveforms sampled at 16 kHz and 3D EMA data. We used 6 sensors fixed in the mid-sagital plane on the lower lip (ll), upper lip (ul), tongue tip (tt), tongue body (tb), tongue dorsum (td) and tongue back dorsum (tbd). Three additional sensors have been used for normalizing the articulatory trajectories with regard to the head movements, one at the bridge of the nose and two located behind each ear. Prior to the features extraction the EMA data axis have been shift and rotated so that each articulatory sample is expressed in the 2D mid-sagital plane. The phonemic segmentation has been obtained by forced aligning French acoustic models onto the acoustic stream.

Features extraction

The silences occurring at the beginning and the end of the records are first discarded as the articulators may move unpredictably. A Linear Predictive Analysis [START_REF] Hermansky | Perceptual linear predictive (plp) analysis of speech[END_REF] is performed on the speech signal using the HTK toolkit [START_REF] Young | The HTK Book (for HTK version 3.3[END_REF]. 12 cepliftered MF-PLPs [START_REF] Woodland | Broadcast news transcription using htk[END_REF] plus the logarithmic energy of the signal comprise the acoustic feature vector extracted from every 25 ms speech frame shifted by 10 ms. The articulatory data are first down sampled to 100 Hz to match the acoustic frameshift. Then, all trajectories are low-pass filtered in order to remove the recording noise using a cut-off frequency of 20 Hz.

Experiments set-up

The memory parameters have been optimized on the development sets. For each corpus, a memory is built for each tracked articulator and along both the horizontal and vertical axis. That is, each memory can be considered as an articulatory special-ized speech recognizer. Each memory provides a phone string recognition hypothesis. The final recognition result is obtained combining these hypotheses using a majority vote and with the constraint that all phone segments within the resulting phone transcription is at least 30 ms long.

Results

Figure 1 shows an example of the memory outputs specialized on the tongue tip dynamics. The acoustic signal to be recognized is the French sentence "juste quelques extrémités de branches gelées" which can be translated by "only few frozen extremities of tree branches". The reference word and phone segmentations are provided at the middle of the figure. The upper and the lower parts, plot the relative index of the nodes along the winning decoding paths, the respective phone string segmentations (with misrecognized phoneme in red) and the articulatory trajectories (in blue) obtained from the articulatory observations of the nodes along the winning paths. For comparison the reference articulatory trajectories are provided with the green curves.

The articulatory specialized phone transcriptions are impressive knowing that no linguistic information (such as phone language model, or word dictionary) is provided and knowing that the memories are based only on few minutes of speech. Only, the acoustic distances computed on the static acoustic parameters and the articulatory dynamics contribute to the transcriptions. As expected, each memory performs differently. Let's consider the /k/ and /s/ which occur at the beginning of the word "extrémités". The places of articulatory constrictions are between the tongue body and the palate for the /k/ and the tongue tip and the palate just behind the upper incisors for the /s/. The realization of these two phones implies a raise of the tongue tip while the tongue tip doesn't move significantly along the front/back axis. The acoustic realization consists in a occlusion followed by a released of the air stream for the /k/ and fricative noise for the /s/. Note now, that the tty specialized memory succeeds in recognizing this phoneme sequence while the ttx fails. The ttx memory has recognized a /f/ which is fricative noise from an acoustic point of view. Regarding the tongue tip, the articulatory realization of a /f/ is close to the realization of a /s/ but the tongue tip doesn't raise in order to facilitate the air flow through the mouth as the constriction is between the upper incisors and the lower lip. Therefore, the movement of the tongue tip along the up/down axis is critical for distinguishing between the /s/ and the /f/. This is a typical example of how the phone recognition results differ across the articulatory specialized memories. Then combining the different recognition hypothesis should result in an improved phone error rate.

Table 2 summarizes the recognition results obtained from each specialized memory as well as the PER averaged over all the memories. The last two rows give PER after the majority vote combination as well as pure acoustic HMM references. For mdem we have trained left-to-right, three states monophone HMMs, each state modeled by a mixture of 8 Gaussian. For MOCHA we didn't train HMMs, instead we use the HMM baseline from [START_REF] Wrench | Continuous speech recognition using articulatory data[END_REF] which match our frontend. Their models consist in left-to-right three states triphone HMMs, each state modeled by a mixture of 2 to 7 Gaussian. In addition they used a bi-gram phone language model. Contrary to the memories, the HMM acoustic feature vectors are supplemented with the first and second time derivatives.

The results show that all the specialized memories give similar PER. The combination is very useful, resulting in a PER reduction ranging between 15% and 20% according to the corpora. Moreover, the PERs resulting from the combination are similar to the HMM baselines.

We compare our PER results with those published by Frankrel and King [START_REF] Frankel | Asr -articulatory speech recognition[END_REF] as our experiment share the same MOCHA corpus and our decoding experiements both take place in the context of phone string recognition. They obtained PERs of 33% in a classification experiment and 44% in recognition experiment while our average PER on the same dataset is 38%. The authors supplemented the acoustic feature vectors with articulatory feature estimated using a recurrent neural network. Contrary to Frankrel and King we do not make explicit use of articulatory measurements, but rather the memory structurally embody the articulatory dynamics. Such a strategy prevent from propagating acoustic-to-articulatory inversion errors to the recognizer.

Conclusions

We have proposed a new episodic memory based approach for speech recognition. The memory is able to combine different episodes it is made of based on their articulatory dynamics. We have shown that the articulatory dynamics strongly influence the recognition results as all articulatory specialized memories perform differently. Combining the articulatory specialized recognition hypotheses results in significant improvements.

Though impressive, regarding the limited size of the training material and the absence of linguistic information, the recognition results are still preliminary and could be significantly improved. We use the simple Euclidean distance for computing the acoustic distances during the decoding and this distance is known not to be robust for speech. The use of a local kernel based distance [START_REF] Wachter | Example Based Continuous Speech Recognition[END_REF] would improve the recognition. Moreover our acoustic frontend consists in only static acoustic features while the first and second derivatives would provide information of the acoustic dynamics. Finally, our majority vote combination is very simple. It seems more efficient to develop a combination procedure based on the notion of critical articulator. The weights of each articulatory specialized memory would be dynamically computed regarding the hypothesized phonemes at time. For a given hypothesized phoneme P memories which are specialized on the critical articulators for the phoneme P should contribute more to the combination.

Figure 1 :

 1 Figure 1: Tongue tip specialized memories recognition hypotheses resulting from the decoding of the French sentence "juste quelques extrémités de branches gelées".

Table 1 :

 1 Synthetic description of the corpora.

	Corpora Sets	Durations	Sentences Phones
		train 16 min 35 sec	368	11179
	fsew	dev	1 min 57 sec	46	1324
		test	2 min 5 sec	46	1457
		train 13 min 59 sec	368	11179
	msak	dev	1 min 41 sec	46	1324
		test	1 min 45 sec	46	1457
		train	8 min 24 sec	319	6355
	mdem	dev	1 min 2 sec	40	817
		test	1 min 3 sec	40	814

Table 2 :

 2 Phone recognition error rates (PER) in %.

		mdem	fsew	msak
	lix	-	41.3	49.0
	liy	-	45.7	44.9
	llx	29.6	42.5	45.9
	lly	27.3	44.1	47.9
	ulx	27.1	43.4	53.2
	uly	26.5	46.4	44.4
	ttx	27.0	41.2	41.9
	tty	26.8	45.7	44.5
	tbx	31.2	43.0	44.8
	tby	25.5	45.9	52.6
	tdx	25.0	42.0	43.1
	tdy	29.6	45.8	45.2
	tbdx	27.3	-	-
	tbdy	29.9	-	-
	vlx	-	43.3	46.5
	vly	-	47.8	53.2
	Average	27.7	44.2	46.9
	Combination	21.8	36.3	39.7
	HMM	22.8	36.0* 38.0*