Tabulate an nD function for access via nearest-value or linear or cubic interpolation.
In other words, the tabulated function can be evaluated using interpolation of order 0 (none), 1 (linear), or 3 (cubic).
The table size and parameter range of each dimension can and must be separately specified.
You can use it anywhere you have an expensive function with multiple parameters with known ranges.
You could use it to build a wavetable synth, for example.
The number of dimensions is deduced from the number of parameters you give, see below.
tabulateNd(C, function, (parameters) ).(val|lin|cub) : _
C: whether to dynamically force the parameter values for each dimension to the ranges specified in parameters: 1 forces the check, 0 deactivates it (constant numerical expression)function: the function we want to tabulate. Can have any number of inputs, but needs to have just one output.(parameters): sizes, ranges and read values.
NOTE: These need to be in brackets, to make them one entity.
If N is the number of dimensions, we need:- N times
S: number of values to store for this dimension (constant numerical expression) - N times
r0: minimum value of this dimension - N times
r1: maximum value of this dimension - N times
x: read value of this dimension
- N times
By providing these parameters, you indirectly specify the number of dimensions; it's the number of parameters divided by 4.
The user facing functions are:
tabulateNd(C, function, S, parameters).val
- Uses the value in the table closest to x.
tabulateNd(C, function, S, parameters).lin
- Evaluates at x using linear interpolation between the closest stored values.
tabulateNd(C, function, S, parameters).cub
- Evaluates at x using cubic interpolation between the closest stored values.
powSin(x,y) = sin(pow(x,y)); // The function we want to tabulate
powSinTable(x,y) = ba.tabulateNd(1, powSin, (sizeX,sizeY, rx0,ry0, rx1,ry1, x,y) ).lin;
sizeX = 512; // table size of the first parameter
sizeY = 512; // of the size of the second parameter
rx0 = 0; // start of the range of the first parameter
ry0 = 0; // start of the range of the second parameter
rx1 = 10; // end of the range of the first parameter
ry1 = 10; // end of the range of the second parameter
x= hslider("x", rx0, rx0, rx1, 0.001):si.smoo;
y= hslider("y", ry0, ry0, ry1, 0.001):si.smoo;
process = powSinTable(x,y), powSin(y,y);
The .val function just outputs the closest stored value.
The .lin and .cub functions interpolate in N dimensions.
To understand what it means to interpolate in N dimensions, here's a quick reminder on the general principle of 2D linear interpolation:
We have a grid of values, and we want to find the value at a point (x, y) within this grid.
We first find the four closest points (A, B, C, D) in the grid surrounding the point (x, y).
Then, we perform linear interpolation in the x-direction between points A and B, and between points C and D.
This gives us two new points E and F.
Finally, we perform linear interpolation in the y-direction between points E and F to get our value.
To implement this in Faust, we need N sequential groups of interpolators, where N is the number of dimensions.
Each group feeds into the next, with the last "group" being a single interpolator, and the group before it containing one interpolator for each input of the group it's feeding.
Some examples:
- Our 2D linear example has two interpolators feeding into one.
- A 3D linear interpolator has four interpolators feeding into two, feeding into one.
- A 2D cubic interpolater has four interpolators feeding into one.
- A 3D cubic interpolator has sixteen interpolators feeding into four, feeding into one.
To understand which values we need to look up, let's consider the 2D linear example again.
The four values going into the first group represent the four closest points (A, B, C, D) mentioned above.
- The first interpolator gets:
- The closest value that is stored (A)
- The next value in the x dimension, keeping y fixed (B)
- The second interpolator gets:
- One step over in the y dimension, keeping x fixed (C)
- One step over in both the x dimension and the y dimension (D)
The outputs of these two interpolators are points E and F.
In other words: the interpolated x values and, respectively, the following y values:
- The closest stored value of the y dimension
- One step forward in the y dimension
The last interpolator takes these two values and interpolates them in the y dimension.
To generalize for N dimensions, linear interpolation:
The first group has 2^(n-1) parallel interpolators interpolating in the first dimension.
The second group has 2^(n-2) parallel interpolators interpolating in the second dimension.
The process continues until the n-th group, which has a single interpolator interpolating in the n-th dimension.
The same principle applies to the cubic interpolation in nD. The only difference is that there would be 4^(n-1) parallel interpolators in the first group, compared to 2^(n-1) for linear interpolation.
This is what the mixers function does.
Besides the values, each interpolator also needs to know the weight of each value in it's output.
Let's call this d, like in ba.interpolate.
It is the same for each group of interpolators, since it correlates to a dimension.
It's value is calculated the similarly to ba.interpolate:
First we prepare a 'float' table read-index for that dimension (id in ba.tabulate):
If the table only had that dimension and it could read a float index, what would it be.
Then we int the float index to get the value we have stored that is closest to, but lower than the input value; the actual index for that dimension.
Our d is the difference between the float index and the actual index.
The ids function calculates the id for each dimension and inside the mixer function they get turned into ds.
The elephant in the room is: how do we get these indexes?
For that we need to know how the values are stored.
We use one big table to store everything.
To understand the concept, let's look at the 2D example again, and then we'll extend it to 3d and the general nD case.
Let's say we have a 2D table with dimensions A and B where:
A has 3 values between 0 and 5 and B has 4 values between 0 and 1.
The 1D array representation of this 2D table will have a size of 3 * 4 = 12.
The values are stored in the following way:
- First 3 values: A is 0, then 3, then 5 while B is at 0.
- Next 3 values: A changes from 0 to 5 while B is at 1/3.
- Next 3 values: A changes from 0 to 5 while B is at 2/3.
- Last 3 values: A changes from 0 to 5 while B is at 1.
For the 3D example, let's extend the 2D example with an additional dimension C having 2 values between 0 and 2.
The total size will be 3 * 4 * 2 = 24.
The values are stored like so:
- First 3 values: A changes from 0 to 5, B is at 0, and C is at 0.
- Next 3 values: A changes from 0 to 5, B is at 1/3, and C is at 0.
- Next 3 values: A changes from 0 to 5, B is at 2/3, and C is at 0.
- Next 3 values: A changes from 0 to 5, B is at 1, and C is at 0.
The last 12 values are the same as the first 12, but with C at 2.
For the general n-dimensional case, we iterate through all dimensions, changing the values of the innermost dimension first, then moving towards the outer dimensions.
To get the float read index (id) corresponding to a particular dimension, we scale the function input value to be between 0 and 1, and multiply it by the size of that dimension minus one.
To understand how we get the readIndexfor .val, let's work trough how we'd do it in our 2D linear example.
For simplicity's sake, the ranges of the inputs to our function are both 0 to 1.
Say we wanted to read the value closest to x=0.5 and y=0, so the id of x is 1 (the second value) and the id of y is 0 (first value).
In this case, the read index is just the id of x, rounded to the nearest integer, just like in ba.tabulate.
If we want to read the value belonging to x=0.5 and y=2/3, things get more complicated.
The id for y is now 2, the third value.
For each step in the y direction, we need to increase the index by 3, the number of values that are stored for x.
So the influence of the y is: the size of x times the rounded id of y.
The final read index is the rounded id of x plus the influence of y.
For the general nD case, we take need to do the same operation N times, each feeding into the next.
This operation is the riN function.
We take four parameters: the size of the dimension before it prevSize, the index of the previous dimension prevIX, the current size sizeX and the current id idX.
riN has 2 outputs, the size, for feeding into the next dimension's prevSize, and the read index feeding into the next dimension's prevIX.
The size is the sizeX times prevSize.
The read index is the rounded idX times prevSize added to the prevIX.
Our final readIndex is the read index output of the last dimension.
The rdtables for the interpolators need a pattern of offsets in each dimension, since we are looking for the read indexes surrounding the point of interest.
These offsets are best explained by looking at the code of tabulate2d, the hardcoded 2D version:
tabulate2d(C,function,sizeX,sizeY, rx0, ry0, rx1, ry1,x,y) =
environment {
size = sizeX*sizeY;
// Maximum X index to access
midX = sizeX-1;
// Maximum Y index to access
midY = sizeY-1;
// Maximum total index to access
mid = size-1;
// Create the table
wf = function(wfX,wfY);
// Prepare the 'float' table read index for X
idX = (x-rx0)/(rx1-rx0)*midX;
// Prepare the 'float' table read index for Y
idY = ((y-ry0)/(ry1-ry0))*midY;
// table creation X:
wfX =
rx0+float(ba.time%sizeX)*(rx1-rx0)
/float(midX);
// table creation Y:
wfY =
ry0+
((float(ba.time-(ba.time%sizeX))
/float(sizeX))
*(ry1-ry0)
)
/float(midY)
;
// Limit the table read index in [0, mid] if C = 1
rid(x,mid, 0) = x;
rid(x,mid, 1) = max(0, min(x, mid));
// Tabulate an unary 'FX' function on a range [rx0, rx1] [ry0, ry1]
val(x,y) =
rdtable(size, wf, readIndex);
readIndex =
rid(
rid(int(idX+0.5),midX, C)
+yOffset
, mid, C);
yOffset =
sizeX*rid(int(idY),midY,C);
// Tabulate an unary 'FX' function over the range [rx0, rx1] [ry0, ry1] with linear interpolation
lin =
it.interpolate_linear(
dy
, it.interpolate_linear(dx,v0,v1)
, it.interpolate_linear(dx,v2,v3))
with {
i0 = rid(int(idX), midX, C)+yOffset;
i1 = i0+1;
i2 = i0+sizeX;
i3 = i1+sizeX;
dx = idX-int(idX);
dy = idY-int(idY);
v0 = rdtable(size, wf, rid(i0, mid, C));
v1 = rdtable(size, wf, rid(i1, mid, C));
v2 = rdtable(size, wf, rid(i2, mid, C));
v3 = rdtable(size, wf, rid(i3, mid, C));
};
// Tabulate an unary 'FX' function over the range [rx0, rx1] [ry0, ry1] with cubic interpolation
cub =
it.interpolate_cubic(
dy
, it.interpolate_cubic(dx,v0,v1,v2,v3)
, it.interpolate_cubic(dx,v4,v5,v6,v7)
, it.interpolate_cubic(dx,v8,v9,v10,v11)
, it.interpolate_cubic(dx,v12,v13,v14,v15)
)
with {
i0 = i4-sizeX;
i1 = i5-sizeX;
i2 = i6-sizeX;
i3 = i7-sizeX;
i4 = i5-1;
i5 = rid(int(idX), midX, C)+yOffset;
i6 = i5+1;
i7 = i6+1;
i8 = i4+sizeX;
i9 = i5+sizeX;
i10 = i6+sizeX;
i11 = i7+sizeX;
i12 = i4+(2*sizeX);
i13 = i5+(2*sizeX);
i14 = i6+(2*sizeX);
i15 = i7+(2*sizeX);
dx = idX-int(idX);
dy = idY-int(idY);
v0 = rdtable(size, wf, rid(i0 , mid, C));
v1 = rdtable(size, wf, rid(i1 , mid, C));
v2 = rdtable(size, wf, rid(i2 , mid, C));
v3 = rdtable(size, wf, rid(i3 , mid, C));
v4 = rdtable(size, wf, rid(i4 , mid, C));
v5 = rdtable(size, wf, rid(i5 , mid, C));
v6 = rdtable(size, wf, rid(i6 , mid, C));
v7 = rdtable(size, wf, rid(i7 , mid, C));
v8 = rdtable(size, wf, rid(i8 , mid, C));
v9 = rdtable(size, wf, rid(i9 , mid, C));
v10 = rdtable(size, wf, rid(i10, mid, C));
v11 = rdtable(size, wf, rid(i11, mid, C));
v12 = rdtable(size, wf, rid(i12, mid, C));
v13 = rdtable(size, wf, rid(i13, mid, C));
v14 = rdtable(size, wf, rid(i14, mid, C));
v15 = rdtable(size, wf, rid(i15, mid, C));
};
};In the intrest of brevity, we'll stop explaining here. If you have any more questions, feel free to open an issue on https://github.com/grame-cncm/faustlibraries/ and tag @magnetophon.