Faust Code for Lagrange Interpolation
The Faust programming language for signal processing [453,450] 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.
// 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);
|
Next Section:
Lagrange Frequency Response Examples
Previous Section:
Maxima Code for Lagrange Interpolation
