Proc. International Seminar on Speech Production, pp. 1-8, 2011.

Generation of a vocal-tract MRI movie based onsparse sampling
Sadao Hiroya 1 , Tatsuya Kitamura 2
1

NTT Communication Science Laboratories, NTT Corporation
3-1 Morinosato-Wakamiya, Atsugi, Kanagawa, Japan
2

Faculty of Intelligence and Informatics, Konan University
8-9-1 Okamoto, Higashinada, Kobe, Japanhiroya.sadao@lab.ntt.co.jp, t-kitamu@konan-u.ac.jp

Abstract. We present a novel technique that can provide a high-quality vocaltract MRI movie during speech production. The method uses MRI vocal-tract
images at the central point of each phoneme and interpolation functions between
adjacent phonemes obtained from electromagnetic articulographic (EMA) data.
It is based on our finding that articulatory parameters are suitable for a sparse
representation. Preliminary results showed that the quality of the obtained
vocal-tract MRI movie is high compared to that for the previous technique.
The method will be useful for constructing a large database of vocal-tract MRI
movies and understanding speech production mechanisms.

1. Introduction
The measurement of articulatory motion during speech production is important for investigating human speech production mechanisms. In the past, X-ray cinematography(Perkell, 1969), an X-ray microbeam system (Kiritani et al., 1975), and an ultrasoundtechnique (Stone et al., 1983) have been used to measure motions of human articulatoryorgans. In the 1990’s, an electromagnetic articulographic (EMA) system was widely used
for speech production research (Perkell et al., 1992). This system can measure the tongue,
lips, incisors, and velum at fairly high rates, but can only track receiver coils in the midsagittal plane attached to the speech articulators. The advantage of this technique is thatthe calculation of velocity and acceleration is more direct with flesh points. Currently, ameasuring system that works in three dimensions is being developed and used (Kaburagi
et al., 2005).
Recently, it has become possible to measure complete two- or three-dimensional
vocal-tract motion during speech production by magnetic resonance imaging (MRI). Oneof the imaging methods is synchronized sampling with an external trigger (Masaki et al.,1999). This can generate a high-quality vocal-tract MRI movie, but the subject has torepeat the same sentence more than 100 times with small trial-to-trial variability. Anotheris a real-time MRI technique (Narayanan et al., 2004), but lower sampling rates are used
as for measuring articulatory motion during speech production. Therefore, a techniquethat can provide a high-quality vocal-tract MRI movie and greatly reduce the number ofsentence repetitions in the MRI scanner is required.

Proc. International Seminar on Speech Production, pp. 1-8, 2011.

We have proposed a method for decomposing EMA-based mid-sagittal articulatory parameters into a set of temporally overlapped event functions and corresponding event vectors using non-negative temporal decomposition (NTD) (Hiroya, 2010) andfound that articulatory motion during speech production can be represented by articulatory positions at the central point of each phoneme and adequate interpolation functions.
In this paper, we present a method for generating a vocal-tract MRI movie basedon sparse sampling: vocal-tract MRI image sequences during speech production can begenerated using MRI images at key frames and interpolation functions obtained by EMA
measurements. To verify this validity, we investigate the articulatory properties obtainedwith EMA and MRI, and compare the proposed method with the previous technique.

2. Non-negative temporal decomposition
The original temporal decomposition (Atal, 1983) approximates the ith parameter yi (t)of time t toŷi (t) =m∑ai,k φk (t),1 ≤ t ≤ T, 1 ≤ i ≤ p,

(1)

k=1

where ai,k is the kth event vector, φk (t) is the kth event function, p is the dimension of
the parameter, T is the length of the parameter sequence, and m is the number of event
functions. Later, Shiraki (2004) and Kim and Oh (1999) assumed that φk (t) is zero fort < tk−1 and t > tk+1 and that the sum of all event functions is one at any time t:ŷi (t) = ai,k φk (t) + ai,k−1 φk−1 (t),tk−1 ≤ t ≤ tk ,where φk (t) + φk−1 (t) = 1. (2)They claimed that the event functions φk (t) should be restricted to the range [0, 1], because
this restriction facilitates parameter modification. The event functions are determined byminimizing the least-squares error between yi (t) and ŷi (t), but the obtained event functions are not restricted to the range [0, 1]. Thus, Kim clipped event functions at the range[0, 1] [i.e. φk (t) = min(1, max(0, φk (t)))], but this resulted in increasing the estimation
errors.
It is difficult to determine the event functions that are restricted to the range [0, 1]
by the least-squares method without a clipping. To overcome this difficulty, we haveproposed a method for decomposition of speech parameters into a set of temporally overlapped event functions φk (t) that are restricted to the range [0, 1] and corresponding event
vectors ai,k using NTD (Hiroya, 2010). The point is that NTD can determine the interpolation functions and the key frames based on a combination of non-negative matrix
factorization (NMF) (Lee and Seung, 1999) and dynamic programming (DP). We consider minimizing the following cost function by the NMF algorithm:
p
tk
m ∑
∑
∑
k=2 t=tk−1 i=1(yi (t) − ai,k φk (t) − ai,k−1 φk−1 (t))2 + α

tk
m ∑
∑

(φk (t) + φk−1 (t) − 1)2 , (3)k=2 t=tk−1

where 1 = t1 < t2 < · · · < tm = T and α is the weight. In line with a previous idea

Proc. International Seminar on Speech Production, pp. 1-8, 2011.

(Virtanen, 2007), we can obtain the multiplicative update rule for the event function:p∑

φk (t) ←ai,k yi (t) + αi=1φk (t)(4)

φk−1 (t)p∑2(ai,k−1 ai,k φk (t) + ai,k−1 φk−1 (t)) + α(φk−1 (t) + φk (t))

(5)

p
∑

(ai,k−1 ai,k φk−1 (t) + a2i,k φk (t)) + α(φk−1 (t) + φk (t))

i=1
p
∑

φk−1 (t) ←ai,k−1 yi (t) + α

i=1

i=1

The distortion d(y(t), ŷ(t)) of the cost function for each interval tk−1 ≤ t ≤ tk
only depends on time tk−1 and tk . Therefore, the event timing
tk = argmint2 ,··· ,tm−1T∑d(y(t), ŷ(t)),

2≤k ≤m−1(6)

t=1that minimizes total distortion for the whole interval 1 ≤ t ≤ T is derived efficiently byutilizing the DP method (Shiraki, 2004). That is, we have(D(tk ) =min

tk−1 ∈Rk−1

D(tk−1 ) +tk∑)
d(y(t), ŷ(t))(7)

t=tk−1

where
Rk−1 = {t|tk−1 − δ ≤ t ≤ tk−1 + δ},(8)

D(tk ) is an accumulated minimal distortion at tk and δ is a search range.
In previous work (Hiroya, 2010), we assumed that event timing corresponds to thecentral point of each phoneme on the time axis and found that the mean errors between themeasured EMA-based mid-sagittal articulatory parameters y(t) and the estimated onesŷ(t) were 0.16 mm. This indicates that articulatory motion during speech production
can be represented by articulatory positions at the central point of each phoneme and byadequate interpolation functions: articulatory parameters are suitable for a sparse representation. This is because the articulatory parameters for a given phoneme are prominentand because the temporal patterns of the articulatory parameters are simple and smooth.
Thus, we had the idea that we would measure vocal-tract MRI images at key frames using
sparse sampling. In the next section, we will show the procedure for the proposed method.

3. Procedure for proposed method
We measured the vertical and horizontal positions of articulators, such as the lips, four
tongue positions, and lower incisor, using a two-dimensional EMA system (CarstensProc. International Seminar on Speech Production, pp. 1-8, 2011.

EMA measurements
Articulatory data
Initial event timings

NTD

Speech sounds

Event timings
Interpolation functions(event functions)

MRI measurements
MRI images at event timings

Generation of a vocal-tract MRI movie

Figure 1. Procedure for proposed method.
AG100). The system can measure the articulators at a sampling rate of 250 Hz and simultaneously record the speech sounds at a sampling rate of 48 kHz.
In NTD, we first set the event vectors as ai,k = yi (tk ) for the initial event timing tk ,which was labeled at the central point of each phoneme on the time axis manually. Then,the event functions and the event timings are determined by DP and the NMF update rulesof Eqs. (4)-(5). Note that event vectors were determined as ai,k = yi (tk ) depending on
the event timings.
Next, we conducted MRI scanning [Siemens MAGNETOM Verio (3T)] for vocaltract motion during speech production. For the MRI scanning, it was necessary that thesame subject in the EMA experiments produce the same sentences as in the EMA experiment through MRI-compatible headphones in order to match the articulatory timing. The
MRI scan was conducted at event timing tk and measured the mid-sagittal section of the
vocal tract.
Using the event functions φk (t) obtained by EMA and vocal-tract MRI images (asevent vectors) bi,j,k (1 ≤ i, j ≤ q) at event timing tk , where q is the pixel length, we
generated a vocal-tract MRI movie x̂i,j (t) = bi,j,k φk (t) + bi,j,k−1 φk−1 (t) (tk−1 ≤ t ≤ tk )at the sampling rate of 250 Hz. The restriction of φk (t) = 1 − φk−1 (t) ∈ [0, 1] would
be beneficial for switching the event vectors of EMA with those of MRI. Shape-based
interpolation (Grevera and Udupa, 1996) between MRI images was introduced in order toobtain sharp contours of the speech articulators.

4. Results
4.1. Comparison between event functions of MRI and EMA
In the proposed method, we need to use the event functions of EMA to generate a vocaltract MRI movie. Thus, we should firstly investigate whether the event functions of EMA
are similar to those of MRI. We created a vocal-tract MRI movie of the Japanese vowel
sequence /aiueo/ based on synchronized sampling with an external trigger and measured
EMA-based articulatory parameters of the same sequence. A Japanese male subject produced the vowel sequence to rhythmic repetitions of a noise-burst train composed of athree-beat rhythm for both MRI and EMA experiments. Sampling rates were 67 and 250
Hz for MRI and EMA, respectively. A 256 × 256 mm2 field of view, 256 × 256-pixelProc. International Seminar on Speech Production, pp. 1-8, 2011.

/a/

/i/

/u/

/e/

/o/

EM A

1
Lips

0.
5

TongueM RI(l
ips)

0

150

200

250

300

350

1
0.
5
0

M RI(whol
e)

Jaw

50

60

70

80

90

100

110

50

60

70

80

90

100

110

1
0.5
0

Figure 2. Event functions of EMA, MRI around the lips of the vocal tractand MRI at the whole articulatory area of /aiueo/.

/a/

/i/

/u/

/e/

/o/

Figure 3. Comparison between upright (EMA) and supine (MRI) postures for /a/ to /o/. White circles represent receiver coil positions of EMA.
image size, and 3-mm slice thickness were used for MRI. All MRI images were converted
to 8-bit gray scale. In NTD, the α was set to 106 and 108 for EMA and MRI, respectively,and δ was around 40 msec. A matrix of the MRI images was converted to a vector for
NTD.
Figure 2 shows the event functions of EMA and MRI. Event functions of MRI
were calculated for two regions: around the lips of the vocal tract, as compared to EMA,
and at the whole articulatory area. Event functions of EMA were similar to those of MRI.
Moreover, there is not a large difference between two event functions of MRI. These
results indicate that we can use event function of EMA to generate a vocal-tract MRI
movie.
The subject spoke in the upright postures for EMA measurements, whereas MRI
measurements require supine postures. However, a previous study (Kitamura et al., 2005)
suggested individual differences in the effects of gravity on the tongue body between
upright and supine postures. As shown in Fig. 3, there is not a large difference in theProc. International Seminar on Speech Production, pp. 1-8, 2011.

/a/

/i/

/u/ /e//o/Speech01.875 (msec)

Dynamic

SparseFigure 4. Scan timings for dynamic and sparse MRI of /aiueo/. Verticalbars indicate scan timings.

/a/

/u/

/b/

Dynamic

SparseFigure 5. Vocal-tract images of dynamic and sparse methods. /a/ and/u/ of /aiueo/, and /b/ of /tabako/.
vocal-tract information between upright (EMA) and supine (MRI) postures of /a/ to /o/for the subject.
4.2. Comparison between dynamic and sparse MRI
We compared the proposed method, called sparse MRI, with a high-speed scanning
method, called dynamic MRI. Both TRs were 100 msec. A 144 × 144-pixel size and10-mm slice thickness were used for MRI. All MRI images were converted to 8-bit grayscale. Dynamic MRI scanned at the rate of 9 Hz and sparse MRI scanned at event timing tk in the mid-sagittal section of vocal tract (Fig. 4). For both scannings, the subject
tried to speak while he listened to his own voice from the EMA experiment through MRIcompatible headphones in order to match the articulation timing. The Japanese vowel sequence /aiueo/ and the word /tabako/ were used. For /tabako/ of sparse MRI, the numberof repetitions required was two because TR was larger than the durations between successive event timings. The image of dynamic MRI at the central point of each phoneme

Proc. International Seminar on Speech Production, pp. 1-8, 2011.

Synchronized

Dynamic

SparseFigure 6. Vocal-tract images between /e/ and /o/ for the synchronizedsampling, dynamic, and sparse methods.
was obtained by linear interpolation between two successive images. Figure 5 shows that
vocal-tract image at the central point of /a/ is the same for both sparse and dynamic MRI
since scan timing of dynamic and sparse MRI was almost the same. However, sparse MRI
is better than dynamic MRI for the lip protrusion for /u/ and lip closure for /b/. This indicates that sparse MRI can measure the most important vocal-tract information for eachphoneme, compared with dynamic MRI.
4.3. Generation of a vocal-tract MRI movie
Figure 6 shows vocal-tract images for the synchronized sampling, dynamic, and sparse
methods in the same timing between /e/ and /o/ of /aiueo/. The image of sparse MRI was
obtained with the proposed method and the image of dynamic MRI was obtained by linear
interpolation between two successive images. We can see that the image of sparse MRI
is more similar to that of synchronized sampling than that of dynamic MRI because theadequate interpolation functions were used for sparse MRI. The quality of the vocal-tract
MRI movie obtained by sparse MRI was high compared to that for dynamic MRI.

5. Conclusions
We presented a method for generating a vocal-tract MRI movie based on sparse sampling.
This is the first study to generate a vocal-tract MRI movie during speech production using
the high temporal resolution of EMA and high spatial resolution of MRI. The method will
be useful for constructing a large database of vocal-tract MRI movies and understanding
speech production mechanisms.

Acknowledgments
The authors thank Drs. H. Gomi of NTT CS Labs and S. Takano of ATR BAIC for many
useful and helpful discussions, and Dr. T. Mochida of NTT CS Labs for helping with the
EMA experiments. This study was partly supported by JSPS KAKENHI (21300071).

References
Atal, B. Efficient coding of LPC parameters by temporal decomposition. In ICASSP,
pages 81–84, 1983.

Proc. International Seminar on Speech Production, pp. 1-8, 2011.

Grevera, G. and Udupa, J. Shape-based interpolation of multidimensional grey-level images. IEEE Med. Imaging, 15(6):881–892, 1996.
Hiroya, S. Non-negative temporal decomposition of speech parameters. In Proc. ICASSP,
pages 5066–5069, 2010.
Kaburagi, T., Wakamiya, K., and Honda, M. Three-dimensional electromagnetic articulography: A measurement principle. J. Acoust. Soc. Am., 118(1):428–443, 2005.
Kim, S.-J. and Oh, Y.-H. Efficient quantisation method for LSF parameters based on
restricted temporal decomposition. Electronics Letters, 35(12):962–964, 1999.
Kiritani, S., Itoh, K., and Fujimura, O. Tongue-pellet tracking by a computer-controlled
X-ray microbeam system. J. Acoust. Soc. Am., 57(6):1516–1520, 1975.
Kitamura, T., Takemoto, H., Honda, K., Shimada, Y., Fujimoto, I., Syakudo, Y., Masaki,S., Kuroda, K., Oku-uchi, N., and Senda, M. Difference in vocal tract shape between
upright and supine postures: Observations by an open-type MRI scanner. Acoust. Sci.
& Tech., 26(5):465–468, 2005.
Lee, D. and Seung, H. Learning the parts of objects by non-negative matrix factorization.
Nature, 401:788–791, 1999.
Masaki, S., Tiede, M., Honda, K., Shimada, Y., Fujimoto, I., Nakamura, Y., and Ninomiya, N. MRI-based speech production study using a synchronized sampling
method. J. Acoust. Soc. Jpn. (E), 20(5):375–379, 1999.
Narayanan, S., Nayak, K., Lee, S., Sethy, A., and Byrd, D. An approach to real-timemagnetic resonance imaging for speech production. J. Acoust. Soc. Am., 115(4):1771–
1776, 2004.
Perkell, J., Cohen, M., Svirsky, M., Mathies, M., Garabieta, I., and Jackson, M. Electromagnetic midsagittal articulometer system for transducing speech articulatory movements. J. Acoust. Soc. Am., 92(6):3078–3096, 1992.
Perkell, J. Physiology of speech production: Results and implications of a quantitative
cineradiographic study. Number 53 in Research Monograph. MIT Press, Cambridge,
1969.
Shiraki, Y. Optimal temporal decomposition for voice morphing preserving ∆ cepstrum.
IEICE Trans. Fundamentals, E87-A(3):577–583, 2004.
Stone, M., Sonies, B., Shawker, T., Weiss, G., and Nadel, L. Analysis of real-time ultrasound images of tongue configuration using a grid-digitizing system. J. Phonetics, 11(3):207–218, 1983.
Virtanen, T. Monaural sound source separation by nonnegative matrix factorization withtemporal continuity and sparseness criteria. IEEE Trans. Audio Speech and Lang. Process., 15(3):1066–1074, 2007.