Implementation of a Variable Fractional Delay with Lagrange Interpolation using Farrow Structure

Question

I am writing a C++ simulation software working in time domain. I generate regularly sampled data, and need them to be delayed, in "real-time", by a variable fractional delay.

This is a pretty common problem, and very close to all the "resampling" issues. After a short bibliographic search, it seems that the most commonly used method is Lagrange Interpolation with a Farrow structure. Based on J. O. Smith, "Interpolated Delay Lines, Ideal Bandlimited Interpolation, and Fractional Delay Filter Design", I used the following implementation:

which includes cascading filters with fixed coefficients $1+z^{-1}$, modulated by cascading factors that depend on the variable delay. This can be seen in the equation

$$ \hat{H}_\Delta(q^{-1}) = \sum_{n=0}^N \hat{h}_\Delta[n](q-1)^n. $$

My implementation in C++ works well for interpolation orders up to 51. Unfortunately, for larger interpolation orders the numerical errors increase rapidly and become dominant. After a quick inspection, it seems that the buffered samples, in the $1+z^{-1}$ cells of increasing index, have an increasing order of magnitude. It helps when I use long double, but this becomes insufficient when the order is of the order of 60: we loose a lot of numerical precision when we combine floating-point values of around $10^{10}$ with values of the order of $1$.

Do you know a variation of the Farrow structure or another implementation for Lagrange interpolation that works well for large orders, i.e. that remains computationally efficient while requiring only reasonable numerical precision in the memory buffers?

Thank you in advance for your time :)

Answer 1

The task at hand is a pretty common one; most symbol timing synchronizers need to do exactly what you describe.

Same goes for many types of resamplers – they need to reconstruct varying fractional inter-sample "snapshots" from the input sample stream to construct a non-rationally related rate output.

So that's where I'd start looking.

In fact, the farrow structure idea takes you pretty far!

Imagine you could have a large bank of linear-phase filters. A linear-phase filter has a constant group delay.

By modulating said filters, you can change the phase slope and hence the group delay. Now, you just pick the filter that fits your needs.

In a polyphase filterbank, such as implemented by the GNU Radio "PFB Arbitrary Resampler" block, exactly that is done: You construct a whole bunch of filters, and for each sample, pick the one that's closest to what you need (or the two closest, and linearly interpolate between these).
Because it's used as a resampler, for every non-1.0 resampling ratio, the next sample will be choosen from a different filter branch, with the appropriate delay.

Now, since you, most of the time, just want to delay, you seldom need the "resampling" capability of that resampler. (Resampling is exactly what happens when you change the delay linearly)

So, instead, you'd have a single filter (if possible, use a filter that you need to have, anyway), and you'd modulate it.

Think about this in frequency domain: assume you had a linear, continuous Fourier transform that you apply to your signal and you multiply that with both the FT of your filter impulse response and $e^{j2\pi\tau f}$ before transforming back to time domain. To adjust delay, you adjust $\tau$. Tadaaah! Variable fractional delay.

Problem certainly is that you're in DSP, so you don't get that linear Fourier transform to implement linear convolution; you get the DFT, but multiplication in discrete frequency domain is equivalent to circular convolution in time domain.

But that's been solved already – fast convolution allows you to apply long filters with little effort (and low numerical problems – the FFT is exceptionally stable, usually even more stable than the direct convolution!). Now, if you go ahead and implement that e.g. via the overlap-save method, all you need to do is make sure that the phase of the complex exponential doesn't "jump" between FFT segments. You can then change (even gracefully!) the $\tau$ between individual convolutions.