Multi-dimensional Interpolation/Extrapolation ============================================= The purpose of this essay is to describe and explain the multi-dimensional interpolation and extrapolation that is described in the verilog-AMS LRM for the $table_model function. The description is done from the perspective of the algorithms and takes into consideration their requirements. Another perspective is that of the user providing the data. This perspective is not described, but it is assumed that it is possible to convert a user specification of multi-dimensional interpolation data to a representation suitable for the algorithms. One-dimensional Interpolation and Extrapolation ----------------------------------------------- The function to be interpolated is described by an ordered collection of (x,y) pairs. In this collection the x values are strictly monotonically increasing, i.e. each x value appears in the collection exactly once. In what follows we call such a collection an (x,y) table. A one-dimensional interpolation may use one of a number of algorithms, including at least the following: - constant value between two x values, i.e. f(x) = yi in the interval xi to xi+1 - linear interpolation - quadratic splines - cubic splines Similarly, several methods for extrapolation must be supported at either end of the range of the data, but in te end there is a finite number of combinations of interpolation and extrapolation that must be supported. One-dimensional interpolation/extrapolation is well understood and is not described further. Multi-dimensional Interpolation and Extrapolation ------------------------------------------------- The Verilog-AMS definition of the $table_model function is based on the concept of isolines. This concept is based on data creation by hierarchical parameter sweeps and leads to a recursive interpolation. In essence, while xi is held constant xi+1 is swept over a certain range and the values recorded. An alternative view, better suited to understand a multi-dimensional interpolation/extrapolation algorith is that of a hierarchy of (x,y) tables, as follows. - For the outermost dimension xn the data is described by the (x,y) table (xn1...xnm,yn1...ynm). The values xni are defined numerically. Conversely, each yni value is described by an (x,y) table at the next lower dimension. - Similarly, the data for an intermediate dimension xj is described by an (x,y) table (xj1...xjk,yj1...yjk), where the xji are defined numerically and the yji are defined by an (x,y) table at the next lower dimension. - The innermost dimension x1 is described by an (x,y) table where both the x and y values are given numerically. - Each (x,y) table may have a different number of entries. - Each (x,y) table is ordered such that the x values are strictly monotonic. The determination of a value f(x1...xn) at coordinates x1...xn now requires a one-dimensional interpolation for the (x,y) table (xn1...xnm,yn1...ynm). The values yn1...ynm are unknown, but a value yni can be determined by a one-dimensional interpolation on the (x,y) table of the next lower dimension. This continues recursively until we reach the innermost dimension, where the recursion ends because there the y values are known numerically. The process is best demonstrated with an example. Consider the 3-dimensional table that in an external file using the Verilog-AMS semantics is defined as follows: x3 x2 x1 f(x1,x2,x3) 1.0 0.0 1.0 0.5 1.0 0.0 2.0 1.0 1.0 0.0 3.0 1.5 1.0 0.0 4.0 2.0 1.0 0.0 5.0 2.5 1.0 0.0 6.0 3.0 1.0 0.5 1.0 1.0 1.0 0.5 3.0 2.0 1.0 0.5 5.0 3.0 1.0 1.0 1.0 1.5 1.0 1.0 2.0 2.0 1.0 1.0 4.0 3.0 2.0 0.0 1.0 1.0 2.0 0.0 3.0 2.0 2.0 0.0 5.0 4.0 2.0 1.0 1.0 2.0 2.0 1.0 2.0 3.0 2.0 1.0 4.0 5.0 2.0 1.0 6.0 6.0 2.0 1.0 8.0 7.0 The data describe one (x,y) table at the outermost level, two (x,y) tables at the next level (because the outermost (x,y) table has two entries), and five (x,y) tables at the innermost level (because the two tables an the middle level have 3 and 2 entries, respectively. Determine the value at x1=3.5, x2=0.25, x3=1.6. For simplicity, use linear interpolation/extrapolation. To do this we must find the numerical values of the outermost (x,y) table. For x3=1 we have an (x,y) table with x2 values 0.0, 0.5, 1.0, but the y values are not known. We descend one more level and find, in turn, the values f(3.5,0.0,1.0)=1.75 f(3.5,0.5,1.0)=2.25 f(3.5,1.0,1.0)=2.75 That is, we have built a one-dimensional interpolation (0.0,1.75),(0.5,2.25), (1.0,2.75). From this we get the value f(1.0,0.25,3.5)=2 Similarly, for x3=2.0: f(3.5,0.0,2.0)=2.5 f(3.5,1.0,2.0)=4.5 The one-dimensional table to interpolate in is (0.0,2.5),(1.0,4.5). From this the value f(2.0,0.25,3.5) is found to be 3.0. Now, we have determined the outermost table completely: (1.0,2.0),(2.0,3.0), and we find the value at (3.5,0.25,1.6) to be 2.6. This can, of course, be extended to an arbitrary number of dimensions, and it works with any one-dimensional interpolation/extrapolation scheme. Data Organization ----------------- In a general purpose programming language the best way to store the data would be by means of a collection of (x,y) tables where the x values are defined numerically and each y value is either defined numerically or by means of a pointer to the (x,y) table from which its value can be determined. Along with each (x,y) table comes the definition of how the table must be interpolated/extrapolated. For VHDL we do not have the option of using pointers (access types require variables, but here we deal with constant data). Also, tha fact that each (x,y) table can have a different size makes it impossible to keep the data as an array of (x,y) tables, unless we allow each (x,y) table to be as large as the largest and thereby allow to waste some space. To work around these issues we can use the following approach: - store the x values in a single vector. Alternatively, store the x values at the innermost dimension in one vector and the x values at the outer dimensions in another vector - store the y values of the innermost dimension in a single vector such that corresponding (x,y) pairs have the same index - create a descriptor for each (x,y) table that records: - the start and end index of the table in the x or y vectors - the dimension index: 1 for innermost, n for outermost dimension in an n-dimensional {x,y) table The descriptors can be stored in a vector whose size is the number of (x,y) tables. For the innermost dimension there is a matching y value for each x value. For each x value at a higher dimension the corresponding y value is determined by an interpolation in a table at the next lower dimension. For performance reasons the association of a table and an x value at a higher dimension should be computed only once and saved along with the table. We establish that there is a single (x,y) table at the outermost dimension. It is possible to compute an n-dimensional interpolation either recursively or iteratively. The data should be stored slightly differently for the two approaches. Recursive algorithm ''''''''''''''''''' The recursive algorithm determines a value at the n-tuple (x1,...xn) by starting at the n-th dimension. Each x value at a dimension > 1 has a corresponding pointer to the table required to determine the corrsponding y value. The pointer is implemented as an index into the vector of (x,y) table descriptors. The following declarations implement this approach: type tdm_table_descriptor is record start_index: integer; -- index of leftmost entry of table end_index: integer; -- index of rightmost entry of table dimension_index: integer; -- index into interpolation/extrapolation -- descriptors. Equals dimension end record; type tdm_table_descriptor_vector is array(natural range <>) of tdm_table_descriptor; type tdm_table_data is record x1: real_vector; -- x values at innermost dimension y1: real_vector; -- corresponding y values xh: real_vector; -- x values at higher dimensions yhi: integer vector;-- vector of indices of each x value to the table -- needed to determine the corresponding y value descr: tdm_table_descriptor_vector; -- table descriptors ip: tdm_interpolation_vector; -- interpolation type for each dimension epl: tdm_extrapolation_vector; -- extrapolation type at lower end -- for each dimension epu: tdm_extrapolation_vector; -- extrapolation type at upper end -- for each dimension end record; As an example, the table used above is stored as follows (the fors column is the array index for easier understanding, it is not stored in tdm_table_data): x1 y1 1 1.0 0.5 2 2.0 1.0 3 3.0 1.5 4 4.0 2.0 5 5.0 2.5 6 6.0 3.0 7 1.0 1.0 8 3.0 2.0 9 5.0 3.0 10 1.0 1.5 11 2.0 2.0 12 4.0 3.0 13 1.0 1.0 14 3.0 2.0 15 5.0 4.0 16 1.0 2.0 17 2.0 3.0 18 4.0 5.0 19 6.0 6.0 20 8.0 7.0 xh yhi 1 0.0 1 2 0.5 2 3 1.0 3 4 0.0 4 5 1.0 5 6 1.0 6 7 2.0 7 descr start_index end_index dimension_index 1 6 1 7 9 1 10 12 1 13 15 1 16 20 1 1 3 2 4 5 2 6 7 3 Interpolation starts at the highest index of descr. Because dimension_index > 1 we know that start_index,end_index apply to (xh,yhi). The table starts at 6 and ends at 7. To determine the corresponding yh: - yhi is an index into descr. Because yhi(6).dimension_index > 1 start_index and end_index describe a table at positions 1...3 in (xh,yhi) - yhi is an index into descr. Because yhi(1).dimension_index == 1 start_index and end_index describe a table at positions 1...6 in (x1,y1) - find y value in from this table, save it in a variable holding the yh values - proceed until all yh values have been determined and y=f(x1...xn) has been found. The example above uses this algorithm. Iterative algorithm ''''''''''''''''''' Because of the isolines approach it is possible to determine the interpolated values iteratively. Interpolation starts at the innermost dimension and ends when the y value at the outermost dimension has been found. The following declarations implement this approach: type tdm_table_descriptor is record start_index: integer; -- index of leftmost entry of table end_index: integer; -- index of rightmost entry of table dimension_index: integer; -- index into interpolation/extrapolation -- descriptors. Equals dimension save_index: integer; -- index into yh vector where to save value -- determined from this table end record; type tdm_table_descriptor_vector is array(natural range <>) of tdm_table_descriptor; type tdm_table_data is record x1: real_vector; -- x values at innermost dimension y1: real_vector; -- corresponding y values xh: real_vector; -- x values at higher dimensions descr: tdm_table_descriptor_vector; -- table descriptors ip: tdm_interpolation_vector; -- interpolation type for each dimension epl: tdm_extrapolation_vector; -- extrapolation type at lower end -- for each dimension epu: tdm_extrapolation_vector; -- extrapolation type at upper end -- for each dimension end record; As an example, the table used above is stored as follows (the fors column is the array index for easier understanding, it is not stored in tdm_table_data): x1 y1 1 1.0 0.5 2 2.0 1.0 3 3.0 1.5 4 4.0 2.0 5 5.0 2.5 6 6.0 3.0 7 1.0 1.0 8 3.0 2.0 9 5.0 3.0 10 1.0 1.5 11 2.0 2.0 12 4.0 3.0 13 1.0 1.0 14 3.0 2.0 15 5.0 4.0 16 1.0 2.0 17 2.0 3.0 18 4.0 5.0 19 6.0 6.0 20 8.0 7.0 xh 1 0.0 2 0.5 3 1.0 4 0.0 5 1.0 6 1.0 7 2.0 descr start_index end_index dimension_index save_index 1 6 1 1 7 9 1 2 10 12 1 3 13 15 1 4 16 20 1 5 1 3 2 6 4 5 2 7 6 7 3 0 Interpolation starts at index 1 of descr and ends at its highest index. This is a symbolic description of the algorithm. Note that the run time function must declare a vector yh whose range is the range of xh. for i from descr'low to descr'high loop if descr(i).dimension_index==1 then interpolate in table x1(descr(i).start_index)..x1(descr(i).end_index) and corresponding y1 values save result at yh(descr(i).save_index else interpolate in table xh(descr(i).start_index)..x1(descr(i).end_index) and corresponding yh values. The yh values for a table have been determined by the time they are needed if descr(i).save_index == 0 then we are done, return value else save result at yh(descr(i).save_index end if end if end loop In the example we interpolate all tables with dimension_index 1 at 3.5, those with dimension_index 2 at 0.25, and the outermost table with dimension_index 3 at 1.6. We obtain: yh 1 1.75 2 2.25 3 2.75 4 1.25 5 4.5 6 2.0 7 3.0 and finally, for descr(8): f(3.5,0.25,1.6)=2.6 Comparison: ''''''''''' Both approaches use the same amount of memory in tdm_table_data, but the iterative approach seems simpler to implement. It also has less overhead because recursion needs to create a new context repeatedly. Note that keeping x1 and xh in separate vectors makes it possible to determine the interpolation polynomials (splines etc.) for the tables at the innermost dimension once as part of the definition of tdm_table_data. If this is done, the tdm_table_data needs to be changed to store the interpolation polynomials instead of x1, y1.