Hi,
PICaxe Basic struggles to implement the more complex mathematical calculations (such as trigonometric or logarithmic functions) because the interpreter is slow and doesn't support signed or fractional (non-integer) numbers. One solution is to use a "lookup table" of results, but normally a table cannot contain every possible value that might be required. Therefore, it is necessary to "interpolate" values between the sample-points which have been stored in the table.
The mathematical function may be considered as a graph of Y versus X coordinates, with a lookup table containing sample points along the curve. The simplest form of interpolation assumes that adjacent points are joined by a straight line, i.e. that the smooth curve is replaced by a "piecewise linear" representation. In this case, the "target" result can be obtained from the coordinates of the start and end points of the appropriate line, then calculating the gradient (slope) of the line and the proportional distance of the target point along that line.
The two subroutines to be described in this thread are part of a family of routines which have a consistent structure as described in the thread here. To minimise the size of the algorithm and the number of (RAM) variables required, a few constraints are placed on the range of the numerical values. The primary limitations here are that the table may contain a maximum of 128 sample points (normally it will be much less) which must be equally-spaced along the X-axis, with 255 intermediate interpolated values available (i.e. defined by a single byte) along each segment between the reference points. Each sample point may have a maximum value in the Y-axis of 16-bits (i.e. two bytes = 65535).
Another constraint is that the interpolating lines must always have a positive (or zero) gradient, i.e. each subsequent Y-value in the table may not be smaller than the previous point. In practice this is not a severe restriction because cyclic trigonometric functions can be considered as separate segments "mirrored" in the X or Y axes, with any negative gradients handled by simply reversing the numerical values along the X-axis. Such a restriction also ensures that the "reverse" lookup (i.e. from y-axis to X-axis) using the same data produces a unique solution (or at worst a contiguous range).
Thus, each interpolation line fits (as the diagonal) across a rectangular "cell" which has a width of 256 units and a maximum height of either 4095 or 65535 units (depending on the algorithm employed). Typically there will be 16 cells (i.e. rectangles with their NE and SW corners touching), to give a total horizontal (x-axis) range of 0 - 4095 and vertical (Y-axis) range of 0 - 65535. These integer values can of course be scaled and/or offset, to suit the numerical range and resolution of the dataset.
These requirements dictate that most cells (and the graph of the whole dataset) will be principally "tall and thin" (e.g. 65536 high by 4096 wide). So, if the dataset converts from a "large" range of numerical values to a smaller range, such as a logarithmic function, then it is necessary to perform a "reverse lookup" function, i.e.. from the Y-axis to the X-axis. This is the reason that two different subroutines are required, one to convert from a point on the X-axis to the corresponding value on the Y-axis and a second to accept a point on the Y-axis and calculate the value along the X-axis. The latter routine is slightly more complex because it needs to "search" through some or all of the (Y) entries in the table, to find the cell which contains the target datapoint.
The lookup tables are "indexed" so that up to 8 different datasets may be used, for example (arc)sine/cos, (arc)tan{2}, exponential/log and NTC resistance/temperature, etc. This is why the number of datapoints per function will be typically limited to around 16, i.e. 16 words of 2 bytes x 8 functions = 256 bytes, or one "page" in the EEPROM (table). The table index number has a value from 0 to 7, but is stored in the top 3 bits of b1 (to avoid conflicts with the flags available in the division subroutine described in the thread linked above). This number is converted by a preliminary lookup table to give a "pointer" address of the first entry in the appropriate table. Typically, the number of cells defined by each table will be a binary power (e.g. 8, 16 or 32), but any value may be used, within the constraints of available memory and scale factors.
So here is the first subroutine to calculate the value on the Y-axis which corresponds to a point on the curve/graph directly above a supplied value along the X-axis. This first example includes a Sine function table for one quadrant, divided into 16 piecewise-linear sections.
The "reverse" functions (arcsin/cos using interpolation from Y to X) are rather more complex because the table needs to be searched and then a "high resolution" division may be needed, so those subroutines will be held over for the next post.
Also, apart for the sine/cos function, a test harness is not included yet, but a very similar subroutine is included in the Sunrise and Sunset calculation in the thread here.
Cheers, Alan.
PICaxe Basic struggles to implement the more complex mathematical calculations (such as trigonometric or logarithmic functions) because the interpreter is slow and doesn't support signed or fractional (non-integer) numbers. One solution is to use a "lookup table" of results, but normally a table cannot contain every possible value that might be required. Therefore, it is necessary to "interpolate" values between the sample-points which have been stored in the table.
The mathematical function may be considered as a graph of Y versus X coordinates, with a lookup table containing sample points along the curve. The simplest form of interpolation assumes that adjacent points are joined by a straight line, i.e. that the smooth curve is replaced by a "piecewise linear" representation. In this case, the "target" result can be obtained from the coordinates of the start and end points of the appropriate line, then calculating the gradient (slope) of the line and the proportional distance of the target point along that line.
The two subroutines to be described in this thread are part of a family of routines which have a consistent structure as described in the thread here. To minimise the size of the algorithm and the number of (RAM) variables required, a few constraints are placed on the range of the numerical values. The primary limitations here are that the table may contain a maximum of 128 sample points (normally it will be much less) which must be equally-spaced along the X-axis, with 255 intermediate interpolated values available (i.e. defined by a single byte) along each segment between the reference points. Each sample point may have a maximum value in the Y-axis of 16-bits (i.e. two bytes = 65535).
Another constraint is that the interpolating lines must always have a positive (or zero) gradient, i.e. each subsequent Y-value in the table may not be smaller than the previous point. In practice this is not a severe restriction because cyclic trigonometric functions can be considered as separate segments "mirrored" in the X or Y axes, with any negative gradients handled by simply reversing the numerical values along the X-axis. Such a restriction also ensures that the "reverse" lookup (i.e. from y-axis to X-axis) using the same data produces a unique solution (or at worst a contiguous range).
Thus, each interpolation line fits (as the diagonal) across a rectangular "cell" which has a width of 256 units and a maximum height of either 4095 or 65535 units (depending on the algorithm employed). Typically there will be 16 cells (i.e. rectangles with their NE and SW corners touching), to give a total horizontal (x-axis) range of 0 - 4095 and vertical (Y-axis) range of 0 - 65535. These integer values can of course be scaled and/or offset, to suit the numerical range and resolution of the dataset.
These requirements dictate that most cells (and the graph of the whole dataset) will be principally "tall and thin" (e.g. 65536 high by 4096 wide). So, if the dataset converts from a "large" range of numerical values to a smaller range, such as a logarithmic function, then it is necessary to perform a "reverse lookup" function, i.e.. from the Y-axis to the X-axis. This is the reason that two different subroutines are required, one to convert from a point on the X-axis to the corresponding value on the Y-axis and a second to accept a point on the Y-axis and calculate the value along the X-axis. The latter routine is slightly more complex because it needs to "search" through some or all of the (Y) entries in the table, to find the cell which contains the target datapoint.
The lookup tables are "indexed" so that up to 8 different datasets may be used, for example (arc)sine/cos, (arc)tan{2}, exponential/log and NTC resistance/temperature, etc. This is why the number of datapoints per function will be typically limited to around 16, i.e. 16 words of 2 bytes x 8 functions = 256 bytes, or one "page" in the EEPROM (table). The table index number has a value from 0 to 7, but is stored in the top 3 bits of b1 (to avoid conflicts with the flags available in the division subroutine described in the thread linked above). This number is converted by a preliminary lookup table to give a "pointer" address of the first entry in the appropriate table. Typically, the number of cells defined by each table will be a binary power (e.g. 8, 16 or 32), but any value may be used, within the constraints of available memory and scale factors.
So here is the first subroutine to calculate the value on the Y-axis which corresponds to a point on the curve/graph directly above a supplied value along the X-axis. This first example includes a Sine function table for one quadrant, divided into 16 piecewise-linear sections.
Code:
picaxe 08m2 ; Or most others
#no_data
symbol TOM = 32 ; Multiplier for Lookup Table number (i.e. flags * 32)
symbol WIDTH = 256 ; Piecewise-Linear horizontal segment length
symbol tempb = b1 ; Temporary (local) and parameter-passing byte
symbol tempw = w1 ; Temporay (loacl) word and parameter-passing word
symbol wx = w2 ; Interpolation X-word and Numerator High word in ddiv
symbol wy = w3 ; Interpolation Y-word and Divisor in ddiv
symbol TABBASE = 0 ; Table of table base addressess (max 8 bytes)
symbol SINBASE = 8 ; Start address of Sine table (length 34 bytes)
symbol SININDEX = 0 * TOM ; Sine/Cos Lookup table pointer
data TABBASE,(SINBASE)
data SINBASE,(0,0,141,12,252,24,45,37,2,49,94,60,39,71,66,81,144,90)
data (1,99,125,106,243,112,83,118,144,122,157,125,117,127,255,127) ; Sine table ($8000 = +1)
cosine: ; Enter with angle in w1 (2^16 = 360 degrees)
w1 = w1 + $4000 ; Add one quadrant and fall through into sine
sine: ; Enter with angle in w1 (2^16 = 360 degrees)
if w1 < $8000 then ; Shift down if in 2nd or 3d quadrants (180 - 360 degs)
posine:
if w1 > $3FFF then ; Reverse 2nd and 4th quadrants
w1 = $8000 - w1
endif
w1 = w1 / 4 ; Max Table X-Axis = 16 * 256 = 4095 (Top 4 bits spare)
b1 = SININDEX ; Select Table and fall into interpXtoY
interpXtoY: ; Retain flags in .0 to .3
wx = w1 ; Frees up b2 and b3
b2 = b1 / TOM + tabbase ; Table address pointer in top 3 bits
read b2,b2 ; Address of start (word) of table
b2 = b2 + b3 + b3 ; Add 2* wx/256 (b3 still contains hi byte of w1) = base of cell
b3 = b2 + 2 ; Point to top of cell
read b3,WORD wy
call readw1
wy = wy - w1 ; Calculate height of cell (wy - wx)
wx = wx / 256 * 65280 + wx ; Horizontal distance across cell (65280 = -1 * 256)
w1 = wx * 256 ** wy + w1 ; Height inside cell (0<wx<256) * gradient + height of cell base
w1 = w1 max $7FFF
else
w1 = w1 - $8000 ; Invert befeore and after positive calculation
call posine
w1 = -w1
endif
return ; From interpolation or inversion routines
readw1: ; PE does not accept READ w1,WORD w1 directly (b2 is pointer)
b3 = b2 + 1 ; Point to high byte
read b2,b2 ; Read low byte
read b3,b3 ; Read high byte
return ; Result in w1Also, apart for the sine/cos function, a test harness is not included yet, but a very similar subroutine is included in the Sunrise and Sunset calculation in the thread here.
Cheers, Alan.
