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

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?

• hm, what kind of data and what kind of rates are we talking about? – Marcus Müller Oct 8 '18 at 13:27
• why do you want to do Lagrange interpolation instead of simply a windowed $\operatorname{sinc}(\cdot)$ interpolator? – robert bristow-johnson Oct 8 '18 at 17:16
• @robertbristow-johnson Allegedly because the Farrow structure is better numerically! It effectively uses Horner's method for calculating the output, which should be better than calculating it directly. Oh, and I used $\hbar$ because that's what the OP used. I'd prefer to stick with the notation they used. YMMV. – Peter K. Oct 8 '18 at 18:11
• well @PeterK. i think i understand why Horner's method of evaluating a polynomial is better, for evaluating polynomials. but i don't understand why limiting one to Lagrange is so good (when you can have a kaiser-windowed sinc()). and if you have a double-wide accumulator (like any decent FIR machine would have), i don't see any numerical error until the final cast to a single-width word for the output. – robert bristow-johnson Oct 8 '18 at 20:47
• no, you compute a table of coefficients in advance. i have found that 512 different phases or fractional delays is enough, you might get away with 256 phases. for each phase you will have something like 16 coefficients (depending on how solid your brick-wall LPF is). maybe 32 coefficients. or in between, but it need not be a power of 2. so, for an arbitrary delay, you compute the two discrete fractional delays closest to your arbitrary delay and linearly interpolate. so that is two FIRs with 16 to 32 terms (depending on how long the FIR is) and a linear interpolation. – robert bristow-johnson Oct 9 '18 at 21:12

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.

• I edited a bit my question, in the hope of making myself more clear: I implemented the Lagrange interpolator using that Farrow structure. And it works well, for "low" orders. However it breaks down for large orders because the numerical precision required in the filter memory buffers increase with the order... and becomes unreasonable for orders of the order of 61. Which I need. – Jean-Baptiste Oct 8 '18 at 16:06