MP8: LPC

In this MP, we will learn how to convert your voice into the voice of a robot. We'll learn how to change the robot's pitch and voice quality, to get slightly different tones of voice.

Here are the contents:

1. Data
2. Linear predictive coding (LPC)
3. LPC synthesis
4. Robot speech
5. Homework

1. Data

To start out with, let's load an arbitrary speech waveform, and create its spectrogram. We'll use librosa to load and process audio, since librosa automatically resamples the file to a desired sampling rate while loading it.

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf
import IPython.display
import librosa

sr=8000
speech, fs = librosa.load(path='speech_waveform.wav', sr=sr)
print('Asked librosa to use sampling rate',sr,'and it provided sampling rate',fs,'with',len(speech),'samples.')
IPython.display.Audio(data=speech,rate=fs)

Asked librosa to use sampling rate 8000 and it provided sampling rate 8000 with 44705 samples.

Most of our analyses will be short-time analyses, e.g., short-time Fourier transform, short-time energy, short-time pitch, and so on. To compute short-time analyses, we need to divide the speech signal into overlapping short frames:

help(librosa.util.frame)

Help on function frame in module librosa.util.utils:

frame(
    x: 'np.ndarray',
    *,
    frame_length: 'int',
    hop_length: 'int',
    axis: 'int' = -1,
    writeable: 'bool' = False,
    subok: 'bool' = False
) -> 'np.ndarray'
    Slice a data array into (overlapping) frames.

    This implementation uses low-level stride manipulation to avoid
    making a copy of the data.  The resulting frame representation
    is a new view of the same input data.

    For example, a one-dimensional input ``x = [0, 1, 2, 3, 4, 5, 6]``
    can be framed with frame length 3 and hop length 2 in two ways.
    The first (``axis=-1``), results in the array ``x_frames``::

        [[0, 2, 4],
         [1, 3, 5],
         [2, 4, 6]]

    where each column ``x_frames[:, i]`` contains a contiguous slice of
    the input ``x[i * hop_length : i * hop_length + frame_length]``.

    The second way (``axis=0``) results in the array ``x_frames``::

        [[0, 1, 2],
         [2, 3, 4],
         [4, 5, 6]]

    where each row ``x_frames[i]`` contains a contiguous slice of the input.

    This generalizes to higher dimensional inputs, as shown in the examples below.
    In general, the framing operation increments by 1 the number of dimensions,
    adding a new "frame axis" either before the framing axis (if ``axis < 0``)
    or after the framing axis (if ``axis >= 0``).

    Parameters
    ----------
    x : np.ndarray
        Array to frame
    frame_length : int > 0 [scalar]
        Length of the frame
    hop_length : int > 0 [scalar]
        Number of steps to advance between frames
    axis : int
        The axis along which to frame.
    writeable : bool
        If ``False``, then the framed view of ``x`` is read-only.
        If ``True``, then the framed view is read-write.  Note that writing to the framed view
        will also write to the input array ``x`` in this case.
    subok : bool
        If True, sub-classes will be passed-through, otherwise the returned array will be
        forced to be a base-class array (default).

    Returns
    -------
    x_frames : np.ndarray [shape=(..., frame_length, N_FRAMES, ...)]
        A framed view of ``x``, for example with ``axis=-1`` (framing on the last dimension)::

            x_frames[..., j] == x[..., j * hop_length : j * hop_length + frame_length]

        If ``axis=0`` (framing on the first dimension), then::

            x_frames[j] = x[j * hop_length : j * hop_length + frame_length]

    Raises
    ------
    ParameterError
        If ``x.shape[axis] < frame_length``, there is not enough data to fill one frame.

        If ``hop_length < 1``, frames cannot advance.

    See Also
    --------
    numpy.lib.stride_tricks.as_strided

    Examples
    --------
    Extract 2048-sample frames from monophonic signal with a hop of 64 samples per frame

    >>> y, sr = librosa.load(librosa.ex('trumpet'))
    >>> frames = librosa.util.frame(y, frame_length=2048, hop_length=64)
    >>> frames
    array([[-1.407e-03, -2.604e-02, ..., -1.795e-05, -8.108e-06],
           [-4.461e-04, -3.721e-02, ..., -1.573e-05, -1.652e-05],
           ...,
           [ 7.960e-02, -2.335e-01, ..., -6.815e-06,  1.266e-05],
           [ 9.568e-02, -1.252e-01, ...,  7.397e-06, -1.921e-05]],
          dtype=float32)
    >>> y.shape
    (117601,)

    >>> frames.shape
    (2048, 1806)

    Or frame along the first axis instead of the last:

    >>> frames = librosa.util.frame(y, frame_length=2048, hop_length=64, axis=0)
    >>> frames.shape
    (1806, 2048)

    Frame a stereo signal:

    >>> y, sr = librosa.load(librosa.ex('trumpet', hq=True), mono=False)
    >>> y.shape
    (2, 117601)
    >>> frames = librosa.util.frame(y, frame_length=2048, hop_length=64)
    (2, 2048, 1806)

    Carve an STFT into fixed-length patches of 32 frames with 50% overlap

    >>> y, sr = librosa.load(librosa.ex('trumpet'))
    >>> S = np.abs(librosa.stft(y))
    >>> S.shape
    (1025, 230)
    >>> S_patch = librosa.util.frame(S, frame_length=32, hop_length=16)
    >>> S_patch.shape
    (1025, 32, 13)
    >>> # The first patch contains the first 32 frames of S
    >>> np.allclose(S_patch[:, :, 0], S[:, :32])
    True
    >>> # The second patch contains frames 16 to 16+32=48, and so on
    >>> np.allclose(S_patch[:, :, 1], S[:, 16:48])
    True

frame_length = int(np.round(fs*0.025))    # 25ms frames
frame_skip = int(np.round(fs*0.01))       # with a 10ms hop between frame starting points
frames = librosa.util.frame(speech, frame_length=frame_length, hop_length=frame_skip, axis=0)
print('There should be %d frames of length %d:'%(1+(len(speech)-frame_length)//frame_skip,frame_length))
print('The computation result includes %d frames of length %d'%(frames.shape[0], frames.shape[1]))

There should be 557 frames of length 200:
The computation result includes 557 frames of length 200

Now, to calculate the spectrogram, we just need to take the FFT of each frame, then compute its log magnitude in decibels.

• To avoid underflow, we'll normalize the spectrogram by its maximum, then threshold at 1e-6.
• We will discard the negative frequencies, which occur in FFT bins frame_length/2 through frame_length.

mstft = np.abs(np.fft.fft(frames))
sgram = 20*np.log10(np.maximum(mstft/np.amax(mstft),1e-6))[:,:int(frame_length/2)]

fig = plt.figure(figsize=(14,4),layout='tight')
ax = fig.subplots(2,1)
ax[0].plot(np.arange(len(speech))/fs,speech)
ax[0].set_title('Speech waveform')
librosa.display.specshow(sgram.transpose(),ax=ax[1],sr=fs,hop_length=frame_skip,x_axis='s',y_axis='hz')
ax[1].set_title('Spectrogram')

Text(0.5, 1.0, 'Spectrogram')

fig = plt.figure(figsize=(14,4),layout='tight')
ax = fig.subplots(2,2)
ax[0,0].plot(np.arange(frame_length)/fs,frames[200])
ax[0,0].set_title('Speech at t=2.0: "wave"')
ax[0,1].plot(np.linspace(0,fs/2,int(frame_length/2)),sgram[200])
ax[0,1].set_title('Spectrum at t=2.0')
ax[1,0].plot(np.arange(frame_length)/fs,frames[230])
ax[1,0].set_title('Speech at t=2.3: "form"')
ax[1,0].set_xlabel('Time (s)')
ax[1,1].plot(np.linspace(0,fs/2,int(frame_length/2)),sgram[230])
ax[1,1].set_title('Spectrum at t=2.3')
ax[1,1].set_xlabel('Frequency (Hz)')

Text(0.5, 0, 'Frequency (Hz)')

2. Linear Predictive Coding (LPC)

Steady-state vowels can be synthesized by passing an excitation signal, $e[n]$, through a series of resonant filters like this:

$$e[n]\rightarrow \fbox{$\stackrel{F_1}{\text{resonator}}$}\rightarrow \fbox{$\stackrel{F_2}{\text{resonator}}$}\rightarrow \fbox{$\stackrel{F_3}{\text{resonator}}$}\rightarrow \fbox{$\stackrel{F_4}{\text{resonator}}$}\rightarrow s[n]$$

• The excitation signal, $e[n]$, has one impulse per pitch period
• The resonators compute feedback, $y[n]=e[n]+a[1]y[n-1]+a[2]y[n-2]$.

Instead of computing resonator coefficients from knowledge about speech production, we will compute them by analyzing a real speech signal. That will be easier if we group all $4\times 2=8$ of the coefficients into a single equation. Just to make sure we have enough modeling power, we'll use 10 feedback loops:

$$s[n]=e[n]-\sum_{k=1}^{10} a[k] s[n-k]$$

$$e[n]\rightarrow \fbox{Order-10 Resonator}\rightarrow s[n]$$

The coefficients $a_k$ are called a linear predictive code (LPC) because they code all of the information about the resonant frequencies of the speech.

• For written exams, you need to know how to compute the LPC coefficients from the matrix equation $\mathbf{r}=\mathbf{R}\mathbf{a}$, where $\mathbf{r}$ and $\mathbf{R}$ are the autocorrelation vector and matrix, and $\mathbf{a}$ is the vector of LPC coefficients.
• For this machine problem, however, you can use the function librosa.lpc to solve the problem for you.
• You should write a function submitted.lpc_analysis that divides the speech into frames, then applies librosa.lpc to compute the coefficients. You can use librosa.util.frame to compute the frames if you wish (that's the easiest method, but other methods are almost as easy).

import importlib, submitted
help(submitted.lpc_analysis)

Help on function lpc_analysis in module submitted:

lpc_analysis(speech, L, S, order)
    Chop speech into frames, and compute LPC coefficients.
    You are allowed to use the following functions, which might be useful:
    librosa.util.frame (axis=0), librosa.lpc (axis=-1).

    @param:
    speech (T,) - input speech signal
    L (1,) - frame length in samples
    S (1,) - number of samples skipped between frame starts

    @return:
    frames (1+(T-L)//S,L) - frames of speech
    coefficients (1+(T-L)//S,order+1)

importlib.reload(submitted)
order = 10
frames, A = submitted.lpc_analysis(speech, frame_length, frame_skip, order)
print('The frames matrix has shape:',frames.shape)
print('The coefficients matrix has shape:',A.shape)

The frames matrix has shape: (557, 200)
The coefficients matrix has shape: (557, 11)

fig = plt.figure(figsize=(14,4),layout='tight')
ax = fig.subplots(2,2)
ax[0,0].stem(np.arange(frame_length)/fs,frames[200])
ax[0,0].set_title('Speech at t=2.0: "wave"')
ax[0,1].stem(A[200])
ax[0,1].set_title('LPC coefficients at t=2.0')
ax[1,0].stem(np.arange(frame_length)/fs,frames[230])
ax[1,0].set_title('Speech at t=2.3: "form"')
ax[1,0].set_xlabel('Time (s)')
ax[1,1].stem(A[230])
ax[1,1].set_title('LPC coefficients at t=2.3')
ax[1,1].set_xlabel('Time (samples)')

Text(0.5, 0, 'Time (samples)')

The LPC coefficients by themselves are not very useful, but they can be used for many useful things. For example, we can find the excitation by subtracting the predicted speech, $\sum_{k=1}^{10}a[k]s[n-k]$, from the original speech:

$$e[n] = s[n] + \sum_{k=1}^{10}a[k] s[n-k]$$

librosa.lpc, by default, sets $a[0]=1$ in each frame, so we can write the above equation as

$$e[n] = \sum_{k=0}^{10}a[k] s[n-k]$$

Note that, in the equation above, the first $10$ samples of $e[n]$ in every frame depend on samples of $s[n]$ that come from the previous frame. To avoid worrying about that problem, we just won't compute the first 10 samples of $e[n]$ in every frame:

importlib.reload(submitted)
help(submitted.lpc_inverse)

Help on function lpc_inverse in module submitted:

lpc_inverse(frames, A)
    Inverse filter speech frames to estimate the excitation signal.

    @param:
    frames (nframes,framelength) - speech frames
    A (nframes,order+1) - LPC coefficients (A[m,0]=1 always)

    @return:
    excitation (nframes,framelength-order) - linear prediction excitation frames

importlib.reload(submitted)
excitation = submitted.lpc_inverse(frames, A)
print('frames matrix has shape:',frames.shape)
print('coefficients matrix has shape:', A.shape)
print('excitation matrix has shape:', excitation.shape)

frames matrix has shape: (557, 200)
coefficients matrix has shape: (557, 11)
excitation matrix has shape: (557, 190)

fig = plt.figure(figsize=(14,4),layout='tight')
ax = fig.subplots(2,2)
ax[0,0].stem(np.arange(frame_length)/fs,frames[200])
ax[0,0].set_title('Speech at t=2.0: "wave"')
ax[0,1].stem(np.arange(frame_length-order)/fs,excitation[200])
ax[0,1].set_title('LPC excitation at t=2.0')
ax[1,0].stem(np.arange(frame_length)/fs,frames[230])
ax[1,0].set_title('Speech at t=2.3: "form"')
ax[1,0].set_xlabel('Time (s)')
ax[1,1].stem(np.arange(frame_length-order)/fs,excitation[230])
ax[1,1].set_title('LPC excitation at t=2.3')
ax[1,1].set_xlabel('Time (seconds)')

Text(0.5, 0, 'Time (seconds)')

Notice two things about the excitation:

1. It's almost an impulse train! This is voiced speech, so the excitation looks a lot like an impulse train.
2. Different frames have different energy.

Let's compute the energy in each frame:

$$G = \sqrt{\frac{1}{N}\sum_{n=0}^{N-1} e^2[n]}$$

importlib.reload(submitted)
help(submitted.lpc_gain)

Help on function lpc_gain in module submitted:

lpc_gain(excitation)
    Calculate the LPC filter gain for each frame.

    @param:
    excitation (nframes,framelength) - inverse-filtered speech signal

    @return:
    gain (nframes,) - gain[m] is the RMS energy of excitation[m,:]

importlib.reload(submitted)
gain = submitted.lpc_gain(excitation)
print('The gain matrix has shape:',gain.shape)

The gain matrix has shape: (557,)

When you calculated the excitation, you might have noticed that it's just the convolution of the speech signal with the LPC filter coefficients (where $a[0]=1$):

$$e[n] = s[n]\ast a[n]$$

For that reason, the excitation spectrum is related to the speech spectrum as $E(\omega)=S(\omega)A(\omega)$, which means that

$$S(\omega)=\frac{E(\omega)}{A(\omega)}$$

If we want to ignore the details of the excitation, and keep only its overall energy level, we can approximate the speech spectrum as:

$$S(\omega)\approx \frac{G}{A(\omega)}$$

nframes, nsamps = frames.shape
lpc_sgram = np.zeros((nframes,int(nsamps/2)))
for frame in range(nframes):
    lpc_sgram[frame,:] = 20*np.log10(gain[frame] / np.abs(np.fft.fft(A[frame,:], n=nsamps)))[:int(nsamps/2)]
    

fig = plt.figure(figsize=(14,6),layout='tight')
ax = fig.subplots(2,1)
librosa.display.specshow(sgram.transpose(),ax=ax[0],sr=fs,hop_length=frame_skip,x_axis='s',y_axis='hz')
ax[0].set_title('Spectrogram, including all details of the excitation')
librosa.display.specshow(lpc_sgram.transpose(),ax=ax[1],sr=fs,hop_length=frame_skip,x_axis='s',y_axis='hz')
ax[1].set_title('LPC Spectrogram, showing only the information captured by $G/A(\omega)$')

Text(0.5, 1.0, 'LPC Spectrogram, showing only the information captured by $G/A(\\omega)$')

3. LPC synthesis

To synthesize speech again, we just need to create an excitation signal, and then pass it through the resonator:

$$e[n]\rightarrow\fbox{Resonant Filter}\rightarrow s[n]$$

$$s[n] = e[n] - \sum_{k=1}^{10} a[k]s[n-k]$$

To get the excitation signal, we need to be careful. When we computed excitation, before, we used overlapping frames, so we need to:

1. Delete the overlapping part from each frame (keep only the last frame_skip samples), then
2. concatenate the frames. We can do this using np.hstack.

e = np.hstack(excitation[:,-frame_skip:])
print('There are %d frames, with a skip of %d samples, total length is %d  (should be %d)'%(nframes,frame_skip,len(e),nframes*frame_skip))

There are 557 frames, with a skip of 80 samples, total length is 44560  (should be 44560)

Now we need to use the synthesis equation. This is a little tricky, because each sample of $s[n]$ is computed using the LPC coefficients from its own frame, but it might depend on samples of $s[n]$ that were computed in the previous frame.

Suppose we use $a[m,k]$ to denote the $k^{\text{th}}$ coefficient in the $m^{\text{th}}$ frame. If $S=$frame_skip, then we can write $m=\text{int}\left(\frac{n}{S}\right)$ to mean that the $n^{\text{th}}$ speech sample is part of the $m^{\text{th}}$ frame. Then the LPC synthesis equation is

$$s[n] = e[n]-\sum_{k=1}^{10} a\left[\text{int}\left(\frac{n}{S}\right),k\right] s[n-k]$$

importlib.reload(submitted)
help(submitted.lpc_synthesis)

Help on function lpc_synthesis in module submitted:

lpc_synthesis(e, A, S)
    Synthesize speech from LPC residual and coefficients.

    @param:
    e (duration,) - excitation signal
    A (nframes,order+1) - linear predictive coefficients from each frames
    S (1,) - frame skip, in samples

    @returns:
    synthesis (duration) - synthetic speech waveform

importlib.reload(submitted)
synthesis = submitted.lpc_synthesis(e, A, frame_skip)
IPython.display.Audio(data=synthesis, rate=fs)

synth_frames = np.array([synthesis[m*frame_skip:m*frame_skip+frame_length] for m in range(int((len(synthesis)-frame_length)/frame_skip)) ])
synth_mstft = np.abs(np.fft.fft(synth_frames))
synth_sgram = 20*np.log10(np.maximum(synth_mstft,1e-6*np.amax(synth_mstft)))[:,:int(frame_length/2)]
fig = plt.figure(figsize=(14,8),layout='tight')
ax = fig.subplots(4,1)
ax[0].plot(np.arange(len(speech))/fs,speech)
ax[0].set_title('Speech waveform')
librosa.display.specshow(sgram.transpose(),ax=ax[1],sr=fs,hop_length=frame_skip,x_axis='s',y_axis='hz')
ax[1].set_title('Spectrogram')
ax[2].plot(np.arange(len(synthesis))/fs,synthesis)
ax[2].set_title('Synthesized speech waveform')
librosa.display.specshow(synth_sgram.transpose(),ax=ax[3],sr=fs,hop_length=frame_skip,x_axis='s',y_axis='hz')
ax[3].set_title('Synthesized Spectrogram')

Text(0.5, 1.0, 'Synthesized Spectrogram')

4. Robot speech

Now that we know how to resynthesize speech, let's modify it so it sounds like a robot.

• Change the excitation so it's a perfect monotone, always exactly $F_0=100$Hz.
• Change the excitation so it's perfectly voiced, with no breathiness.

We can do this by setting $e_{\text{robot}}[n]=G p[n]$, where $G$ is the excitation gain in each frame, and $p[n]$ is an impulse train with a period of $T_0=100/F_s$:

$$p[n]=\left\{\begin{array}{ll}1 & n=\text{integer multiple of}~T_0\\0&\text{otherwise}\end{array}\right.$$

p, e_robot = np.zeros(len(e)), np.zeros(len(e))

p[::int(fs/100)] = 1

for n in range(len(p)):
    e_robot[n] = gain[int(n/frame_skip)] * p[n]

s_robot = submitted.lpc_synthesis(e_robot, A, frame_skip)

IPython.display.Audio(data=s_robot, rate=fs)

Let's make it a female robot! We can do that by setting $F_0$ to a value that sounds female, for example, $F_0=200$Hz.

p, e_robot = np.zeros(len(e)), np.zeros(len(e))

p[::int(fs/200)] = 1

for n in range(len(p)):
    e_robot[n] = gain[int(n/frame_skip)] * p[n]

s_robot = submitted.lpc_synthesis(e_robot, A, frame_skip)

IPython.display.Audio(data=s_robot, rate=fs)

As you can imagine, there are lots of other modifications we can make. For example, let's make the pitch frequency rise for the first half, then fall for the second half:

$$F_0(t) =\min\left(200+40t,400-40t\right)$$

t = np.arange(len(e))/fs
F0 = np.minimum(200+40*t,400-40*t)
T0 = fs/F0

fig = plt.figure(figsize=(14,4),layout='tight')
ax = fig.subplots(2,1)
ax[0].plot(np.arange(len(F0))/fs,F0)
ax[0].set_title('F0 in Hertz')
ax[1].plot(np.arange(len(T0))/fs,T0)
ax[1].set_title('T0 in samples')

Text(0.5, 1.0, 'T0 in samples')

e_robot = np.zeros(len(e))
n = 0
while n < len(e_robot):
    e_robot[n] = gain[int(n/frame_skip)]
    n += int(np.round(T0[n]))
fig = plt.figure(figsize=(14,6),layout='tight')
ax = fig.subplots(3,1)
ax[0].plot(np.arange(len(e_robot))/fs,e_robot)
ax[0].set_title('Robot excitation with pitch frequency that rises then falls')
ax[1].plot(np.arange(frame_length)/fs,e_robot[2*fs:2*fs+frame_length])
ax[1].set_title('Robot excitation in the frame at t=2.0s')
ax[2].plot(np.arange(frame_length)/fs,e_robot[int(2.3*fs):int(2.3*fs)+frame_length])
ax[2].set_title('Robot excitation in the frame at t=2.3s.  Notice: The pitch period is shorter here!')

Text(0.5, 1.0, 'Robot excitation in the frame at t=2.3s.  Notice: The pitch period is shorter here!')

s_robot = submitted.lpc_synthesis(e_robot, A, frame_skip)

IPython.display.Audio(data=s_robot, rate=fs)

Homework

Once you have finished writing the functions in submitted.py, test them using grade.py before submitting your code to the autograder:

!python grade.py

....
----------------------------------------------------------------------
Ran 4 tests in 0.601s

OK