GOLF: A Singing Voice Synthesiser with Glottal Flow Wavetables and LPC Filters

2.1 Glottal source‑filter model

In the source‑filter model, we have the following simplified voice production model:

where represents the periodic vibration from the vocal folds, represents random components of the glottal source, represents the vocal‑tract filter and represents the radiation at the lips (Degottex, 2010). Since this formulation is linear, the radiation filter and the glottal pulse can be merged into a single source called the radiated glottal pulse. If we assume is a first‑order differentiator (Lu and Smith III, 2000), then is the derivative of the glottal pulse, which can be described by the LF model (Fant et al., 1985), a four‑parameter model of glottal flow. is usually modelled with a linear prediction (LP) filter.

2.2 Linear predictive coding

Linear predictive coding (LPC) assumes that the current voice sample can be predicted from a finite number of previous samples by a linear combination with residuals :

where are the linear prediction coefficients. In this work, we only evaluate (so the audio length is ) and assume that, when , is always zeros. In our case, this is the same as filtering the residuals, equivalent to the glottal source, with an ‑order all‑pole filter. This filter has an infinite impulse response (IIR). We can use the LPC filter to represent the response of the vocal tract if the vocal tract is approximated by a series of cylindrical tubes with varying diameters (Markel and Gray, 1976), providing a physical interpretation.

Using LPC for neural audio synthesis is not new (Oh et al., 2020; Subramani et al., 2022; Valin and Skoglund, 2019), and works have been conducted to incorporate IIR filters and train them jointly with deep learning models (Bhattacharya et al., 2020; Colonel et al., 2022; Kim et al., 2022; Kuznetsov et al., 2020; Nercessian, 2020; Nercessian et al., 2021; Steinmetz et al., 2022). The difficulty of training IIR in a deep learning framework (e.g., PyTorch) using Equation 2 is that its computation is recursive, i.e., the output at each step depends on the previous results, and to make the calculation differentiable, separated tensors are allocated in each step. This results in significant computational overhead for allocating fragmented memory, which can lead to performance issues, especially for long sequences.

One way to mitigate this issue is to pre‑allocate a shared tensor before computation as a buffer and reuse it for each step. However, in‑place modification is mostly not differentiable in these frameworks.¹ Some studies sidestep the recursion by approximating IIR in the frequency domain using discrete Fourier transform (DFT). This is usually called frequency sampling (FS; Colonel et al., 2022; Kim et al., 2022; Nercessian, 2020; Nercessian et al., 2021; Oh et al., 2020; Steinmetz et al., 2022; Wright and Valimaki, 2022). The resulting finite impulse response (FIR) is a folded version of the original IIR, introducing errors called time‑domain aliasing (Smith, 2024). Directly multiplying the DFT of the input with the transfer function of the filter (Oh et al., 2020) does not equal its time‑domain convolution but a circular convolution with the folded FIR. The aliasing effect cannot be avoided but only mitigated by increasing the DFT resolution, which also increases the filter’s computational cost.

Moreover, the IIRs used in practice are usually low order, such as the classic biquad filters (Colonel et al., 2022; Nercessian, 2020; Nercessian et al., 2021) and one‑pole smoothing filters (Wright and Valimaki, 2022). In these cases, filtering the signal directly is faster than using the FS method, especially when the sequence is long.

3.1 Glottal flow wavetables

We adopt the transformed‑LF model (Fant, 1995) for generating glottal pulses. This model re‑parametrises the LF model (Fant et al., 1985) using just one parameter, , which was found to correspond to the perceived vocal effort well and covers a wide range of different glottal flow shapes.

To save time computing the pulses, we pre‑compute a set of the glottal flow and synthesise the waveform using wavetable synthesis. Thus, the computation is limited to indexing the memory plus naive interpolation. We sample values of with equal spacing inside following the value range suggested by Degottex (2010). We calculate the flow derivative function in continuous time for each sampled and then sample points in one period to get its discrete version. The details for calculating are given in Gobl (2017). By stacking these sampled glottal flows, we build wavetables , with each row containing one period of a sampled glottal pulse (Figure 1). The rows are sorted based on .

Figure 1

The model generates glottal pulses by linearly interpolating the two axes. The encoder network first predicts instantaneous frequency and the fractional index for . We then use the instantaneous phase to interpolate the waveform as:

where . The bilinear interpolation maps and as the row/column indexes to the matrix and returns the interpolated value using the four nearest neighbours. The wavetables are fixed in our case, contrary to Shan et al. (2022), and we only pick one wavetable at a time, not a weighted sum. We normalise each row of to have unit energy and align their negative peaks to the same position.

3.2 Frame‑Wise LPC synthesis

Time‑varying LPC synthesis is usually done by linearly interpolating the LPC coefficients to the audio resolution and filtering sample by sample. This is not parallelisable and slows down the training process. As an alternative, we approximate LPC synthesis by treating each frame independently and using overlap‑add:

where are the filter coefficients at the frame, and are the windowing functions, is the gain and is the hop size. is fixed to the square window. In this way, the computation can be parallelised on the frame level.

3.3 LPC coefficients parametrisation

For the LPC filter to be stable, all its poles must lie inside the unit circle on the complex plane. Stability can be guaranteed using robust representations, such as reflection coefficients (Subramani et al., 2022). The representation we choose in this work is cascaded ‑order IIR filters, and we solve the stability issue by ensuring all the ‑order filters are stable. We use the coefficient representation from Nercessian et al. (2021) to parametrise the IIR filter’s coefficients from the encoder’s outputs and cascade filters together to form an ‑order LPC filter . In practice, we multiply the filter polynomials first and perform one function call of the ‑order LPC filter, as it runs faster than sequentially applying each filter section.

3.4 Unvoiced gating

The encoder’s prediction of the instantaneous frequency is always non‑zero, keeping the oscillator working. Without constraint, the model would utilise these harmonics in the unvoiced region, creating buzzing artefacts (Wu et al., 2022). We propose to mitigate this problem by jointly training the model to predict the voiced/unvoiced probabilities as and feeding the gated frequency to the oscillator instead.

4.1 Gradients to the coefficients

Taking the derivatives of Equation 2 with respect to , we get:

does not depend on , so the initial conditions . is simply a matrix multiplication of . For , we use the fact that so we can get all the entries of using just . The details of computing are given in the Appendix.

4.2 Gradients to the input

To get the derivatives for input , let us first re‑write Equation 2 in the following convolutional form:

where . In this way, there is no output on the right‑hand side, so we can derive the gradients more easily. We assume , so we only need to evaluate in Equation 7. This assumption also implies that is always true, an assumption we stated in Section 2.2. is not practically computable, so we only use it for demonstration. The derivative of loss with respect to depends on future samples , which is:

By swapping the variables and and considering the equivalence of Equations 2 and 7, Equation 8 can be simplified to:

where and denotes function composition. Equation 9 shows that we can get the derivatives by just filtering with the same filter, but running in the reverse direction. The initial conditions are zeros, as the loss does not depend on them.

In conclusion, backpropagation through an IIR filter consists of two passes of the same filter and one matrix multiplication. The computation of the IIR does not need to be implemented in an automatic differentiation framework since we view an IIR as a single differentiable function, not a batch of differentiable multiplications and additions, avoiding the computational overhead issue we discussed in Section 2.2, and the filter can be highly optimised.

Interestingly, since the differentiable IIR calls itself for backpropagation, the backpropagation computational graph is differentiable, enabling higher‑order gradients. However, this goes beyond the paper’s scope but is worth exploring in the future.

5.1 Dataset

We test GOLF as a neural vocoder on the MPop600 dataset (Chu et al., 2020), a high‑quality Mandarin singing dataset featuring nearly 600 singing recordings with aligned lyrics sung by four singers. We use the audio recordings from the f1 (female) and m1 (male) singers. For each singer, we select the first three recordings as the test set, the following 27 as the validation set, and the rest as training data (around 3 hours). All the recordings are downsampled to 24 kHz, following the same procedure as Wu et al. (2022). The vocoder feature is the log mel spectrogram. We compute the feature with a window size of 1024 and 80 mel‑frequency bins and set the hop size to 120 (equal to 200 Hz frame rate). We normalise the feature to fall between zero and one and slice the training data into 2‑second excerpts with 1.5 seconds of overlap.

5.2 Model details

We adapt the encoder from SawSing with some modifications. The architecture details are shown in Figure 2. The number of parameters is 0.7 M. The synthesis parameters are predicted at the same rate as the mel spectrogram (200 Hz). The wavetable fractional index is predicted at 20 Hz to prevent modulation artefacts from rapid wavetable switching. We set and for dense sampling and accurate interpolation. For a sample rate of 24 kHz, the wavetables faithfully represent the glottal flow down to 12 Hz, which is much lower than the minimum f0 of the dataset (greater than 80 Hz). All the parameters except and are linearly upsampled to the audio rate. and are the Gaussian noise’s gain and filter coefficients. We use the Hanning window for and set for both LPC filters. We use 75% overlap for frame‑wise LPC.

Figure 2

We implemented the IIR in C++ with multi‑threading and CUDA for batch processing. The differentiability is achieved by registering the above backward computation and binding the optimised code in PyTorch. We contributed the implementation to TorchAudio (Hwang et al., 2023) in 2021 as part of the function torchaudio.functional.lfilter.² The initial implementation was developed independently in parallel with dynoNet, a similar tool by Forgione and Piga (2021). Our implementation provides a naive CUDA kernel, which enables it to be run on NVIDIA GPUs (the most common hardware for deep learning) for convenience. In comparison, dynoNet uses SciPy’s lfilter function as the backend and does not support its running on GPU.

We compare GOLF with three DDSP‑based baselines using the same encoder (without the downsampler for ) to predict their synthesis parameters. The size of the last linear layer is adjusted to match the number of parameters. The first two are the original DDSP (Engel et al., 2020) and SawSing (Wu et al., 2022). We set their noise filter length to 80 and harmonic filter length to 256 for SawSing. The third model is PUlse‑train LPC Filter (PULF), which is similar to GOLF but replaces the glottal flow wavetables with a band‑limited pulse train (Nercessian, 2023) using additive synthesis, while for the harmonic source is increased to 26 to accommodate the glottal pulse response. We set the number of oscillating sinusoids to be more than 150 for all baselines. We do not compare GOLF with Nercessian (2023) and Yoshimura et al. (2023b), because their architectures are based on the original source‑filter model and use additional post‑nets to enhance the voice, which makes direct comparison difficult.

5.3 Training configurations

We train separate models for each singer, resulting in eight models. The loss function is the summation of the same multi‑resolution short‑time Fourier transform (STFT) loss (MSSTFT) and f0 loss used in SawSing with the fast Fourier transform (FFT) sizes set to , plus a binary cross‑entropy loss on voiced/unvoiced prediction. Our MSSTFT follows the original DDSP (Engel et al., 2020) implementation, where the L1 distance is calculated between the magnitude of the ground truth and the predicted STFTs both in linear and log scale. Similar to Wu et al. (2022), we stop the gradients from the harmonic source to the f0s and voiced decisions to stabilise the training.

Although this paper uses reconstruction losses without adversarial training, which might result in an overly smooth spectral response, we empirically found this effect is not significant. The reason is most likely that, in a singer‑dependent setting, the compressed information of the mel spectrogram is enough to capture the singer’s voice characteristics. We verified this by training a speaker‑independent GOLF on a multi‑speaker dataset using almost the same configurations and found that the smoothing effect became more noticeable.

We used Adam (Kingma and Ba, 2015) for all optimisations. For DDSP and SawSing, the batch size and learning rate are set to 32 and 0.0005; for GOLF and PULF, the numbers were set to 64 and 0.0001 to reduce training instability. We used the ground truth f0s [extracted by WORLD (Morise et al., 2016)] for the harmonic oscillator of PULF during training due to stability issues. We trained all models for 800,000 steps to reach sufficient convergence (based on several early‑stage test runs). We picked the checkpoint with the lowest validation loss as the final model.

6.1 Objective evaluation

The objective metrics we choose are MSSTFT, mean absolute error (MAE) in f0 and Fréchet audio distance (FAD; Kilgour et al., 2019) on the predicted singing of the test set.³ Table 1 shows that DDSP has the lowest MSSTFT and f0 errors, while SawSing reaches the lowest FAD. GOLF and PULF show comparable results in f0 errors to other baselines. We report the memory usage measured on an Ubuntu 20.04 LTS machine with an i5‑4790k processor and an NVIDIA GeForce RTX 3070 GPU when training these models and their real‑time factor (RTF) in Table 2. GOLF requires approximately 35% of the memory needed by the other models for training, and demonstrates significantly faster performance, particularly on a central processing unit (CPU).

Given the evaluation results and the fact that GOLF’s number of synthesis parameters (49) is roughly four to six times less than DDSP (233) and SawSing (338), it is clear that GOLF’s synthesis parameters are more compact representations of the singing voice. Comparing the differences between GOLF and PULF in Table 2, we can see that wavetable synthesis is the factor yielding significant speedup. Other baselines synthesise band‑limited harmonic sources with many sinusoids oscillating simultaneously, thus increasing the computational cost.

6.2 Subjective evaluation

We conducted an online multiple stimuli with hidden reference and anchor (MUSHRA; Series, 2014) listening test using Go Listen (Barry et al., 2021). We picked one short clip from each test set recording, resulting in six clips with durations ranging from 6 to 11 seconds. Thus, the test is divided into six examples, each consisting of one ground truth clip and four synthesised results from different models. Their loudness was normalised to −16 dB LUFS using pyloudnorm (Steinmetz and Reiss, 2021). The order of the examples and the stimulus were randomised for each subject. Each subject was asked to rate these stimuli’s quality using a score ranging from 0 to 100. We collected responses from 33 anonymous participants. We excluded one response, as the participant did not indicate the use of headphones. We normalised scores from one to five and report the mean score in Figure 3. DDSP has the highest scores overall, and a Wilcoxon signed‑rank test shows that it is not significantly different from the f1 ground truth (p = 0.168). We apply the one‑sided Wilcoxon test on GOLF and PULF to compare them with SawSing. The results show that GOLF outperforms SawSing (p < 0.0001), and PULF performs better than SawSing on m1 (p < 0.022).

Figure 3

PULF’s MUSHRA score is much worse on f1, with noticeable artefacts in the unvoiced region and random components of the voice. After investigation, we found the noise gains predicted by PULF fluctuating at high speed, producing a modulation effect. This behaviour is also found in PULF trained on m1 and GOLF, even on the harmonic gains . Still, the degree of fluctuation is small and barely noticeable in the test samples. This amplitude modulation effect cannot be observed in MSSTFT spectrograms; thus, our training loss does not capture it. We will discuss the possible reason for this phenomenon in Section 7.5.

In addition, SawSing produced low scores for both singers because of the buzzing artefacts in the unvoiced region. This was pointed out in the original paper (Wu et al., 2022). Although unvoiced gating (Section 3.4) reduces this problem to a large degree, human ears are susceptible to this effect. This could be an inherent problem when using a sawtooth as the harmonic source.

7.1 Waveform‑Matching experiment

GOLF and PULF are the only two models that introduce non‑linear phase response because of IIR filtering. We are interested in whether their synthesised waveform matches the ground truth better than the other models. The filters in DDSP and SawSing are all zero‑phase, and the initial phases of the sinusoidal oscillators are fixed to zero. We hypothesise that, if a model can more faithfully represent the human voice mechanism, its waveform should be more similar to the reference recordings, even if its training objective does not consider phase information. We emphasise that this analysis does not target human perception but the phase reconstruction ability of the models. Humans cannot perceive the absolute frequency phase, but accurate reconstruction could be important in sound matching and mixing use cases.

We perform the test on one of the samples from m1 that we use in the subjective evaluation. We created a new parameter called phase offset sampled at 20 Hz. We linearly upsampled phase offset and added it to the instantaneous phase before feeding it to the oscillator, introducing a slowly varying phase shift. We optimised this parameter by minimising the predicted waveform’s L2 loss to the ground truth using Adam with a 0.001 learning rate and 1000 steps. We wrapped the differences between the points of phase offset during optimisation to [−0.5, 0.5]. We ran this optimisation five times for each model. Each time, the phase offset was initialised randomly. We report the minimum and maximum final losses from these trials. Table 3 shows that the lowest losses GOLF and PULF can reach are much smaller than the others, with GOLF having the smallest of all.

The L2 loss shown in Table 3 demonstrates that GOLF matches phase‑related characteristics more accurately than other models. Figure 4 shows that GOLF produces the most similar waveform to the ground truth. Other baselines’ waveforms are similar because they use the same additive synthesiser. Optimising the oscillators’ initial phases can reduce their L2 loss, but this does not account for time‑varying wave shapes.

Figure 4

7.2 Spectral analysis

To further understand the synthesised components, we compute the harmonics () and noise () spectrograms of each vocoder on the test set and take the average over the frames. The results are shown in Figure 5.

Figure 5

The first difference we notice is that, for GOLF and PULF, there is a substantial cut in the frequencies from 6 kHz to 10 kHz (except PULF on m1). The possible reason is that the voice usually has less energy in these high frequencies. As we factorise GOLF’s LPC filters into 11 ‑order filters (in PULF, the number is 13), to utilise the limited number of filters better, the model might discard the low energy frequencies and focus on the more important ones. This claim is supported by the differences in harmonic spectral envelopes between GOLF and PULF on m1, as PULF has two more filter sections than GOLF and does not cut these frequencies. In addition, as males’ fundamental frequencies are lower than females’, the intervals between the harmonics are narrower and less sparse, giving more energy to learning the filters in this range. This reasoning is based on a similar finding in our recent work Yu and Fazekas (2024), where we found that, if the signal to the GOLF filter has a more even energy distribution across the frequencies, the filter will learn better.

Another observation is that the harmonics from both PULF and GOLF show strong energies in frequencies above 10 kHz. In contrast, their noise components in these frequencies are much weaker than the other two baselines. We are still unsure why the harmonics filter acts like the noise filter near the Nyquist frequency. A possible cause may be the energies of high‑order partials, which vary more rapidly than the fundamental frequency. In very high frequencies, the contours of these partials are noise‑like, so the LPC filters tend to model them as noise. This is confirmed after a closer look at the harmonics spectrogram of GOLF and PULF, where we cannot see a clear harmonic structure in the filtered high frequencies.

We hypothesise that this effect is not observed in DDSP and SawSing due to the parametrisation we used in GOLF and PULF. For example, DDSP has to take the varying f0 into account when predicting corresponding harmonic gains to boost a specific frequency range. On the contrary, GOLF only needs to predict fixed coefficients in one of the filter sections to achieve the same effect. The in SawSing is an FIR filter, so it can model both the peaks and valleys in the spectrum, while the LPC filters are all‑pole and can only model the peaks (resonances). This results in high peak gains in the high frequencies we are discussing. The windowing in our frame‑wise LPC could also encourage this effect, as it smears the spectral peaks and makes them flatter. We will analyse the filter sections of the LPC filters in Section 7.4 to provide more insight.

7.3 Separability of components

The harmonics of SawSing are flatter and have more energy in high frequencies than the other models, as shown in Figure 5. After inspecting their spectrograms, we also found that SawSing harmonics tend to contribute more consonant noises, indicating more overlap between the harmonics and noise components. A more apparent separation between the harmonics and noise components could help analyse and manipulate the voice signal more effectively. To examine this, we compute the cosine distance between the harmonics and noise spectrograms of each model on the test set as:

where and are the harmonics and noise magnitude spectrograms, respectively. We only consider the frequency bins below 6 kHz to exclude the effects we discussed in Section 7.2.

The results are shown in Table 4. We can see that SawSing has the lowest cosine distance, indicating the least separability between the harmonics and noise. GOLF and PULF separate the components better than DDSP and SawSing, suggesting that the LPC filters are more suitable in the harmonic‑plus‑noise voice model.

7.4 Filter sections statistics

As we factorise the LPC filter into cascaded ‑order filters (Section 3.3), we can decompose the spectra further and analyse the behaviour of each filter section. We factorise every ‑order filter section into two poles . If the poles form a complex conjugate pair (), then the filter is a peak filter. For real poles, this equals cascading two first‑order high/low pass filters (Smith, 2024).

The number of occurrences of real‑pole filters can be found in Figure 6. We see that, most of the time, the poles are complex conjugate pairs (peak filters). The noise filters have more real poles than the harmonic filters. This is reasonable, as the harmonics contain more resonant structures such as formants, while the random components are more broadband. The only exception is PULF for the m1 singer, having two filters with significantly more real poles. This phenomenon may be attributed to the flat spectrum of pulse trains, which necessitates the model’s use of real poles (high/low pass filters) to adjust the spectral tilt.

Figure 6

Since GOLF and PULF filter sections most often predict conjugate poles, the following analysis will focus only on peak filters. We plot the statistics of (1) the resonant frequencies and (2) the resonant gains of the peak filters in Figure 7. The resonant frequency is computed as . Although the actual peak frequency slightly deviates from this due to the existence of the conjugate pole, this is a good approximation (Smith, 6 30). The resonant gain is .

Figure 7

From Figure 7, we see that around two to three extremely high peak filters are used above 10 kHz for the harmonic filters, which aligns with the high energy frequencies we discussed in Section 7.2. Both models learnt to assign distinct filters to formants, and their frequency statistics are aligned with each other, at least for (the one below 1 kHz) and (the one between 1 kHz and 2 kHz). We see, in the male singer m1, there are more aligned filters, possibly belonging to (around 3 kHz and above 20 dB) and (around 4 kHz). This is likely because the intervals between the harmonic partials are narrower for the male voice, so the energy distribution is denser, making it easier to learn the formants, which is the same reason we discussed in Section 7.2. Due to broadband noise, the noise filter statistics between the two models tend to be more similar than those of the harmonic filters.

7.5 Interleaving behaviour of resonant frequencies

Although the time‑varying behaviour of the resonant frequencies is not constrained, we found that the learnt formants generally correspond to the filters one‑to‑one. In other words, the models learnt to implicitly assign the same filter to the same formant across the frames, which could be used for formant tracking. We believe this is because the formants have a high signal‑to‑noise ratio and distinct structures, which make it easier for the model to learn.

Other filter sections do not always have this one‑to‑one mapping and sometimes have different resonant frequency trajectories interleaved inside one filter section. To demonstrate this, we pick three filter sections from the harmonic filters of GOLF on the f1 singer and plot the trajectories of the poles in the complex plane in Figure 8. Filter 1 has three separated groups of poles, and in this short 20‑frame slice, it jumps between them. This effect is not observable in the filtered spectrogram with only filter 1, as the overlap‑add of frame‑wise LPC smooths these abrupt changes. We suspect this is the same cause resulting in fluctuating noise gains discussed in Section 8. The fluctuation pattern we see is the gain oscillates with a cycle of two frames, which can also be seen as jumping between two groups of gain values. The MSSTFT loss cannot capture this change as the smallest FFT size in our loss is 512, which can cover at least five synthesis frames (120 4 + 1 = 481), and thus is large enough to smooth this effect in the spectrogram. We test this idea using an FFT size of 128 to calculate the spectrograms of PULF’s noise component. The resulting spectrograms show a perceivable periodic pattern, thus verifying our hypothesis.

Figure 8

7.6 Statistics of

To analyse how GOLF uses to synthesise the glottal flow, we compute the minimum, maximum and average values of [inversely mapped from ] on the frames when and plot the synthesised glottal flows using these statistics in Figure 9.

Figure 9

From the plots, we see that there is no significant difference between the two singers. The models learn to use higher values (which result in a fast‑decaying spectrum) more and barely use an below 1, which corresponds to a flatter spectrum. This is expected, as our training objective is a spectral loss, so GOLF likely learns a response that matches the shape of the spectrum on average (see the similar decaying structure of in Figure 5). If the time‑domain loss is used, the model might learn a very different distribution of , as the target becomes learning the most similar waveform instead of its magnitude response.