The function ATN returns the angle that has the given X and Y components. The function HYP gives the length of the hypotenuse of the right triangle formed by X and Y (the magnitude of the vector). Taken together these operators allow you to convert from rectangular to polar coordinates. These were not documented in the early versions of the Stamp manual, but recently they have become official and more information can be found in the online help within the Stamp Editor, version 2.0 beta 1.
The input is -127 to +127 for X and Y, and for the output, angles around the circle go from 0 to 256 brads (instead of 0-360 degrees or 0 to 2pi radians). You can easily convert brads to degrees with a */360 operation as shown below. The Stamp measures angles counterclockwise from the positive X axis. The output magnitude is in the range of 0 to 179. Keep in mind that the resolution is only 8 bits, if that.
'{$STAMP BS2}
' Calculate inverse tangent (arctangent) and magnitude of a vector.
X con 100 ' in the range -127 to +127
Y con -100 ' in the range -127 to +127
angle var word
magnitude var word
angle = (X ATN Y) */ 360 ' order of X and Y matters. */360 converts from brads to degrees
debug ? angle ' gives 315 degrees, 45 degrees less than a full circle.
magnitude = X HYP Y ' the order of X and Y does not matter for HYP
debug ? magnitude ' gives 141 (looks like the square root of 2, huh? SQR(100^2 +100^2))
end
Incidentally, undocumented features can also turn up as "bugs". You might try to use "ATN" or "HYP" as names for variables in one of your own programs and be surprised to get an error message, "name already defined". You look in the "reserved words" list, and those names are not there. But they are in fact reserved by the IDE.
The BS2 has SIN and COS functions, but you have to understand how they work in order to make them work for you. Tnis note shows a few examples, including a way to calculate the inverse SIN and COS functions.
X = COS A ' in normal trig: 0<=A<360 degrees, -1<=X<=+1
' or 0<=A<=2*pi radians
' but...
' in stamp trig: 0<=A<255 brads, -127<=X<=+127
Note the scaling that the stamp uses. A full circle is 0-256 brads, instead of the familiar 0-360 degrees or 0 to 2*pi radians. The term "brad" stands for "binary radian". And the cosine of the angle goes from -127 to +127, to represent cosine values from -1 to +1. You can think of the range of X as a fraction that goes from -127/128 up to 127/128, always with a denominator of 128. Even though X covers only a small range in Stamp math, it must always be defined as a word variable, in order to represent negative numbers correctly. You might ask why the denominator is 128 instead of 127. The answer is, it doesn't make much difference, given the accuracy of the function in the stamp.
This is not a high accuracy function. When the angle goes from 0 up to 64 (one-quarter of a circle), the cosine goes from 127 down to zero. So the resolution is only about one part in 64. At small angles, near A=0, the cosine function changes very slowly. On the BASIC Stamp, the angles A=0, A=1, A=2, and A= 3 all produce the same result, X=127.
To convert from brads to degrees, multiply by 45 and then divide by 32 (45/32 = 360/256 simplified). Or simply use the */ operator with 360 on the right hand side. (The */ operator does an implicit division by 256.) Example, to convert 256 brads to degrees:
degrees = 256 * 45 / 32 ' gives 360 degrees
or...
degrees = 256 */ 360 = 360 ' also gives 360 degrees
The two have the same accuracy, but the */ form runs faster and takes less memory. To convert the other way, from degrees to brads, multiply degrees times 32 and divide by 45, or use */ or ** with an approximation of the fraction 256/360. Details on */ and ** elsewhere. For example, to convert 157 brads to degrees:
brads = 157 * 32/45 ' gives 111 brads
or...
brads = 157 ** 46603 ' gives 111 brads (46603/65536~=32/45)
or...
brads = 157 */ 182 ' gives 111 brads (182/256~=32/45)
And sometimes you want the value of X to span a more convenient decimal range, say from -100 to 100 instead of -127 to 127. Use the conversions following conversions. X1 is a value between -100 and +100, that represents fraction -100/100 up to +100/100. X is a value between -127 and +127 that represents the fraction -127/128 up to +127/128. The process of changing a denominator like that is called renormalizing.
X1 = abs X * 25 / 32 * (-2*X.bit15+1) ' renormalization X/128 to X1/100
or...
X1 = abs X */ 200 * (-2*X.bit15+1) ' 25/32=200/256 using *? operator
and to go the other way, from ...
X = abs X1 * 2 ** 41943 * (-2*X.bit15+1) ' renormalization X1/100 to X/128
' 32/25~=2*41939/65536
or...
X = abs X1 */ 983 / 3 * (-2*X.bit15+1)
' 32/25 ~= (983/256)/3
The trick here is that the values of X and X1 can be negative; the BASIC stamp cannot perform division on negative numbers. The procedure is to take the absolute value of X or X1, then multiply and divide, and then restore the sign based on the sign bit of the original number.
Here is a simple demo program that illustrates the use of the SIN and COS on the BASIC Stamp. It steps through all the angles from 0 to 360. The angle is printed in degrees and the cosine and sine of the angle are printed as decimal numbers in the range of -1.00 to +1.00. Then for good measure it calculates the tangent of the angle and prints that too, as a number from 0.00 up to 65535 (infinity). Recall that the tangent is equal to the sine divided by the cosine.
for A1=0 to 360
A=A1 */ 182 ' <--convert degrees to brads
X=COS A
Y=SIN A
X1 = abs X */ 200 * (-2*X.bit15+1) ' <-- renormalize.
Y1 = abs Y */ 200 * (-2*Y.bit15+1)
' now have the values in +/- decimal form to display in -0.xx form
debug "COS ",DEC A1,"= ",X1.bit15*13+32,"0.",DEC2 X1," "
debug "SIN ",DEC A1,"= ",Y1.bit15*13+32,"0.",DEC2 Y1," "
Z=abs Y * 100 / abs X * (X.bit15^Y.bit15*-2+1) ' calculate the tangent.
if X<>0 then okay
Z=65535
okay:
debug "TAN ",DEC A1,"= ",Z.bit15*13+32,DEC abs Z/100,".",DEC2 Z,cr
pause 100 ' slow it down so I can see it
next
The calculation of the tangent involves taking the absolute value of both numbers, and then restoring the sign at the end. The sign of the result is positive if both X and Y are positive or if they both are negative. Note how the decimal value is split up for display with the decimal point. It takes some additional tricks to display the sign and the decimal point correctly, as detailed elsewhere.
Sometimes you need the inverse trig funcitons, which come up in a lot in computer work with motion control and data analysis and in astronomy and time calculations. The inverse COS is sometimes written ARCCOS, or COS-1. ARCCOS X is the angle A that would give the value X=COS A.
The cosine function is periodic, so that many values of A will give X=COS A. So for convenience the value of A is restricted to 0<A<180 degrees (pi radians, 128 brads).
A = COS-1 X ' in normal trig, -1<=X<=1, 0<=A<=180 degrees or 0<=A<=pi
' but...
' in stamp trig -127<=X<=127, 0<=A<=128 brads
Similarly, SIN-1 Y is the angle A that would give Y = SIN A. The angle A is restricted to the values -90 to +90 degrees (+/-pi/2 radians, 64 to 192 brad). The BASIC Stamp doesn't like negative angles, so it will use 64 to 192 brads.
For example, in averaging vector wind direction, the sine and cosine components of the wind vector are separated out and accumulated. Then at the end, the totals for the x and y directions are divided by number of samples. The value (X2 + Y2) 1/2 is the average vector magnitude of the windspeed, and arctan Y/X is the average vector direction.
The quick trick to calculate COS-1 is to use the stamp's built-in cos function and loop to find the angle that will come closest to the given value of X. The following calculation restricts the domain even further, to 0<=X<=1 (0<=X<=127 for the stamp scaling), and the range to 0<=Y<=90 degrees (0<=Y<=64 brads for stamp scaling). These are angles in the first quadrant. Trigonometric identities can be used to calculate the answer for the second quadrant.
X var byte ' input value 0 to 99 (0.00 to 0.99)
A var byte ' angle in brads, 0<=A<=64.
Z var byte ' helper variable
for X = 0 to 99 ' demo, just step through the values of X
' calculate the arctan here:
Z = X */ 983 / 3 ' renormalize to 127
A=63-(Z/2) ' first approx, to minimize iterations
loop:
if COS A <= Z then done
A = A+1
goto loop
done:
A = A */ 360 ' convert brads to degrees using */ operator
' now have the angle in decimal degrees, 0->90
debug "arccos 0.", dec2 X," = ",dec A
pause 100 ' slow it down so I can see it
next ' do some more
The core of the routine to calculate arctan X is shown in red. It includes a step to renormalize the decimal fraction into a /127 fraction, and then at the end a step to renormalize from brads to degrees. This is neither blindingly fast nor NIST accurate--The initial angle=63-(X/2)' is an initial approximation to keep the number of iterations at a minimum (13 iterations worst case).
The above covers X=100 to X=0, to give zero to 90 degrees. How to find the other half, from X=0 to X=-100, to give 90 to 180 degrees? Use the identity:
arccos (-X) = 180 - arcos (X)
So if you had to find COS-1 (-37), you would use the above algorithm to find A=COS-1 (+37), and then calculate COS-1 (-37) = 180-A degrees.
To find SIN-1 Y, it is possible to develop an routine very similar to the above, and to use the relation
arcsin (-X) = -arsin (X)
or, if you need to find sin-1 Y, then find cos-1 Y instead, and use the identity
sin-1 Y = 90 - cos-1 Y
This is simply a restatement of the well-know fact that the two opposite angles in a right triangle add up to 90 degrees.
The BS2 doesn't offer a TAN function. But recall from high school trig that if tan A = a/b, then cos A = b/(a2+b2)1/2. If you need to find A = TAN-1 Z, then you can first calculate what the cosine of that angle would be, and then calculate the angle itself. For example, to find TAN-1 3, you can instead find COS-1 {1/(32+1)1/2} ...
TAN-1 3 = COS-1 {1/(32+1)1/2} = COS-1 (1/101/2)In stampese with angles in brads and normalization to 128, this becomes:
COS-1 40/128
which gives an angle close to the correct value, 71 degrees, from the above algorithm.
The Y=SQR X function of the BS2 takes a 16-bit argument X and returns an 8-bit result Y. This results in low resolution and low accuracy. E.g. the result is Y=10 for all X=100 to X=120, and then it is 11 for all values of X from 121 to 143. Suppose you want to find the square root of 115. The stamp will tell you that the square root is 10. But if 115 is 3/4 of the way between 110 and 120. So the square root is approximately 3/4 of the way between 10 and 11, that is, about 10.75. This is linear interpolation. Of course, the stamp can't use the decimal point, but it can give you an answer of 1075, where it is understood that the decimal point is two places over. The square root of 65536 is 256, but we can safely move the decimal place two over to become 25600, which represents 256.00. The interpolation formula just takes one line of BS2 code:
' To calculate square root on BS2
' using interpolation for greater resolution
' best for larger values of X
X var word
Y0 var byte
Y var word
for x=0 to 65535
Y0 = sqr X
Y = X-(Y0*Y0)*100+Y0/(2*Y0+1)+(Y0*100)
debug dec X,tab,dec Y0,tab,dec Y/100,".",dec2 Y,CR
pause 50
next
The demo program steps through all possible values of X, and shows both the simple value Y0 that comes from the BS2's built-in SQR function, and the interpolated value to two decimal places. The interpolation is more accurate at higher values of X, where the function is "flatter", that is, naturally more of a straight line.
Here is another way to extend the accuracy of the square root, by continuation. This uses an algorithm called the "high school algorithm". The method is similar to long division. When the BS2 takes the square root of a given number, it finds the largest whole number that, when squared, is still less than the given number.
y = sqr x ' this is what the BS2 calculates
'y*y <x which squared is less than x '(y+1)*(y+1)>x but y+1 squared is greater than x
'R = x - (y * y) This is the "remainder"
The following:starts with of this remainder as well as the Stamp's own SQR result, to extend the result to 14 bits.
sqr64:
y=sqr x ' BASIC Stamp does its thing to 8 bit result
R=x-(y*y) ' remainder to start
y=y*2
for ix=5 to 0 ' continuation, using the "high school" algorithm.
R=(R<<2)
y=y&~1<<1|1 ' how do you like this for a trick formula?!
if y>R then sqr64b
R=R-y
y=y|2
sqr64b
next
y=y*/200 ' calculate to two extra decimals
debug dec y/100,".",dec2 y,cr
return
The logarithm is important for calculations involving decibels, power levels and noise figures. Also it is important for some kinds of data processing where the signal ranges over a very wide range of levels. For example, light intensity and sound level can range over many orders of magnitude. Our sensory apparatus, our eyes and ears, are intrinsically responsive to the logarithm of the signal. Some kinds of sensors, such as thermistors, have a logarithmic response, so the best equation to derive the temperature from the resistance requires the log function. Many problems involving first order differential equations such as radioactive decay involve logs and exponentials.
What to do on the BS2? It has neither a log nor an exp function, or does it?
The BS2 does have a built-in log of sorts, the NCD function. The BASIC Stamp manual explains NCD x in terms of the highest bit set in the variable x. For example, the number 4097 in binary is, %1000000000001, and the thirteenth bit is the highest one set, so NCD 4097=13.
Also, 2^13=8192 is the exact power of two that is greater than (x=4097), and 2^12=4096 is the power of two that is less than x. The function (NCD x - 1) is the integer part, or characteristic, of the base 2 logarithm of X.
Let us suppose you need to make a VU meter for sound level that will be displayed as a bar graph on a 8 segment display, and that sound level from an analog to digital converter ranges from 0--256. Then the value to put on the display is simply
y = NCD x ' value that ranges from 0 to 8
' as x ranges from 0 to 255.
The first segment comes on for a sound level of 1, the second segment for sound levels from 2->3, and so on up to the eigth segment for sound levels from 128 ->255
The inverse of the NCD operator is the DCD operator. It in effect raises 2 to a certain power. For example, y = DCD 7 returns 128, because 128 is the seventh power of two. This is kind of a crude exponential function. The Stamp accepts values from x=0 to x=15, and returns the corresponding 16 values y=1, y=2, y=4, y=8 and so on up to y=32768. Pretty crude, when viewed as an exponential.
What if instead of just 8 segments, you want to make a VU meter to display sound level on a 100 segment bar graph on an LCD. The bitlog function is a simple way to expand the scale.
The graph below shows the values of x={1,2,4,8,16,32}, with the corresponding values of NCD x - 1 = {0,1,2,3,4,5} on the left hand y axis. The graph shows straight lines connecting the exact binary powers. The right hand y axis shows values that are four times the values on the left hand y axis. There are small dots on the straight lines that connect the main black dots. Those are going to be our points of interpolation. given by the bitlog function.
Okay,
it is a form of linear interpolation. The real logarithm connects the
major black dots (which are called the characteristic, or integer
part of the logarithm) with a curve (called the mantissa of the
logarithm). The curve "C" is shown in the inset to the right,
and the linear interpolation is the points along the straight line
that connects a and b. Not at all perfect, but the
linear approximation is an improvement for things like a VU meter.
The small figure to the right shows two of the binary powers, a
and b connected by the curve C. If we don't use the
interpolation, the resulting curve is just one big step,
S.
Lets say we want to interpolate 4 points. The formula to calculate the values is shown as in inset in the above graph. The first line of the formula calculates the NCD-1 characteristic, limiting the value to >0 The second line multiplies that preliminary result by 4, and then adds in the interpolation-- a value of 0,1,2 or 3 taken from the two bits next down from the most significant bit in the original value. Thus the name, "bitlog". The formula is just a tricky way of extracting the value of those two bits. All it is doing is that interpolation. The business with the min is there simply to handle values of x less than four, for which we cannot have four interpolated points, and to handle the value x=0.
Here is the formula for 8 (23) interpolated points:
y=NCD x min 1 - 1
y=y*8+(x- dcd y >>(y min 3 - 3)
If the input value of x can vary from 1 to 4095, then the value of the bitlog as defined by this formula varies from 0 to 96, just right for display on a 100 segment LCD bar graph. Here is the computation when you substitute x=4096...
y=NCD 4096 ' results in y=12, then...
y=12 min 1 - 1 ' results in y=11, then...
y=y*8 + (4096-dcd 11 >> (11 min 3 - 3))
^^^ ^^^^^^ ^^^^^^^^^^^^
2048
^^^^^^^^^^
2048 >> 8
^^^^^^^^^^^^^^^^
88 + 8
= 96
Here is the general formula, where there are to be 2^n interpolated points:
y=ncd x min 1 -1 : y=y*dcd n+((x-dcd y)>>(y min n-n))
For example,with n=8, there would be 256=28 interpolated points between each exact power of two.
y=NCD x min 1 - 1
y= y*256+((x-dcd y)>>(y min 8 - 8))
Entering with x=10000, the formula comes up with y=3384. Divide that by 256 on a calculator to get 13.218. Compare that to the real value, lg 10000 = 13.276. (Lg stands for the logarithm base 2) The bitlog function is a crude linear interpolation, but it gives a better approximation than just 13.000.
To convert to base 10, if you must, multiply the lg value times a constant factor, using the ** operator. (factor = 100/256*0.30103*65536=7706).
y10 = y ** 7706
This returns a value of y10=397 when x=10000. Compare that to the true value, 100*log 10000 = 400.
The inverse of the lg operator is the power of 2, as on the Stamp the inverse of NCD is DCD. Here is the inverse of the bitlog function. The value n in the formula is, as above, the number of steps of interpolation:
x=dcd (y/dcd n)
x=(x>>n)*(y-x)+x
There is a nice discussion of the bitlog function and its inverse in
Jack Crenshaw' Math Toolkit for Real-Time Programming, (c) 2000, CMP Books, ISBN: 1-929629-09-5
The next chore is to find a better approximation to the curve between the exact powers of two. The operator (NCD X - 1) returns the integer part of the base 2 logarithm of X. This integer part is usually called the "characteristic". The fractional part is called the "mantissa". For example 2^5.75 = 57.8368, and the other way of writing this is, log257.8368 = 5.75. The characteristic is 5, the mantissa is 0.75. On the Stamp with integer math, we have to move the decimal point over to come up with a value of 575 = 100* log2 58.
The following program uses the NCD function built into the BS2 to find the characteristic, and then looks up the mantissa in a table with interpolation. The table here consists of five word-size constants, but the program makes it easy to change the number of table entries where more memory is available, such as on the BS2SX, or in external eeprom. The values to use in a 66 byte table are listed later on.
It is important to understand that we do not need a separate table for every octave between powers of two. the same table works for every octave. It is important to understand that the logarithm has a property of periodicity. The shape of the mantissa between 1 and 2 is in all respects congruent to the shape it has between 2 and 4, or between 4 and 8, or between any two numbers that are multiple of two apart. The difference is only a question of scale. In this way the exponential functions are related to the trigonometric functions, one of the deep notions of mathematics.
' LOG2T.BS2
' Tracy Allen
' Example of how to determine log N
' Integer part (characteristic) comes from BS2 function, NCD.
' Fractional part (mantissa) from interpolation in a table in eeprom
' Could use external memory for the table if available.
' Conversion to other bases using constant basex:
' 65536*Log(2)= 65536*0.30103=19728. (19728.3)
' 65536*Ln(2)=65536*0.69315=45426 (45426.2784)
Ld con 2 'number of table entries is (2^LD+1)
basex con 65535 ' 65535 for base 2, 19728 for base 10, 45426 for base e
Ldata data word 0, word 322, word 585, word 807, word 1000
' x 0.000 1.250 1.500 1.750 2.000
'log2(x) 0.000 .3219 .5850 .8074 1.000
N var word ' input to log2 routine
M var word ' output and intermediate variable
Q var word ' lower estimate, (and log base e)
Q0 var Q.byte0 ' for table lookup, implied array
P var word ' upper estimate, (and log base 10)
L var nib ' integer part of log2(N)
K var nib ' adjusted for interpolation
J var nib ' fractional catagory of log2(N)
I var nib ' index into table
demo: ' just pick a number N and find its log
N=48756
debug "N=",dec N,CR ' show the number
gosub log2 ' find lg2 base 2 and display:
debug "lg N=",dec M/1000,".",dec3 M,CR
P = M ** 19728 ' convert to log base 10 & display:
debug "log N=",dec P/1000,".",dec3 P,CR
Q = M ** 45426 ' convert to ln base e and display:
debug "ln N=",dec Q/1000,".",dec3 Q,CR
end
log2: ' input N, returns M = log2(N)
L = NCD N - 1 ' L=integer part of log2(N)
M = DCD L ' M < N < 2M between powers of two
K = L min Ld - Ld ' restrict value for interpolation
J = N - M >> K ' fractional catagory, 2^Ld catagories
for I=0 to 3 ' read values from table
read J*2+I,Q0(I) ' Q < frac. part of log2(N) < P
next
' now interpolate
' note 32 bit math, divide by M is >> K
M = N-M-(M>>Ld*J)**(P-Q)<<(16-K)+(N-M-(M>>Ld*J)*(P-Q)>>K)+Q+(L*1000)
return
Values for Ld=5, 32 catagories, 66 byte table for mantissa: 0,444, 875, 1293, 1699, 2095, 2479, 2854, 3219, 3576, 3923, 4262, 4594, 4919, 5236, 5546, 5850, 6147, 6439, 6724, 7004, 7279, 7549, 7814, 8074, 8329, 8580, 8826, 9069, 9307, 9542, 9773, 10000. The entries in the Ld=2 table that is implemented in the above demo program are underlined. The table can be expanded even more if you have enough memory available, say, on a BS2SX or BS2e.
Note that the math can be simplified if Ld>6, a 258 byte table.
M = N-M-(M>>Ld*J)*(P-Q)>>K+Q+(L*1000) ' if Ld>6
This is due to considerations of integer overflow. With a larger table, the differences between neighboring values in the table become smaller, and at Ld>6 become small enough so that the double precision math is no longer necessary. Of course, the more elements there are in the table, the more accurate the result.
What about base 10 or base e? Once you have base 2, then the others can be calculated simply by multiplying the base 2 value times a constant, 0.30103 for base 10, or 0.69315 for base e. This is accomplished in the above program by using the ** operator in its role as a fractional multiplier.
The above table lookup is fast. However, the table does take up a lot of eeprom for the best accuracy. The following is an implementation of exercise 24 in chapter 1.2.2 of Donald Knuth's, "Art of Computer Programming". It computes the logarithm by successive approximation. The BASIC Stamp code for the actual calculation is only 9 or 10 lines long. It takes longer to explain it than to do it!
The first program given below shows how to calculate the manitssa. The algorithm assumes that the number to be converted lies between 1 and 2, so the logarithm (base 2) will lie between 0 and 1. Since the BASIC Stamp can't do fractional parts from 1 to 2 directly, we have to use a trick. We will use integers from 32768 to 65535 to represent the fractions from 1=32768/32768 to 2 =65536/32768. That is 32768 fractional parts. The logarithms of numbers from 1 to 2 range from lg1=0 to lg 2=1. In our system with 32768 as "unity", we have 0=0/32768 to 1=32768/32768. The algorithm below works by examining the square of the number at each step, and always adjusting it by a division by 2 so that it stays in the interval [1,2). The core of the BASIC Stamp code is highlighted in red color.
' LOGRTHM.BS2
' logarithm base 2 by calculation
' y = lg n
' for 1<=n<2 0<=y<1
' represented on the Stamp as fractions
' 32768/32768<=n<=65535/32768
' 0/32768<=y<=32767/32768
' Tracy Allen
' implements Knuth, ArtComProg ch. 1.2.2, exercise 24.
'
x var word ' input, a number between 1(32768) and 2(65536)
xf var x.bit15 ' the most significant bit of x
x2 var word ' auxiliary variable, high word of x^2
x2f var x2.bit15 ' msb of x2
lgx var word ' result, log base 2 of x
lgx0 var lgx.bit0 ' lowest bit of lgx, for bit addressing
bitk var bit ' auxiliary bit for calculation
k var nib ' index for steps of approximation
x=39457 ' example x=39457/32768 = 1.204132
debug dec x**20000,cr ' prints 12041 (decimal conversion)
lgx=0 ' initialize logarithm
for k=14 to 0 ' 15 bit result
x2=x**x ' square of x, high word
x=x*x ' low word
lgx0(k)=x2f ' this bit is 1 if x^2>=2
bitk=~x2f ' complement it for calculation
x=x2<<bitk+(bitk&xf) ' adjusted value of x
next ' next bit
debug dec lgx,32,cr, lgx**20000
' print the result
' as fraction e.g. 8781/32768
' as decimal value e.g.
' log2 x = 2680/10000=0.2680
debug dec lgx**60206,cr,dec lgx**13863
' log10 x = 8067/100000 =0.08067
' ln x = 1857/10000=0.1857
Notes:
The situation commonly encountered with the BASIC Stamp is that a logarithm is needed for an integer, N, that lies between 1 and 65535. These numbers range over 16 powers of 2. We are simply interested in lg N. (lg is shorthand for log base 2, and we can convert to any other base by means of a simple ** multiplication). The characteristic (base 2) is found by using the NCD operator. Then the algorithm above is applied to resolve the mantissa. An example may help to understand what is going on. Say you need the logarithm of 3456, taken as a raw reading from an analog to digital converter. The characteristic is 11, that is, 211=2048 is the highest power of 2 that is less than 3456. In fact, from fundamentals,
3456 = 211+x = 2048 * 2x
or to write this a different way,
lg 3456 = 11+x
characteristic=11
mantissa=x 0<=x<1, value to be determined
To find the mantissa, rewrite the equation as,
2x = 3456/2048, or x = lg 3456/2048
Now 3456/2048 is our number between 1 and 2. The numerator and denominator can both be multiplied times 16, which gives the fraction, 55296/32768. This is the condition of the algorithm just above, where 32768 represents unity. The following algorithm first calculates the characteristic by using the NCD operator, and then normalizes the denominator to 32768, and then applies the above iterative algorithm to calculate the mantissa. The core of the BASIC Stamp code is highlighted in red color.
' LOGRTHM2.BS2
' logarithm base 2 by calculation
' Tracy Allen
' implements Knuth, ArtComProg ch. 1.2.2, exercise 24.
' characteristic and mantissa in base 2, e and 10
y var word ' to hold an input number from 1 to 65535
x var word ' for processing the number
xf var x.bit15 ' high bit of x, note alias
x2 var word ' for squaring the number
x2f var x2.bit15 ' high bit of x2, note alias
lgx var word ' will be the lg (base 2) of y, the mantissa
lgx0 var lgx.bit0 ' lowest bit of lgx, for bit addressing
bitk var bit ' temporary bit
k var nib ' loop and array index
cc var nib ' characteristic of the lg
lg var word ' to hold the log base 2
log var word ' to hold the log base 10
ln var word ' to hold the log base e
pick:
debug cr,"Enter a number from 1 to 65536:"
serin 16,$54,[DEC y] ' get the number
cc=ncd y - 1 ' find the characteristic
x=y << (15-cc) ' adjust for a denominator of 32768
'debug cr,"y=",dec dcd cc," * 1.",dec4 x**20000
' optionally, show the decompostion
lgx=0 ' initialize accumulator
for k=14 to 0 ' 15 steps of precision
x2=x**x ' high byte of x squared
lgx0(k)=x2f ' high bit of x squared is this bit of log.
bitk=~x2f ' complement of that bit
x=x2<<bitk+(bitk&xf)' adjust x
next ' repeat
' lgx now holds the mantissa, log base 2 of x
' cc holds the characteristic
' show it:
debug cr,"lg x=",dec cc,".",dec4 lgx**20000,cr ' 0 to 15.9999
' combine it into one 16 bit word (but lose one digit!):
lg=cc*1000+(lgx**20000/10)
debug "lg x=",dec lg/1000,".",dec3 lg ' lg from 0 to 15.999
' convert it to log base 10:
log=lg**19728
debug "log x=",dec log/1000,".",dec3 log ' log from 0 to 4.816
' or to get one more digit resolution in log base 10:
' still 0.30103 * log2 x
' cc*10000*0.30103, in stampese is cc*10000**19728, or cc*4000**49321
' to avoid overflow in the first cc*4000, and to get best resolution in the 49321.
' First term **6021 because lgx as it comes out of the loop has an implied denominator
' of 32768. Multiply x2 to get denominator of 216 for **, with 0.30103 * 10000
log=lgx**6021+(cc*4000**49321)
debug "log x=",dec log/1000,".",dec4 log ' log from 0 to 4.8164
' convert log base 2 to log base e:
ln=lg**45426
debug "ln x=",dec ln/1000,".",dec3 ln ' ln from 0 to 11.089
pause 1000
goto pick ' get another number to try
At the end the program demo combines the characteristic and the mantissa into one 16-bit quantity for display, and converts to other bases as examples of how to manipulate the logarithm.
Sometimes the best way to condition the output of a nonlinear sensor or process is by means of a lookup table. The question usually arises of how many entries are needed in the table. For example, if a 12-bit analog to digital converter provides the data, then you might need a table with 4096 elements in order to translate each reading into its corresponding output engineering or scientific units. However, this is seldom necessary. The curve may be smooth enough that a few points may be chosen to represent the whole curve, and intermediate points are constructed on the fly by interpolation. This might be the calibration curve for a thermistor, or an altitude vs pressure curve, or the control curve for a process variable etc. The table holds the output values corresponding to a set of input values.
There are various kinds of interpolation that take into consideration more or less the local curvature of the line segments that the table represents. Linear interpolation is the simplest. It takes the pairs of table entries as end points of straight line segments. A computed output value will be an equal proportion along the output segemnt as the input value is along between the input end points. That is what I am going to talk about here.
For example, if the input value x lies 1/3 of the distance between two of the input table entries, then the output value computed will lie 1/3of the way between the corresponding two output table entries.
In the table below, the input values start at zero and increases evenly in jumps of 128. On the BASIC Stamp it is convenient to use binary intervals, to make the calculations easier. The output values also start at zero, but they do not increase in even jumps. In general, an interpolation table does not have to start at zero, nor do the values have to be evenly spaced. But let's take a closer look at this particular table and linear interpolation.
|
input, Xi |
output, Yi |
|
|
0 |
0 |
|
|
128 |
444 |
|
|
256 |
875 |
|
|
384 |
1293 |
|
|
512 |
1699 |
|
|
640 |
2095 |
|
|
668 |
2479 |
|
|
796 |
2854 |
On to Stamp code, an example. The result of an analog to digital conversion will be held in a variable, RESULT, and its value may be an integer from 0 to 4095 counts. Give the table 33 entries, and 32 intervals. The zeroth interval is from 0 to 127, the first from 128 to to 255, ... , up to the 31st 3968 to 4095. The width of each catagory is 128. The reading from the A/D converter will fit into one of 32 catagories (4096/128=32). Which catagory to use is to be found by dividing the current A/D reading by 128. The remainder after dividing is the position within the catagory, for interpolation. For example, if the A/D reading happens to be 427, that divided by 128 gives 3 with a remainder of 43. The number is in the third catagory, and it lies 43/128 of the way between the end points. The algorithm picks the corresponding y endpoint entries from the table (1293 and 1699 in this case), and then calculates the difference between them (1699-1293=406), takes 43/128 of that (136), and adds to the starting value of the catagory (1293+136=1429). That is the result returned by the subroutine.
' Interpol.BS2
' Tracy Allen
' Example of how to interpolate in a table
table data word 0,word 444,word 875,word 1293,word 1699,word 2095
data word 2479,word 2854,word 3219,word 3576,word 3923
data word 4262,word 4594,word 4919,word 5236,word 5546
data word 5850,word 6147,word 6439,word 6724,word 7004
data word 7279,word 7549,word 7814,word 8074,word 8329
data word 8580,word 8826,word 9069,word 9307,word 9542
data word 9773,word 10000.
X var word ' this will be the input variable
Z con 4096 ' the maximum value of X, 0<=X<Z
' make this an exact power of 2.
L con 32 ' the number of intervals in the table
K con Z/L ' the width of catagories of X
M con 65536/K ' for the interpolation formula
J var byte ' Jth table entry, 0<=J<=32
F var word ' the remainder of X in the catagory X=J*K+F
Q var word ' for lower bound from table read, Q<=X
Q0 var Q.byte0 ' for table lookup, array of bytes
P var word ' for upper bound from table read, P>X
I var nib ' index into table for lowbyte highbyte
Y var word ' this will be the output variable
X=427 ' for example, an input datum
J=X/K ' choose the catagory within the table
F=X//K ' remainder
for I=0 to 3 ' read values that bracket the catagory.
read J*2+I+table,Q0(I) ' Q and P are contiguous 4 bytes
next
Y=P-Q*F/K+Q ' simple interpolation, works only when P-Q and F are small.
' Y=(P-Q)**(M*F) +Q ' better alternative
' Y=(P-Q)**(F<<(17-NCD K))+Q ' another alternative, using fast shifts
debug DEC J,tab,DEC Q,tab,DEC P,tab,dec Y,CR
next
The interpolation formula needs some commentary. The basic interpolation formula is
Y = Q + ((P-Q) * F/K) ' e.g. 1293+((1699-1293)*43/128)this is the the same as the formula in the program, just rearranged.
Y=P-Q*F/K+Q
This simple formula will work directly when both P-Q and F are small, so long as their product is less than 65536. In the example, 408*43 is only 17544. It turns out that in this table, the value of F will never be greater than 127, and P-Q will never be greater than 444, so the product will never be greater than 56388. So the simple formula will be okay. But it is something you have to consider. The formula would not work for example if the P-Q value was 4444 instead of 444.
The alternative formula takes advantage of the ** operator: F/K is approximated by a new fraction that has a denominator of 65536. The numerator comes from solving the equation:
Numerator/65536 = F/K
Numerator = (65536/K)*F
The factor 65536/K may be computed in advance, and it is 512 in this example, where K=128. Note that the product of M*F can never be greater than 65536, because F is always less than K. (Another alternative formula accomplishes the task by implementing the factor (65536/K) as a binary shift. That works only when the catagories are binary widths, but it is fast and easy to transfer to a PIC. When the catagories are 128 wide, then 17- NCD K is equal to 9. A shift 9 places left is the same as multiply by 512.
Another caveat is that the value of P-Q cannot be negative, because division does not work correctly on the Stamp for negative numbers. The above program uses data values that increase as the input increases. If the successive values are going to be decreasing, or if the curve to be approximated has both increasing and decreasing segments, the formula must be adapted or generalized. Here is an interpolation formula that can handle both positive and negative values. The absolute value is separated out for treatment, and then the sign is restored.
sign var Y.bit15
Y=P-Q
Y= -sign^(abs Y**(M*F))+sign +Q ' can handle both + and -
If the table entries are bytes instead of words, reading them from the table becomes easier: Instead of the for-next loop, it is simply,
read J+table,Q
read J+1+table,P
Here is an example of a thermistor linearization, where the thermistor is part of a simple voltage divider connected to an A/D converter:
;----------- ADC input
15k |
5.12V ---/\/\---o---/\/\--- Vss
reference thermistor
This particular thermistor obeys the simplified Steinhart-Hart equation
T=(B/(ln(R)-A))-C
where A=-6.913, B=5078, C=42.24 (from the manufacturer's data sheet), and R is the sensor resistance in ohms and T is the temperature in degrees Celsius. The code produced by the ADC directly related to the thermistor resistance. (5.12 volts gives a code of 4096)
code = 4096 * R/(15000+R) or R = 15000*code/(4096-code)
This can be substituted for R in the previous equation:
T=(B/(ln(15000/(4096/code-1))-A))-C
This can then be plugged into a spreadsheet program like Excel™ to develop a table of values of temperature corresponding to different equally spaced codes.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
The final column contains the Kelvin values, times 100, which is what I will use for the table entries. Kelvin has the advantage of always being positive. Note that this program uses adapts the simple interpolation so that it can deal with a strictly decreasing function of the input code.
Here is the program:
' {$stamp bs2e}
' Tracy Allen, http://www.emesystems.com
table data word 44999,word 34482,word 32450,word 31320,word 30532,word 29923
data word 29423,word 28995,word 28617,word 28277,word 27966
data word 27675,word 27402,word 27141,word 26891,word 26647
data word 26409
X var word ' this will be the input variable
Z con 2048 ' the maximum value of X, 0<=X<Z ' make this an exact power of 2.
L con 16 ' the number of intervals in the table
K con 128 ' Z/L the width of catagories of X
M con 512 ' 65536/K for the interpolation formula
J var byte ' Jth table entry, 0<=J<=32
F var word ' the remainder of X in the catagory X=J*K+F
Q var word ' for lower bound from table read, Q<=X
Q0 var Q.byte0 ' for table lookup, array of bytes
P var word ' for upper bound from table read, P>X
I var nib ' index into table for lowbyte highbyte
Y var word ' this will be the output variable
degC var word
' Interpol.BS2e
' Tracy Allen info@emesystems.com
' Unidata Red thermistor, 15kohm termination to 5.120 volts
thermistor:
gosub ADread ' not shown, returns wx
debug dec wx,tab
x=wx max 2047
J=X/K ' choose the catagory within the table
F=X//K ' remainder
for I=0 to 3 ' read values that bracket the catagory.
read J*2+I+table,Q0(I) ' Q and P are contiguous 4 bytes
next
Y = Q - ((Q-P)**(M*F)) ' interpolate
degC = Y-27315 ' convert Kelvin*100 to degrees C *100
debug DEC J,tab,DEC Q,tab,DEC P,tab,dec Y,tab, rep "-"\degC.bit15,dec abs degC/100, ".",dec2 abs degC,CR
pause 1000 ' show values including degrees C, then pause
goto thermistorThe program displays the ADC code returned, and the endpoints of the temperature interval, and the result of the interpolation in Kelvin*100 and in degrees Celsius. Even with only the 16 catagories, the interpolation is accurate (compared to the Steinhart-Hart formula) within 0.1 degree Celsius.
<top> <index>
<home> 
< mailto:info@emesystems.com
>