Faust Code for Lagrange Interpolation

The Faust programming language for signal processing [455,453] includes support for Lagrange fractional-delay filtering, up to order five, in the library file filter.lib. For example, the fourth-order case is listed below:

// fourth-order (quartic) case, delay d in [1.5,2.5]
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 {
  o = 1.49999;
  dmo = d - o; // assumed nonnegative
  id = int(dmo);
  fd = o + frac(dmo);
  fdm1 = fd-1;
  fdm2 = fd-2;
  fdm3 = fd-3;
  fdm4 = fd-4;
};

An example calling program is shown in Fig.4.12.

Figure: Faust program tlagrange.dsp used to generate Figures 4.13 through 4.16.

 

// tlagrange.dsp - test Lagrange interpolation in Faust

import("filter.lib");

N = 16; % Allocated delay-line length

% Compare various orders:
D = 5.4;
process = 1-1' <: fdelay1(N,D),
                  fdelay2(N,D),
                  fdelay3(N,D),
                  fdelay4(N,D),
                  fdelay5(N,D);
// To see results:
// [in a shell]:
//   faust2octave tlagrange.dsp
// [at the Octave command prompt]:
//   plot(db(fft(faustout,1024)(1:512,:)));

// Alternate example for testing a range of 4th-order cases
// (change name to "process" and rename "process" above):
process2  = 1-1' <: fdelay4(N, 1.5),
                    fdelay4(N, 1.6),
                    fdelay4(N, 1.7),
                    fdelay4(N, 1.8),
                    fdelay4(N, 1.9),
                    fdelay4(N, 2.0),
                    fdelay4(N, 2.1),
                    fdelay4(N, 2.2),
                    fdelay4(N, 2.3),
                    fdelay4(N, 2.4),
                    fdelay4(N, 2.499),
                    fdelay4(N, 2.5);