A new cosine series antialiasing function and its application to aliasing-free glottal source models for speech and singing synthesis

TL;DR

This paper introduces a novel cosine-series-based antialiasing filter for generating aliasing-free excitation signals, enhancing speech and singing synthesis models with improved accuracy and flexibility.

Contribution

It presents a new antialiasing filter design procedure and applies it to improve existing speech source models, with open-source MATLAB implementations.

Findings

Successfully generated aliasing-free excitation signals

Enhanced the reliability of speech synthesis models

Provided open-source tools for speech science

Abstract

We formulated and implemented a procedure to generate aliasing-free excitation source signals. It uses a new antialiasing filter in the continuous time domain followed by an IIR digital filter for response equalization. We introduced a cosine-series-based general design procedure for the new antialiasing function. We applied this new procedure to implement the antialiased Fujisaki-Ljungqvist model. We also applied it to revise our previous implementation of the antialiased Fant-Liljencrants model. A combination of these signals and a lattice implementation of the time varying vocal tract model provides a reliable and flexible basis to test fo extractors and source aperiodicity analysis methods. MATLAB implementations of these antialiased excitation source models are available as part of our open source tools for speech science.

Equations

$$E(t)$$

$$\alpha$$

$$p(\tau )$$

$$={\mathbf{p}}^T{\mathbf{u}}_{\tau }-{\mathbf{p}}^T\mathbf{B}{\mathbf{u}}_{\tau -1},$$

$$\mathbf{B}=$$

$$\mathbf{p}=$$

$$p_h(\tau )$$

$$\displaystyle\!\!\!\!\!\mathbf{h}_{\tau}\!=\!\!\left[\!\!\!\!\begin{array}[]{l}h(\tau)\!\ast\!\tau^{0}u(\tau)\\ h(\tau)\!\ast\!\tau^{1}u(\tau)\\ \vdots\\ h(\tau)\ast\tau^{n}u(\tau)\end{array}\!\!\!\!\right]\!\!,\mathbf{h}_{\tau-1}\!=\!\!\left[\!\!\!\begin{array}[]{l}h(\tau)\!\ast\!(\tau\!-\!1)^{0}u(\tau\!-\!1)\\ h(\tau)\!\ast\!(\tau\!-\!1)^{1}u(\tau\!-\!1)\\ \vdots\\ h(\tau)\!\ast\!(\tau\!-\!1)^{n}u(\tau\!-\!1)\end{array}\!\!\!\right]\!\!.$$

$$h(t)$$

$$\displaystyle\mathbf{h}_{\tau}=\left\{\begin{array}[]{ ll}0,&t\leq-t_{w}\\ \mathbf{C}\mathbf{c}_{t}+\mathbf{S}\mathbf{s}_{t}+\mathbf{U}\mathbf{t}_{n+1},&-t_{w}<t\leq t_{w}\\ \mathbf{V}\mathbf{t}_{n},&t_{w}<t\end{array}\right.$$

$${\mathbf{h}}_{\tau -1}=\mathbf{C}{\mathbf{c}}_{t-1}+\mathbf{S}{\mathbf{s}}_{t-1}+\mathbf{U}{\mathbf{d}}_{n+1},1-t_w<t\le 1+t_w,$$

$$({\mathbf{c}}_t{)}_k=\cos \left(\frac{k\pi t}{t_w}\right),({\mathbf{s}}_t{)}_k=\sin \left(\frac{k\pi t}{t_w}\right),({\mathbf{t}}_n{)}_k=t^{k-1},$$

$$({\mathbf{c}}_{t-1}{)}_k$$

$$({\mathbf{d}}_{n+1}{)}_k$$

$$\mathbf{C}{\mathbf{c}}_t+\mathbf{S}{\mathbf{s}}_t+\mathbf{U}{\mathbf{t}}_{n+1}=0$$

$$\mathbf{C}{\mathbf{c}}_t+\mathbf{S}{\mathbf{s}}_t+\mathbf{U}{\mathbf{t}}_{n+1}=\mathbf{V}{\mathbf{t}}_n$$

$$\displaystyle\begin{array}[]{l}\mathbf{C}_{0,k}=0,\ \ 1\leq k\leq m,\\[4.0pt] \displaystyle\mathbf{S}_{0,k}=\left(\frac{t_{w}}{k\pi}\right)h_{k},\ \ 1\leq k\leq m,\\[4.0pt] \mathbf{U}_{0,1}=h_{0},\\[4.0pt] \mathbf{U}_{0,k}=0,\ \ 1<k\leq n+1,\\[4.0pt] \mathbf{V}_{0,0}=1,\ \ \mbox{(note that this is missing in \cite[cite]{[\@@bibref{}{milenkovic1993voice}{}{}]})},\\[4.0pt] \mathbf{V}_{0,k}=0,\ \ 1\leq k\leq n,\end{array}$$

$$\displaystyle\!\!\!\!\!\!\!\begin{array}[]{l}\displaystyle\mathbf{C}_{r,k}=-\left(\frac{rt_{w}}{k\pi}\right)\mathbf{S}_{r-1,k},\ \ 1\leq k\leq m,\\[8.0pt] \displaystyle\mathbf{S}_{r,k}=\left(\frac{rt_{w}}{k\pi}\right)\mathbf{C}_{r-1,k},\ \ 1\leq k\leq m,\\[8.0pt] \displaystyle\mathbf{U}_{r,k}=\left(\frac{r}{k}\right)\mathbf{U}_{r-1,k-1},\ \ 1\leq k\leq n+1,\\[5.0pt] \displaystyle\mathbf{U}_{r,0}=-\sum_{k=1}^{n+1}(-t_{w})^{k}\mathbf{U}_{r,k}-\sum_{k=1}^{m}(-1)^{k}\mathbf{C}_{r,k},\\[8.0pt] \displaystyle\mathbf{V}_{r,k}=\left(\frac{r}{k}\right)\mathbf{V}_{r-1,k-1},\ \ 1\leq k\leq n,\\[6.0pt] \displaystyle\mathbf{V}_{r,0}=\!\!\sum_{k=0}^{n+1}(t_{w})^{k}\mathbf{U}_{r,k}\!+\!\!\sum_{k=1}^{m}(-1)^{k}\mathbf{C}_{r,k}\!-\!\!\sum_{k=1}^{n}(t_{w})^{k}\mathbf{V}_{r,k},\end{array}$$

$$p_h(t)$$

$$\displaystyle\begin{array}[]{l}\mathcal{Q}_{1}(t)=\{-t_{w}<t\leq t_{w}\wedge t\leq 1-t_{w}\}\\ \mathcal{Q}_{2}(t)=\{t_{w}<t\leq 1-t_{w}\}\\ \mathcal{Q}_{3}(t)=\{1-t_{w}<t\leq t_{w}\}\\ \mathcal{Q}_{4}(t)=\{1-t_{w}<t\leq 1+t_{w}\wedge t_{w}<t\}\end{array},$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\begin{array}[]{llll}\mathbf{c}_{0}\!=\mathbf{p}^{T}\mathbf{C},&\!\!\mathbf{s}_{0}\!=\mathbf{p}^{T}\mathbf{S},&\!\!\mathbf{u}_{0}\!=\mathbf{p}^{T}\mathbf{U},&\!\!\!\mathbf{v}\!=\mathbf{p}^{T}\mathbf{V},\\ \mathbf{c}_{1}\!=\mathbf{p}^{T}\mathbf{B}\mathbf{C},&\!\!\mathbf{s}_{1}\!=\mathbf{p}^{T}\mathbf{B}\mathbf{S},&\!\!\mathbf{u}_{1}\!=\mathbf{p}^{T}\mathbf{B}\mathbf{U}&\end{array}$$

$${\psi }_k(t)$$

$$h(0)$$

$$h(\pm t_w)$$

$${\frac{d^{2p}h(t)}{d{t}^{2p}}|}_{t=\pm t_w}$$

$$\mathbf{q}=\mathrm{R}\mathbf{g},$$

$$\mathrm{R}$$

$$\mathbf{q}$$

$$\mathbf{g}$$

$$\{h_k{\}}_{k=0}^4$$

$$0.0460129686,0.0036907422\}$$

$$\{h_k{\}}_{k=0}^5$$

$$0.0726831633,0.0125124215,0.0007833203\},$$

$$p_h(t)$$

Full text

Index Terms: antialiasing, glottal source, piece-wise polynomial, piece-wise exponential, cosine series

1 Introduction

Voice quality plays important roles in speech communication, especially in para- and non-linguistic aspects. To test such aspects of speech communication, it is important to use relevant test stimuli that sound natural to listeners and that at the same time have to be precisely determined and controlled. In addition, it is desirable for the stimuli to be easy to interpret in terms of voice production, as well as auditory perception.

Source filter models with source filter interaction[1] and glottal excitation models[2, 3, 4, 5, 6, 7, 8, 9] may provide practical and useful tools. However, because glottal excitation comprises several types of discontinuity, aliasing introduces spurious signals that interfere with reliable subjective tests. This paper introduces a systematic procedure to eliminate the aliasing problem by deriving a closed form representation of the antialiased excitation. We also introduce a new set of antialiasing functions using a cosine series to keep the level of spurious signals around the Nyquist frequency low. We will make them accessible by providing MATLAB implementations as well as interactive GUI-based tools[10].

2 Background and related work

One of the authors has been developing a speech analysis, modification and resynthesis framework and related tools[11, 12, 13]. They are based on $f_{\rm o}$††We use “$f_{\rm o}$” instead of “F0” based on Reference[14]. -adaptive procedures, which require reliable and precise extraction of $f_{\rm o}$ and aperiodicity information. Development of such $f_{\rm o}$ extractors requires dependable ground truth. The aliasing-free L–F model (Fant–Liljencrants model) provided the ground truth for developing our new source information analysis framework[15].

A search for aliasing-free glottal source models produced one reference[16], which is not directly applicable to antialiasing of the L–F model for two reasons. First, it provides a procedure to antialias a piece-wise polynomial function, whereas the L–F model is a piece-wise exponential function. We had to derive the closed-form representation of the antialiased L–F model by ourselves[17]. Second, the reference by Milenkovic[16] had some typos, some equations were missing and a sample implementation was not available. In this article, we repaired the procedure and developed executable MATLAB functions. Retrospectively we found BLIT (band-limited impulse train)-based methods[18, 19].We found that literatures on digital representation of analog musical signals also provide aliasing reduction methods[20, 21, 22]. Our formulation with a new cosine series is more flexible and provides better aliasing suppression.

There are two important reasons for deriving an aliasing-free Fujisaki–Ljungqvist model. First reason is the relevance of the model. We found that the L–F model does not necessarily model the actual glottal source behavior, especially extreme voices[23]. Such voices sometimes consist of stronger discontinuities than the L–F model provides. The Fujisaki–Ljungqvist model provides several different levels of discontinuity and is reported to fit actual speech samples better[6].

The second and the most important reason is that it enables to develop a general procedure for antialiasing glottal source models. The Fujisaki–Ljungqvist model is a piece-wise polynomial, whereas the L–F model is a piece-wise exponential. The other popular glottal source models can be represented using both or one of these representations[2, 3, 4, 5, 6, 7, 8, 9]. Developing a procedure to make the Fujisaki–Ljungqvist model aliasing-free provides the necessary means of attaining this goal.

3 Fujisaki–Ljungqvist model

The excitation signal $E(t)$ of the Fujisaki–Ljungqvist model is defined by the following equation in the continuous time domain[5, 6]:

[TABLE]

where

[TABLE]

Figure 1 shows the Fujisaki–Ljungqvist model parameters and waveforms. The excitation signal $E(t)$ is a time derivative of the glottal flow $U_{g}(t)$. When the fundamental period $T$ is normalized to one, six model parameters determine the shape. They are $A,B,C,R,F,$ and $D$.

The Fujisaki–Ljungqvist model is a piece-wise polynomial function. Because filtering for antialiasing is a linear operation, antialiasing of the Fujisaki–Ljungqvist model is solved by adding each antialiased polynomial pulse $p(\tau)$, which is defined in $(0\leq\tau\leq 1)$ using relevant scaling:

[TABLE]

where the scaling is given by $\tau=t/T$.

4 Antialiasing a polynomial pulse

We use the framework proposed by reference[16]††Note that we removed the model weighting coefficient $w_{m}$ in the original reference[16]. It significantly simplifies the following discussion. . It starts from a matrix representation by introducing step functions and antialiasing in the continuous time domain, followed by equalization in the discrete time domain. The simplification and this equalization, as well as their application to the Fujisaki–Ljungqvist model, are our contribution.

This section fixes typos, adds missing equations and clarifies and adds missing conditions in the original reference[16]. All equations are defined in the continuous time domain.

4.1 Matrix form

Let us start with the introduction of the two step functions $u(\tau)$ and $u(\tau-1)$, which yields

[TABLE]

where $\mathbf{B}$ is a lower triangular matrix defined as follows:

[TABLE]

where the coefficient $b_{nk}=n!/(k!(n-k)!)$ is a binomial coefficient, and $\mathbf{p}$, $\mathbf{u}_{\tau}$, and $\mathbf{u}_{\tau-1}$ are vectors defined as follows:

[TABLE]

4.2 Antialiasing using a cosine series

The antialiased filtered pulse $p_{h}(\tau)$ has the following form:

[TABLE]

where $h\ast u$ represents the convolution of $h$ and $u$. Convolution is a linear operation and yields the following:

[TABLE]

where $\mathbf{h}_{\tau}$ and $\mathbf{h}_{\tau-1}$ represent the following vectors:

[TABLE]

Equation 30 provides the antialiased signal to be discretized. We use the following cosine series for antialiasing the filter response $h(t)$:

[TABLE]

where $m$ represents the exponent of the highest order term and $t_{w}$ represents the length of the filter, which is determined according to the sampling frequency.

Assigning specific values to the coefficients $h_{k}$ determines $\mathbf{h}_{\tau}$ and $\mathbf{h}_{\tau-1}$. We derive the closed-form representation of Equation 29 and then discretize it. Equations 44 and 45 determine $\mathbf{h}_{\tau}$ and $\mathbf{h}_{\tau-1}$, respectively. The next step is to determine the matrices $\mathbf{C},\mathbf{S},\mathbf{U}$, and $\mathbf{V}$, with constant coefficients.

[TABLE]

where each element of the vectors $\mathbf{c}_{t}$, $\mathbf{s}_{t}$, $\mathbf{t}_{n+1}$, and $\mathbf{t}_{n}$ is defined as follows, where $(\mathbf{c}_{t})_{k}$ represents the $k$-th element of $\mathbf{c}_{t}$:

[TABLE]

and vectors $\mathbf{c}_{t-1}$, $\mathbf{s}_{t-1}$, and $\mathbf{d}_{n+1}$ are defined thus:

[TABLE]

4.3 Recursive determination of coefficients

The recursion uses the following continuity constraints and initial values. The continuity gives the following:

[TABLE]

The initial condition is for the first row of the matrices, where $m$ represents the highest exponent of the selected cosine series antialiasing function. In the following equation, $\mathbf{C}_{r,k}$ represents the element of the $r$-th row and $k$-th column:

[TABLE]

For $r=1$ through $r=n$ the following recursion determines the coefficients.

[TABLE]

where the fourth and sixth lines were placed in a confusing manner in Reference[16]. Further, the reference has a typo in Eq. 56, which defines $\mathbf{U}_{0,1}$.

4.4 Closed form representation

The following equation provides the antialiased polynomial value at a given $t$:

[TABLE]

where $\wedge$ represents logical AND. Note that we have to refine the conditions given in Reference[16] to make this procedure to work properly. The coefficient vectors are defined as follows:

[TABLE]

5 Antialiasing function

We introduce two new cosine series antialiasing functions here. Because of the infinite frequency range of the glottal excitation models, which have discontinuities, commonly used time windowing functions with 6 dB/oct sidelobe decay[24, 25, 26, 27] introduce significant spurious due to aliasing. In our previous derivation, we used one of Nuttall’s windows[28] for the antialiasing function. The maximum side lobe level of the window is $-82.60$ dB and the decay speed of side lobes is 30 dB/oct. This decay rate is not steep enough to suppress spurious components around $f_{\rm o}$, and the sidelobe level is not low enough to suppress spurious components around the Nyquist frequency. The other windows listed in [28] cannot solve both issues at the same time.

5.1 Design procedure

Using s similar process as that used in [28], we design a new set of windows to satisfy both the sidelobe level and the decay conditions. We use a cosine series to design the antialiasing function. The elements of the cosine series have the following form:

[TABLE]

For the designed function to behave properly for antialiasing, the following conditions have to be satisfied. Note that the odd ordered derivatives are always zero for $t=\pm t_{w}$.

Sum of coefficients should be one

This determines the height of the function at the origin. For simplicity, we set it to one.

Level at the end point should be zero

The function should be continuous at the end point.

Derivatives at the end point should be zero

Depending on the required slope, the derivatives at the end point have to be equal to zero. For the decay rate $6+(12\times P)$ dB/oct, derivatives up to the order $2P$ are zero.

These conditions are summarized by the following equations.

[TABLE]

Once the desired decay is decided, it provides $P+2$ conditions. If the number of coefficients of the cosine series is equal to $P+2$, there is no room for adjustment. By adding one adjustable coefficient $q_{0}$, we can control the sidelobe level. It yields the following equation:

[TABLE]

where

[TABLE]

The solution $\mathbf{g}$ is given by

[TABLE]

where the elements of $\mathbf{g}$ provide the coefficients of the cosine series. Designing an antialiasing function means tuning the parameter $q_{0}$ to minimize the target cost. This time the target is the level of the maximum sidelobe level. For a 42 dB/oct decay, a five-term cosine series is designed. For a 54 dB/oct decay, a six-term cosine series is designed.

Numerical optimization yielded the following coefficients:††The coefficients are rounded to ten digits to the right of the decimal point. The best coefficients were selected from rounded numbers.

[TABLE]

Substituting these coefficients in Eq. 77 provides an antialiased Fujisaki–Ljungqvist model. Half of the window length $t_{w}$ is set to make the first zero of the frequency domain representation coincide with half of the sampling frequency $f_{s}$.

Figure 2 shows the gain of these new functions and the Nuttall window that was used for the aliasing-free L–F model[17]. The five-term function has a $-99.23$ dB maximum sidelobe level with a 42 dB/oct decay. The six-term function has a $-114.24$ dB maximum sidelobe level with a 54 dB/oct decay. We decided to use the six-term function afterward.

5.2 Discrete time domain: equalization

The final stage, which is processed in the discrete time domain, is equalization.

The designed antialiasing functions introduce severe attenuation around the Nyquist frequency. In our antialiased L–F model, an FIR equalizer was designed to compensate for this attenuation[17]. We used a simple IIR filter with six poles in this Fujisaki–Ljungqvist model and found that it equalizes this attenuation effectively. A sample implementation produces equalized gain deviations from the FIR version within $\pm 0.2$ dB from 0 to 16 kHz for 44100 Hz sampling.

6 Application to glottal source models

The antialiased Fujisaki–Ljungqvist model output is obtained by the calculation of each antialiased polynomial and sum together. Discretization is the calculation of the antialiased Fujisaki–Ljungqvist model value at each sampling instance. Finally, applying the discrete time IIR equalizer to the discretized samples provides the discrete signal of the antialiased Fujisaki–Ljungqvist model.

6.1 Examples: Fujisaki–Ljungqvist model performance

We implemented this procedure using MATLAB and prepared high-level APIs. One function generates an excitation source signal using the given $f_{\rm o}$ trajectory and the time-varying Fujisaki–Ljungqvist model parameter set $A(t),B(t),C(t),R(t),F(t)$ and $D(t)$.

Figure 3 shows the power spectra of the Fujisaki–Ljungqvist model outputs. Direct discretization generates aliasing noise around $-60$ dB from the peak harmonics level. The noise level of the final equalization around the fundamental component is about $-180$ dB when using the six-term proposed function. Antialiasing using the Nuttall-11 window introduces approximately 20 dB higher spurious levels. (Note that the noise level using nuttallwin of MATLAB as antialiasing is approximately $-120$ dB. )

6.2 Revision of antialiased L–F model

We revised our previous antialiased L–F model[17] using the new six-term cosine series and the six-pole IIR equalizer. We also reformulated the algorithm using the element function $\varphi_{k}(t)=I_{1}+I_{2}+I_{3}$ as a building block to define the antialiased and normalized complex exponential pulse $p_{h}(t)$:

[TABLE]

where $\{h_{k}\}_{k=0}^{6}$ is given by Eq. 97. The factor $2t_{w}h_{0}$ is for the gain normalization. For each $k>0$, the explicit form is as follows:

[TABLE]

where $\alpha=\pi/t_{w}$. The functions $I_{1}$, $I_{2}$, and $I_{3}$ are defined in $(-t_{w},t_{w}]$, $(t_{w},1+t_{w})$, and $[1-t_{w},1+t_{w})$ ,respectively. They are 0 outside. Note that the second interval overlaps with the third one. The constant $\beta$ can be a complex number depending on the piece. For example, the first piece of the L–F model, $\beta$ is a complex number, and $\Im[p_{h}(t)]$ provides the antialiased result. For $k=0$, the zeroth-order antialiased polynomial pulse is applicable.

We conducted a set of tests using this revised L-F model. Similar to Fig. 3, the noise level around the fundamental component was about $-180$ dB from the peak harmonics level. This revised model is also available as a set of open access MATLAB functions.

7 Conclusions

We formulated and implemented a procedure to generate aliasing-free glottal source model output. We antialiased the Fujisaki–Ljungqvist model using a newly designed cosine series antialiasing function, followed by an IIR digital equalizer. We also revised our antialiased L–F model using the six-term cosine series and the IIR equalizer. The proposed procedure is general enough to be applicable to other glottal source models and any signal models consisting of polynomial and complex exponential segments. These antialiased models are available as open access MATALB procedures with interactive GUI tools for education and research in speech science[10]. The antialiased glottal excitation signals also provide a reliable and flexible means to test $f_{\rm o}$ extractors and source aperiodicity analysis procedures.

8 Acknowledgments

This work was supported by JSPS KAKENHI Grant Numbers JP15H03207, JP15H02726 and JP16K12464.

Appendix A Supplement

The main text was accepted for publication in Proc. of Interspeech2017. This appendix provides materials that were dropped owing to paucity of space.

A.1 Test signals

We used the frequency-modulated $f_{\rm o}$ trajectory $f_{0i}(t)$ for generating the figures shown in this article. The following equation provides the details:

[TABLE]

where $f_{d}$ determines the depth of the vibrato in terms of the musical cent, and $f_{m}$ determines the rate of the vibrato. We used $f_{d}=6$ (cent) and $f_{m}=5.2$ (Hz) to generate the figures. The Fujisaki–Ljungqvist model parameters $\{A,B,C,R,F,D\}$ were set as follows:

[TABLE]

We also generated test signals using the antialiased L–F model. The same frequency-modulated $f_{\rm o}$ trajectory was used. The L–F model parameters $\{t_{p},t_{e},t_{a},t_{c}\}$ were set as follows:

[TABLE]

where the specific values are the average value of the modal voices reported in Reference[29].

A.2 Spectrogram of Fujisaki–Ljungqvist model

Figure 4 shows the spectrogram of the generated excitation source signals using direct discretization, the Nuttall-11 window, and the proposed six-term cosine series. The spectrum slices in Fig. 3 are sampled from these spectrograms.

We used a self-convolved version of a Nuttall’s window for calculating these spectrograms. It is the 12th item of Table II of Reference[28]. This self-convolved window has the maximum sidelobe level at $-186.64$ dB and a 36 dB/oct decay rate. The low sidelobe level and the steep decay rate of this window allow us to inspect low-level spurious of the proposed procedure. The window length and frame shift were 40 ms and 2 ms, respectively.

A.3 IIR equalizer design

The target equalizer shape is designed to equalize attenuation up to 68 dB for six-term and 58 dB for five-term cosine series, respectively. The FFT buffer length was 32,768, and the length of the FIR response was 161 taps. The original equalizer shape, represented as an absolute spectrum, was converted to the FIR response and truncated using one of the Nuttall windows[28]. The 12th item of Table II of the reference was used. It has a maximum sidelobe level of $-93.32$ dB and a decay rate of 18 dB/oct.

The IIR equalizer was designed using the autocorrelation coefficients of this truncated FIR response, by applying LPC analysis.

Figure 5 shows the frequency response of each equalizer and the pole locations. The poles are not close to the unit circle, indicating the numerical stability of the equalizers.

A.4 MATLAB functions

The following MATLAB functions are prepared to generate the excitation source signals.

**antiAliasedPolynomialSegmentR **

Generate a time normalized segment of an antialiased polynomial pulse.

antialiasedFLmodelSingleR

Generate one cycle of an antialiased F-L model excitation signal.

AAFjLjmodelFrom$f_{\rm o}$TrajectoryR

Generate an antialiased F-L model excitation signal using the given $f_{\rm o}$ trajectory and constant F-L model parameters.

AAFjLjmodelFrom$f_{\rm o}$TrajectoryTVR

Generate an antialiased F-L model excitation signal using the given $f_{\rm o}$ trajectory and time varying F-L model parameters.

We used AAFjLjmodelFrom$f_{\rm o}$TrajectoryR to generate the test signal used to draw Figs. 3 and 4.

Figure 6 shows an example of an excitation signal with time varying Fujisaki–Ljungqvist model parameters. The total airflow of each pitch cycle is kept constant.

A.5 Application to aliasing-free L–F model

Figure 7 shows the L–F model waveform with time parameters to define it. There is a typo in the L–F model definition given on page 6 of the original reference[4], while other descriptions in the reference are correct. The fixed equation that defines the L-F model is as follows:

[TABLE]

where $E(t)$ is defined as the time derivative of the glottal airflow $U_{g}(t)$. It is convenient to normalize the time axis by $T_{0}$ and the amplitude by $E_{e}$ without loss of generality. The coefficients that can then be determined from the design parameters are $E_{0}/E_{e},\alpha,\omega_{g}$, and $\beta$. Because the airflow is zero while the vocal fold is closed, the following constraint holds:

[TABLE]

The following steps provide the parameter values. First, substitute $t=t_{c}$ in Eq. 106. Solving it yields $\beta$. Then, use Eq. 107 and $\omega_{g}=t_{p}/\pi$ to determine $\alpha$. We used the numerical optimization function fzero of MATLAB to implement these steps.

A.5.1 Spectrum slice and spectrogram

Figure 8 shows spectrum slice and spectrogram of a generated test signal. The same frequency-modulated $f_{\rm o}$ trajectory was used and the details are given in Appendix A.1