​ ​ ​ ​ "General Purpose" Interpolation subroutines
Results 1 to 4 of 4

Thread: "General Purpose" Interpolation subroutines

  1. #1
    Senior Member
    Join Date
    Feb 2012
    Location
    London
    Posts
    2,872

    Default "General Purpose" Interpolation subroutines

    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.

    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 w1
    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.

  2. #2
    Senior Member
    Join Date
    Feb 2012
    Location
    London
    Posts
    2,872

    Default

    Hi,

    A recent topic in the Active Forum reminded me that I didn't "finish" this thread, so here are a few more details. Not yet the promised "reverse lookup" (X from Y) algorithm, but code for ARCTAN and an example of how the subroutine can use multiple tables (e.g. for SIN and ATAN).

    Unlike SIN and COS, which have a maximum range from zero to 1.0 (or down to -1.0 if angles outside the "first quadrant" are considered), the TAN function spans all the way from zero up to "infinity" (and the negative equivalent in other quadrants) which makes it difficult to create a complete integer-based lookup table. However, angles between 0 and 45 degrees do convert to a range from 0 to 1.0 (i.e. the same range as for SIN and COS). Then, angles between 45 and 90 degrees can be calculated from the inverse representation (i.e. by dividing unity by the complementary angle). Thus, {A}TAN is best implemented by using 8 octants, each spanning 45 degrees, whilst SIN and COS need only 4 quadrants of 90 degrees each.

    Normally, ATAN would be considered the "reverse" (or inverse) function, but ATAN (or more specifically ATAN2 which can receive two independent directional vectors) is generally more useful than TAN. The gradient of the lookup curve (table) doesn't change greatly over the range of 0 to 45 degrees (i.e. the X and Y axes can have comparable resolutions), so it's practical to create a table using ATAN for the forward (X to Y) conversion direction. TAN can then be calculated by the reverse lookup algorithm, or from the knowledge that TAN(x) = SIN(x) divided by COS(x).

    Accuracy:


    The thread linked at the top of this post paid particular attention to the accuracy of the ATAN function, so let's consider what compromises are necessary or worthwhile here. The available PICAxe table space is 256 bytes (in either EEPROM or TABLE memory), i.e. an average of 16 words per table if all 8 potential tables were implemented. However, the two tables below can handle all the "Trig" functions, so probably a total of 3 or 4 tables is as many as might be needed in practice, so a larger number of points per table can be considered.

    Three basic factors determine the "accuracy" of the results: the "resolution" (i.e. numerical range) of the tables' "output" values (Y), the resolution of the input values (X) and the error due to the degree of "curvature" of the table values (i.e. how far the actual function deviates from a straight line joining each pair of reference points). Also the "scaling" of the input/output values (e.g. representing 1.0 by an integer value of 10,000 , or 36,000 to represent 360 degrees, etc.) may affect the resolution/accuracy. But here, using the full range of a Word (i.e. 0 - 65535) to represent one revolution (360 degrees) and the range -1.0 to +1.0, makes the best use of the available number system.

    The table output values (Y) are stored in full Words (2 bytes), so are unlikely to limit the accuracy/resolution any more than any subsequent numerical processing. However, the chosen subroutine limits the number of interpolation points to 256 within each "line segment", so the use of 16 lines (17 lookup points) limits the "input" resolution to 4096 values. This is the limiting factor for the fairly linear {A}TAN function, so it may be worthwhile to increase the table to 32 segments (33 words, or 66 bytes). It's probably unnecessary to go any further because the whole table represents only one octant (45 degrees) so the effective output value also has a nominal resolution of 8192 points. For compatibility, the Y values have been scaled to a value of 32768 = Unity, but the "error" value at the centre of the "worst" line is no more than 2 in the Least Significant digit, so becomes negligible when scaled to a full revolution.

    The SIN function has considerably more curvature so the worst interpolation error is about 5 in the LS digit, when using a 32 section lookup table. Therefore, it may be worthwhile to increase the table to 64 segments (130 bytes) to reduce the worst interpolation errors to around 1 in the LS digit. However, even the 5 digits error (for 32 lines) scales to less than 0.0002 (i.e. 5 / 32768) , or a little more if alternate data values are simply deleted from the table (because the remaining points no longer have a symmetrical best fit to the maximum errors).


    So here are sample subroutines for SIN, COS and ARCTAN with a test harness and conversion factors for scaling to integer values of 10,000 (i.e. 4 decimal places) and 36,000 (i.e. 360 degrees with two decimal places). Also the octant-handling for ATAN2 is shown, but adding a test harness and division routine for this will exceed the 10,000 character forum limit.

    Code:
    ; SIN,COS and ARCTAN functions,	AllyCat,	March 2018  
    #picaxe 08m2		; And most others
    
    symbol TOM = 32				; Multiplier for Lookup Table number (i.e. flags * 32)
    symbol tempb = b1			; Temporary (local) and parameter-passing Byte 
    symbol tempw = w1			; Temporary (local) 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 NEGX = 32768
    symbol degs = w9
    symbol tan = w10
    
    symbol TABBASE = 0			; Table of table base addressess (max 8 bytes)
    symbol SINBASE = 8			; Start address of Sine table (length 130 bytes)
    symbol SINTAB = 0 * TOM			; Sine/Cos Lookup table pointer
    symbol ATANBASE = 138
    symbol ATANTAB = 1 * TOM
    
    data TABBASE,(SINBASE,ATANBASE)
    
    data SINBASE,(0,0,36,3,72,6,107,9,140,12,171,15,200,18,226,21,249,24)	; 64 segment SINE table
    data (12,28,26,31,36,34,41,37,39,40,32,43,18,46,252,48)
    data (223,51,187,54,142,57,88,60,24,63,207,65,124,68,30,71)
    data (181,73,65,76,193,78,53,81,156,83,247,85,68,88,131,90)
    data (181,92,216,94,237,96,243,98,234,100,209,102,168,104,111,106)
    data (37,108,203,109,96,111,228,112,86,114,183,115,6,117,67,118)
    data (109,119,134,120,139,121,126,122,94,123,43,124,229,124,140,125)
    data (31,126,158,126,11,127,99,127,168,127,218,127,247,127,2,128)
    
    data ATANBASE,(0,0,23,5,44,10,61,15,69,20,68,25,54,30,26,35,238,39)		; 32 segment ATAN table
    data (176,44,94,49,247,53,122,58,230,62,56,67,114,71,146,75)
    data (152,79,132,83,85,87,12,91,168,94,42,98,147,101,225,104)
    data (24,108,53,111,59,114,41,117,255,119,192,122,107,125,1,128)
    
    test:  
    	for tan = 0 to 10
    		w1 = tan * 1000 ** 53687	; Scale to unity = 8192 (32 line segments)
    		b1 = 0				; Set the ATAN2 flags to first quadrant
    		call arctan
    		w1 = w1 ** 36000		; Scale to hundredths of a degree
    		sertxd(cr,lf,"ArcTan ")
    		if tan = 10 then
    			sertxd("1.0")
    		else
    			sertxd("0.",#tan)
    		endif
    		b1 = w1 // 100	
    		w1 = w1 / 100
    		sertxd(" = ",#w1,".")
    		call showlz
    		sertxd(" degrees")
    	next
    	for degs = 0 to 360 step 15
    		w1 = degs * 100			; full scale = 36000 (hundredths of degrees)
    		w1 = w1 ** 53769 + w1	; scale up to 65536 = 360 degrees
    		call sine
    		sertxd(cr,lf,"Sin ",#degs," = ")
    		call showdec
    		w1 = degs * 100
    		w1 = w1 ** 53769 + w1	
    		call cosine
    		sertxd(" Cos = ")
    		call showdec
    	next
    	end
    	
    showlz:			; Reinsert a suppressed leading zero
    	if b1 < 10 then
    		sertxd("0")
    	endif	
    	sertxd(#b1)
    return
    	
    showdec:			; Show w1 signed with 4 decimal digits
    	if w1 => NEGX then 
    		sertxd("-")
    		w1 = -w1
    	endif	
    	w1 = w1 ** 20001			; Convert 32768 to 4 decimal digits
    	b1 = w1 / 10000	
    	w1 = w1 // 10000
    	sertxd(#b1,".")
    	b1 = w1 / 100
    	call showlz
    	b1 = w1 // 100
    	goto showlz 	; Then use its return
    
    atan2:			; ARCTAN2 subroutine (signed vectors in wx and wy) 
    	b1 = 0 				; First octant (0 - $1FFF)	
    	if wx >= NEGX then 		; Value is negative
    		b1 = 7
    		wx = -wx			; Convert negative to positive
    	endif					; WX is (now) positive
    	if wy >= NEGX then 		; Value is negative
    		b1 = b1 xor 3
    		wy = -wy			
    	endif					; WY is (now) positive
    	if wx > wy then 			; WX/WY is not a proper fraction ( <1 ) so swap
    		b1 = b1 xor 1			; Mark that WX > WY (Octant 1,2 etc)
    		w1 = wx			; )
    		wx = wy			; }Exchange WX and WY (~5 code bytes less than SWAP)
    		wy = w1			; )
    	endif					; Divisor larger than numerator so "fractional" result
    
    ;** Call Double-Word Division subroutine (w1 = wx:00 / wy) , then fall into arctan **
    
    arctan:				; w1 = 8192 represents unity (for 32 segment table)
    	b1 = b1 + ATANTAB		; Select Table
    	call interpXtoY
    	w1 = w1 + 2 / 4			; Full scale for one octant is 8192
    	if bit8 = 1 then			; LSB of b1
    		w1 = $2000 - w1		; Reverse the sector from 90 to 45 degrees
    	endif
    	w1 = b1 & 7 * $2000 + w1		; Add the other octants
    return
    
    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)	
    possine:	
    		if w1 > $3FFF then		; Reverse 2nd and 4th quadrants
    			w1 = $8000 - w1
    		endif
    		b1 = SINTAB			; Select Table 	
    		call interpXtoY
    	else
    		w1 = w1 - $8000			; Invert before and after positive calculation
    		call possine
    		w1 = -w1		
    	endif 		
    return
    
    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
    	b3 = b2 + 1				; Point to high byte
    	read b2,b2				; Read low byte
    	read b3,b3				; Read high byte
    	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
    return
    Cheers, Alan.
    Last edited by AllyCat; 29-03-2018 at 00:01.

  3. #3
    Senior Member
    Join Date
    Feb 2012
    Location
    London
    Posts
    2,872

    Default

    Hi,

    Here is the subroutine to perform interpolation in the "reverse" Lookup direction, i.e. from Y values to X values. It can use the same table(s) as above to perform "inverse" functions such as ARCSIN and ARCCOS, but also may be needed for functions with "asymmetric" scales such as Square Root or Logarithms. However, it is not the overall scales which are significant (because they can be adjusted by simple multiplication or division) but the change in the gradient of the curve (or transfer characteristic) over the range of the table.

    For example with a square-law function, when the X-value changes from 0 to 1, the Y value also increases from 0 to 1 so the average gradient is 1. However, when X changes from 10 to 11, Y changes from 100 to 121 so the gradient is 21, and when X changes from 100 to 101, Y changes from 10000 to 10201 with a gradient of 201. Typically the Y-scale can be from 0 to 65535, whilst the X-scale with 16 line segments (each having 256 sub-points) would be from zero to 4095. Hence it's sensible to create a "SQUARES" table for the forward direction (even though the basic calculation doesn't need to use a lookup table) so that the larger range of values appears along the Y-axis.

    The following program demonstrates the subroutine in a simple test harness. It uses only one SINE table, but with the choice of two resolutions (number of line segments), smaller than in the previous post to speed the table search and to show how the scaling needs to change. To select the lookup data from multiple tables, a subroutine is included, but this could be a Macro, or even replaced by a fixed constant assignment, if only one table is required.

    The main data values are passed to and from the subroutine in W1 and may need to be scaled before and/or afterwards, depending on the number of table segments and the User-Interface parameters (degrees, radians, etc.). Also, B1 contains flags which are combined to indicate the lookup table number and the Sign, Quadrant or Octant values of trigonometric functions, etc.. The table address is recalculated after the interpolation stage to avoid the need to store it in an additional variable.

    Code:
    ; Interpolation subroutine for "reverse" lookup through tables (X from Y) , AllyCat , April 2018
    #picaxe 08m2					; Or most others
    #define fastdiv
    #define sinseg32
    
    symbol NEGX = 32768				; Negative numbers boundary
    symbol TAM = 32					; Lookup Table number Multiplier
    symbol tempb = b1				; Temporary (local) and parameter-passing byte 
    symbol tempw = w1				; Temporary (local) and parameter passing word
    symbol wx = w2					; Interpolator X-word
    symbol wy = w3					; Interpolator Y-word
    symbol tabbase = 0				; Table of table bases (max 8 pointer bytes)
    symbol sinebase = 8				; Start address of Sine table (length 66 bytes)
    symbol nextbase = sinebase + 66
    symbol sinindex = TAM * 0			; Sine/Cos Lookup table pointer flags
    data tabbase,(sinebase,nextbase)		; Pointers to base of tables
    
    #ifdef sinseg32
    symbol SINSC = 2								; Scale multiplier for 360 degrees = 65536
    data sinebase,(0,0,72,6,141,12,201,18,250,24,27,31,41,37,33,43)			; 32 Segment Sine table
    data (254,48,188,54,90,60,209,65,32,71,67,76,55,81,250,85,133,90)		; With Symmetrical errors
    data (219,94,246,98,212,102,114,106,206,109,232,112,187,115,71,118)
    data (138,120,130,122,47,124,143,125,162,126,103,127,222,127,255,127)
    #else										; 16 line segments
    symbol SINSC = 4								; Scale multiplier
    data sinebase,(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)
    #endif
    
    demo:
    	for b0 = 0 to 100 step 5				; Hundredths of unity
    		w1 = b0 * 100 ** 53687 * 4 max 32767		; Scale Unity = 32768 
    		call arcsin
    		w1 = w1 ** 36000				; Scale to hundredths of a degree
    		b1 = w1 // 100					; Hundredths
    		w1 = w1 / 100					; Unit degrees
    		sertxd(cr,lf,#b0,"% = ",#w1,".")
    		w1 = b1 / 10
    		b1 = b1 // 10
    		sertxd(#w1,#b1," degrees")	
    	next b0
    end
    
    arcsin:
    	if w1 >= NEGX then 				; Value is negative
    		w1 = -w1				; Convert to a positive value
    		call parcsin
    		w1 = -w1				; Shift result by two quadrants
    	else						; Positive calculation
    parcsin: 			; Calculate ArcSin of a positive number (in w1)
    		b1 = sinindex 				; Select SINE table
    		call interpYtoX
    		w1 = w1 * SINSC				; Scale dependent on number of segments
    	endif
    return
    
    arccos:
    	call arcsin
    	w1 = $4000 - w1					; Cos is shifted one quadrant (90 degrees)
    return
    
    getbase:						; Load b2 with base address of the required table
    	b2 = b1 / TAM + tabbase				; b1 contains flags for table number
    	read b2,b2
    return
    
    interpYtoX:						; Enter with table number (*32) and flags in b1 (retained)
    	wx = w1						; Free up b2 and b3 to use separately
    	call getbase					; Load b2 with the base address of SINE table
    	do
    		b2 = b2 + 2				; Move to RHS of (next) cell
    		read b2,word wy				; Read Y-value at top of cell (NE corner)
    	loop until wy => wx				; Repeat until target point is within the cell  	
    	dec b2						; Point to High byte of LH corner (SW) of cell 
    	read b2,b3					; Read High byte of base of cell
    	wx = -b3 * 256 + wx   				; Subtract high byte of base of cell
    	wy = -b3 * 256 + wy 
    	b3 = b2 - 1					; Pointer to low byte
    	read b3,b2					; Read Low byte of base of cell
    	wx = wx - b2					; Subtract low byte of base of cell
    	wy = wy - b2
    	call getbase					; Read start of table again
    	b3 = b3 - b2 				; Calculate cell number * 2 (or divide by 2 to load directly in w1)
    
    #ifdef fastdiv						; Fast version if cell height is restricted to 4095
    	wx = wx * 16					; Multiply numerator by 16
    	wy = wy / 16					; Divide divisor by 16	
    #else							; Scale numerator and divisor for best resolution
    	for b2 = 0 to 7				 	; Multiply numerator/divisor by 256 to give result = 0 - 255
    		if wx < 32768 then 			; WX can be doubled
    			wx = wx + wx			; Multiply numerator by 2
    		else			
    			wy = wy / 2			; Divide divisor by 2	
    		endif	
    	next 	
    #endif				
    	wx = wx / wy max 255				; Horizontal distance within cell  (max = 255)
    	w1 = b3 * 128 + wx				; Add the value of the LH edge of cell
    return
    Details of the Subroutine's Operation:

    The reverse subroutine is rather more complex than the forward one for several reasons: Firstly, the algorithm needs to search through the list of Y-values to find the relevant "cell" (because they are not of equal height), whilst the constant-width of the cells in the X-direction (256 points) made calculation of the cell number particularly easy. The search process is not difficult, but it does consume more time and also it's important to ensure that the table data is constructed such that the search cannot fail (loop endlessly) because the subroutine does not set any timeout (overflow) limits.

    But the main complication is the calculation of the interpolation points within each cell. In the forward direction, the "base line" (X) was fixed at 256 points, so calculating the ratio of Y/X (i.e. the line gradient) was easy because division by 256 (a byte) in binary is simple. However, in the reverse direction, the divisor (for X/Y) is not only a variable but can be very large, potentially up to most of a full Word (two bytes):

    The "double-Word" division subroutine that I have previously documented can perform the required function but it has disadvantages; It is relatively slow and uses the PICaxe b1 variable which (in this family of subroutines) already contains various table and maths flags. Of course another variable (or two) could be introduced, but I have tried to keep the "RAM footprint" of these routines as small as possible. The introduction of another variable is also avoided by reading and subtracting the table's lower word one byte at a time, but this doesn't affect the executed code very much because the WORD qualifier is basically a predefined Macro anyway.

    Thus, two alternative (conditionally compiled) division routines are included. The simpler assumes that the height of each cell will never exceed a value of 4095, so it is safe to multiply by 16 without causing an overflow. At the same time the divisor is divided by 16, which can't cause an overflow, but might lead to some loss of accuracy. However, these conditions cannot occur arbitrarily or "unexpectedly" because they are determined by the data in the lookup table. So it is a valid simplification provided that there are sufficient cells and/or the scale of the Y-axis is restricted. For example, a SINE curve divided into 16 cells, with unity represented as 32768 in the Y-direction, should stay just within bounds (the "tallest" cell is just less than 2048 * PI / 2).

    The more "intelligent" method also employs an intrinsic multiplication by 256, but this is done by successively doubling the numerator when it can't cause a word overflow, or otherwise by halving the divisor. Then, the division can be performed by the PICaxe Basic interpreter quickly, with reasonable accuracy. This method is likely to be needed with exponential or "power" function tables, or if you're not confident about analysing the structure of the tables.


    Still to come: A more complete test harness to show all the Trigonometric expressions.

    Cheers, Alan.
    Last edited by AllyCat; 07-04-2018 at 16:00.

  4. #4
    Senior Member
    Join Date
    Feb 2012
    Location
    London
    Posts
    2,872

    Default

    Hi,

    Finally, a demonstration and test harness for all the Trigonometric functions above, in 5 degree increments. It can be run in the simulator but might take half an hour or more. At 900+ bytes, it's not intended as a formal Code Snippet, but shows how the modules and tables' scale/resolution work together.

    Two methods for TAN are shown, one dividing Sin / Cos (which needs one additional storage word), the other using the ATAN table in reverse. The tables have only 16 segments because the program is very close to the forum's 10,000 character posting limit.

    Code:
    ; General Purpose Interpolation routines with lookup tables for Trig Functions;
    ; Angles in Trig functions use One revolution = 65536 (16 bits)
    ; AllyCat , July 2017, updated April 2018
    
    #picaxe 20x2        		; Or any M2, except to compare embedded X2 SIN function
    #define fastdiv			; Delete if InterpYtoX has steep gradients
    
    symbol NEGX = 32768		; Boundary of Negative numbers (Twos complement)
    symbol TAM = 32			; Lookup Table number Multiplier
    symbol tempb = b1		; Temporary (local) and Parameter-passing byte 
    symbol tempw = w1		; Temporary and Parameter-passing word
    symbol wx 	= w2		; Interpolator X-word and Division Extension word
    symbol wy  	= w3		; Interpolator Y-word and Divisor
    symbol carry = bit31		; Top bit of w1
    symbol tabbase = 0		; Table of table bases (max 8 bytes)
    symbol atanbase = tabbase + 8	; Address of base of ATAN lookup table (length 34 bytes)
    symbol sinebase = atanbase + 34	; Start address of Sine table (length 34 bytes)
    symbol tanindex = TAM * 0		; ArcTan Lookup table pointer
    symbol sinindex = TAM * 1		; Sine/Cos Lookup table pointer
    
    data tabbase,(atanbase,sinebase)	; Pointer to bases of tables
    
    ; 16 line segments
    symbol TANSC = 8		; Scale multiplier
    data atanbase,(0,0,45,10,70,20,56,30,241,39,98,49,126,58,61,67,151,75)	; "Balanced" error curves
    data (137,83,17,91,47,98,231,104,57,111,45,117,197,122,5,128)			; ATan table
    
    symbol SINSC = 4			; Scale multiplier
    data sinebase,(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 ($7FFF = 1)
    
    ; TEST CODE
    
    	sertxd(cr,lf,"DEGs X2%_SIN  COS ASIN {Sin/Cos} ATAN2 {TAN}")
    
    	for b0 = 0 to 72 				; One revolution in 5 degree increments
    		wx = b0 * 5			; Degrees
    		w1 = sin wx				; X2 only
    		sertxd(cr,lf,#wx," ") 
    		if w1 > 127 then			; It's a negative value
    			sertxd("-")
    			w1 = w1 - 128
    		endif
    		sertxd(#w1,"% ")
    		w1 = b0 ** 14563		; 910.222  => 14563/65536 + 910
    		w1 = b0 * 910 + w1		; Next 5 degree step
    		w12 = w1				; Copy to use with cosine
    		call sine
    		w10 = w1				; Store sine result
    		call showS4dp
    		w1 = w12
    		call cosine
    		w11 = w1				; Store cosine result
    		call showS4dp
    		w1 = w10
    		call arcsin
    		call showdeg			; Report the angle
    		wx = w10				; Sine of initial angle
    		wy = w11				; Cosine of initial angle
    		call calcTAN			; TAN = SIN / COS
    		wx = w10				; Sine of initial angle
    		wy = w11				; Cosine of initial angle		
    		call arctan2
    		call showdeg
    		w1 = w12
    		call TAN				; Report the TAN value
    	next
    end
    
    TAN:		; Using ATAN table in reverse
    	b1 = w1 / 8192 + tanindex		; Set the octant flags and table pointer
    	w1 = w1 // 8192 * 4 max 32767		; Mask and scale		Octant		0 1 2 3 4 5 6 7
    	if bit8 = 1 then			; Mirror about 45 degrees	Mirror		0 1 0 1 0 1 0 1 bit8
    		w1 = 32767 - w1			; 32768 overflows ddiv		Negative	0 0 1 1 0 0 1 1 bit9
    	endif					;				Complement	0 1 1 0 0 1 1 0 bit8xorbit9
    	sertxd(" {")
    	if bit9 = 1 then
    		sertxd("-")
    	endif	
    	call interpYtoX									
    	w1 = w1 * 8
    	if bit8 = bit9 then
    		sertxd("0.")
    	else						; Calculate 1 / w1	
    		wy = w1 / 2			; Divisor (15 steps per loop)
    		w1 = 32768
    		wx = 0
    		b1 = b1 & 7			; Set division to 15 bits
    		call ddiv
    		sertxd(#w1,".")
    		w1 = 0				; Remainder in wx
    		b1 = b1 & 7
    		call ddiv				; Fractional part
    	endif
    	w1 = w1 ** 20000			; Scale 10,000 = unity
    	call show4dp
    	sertxd("}")			 	
    return	
    
    calcTAN:			; Using SIN/COS
    	sertxd(" {")
    	b1 = 0
    	if wx => NEGX then
    		wx = -wx
    		b1 = 2
    	endif
    	if wy => NEGX then
    		wy = -wy
    		b1 = b1 xor 2
    	endif	
    	w1 = wx + wx
    	wx = 0
    	call ddiv				; Integer result in w1, remainder in wx
    showTAN:
    	if bit9 = 1 then 
    		sertxd("-")
    	endif
    	sertxd(#w1,".")			; Integer part
    	b1 = b1 & 7
    	w1 = 0
    	call ddiv				; Fractional part
    	w1 = w1 ** 20000
    	call show4dp
    	sertxd("}")
    return		
    
    showS4dp:					; Report (fractional) w1, signed to 4 digits (i.e. * 10,000)
    	sertxd(" ") 
    	if w1 > NEGX then			
    		sertxd("-")
    		w1 = -w1
    	endif
    	w1 = w1 + 2 ** 20000 		; Convert to 4 decimal places
    show4dp:					; Unsigned
    	if w1 > 999 then nosupz
    	sertxd("0")
    	if w1 > 99 then nosupz 
    	sertxd("0")
    	if w1 > 9 then nosupz 
    	sertxd("0")	
    nosupz:
    	sertxd(#w1)
    	return
    
    showdeg:
    	w1 = w1 ** 36000		; Convert to hundredths of a degree
    	b1 = w1 // 100			; Hundredths
    	w1 = w1 / 100			; Units
    	sertxd(" ",#w1,".")
    	if b1 < 10 then 
    		sertxd("0")			; Don't suppress leading zero in decimal
    	endif
    	sertxd(#b1)
    	return
    
    cosine:
    	w1 = w1 + $4000			; Now fall through into sine 
    sine:						; Enter with angle in w1 (2^16 = 360 degrees)
    	if w1 < NEGX then posine	; Shift down if 180 - 360 degrees	
    		w1 = w1 - $8000
    		call posine
    		w1 = -w1
    	return
    posine:
    	if w1 > $3FFF then		; Reverse 2nd and 4th quadrants
    		w1 = $8000 - w1
    	endif
    	b1 = sinindex 			; Select Table + first quadrant (0)
    	w1 = w1 / SINSC			; b1=Table number *32 (can retain flags in b1.0 to b1.4)
    	goto interpXtoY			; Then use its return
    
    interpXtoY:	 				; b1=Table number *32 (can retain flags in b1.0 to b1.4)
    	wx = w1				; Free up b2 and b3. Note value in b3 used later as hi byte of w1  
    	call getbase			; Get pointer to base of table into b2 
    	b2 = b2 + b3 + b3		; Add 2* cellnumber (b3 is hi byte of w1) = SW corner of cell
    	b3 = b2 + 2				; Point to RHS of cell
    	read b3,WORD wy			; Read top (NE corner) of cell
    	b3 = b2 + 1				; Point to high byte
    	read b2,b2				; Read low byte of height of base of cell (for SW corner)
    	read b3,b3				; Read high byte to form whole word in w1 (= height of baseline)
    	wy = wy - w1 			; Calculate height of cell (range 0-65535)
    	wx = wx / 256 * 65280 + wx	; Subtract Cell number * 256 from wx (65280 is -1 * 256) to give: 
    	w1 = wx * 256 ** wy + w1	; Hor.Dist.inside cell (0-255) * gradient + height of cell base
    	w1 = w1 max $7FFF 		; Avoid overflow into negative flag (probably unnecessary)
    return 
    
    arccos:
    	call arcsin
    	w1 = $4000 - w1			; Cos is shifted one quadrant (90 degrees)
    return
    
    arcsin:
    	if w1 >= NEGX then 		; Value is negative
    		w1 = -w1			; Convert to a positive value
    		call parcsin
    		w1 = -w1			; Shift result by two quadrants
    	else					; Positive calculation
    parcsin:
    		b1 = sinindex 		; Select Sine Table
    		call interpYtoX
    		w1 = w1 * SINSC	
    	endif
    return
    
    getbase:					; Use flags in b1 to fetch base address of required table into b2
    	b2 = b1 / TAM + tabbase
    	read b2,b2
    return
    
    interpYtoX:					; Enter with table number (*32) and flags in b1 (retained)
    	wx = w1				; Free up b2 and b3 to use separately
    	call getbase
    	do
    		b2 = b2 + 2			; Move to RHS of (next) cell
    		read b2,word wy		; Read Y-value at top of cell (NE corner)
    	loop until wy => wx		; Repeat until target point is within cell
    	dec b2				; Back to high byte of LHS
    	read b2,b3
    	wx = -b3 * 256 + wx   		; Subtract high byte of base of cell
    	wy = -b3 * 256 + wy
    	b3 = b2 - 1				; Pointer to low byte
    	read b3,b2
    	wx = wx - b2			; Subtract low byte of base of cell
    	wy = wy - b2
    	call getbase
    	b3 = b3 - b2 			; Cell number * 2
    
    #ifdef fastdiv				; Fast version if cell height is restricted to 4095
    	wx = wx * 16			; Multiply numerator by 16
    	wy = wy / 16			; Divide divisor by 16
    #else						; Scale numerator and divisor for best resolution
    	for b2 = 0 to 7			; Multiply numerator/divisor by 256 to give result = 0 - 255
    		if wx < 32768 then 	; WX can be doubled
    			wx = wx + wx	; Multiply numerator by 2
    		else
    			wy = wy / 2	; Divide divisor by 2
    		endif
    	next 	
    #endif				
    	wx = wx / wy max 255		; Horizontal distance within cell  (max = 255)
    	w1 = b3 * 128 + wx		; Add the value of the LH edge of cell
    return
    
    ddiv:						; Divide double-word numerator by 15-bit divisor
    	for b1 = b1 to 239 step 16	; Retain starting value and repeat for each bit position
    		wx = wx + wx + carry	; Rotate numerator and result left (top bit lost)
    		w1 = w1 + w1		; Shift done
    		if wx => wy then 		; Skip if can't subtract
    			w1 = w1 + 1	; Add the flag to result (in low word)
    			wx = wx - wy	; Subtract		
    		endif
    	next b1
    	b1 = b1 & 15			; Clear the counter
    return
     
    arctan2:
    	b1 = 48 		; Clear three lowest flag bits for first octant and select 12-bit ddiv
    	if wx >= NEGX then		; Value is negative
    		b1 = b1 xor 7		; Set the three flags
    		wx = -wx			; Convert negative to positive
    	endif					; wx is (now) positive
    	if wy >= NEGX then 		; Value is negative
    		b1 = b1 xor 3
    		wy = -wy
    	endif					; wy is (now) positive
    	if wx > wy then 			; wx/wy is not a proper fraction (< 1) so swap
    		b1 = b1 xor 1		; Mark that wx > wy (Octant 1,2 etc)
    		w1 = wx			; )
    		wx = wy			; }-Exchange wx and wy (~5 code bytes fewer than SWAP)
    		wy = w1			; )
    	endif					; Divisor is now larger than numerator so "fractional" result
    	w1 = 0
    	call ddiv				; Result in w1=wx/wy, Can now use wx,wy as temp variables
    arctan:					; Enter with positive 12-bit value in w1 (and 3 flags set)
    	b1 = b1 & 7 + tanindex
    ; Scale here if needed  	
    	call interpXtoY			; Fullscale = 4096 for one octant, returns max $8000
    	w1 = w1 + 2 / 4			; Round up and scale down
    	if bit8 = 1 then			; LSB of b1
    		w1 = $2000 - w1		; Reverse the sector from 90 to 45 degrees
    	endif
    	w1 = b1 & 7 * $2000 + w1	; Add the appropriate octant number (45 degrees = 4096)
    return
    Cheers, Alan.

Bookmarks

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •