Tuning by Lagrange Interpolation

Researchers at the Helsinki University of Technology have historically used Lagrange interpolation for digital-waveguide fine-tuning [5,18,6]. This has the advantage of being robust under rapidly time-varying conditions, as opposed to allpass interpolation which can exhibit artifacts when the delay changes too rapidly. However, Lagrange interpolation, like all FIR filters except pure delays, has the disadvantage of introducing some gain error in the string feedback loop, unlike allpass interpolation.

As discussed further in [14], the closed-form expression for Lagrange-interpolation coefficients is

$$h_\Delta(n) = \prod_{\stackrel{k=0}{k\neq n}}^N \frac{\Delta-k}{n-k}, \quad n=0,1,2,\ldots,N$$

where $N$ is the length of the interpolator and $\Delta$ is the desired delay, centered about $(N-1)/2$. Note that $N=1$ corresponds to linear interpolation. For audio delay-line interpolation, the fourth-order case $N=4$ is often sufficient (and sometimes overkill).

For the same reasons discussed in §3.6.1, fourth-order Lagrange interpolation of a delay line is best viewed as a five-tap interpolating read of the delay line at the desired fractional delay. In normal object-oriented languages, an interpolating read is naturally implemented internal to the interpolating delay-line object. In Faust, the same effect is obtained by specifying five delay lines, as shown in Fig.8. While it appears that five delay lines are needed in the Faust implementation, only one is actually used in the generated C++ code, thanks to compiler optimizations. Faust implementations of Lagrange interpolation of orders 1 through 4 may be found in the file filter.lib distributed with Faust (starting with version 0.9.9.3).

Figure 8: Faust implementation of fourth-order Lagrange interpolation.

 

fdelay4(n,d,x) = delay(n,id,x) * fdm1*fdm2*fdm3*fdm4/24 
               + delay(n,id+1,x) * (0-fd*fdm2*fdm3*fdm4)/6
               + delay(n,id+2,x) * fd*fdm1*fdm3*fdm4/4
               + delay(n,id+3,x) * (0-fd*fdm1*fdm2*fdm4)/6
               + delay(n,id+4,x) * fd*fdm1*fdm2*fdm3/24
with {
  id = int(d-1);
  fd = 1+frac(d);
  fdm1 = fd-1;
  fdm2 = fd-2;
  fdm3 = fd-3;
  fdm4 = fd-4;
};