Figure 1: Discrete-Time Speech Production Model

From a signal processing standpoint, it is very useful to think of speech production in terms of a model, as in Figure 1. The model shown is the simplest of its kind, but it includes all the principal components. The excitations for voiced and unvoiced speech are represented by an impulse train and white noise generator, respectively. The pitch of voiced speech is controlled by the spacing between impulses, TpTp, and the amplitude (volume) of the excitation is controlled by the gain factor GG.

As the acoustical excitation travels from its source (vocal cords, or a constriction), the shape of the vocal tract alters the spectral content of the signal. The most prominent effect is the formation of resonances, which intensifies the signal energy at certain frequencies (called formants). As we learned in the Digital Filter Design lab, the amplification of certain frequencies may be achieved with a linear filter by an appropriate placement of poles in the transfer function. This is why the filter in our speech model utilizes an all-pole LTI filter. A more accurate model might include a few zeros in the transfer function, but if the order of the filter is chosen appropriately, the all-pole model is sufficient. The primary reason for using the all-pole model is the distinct computational advantage in calculating the filter coefficients, as will be discussed shortly.

Recall that the transfer function of an all-pole filter has the form

where PP is the order of the filter. This is an IIR filter that may be implemented with a recursive difference equation. With the input G·x(n)G·x(n), the speech signal s(n)s(n) may be written as

Keep in mind that the filter coefficients will change continuously as the shape of the vocal tract changes, but speech segments of an appropriately small length may be approximated by a time-invariant model.

This speech model is used in a variety of speech processing applications, including methods of speech recognition, speech coding for transmission, and speech synthesis. Each of these applications of the model involves dividing the speech signal into short segments, over which the filter coefficients are almost constant. For example, in speech transmission the bit rate can be significantly reduced by dividing the signal up into segments, computing and sending the model parameters for each segment (filter coefficients, gain, etc.), and re-synthesizing the signal at the receiving end, using a model similar to Figure 1. Most telephone systems use some form of this approach. Another example is speech recognition. Most recognition methods involve comparisons between short segments of the speech signals, and the filter coefficients of this model are often used in computing the “difference" between segments.

Download the file coeff.mat for the following section.

Download the file coeff.mat and load it into the Matlab workspace using the load command. This will load three sets of filter coefficients: A1A1, A2A2, and A3A3 for the vocal tract model in Equation 1 and Equation 2. Each vector contains coefficients {a1,a2,⋯,a15}{a1,a2,⋯,a15} for an all-pole filter of order 15.

We will now synthesize voiced speech segments for each of these sets of coefficients. First write a Matlab function x=exciteV(N,Np) which creates a length NN excitation for voiced speech, with a pitch period of Np samples. The output vector xx should contain a discrete-time impulse train with period Np (e.g. [1 0 0 ⋯⋯ 0 1 0 0 ⋯⋯]).

Assuming a sampling frequency of 8 kHz (0.125 ms/sample), create a 40 millisecond-long excitation with a pitch period of 8 ms, and filter it using Equation 2 for each set of coefficients. For this, you may use the command

s = filter(1,[1 -A],x)

where AA is the row vector of filter coefficients (see Matlab's help on filter for details). Plot each of the three filtered signals. Use subplot() and orient tall to place them in the same figure.

We will now compute the frequency response of each of these filters. The frequency response may be obtained by evaluating Equation 1 at points along z=ejωz=ejω. Matlab will compute this with the command [H,W]=freqz(1,[1 -A],512) , where AA is the vector of coefficients. Plot the magnitude of each response versus frequency in Hertz. Use subplot() and orient tall to plot them in the same figure.

The location of the peaks in the spectrum correspond to the formant frequencies. For each vowel signal, estimate the first three formants (in Hz) and list them in the figure.

Now generate the three signals again, but use an excitation which is 1-2 seconds long. Listen to the filtered signals using soundsc. Can you hear qualitative differences in the signals? Can you identify the vowel sounds?

INLAB REPORT

Hand in the following:

Suppose we have a discrete-time random process {...,S-1,S0,S1,S2,...}{...,S-1,S0,S1,S2,...} whose elements have some degree of correlation. The goal of forward linear prediction is to predict the sample SnSn using a linear combination of the previous PP samples.

PP is called the order of the predictor. We may represent the error of predicting SnSn by a random sequence enen.

An optimal set of prediction coefficients akak for Equation 4 may be determined by minimizing the mean-square error E[en2]E[en2]. Note that since the error is generally a function of nn, the prediction coefficients will also be functions of nn. To simplify notation, let us first define the following column vectors.

Then,

The second and third terms of Equation 7 may be written in terms of the autocorrelation sequence rSS(k,l)rSS(k,l).

Substituting into Equation 7, the mean-square error may be written as

Note that while aa and rSrS are vectors, and RSRS is a matrix, the expression in Equation 10 is still a scalar quantity.

To find the optimal akak coefficients, which we will call a^a^, we differentiate Equation 10 with respect to the vector aa (compute the gradient), and set it equal to the zero vector.

Solving,

The vector equation in Equation 12 is a system of PP scalar linear equations, which may be solved by inverting the matrix RSRS.

Note from Equation 8 and Equation 9 that rSrS and RSRS are generally functions of nn. However, if SnSn is wide-sense stationary, the autocorrelation function is only dependent on the difference between the two indices, rSS(k,l)=rSS(|k-l|)rSS(k,l)=rSS(|k-l|). Then RSRS and rSrS are no longer dependent on nn, and may be written as follows.

Therefore, if SnSn is wide-sense stationary, the optimal akak coefficients do not depend on nn. In this case, it is also important to note that RSRS is a Toeplitz (constant along diagonals) and symmetric matrix, which allows Equation 12 to be solved efficiently using the Levinson-Durbin algorithm (see [1]). This property is essential for many real-time applications of linear prediction.

An important question has yet to be addressed. The solution in Equation 12 to the linear prediction problem depends entirely on the autocorrelation sequence. How do we estimate the autocorrelation of a speech signal? Recall that the applications to which we are applying LPC involve dividing the speech signal up into short segments and computing the filter coefficients for each segment. Therefore we need to consider the problem of estimating the autocorrelation for a short segment of the signal. In LPC, the following "biased" autocorrelation estimate is often used.

Here we are assuming we have a length NN segment which starts at n=0n=0. Note that this is the single-parameter form of the autocorrelation sequence, so that the forms in Equation 13 and Equation 14 may be used for rSrS and RSRS.

Download the file test.mat for this exercise.

Write a function coef=mylpc(x,P) which will compute the order-P LPC coefficients for the column vector xx, using the autocorrelation method (“lpc" is a built-in Matlab function, so use the name mylpc). Consider the input vector xx as a speech segment, in other words do not divide it up into pieces. The output vector coef should be a column vector containing the PP coefficients {a^1,a^2,⋯,a^P}{a^1,a^2,⋯,a^P}. In your function you should do the following:

To test your function, download the file test.mat, and load it into Matlab. This file contains two vectors: a signal xx and its order-15 LPC coefficients aa. Use your function to compute the order-15 LPC coefficients of xx, and compare the result to the vector aa.