**Problem Statement**

I am building a profiling temperature controller for a small oven that
is used for enameling and the preparation of investment casting molds.
One of the necessary details is a way to read a thermocouple that is to
indicate temperature in degrees F and be suitable for use in a control
loop. Thermocouples are only slightly nonlinear. Nevertheless, a simple
way to linearize them also works well with functions that have much greater
nonlinearity, and I present it here.

Function approximation is often done by expanding a polynomial. The polynomial can require many terms, even if the function is only modestly nonlinear, and determining the best coefficients can be time consuming. (Thermocouple polynomials are typically ninth order.) Evaluating the polynomial at run time may take too long, especially on the small, slow processors used in many embedded systems. For these systems, a good routine will execute in few cycles using fixed-point arithmetic, and have adequate accuracy and resolution for the job at hand. It may be important that the resolution be greater than the accuracy. Control systems usually need to differentiate the approximated result, and smooth differentiation requires high resolution.

Polynomial approximations require more terms as the range of the function increases. The technique described divides the function into segments small enough so that each segment is adequately characterized by a parabola. For each segment, we must calculate ax2 + bx + c. The programmer's task is to determine the proper form of x, and the values of a and b. Errors can be minimized by proper choice of the end points of the segments and the internal value at which the error becomes zero. The method I use here is devoid of subtlety; I simply use some of the leading bits - in this case, three - to define the segment, and construe the rest as a fraction between zero and one. I make the end points exact, and force the error to zero also at the midpoint of the segment. With such a "tame" curve as a thermocouple's, more sophistication gives no better results, not even for continuity of slope. To achieve smooth control in the intended application, I want to read to the nearest ¼ F. I therefore calculate temperature times four.

**Design Method**

My measurements come from a 12-bit converter, making the full range
4096 counts. This corresponds to 50 mv., given the gain of the converter
system, and represents 2250 degrees F. After dividing to identify one of
eight segments, the remainder is up to 511. The variable "x" now takes
the form [remainder/512], and "c" is evidently the temperature of the beginning
if the segment. It is fairly easy to see what "a" and "b" must be. The
temperature at the end of the segment is a + b + c, and the middle temperature
is a/4 + b/2 + c. The computing equation is x(ax + b) + c, and the multiplications
by x require normalizing divisions. Raw thermocouple data are in tables
with one degree increments that give the signal to the nearest microvolt.
I use a spreadsheet to interpolate the temperature values from the table
and to calculate the coefficients. They are then used in a Forth table.
Tables 1 and 2 are examples of the work. The first column is the count
from the converter. The second, the corresponding millivolts. The next
three columns are read from the table and entered by hand; highest temperature
not exceeding the millivolt column, millivolts at that temperature, and
millivolts one degree higher. The next columns are the interpolated temperature,
four times that, rounded to an integer (eight times for Celsius), and the
calculated coefficients. (The column marked "4T (c)" shows temperatures
at the end and midpoint of segments. These are all needed for subsequent
calculations, but only the endpoint values are coefficients.) Working code
is also shown.

**Implementation Details**

Since the error is forced to zero (within the accuracy of the coefficients)
at the ends and at the middle of each segment, the obvious places to look
for errors are the one- and three-quarter points. I have found no error
exceeding ¼ degree, the best that could be expected. Since it is
unlikely that any given thermocouple will give a reading closer than two
degrees of the reference value, the approximation is clearly better than
necessary. It might seem that four segments would be adequate. There is
good reason to retain eight, and little incentive to reduce the number.
(Naturally, it is desirable to make the number of segments a power of two.)
The computation time would be the same in either case, and only 12 cells
would be saved. However, the effect on the computation would be drastic.
The maximum x (before normalization) would double, "b" would double, and
"a" would quadruple. It would not then be possible to control round-off.
Notice the rounding step in the second line of "interpolate", adding 256
to the product: 2@ ROT * 256 + 512 / + . That keeps the error from being
one sided over the range. In order to do that, */ cannot be used, so the
multiplication must be kept in bounds. With only four segments, that couldn't
be done in single precision. The net result would be going from unnecessarily
good to unacceptably poor. Without additional tricks, there is nothing
in between. Such tricks aren't warranted here.

We could get by with four segments if the precision were limited to one degree; fine for display and adequate for proportional control, but skimpy for the derivative. However, the raw converter data could be used for that. The sensitivity of the thermocouple varies between 20 and 24 microvolts (1.6 to 2 counts) per degree over the range, a variation of ± 10 percent. In some cases, the gain variation in the derivative might actually be less objectionable than the inevitable staircasing of the linearized value. Clearly, we can't read to ¼ degree with a 12-bit converter. What we can do is, given a count, report accurately the temperature that would produce it. Table Two shows the (unimplemented) calculated coefficients to return 8x Celsius temperature directly. I have no reason to believe that its performance would be inferior.

There are a few variations that make it possible to use this two-step interpolation method with more difficult functions. Increasing the number of segments is the most obvious. There is much less to be gained by moving the end points off the true curve than with segmented (piecewise) linear interpolation, but moving the internal point of no error from the midpoint toward a region of greater curvature can sometimes halve the maximum error in the segment. That complicates the computation of the a's and b's, but not inordinately. The endpoint remains a + b + c, but the internal point of no error becomes a/n2 + b/n + c, in which n is the point in the segment where the error is to be removed, expressed as a fraction of the segment size. Segmented third- order interpolation would probably handle the most difficult cases encountered in practice.

Common practice might place "interpolate" as a `DOES>` in "fahrenheit.
Separating them allows more than one table to use the same code, provided
that the segment sizes are the same. There are two reasons why "Interpolate"
does not separate the argument into index and remainder: the address of
the table is already on the stack when it begins, complicating the stack,
and there may be a need to adjust the index before using it for tables
which do not start from zero.

**Another Example**

The rough and ready sine/cosine generator shown in Listing Two, with
the calculations in Table Three, is another example of what this interpolation
method can do. It has roughly slide-rule accuracy, enough for many purposes.
As written, the routines work for angles of any size, positive or negative.
In the application described, this generality is unnecessary. It is merely
a side effect of the 8181 AND needed for cosine to work, and the cyclic
nature of the function.

I have a two-phase incremental rotation encoder with 2048 pulses per turn on each phase. Such encoders can produce errors if they move (or vibrate) around a single transition, and my (patented) circuit to prevent that automatically provides double or quadruple resolution. It is not magic: the four states of the two phases already contain the extra information. This is not the place for hardware discussion, but I will be happy to respond privately to any who want to know how to do this with two XOR gates in front of the a counter, or with software.

Some day, this encoder may be used in a robot arm, where trigonometry would be needed to calculate the hand position. It will be convenient to get sines and cosines directly, rather than to determine the quadrant as a preliminary. I therefor generate the values over a full turn, with the count of 8192 representing 360 degrees. The values returned are 512 times actual, allowing 9 bits of precision, about three decimal digits; that is as much as can be had without modifying "interpolate". "Interpolate" also dictates segments of 512 counts spanning 22.5 degrees, so that one turn requires sixteen of them.

The computed table entries make the error zero at the center of the segment, and no significant improvement seems possible, at least as far as I have investigated The "corrected" table entries in the listing slightly increase the average error over the entire segment, but provide better continuity of slope at the peaks.

**Summary**

This is a useful (but not earthshaking) way to produce arbitrary functions
by segmented second- order interpolation. The accuracy is modest, but enough
for many applications, and can approach all that can be expected from integer
calculation. Computation time is much less than for ordinary polynomial
expansions of the same accuracy (which usually need many more terms), especially
on processors without cell-wide hardware multipliers. The examples shown
are for 16- bit systems. Of course, 32-bit systems can directly extend
the method to much higher precision, but they need more segments to attain
it. I have shown that even with 16 bits, the method provides as much accuracy
as is useful for thermometry, and accurate enough trigonometry for most
control applications.

Acknowledgments: Sorry!

The idea of making "interpolate" a separate word came from reading
Julian V. Noble's recent Finite State Machine article. I invented the rest
about 15 years ago out of necessity, for use on a 12 Mhz 8086 that would
not otherwise have been fast enough despite its built-in assembly-coded
polynomial evaluator. I need to use it again, and polished it out of pride
of craftsmanship. It has likely been done many times by many others, but
I don't know who or when. Priority is hereby ceded to all who wish to claim
it.

`CREATE fahrenheit ( -- adr ) -8 , 1108 , 128 ,`
`
-26 , 1137 , 1228 ,`
`( 4 times actual temperature ) -14 , 1083 , 2339 ,`
`
-2 , 1059 , 3408 ,`
`
14 , 1055 , 4465 ,`
`
20 , 1086 , 5534 ,`
`
26 , 1125 , 6640 ,`
`
36 , 1174 , 7791 ,`

`: >temp ( n -- 4*temperature )`
` 512 /MOD DUP 0 8 WITHIN`
` IF fahrenheit interpolate`
` ELSE ABORT" Out of range." \ Real
code will shut down.`
` THEN ;`

`\ Words for testing:`

`: mv 4096 25 */ 1+ 2/ ;`

`: c >temp 4 /mod 1 .R ASCII . EMIT 25 * . ;`

`: t mv c ;`

`: sin ( n -- 512*sine ) 8191 AND`
` 512 /MOD sine-table interpolate
;`

`: cos ( n -- 512*cosine ) 2048 + sin ;`

**
Table 1. Fahrenheit Coefficient Calculation**

`Count Millivolts T lower mV lower mV upper
T 4T(c) a b`
`----------------------------------------------------------------------`
` 0 0.000
32 0.000 0.022
32.00 128 -8 1108`
` 256 3.125
169 3.104 3.127
169.91 680`
` 512 6.250
307 6.249 6.271
307.05 1228 -26 1137`
` 768 9.375
447 9.363 9.385
447.55 1790`
`1024 12.500 584
12.484 12.505 584.76 2339
-14 1083`
`1280 15.625 719
15.622 15.646 719.13 2877`
`1536 18.750 852
18.749 18.772 852.04 3408
-2 1059`
`1792 21.875 984
21.872 21.895 984.13 3937`
`2048 25.000 1116
24.996 25.020 1116.17 4465
14 1055`
`2304 28.125 1249
28.124 28.148 1249.04 4996`
`2560 31.250 1383
31.237 31.260 1383.57 5534
20 1086`
`2816 34.375 1520
34.366 34.389 1520.39 6082`
`3072 37.500 1659
34.480 37.502 1660.00 6640
26 1125`
`3328 40.625 1802
40.619 40.640 1802.29 7209`
`3584 43.750 1947
43.734 43.756 1947.73 7791
36 1174`
`3840 46.875 2096
46.861 46.881 2096.70 8387`
`4096 50.000 2250
49.996 50.016 2250.20 9001`
` ``
`

**
Table 2. Celsius Coefficient Calculation**

`Count Millivolts T lower mV lower
mV upper T 8T(c) a b`
`---------------------------------------------------------------------`
` 0 0.000
0 0.000
0.050 0.00 0 -8 1230`
` 256 3.125
76 3.100
3.141 76.61 613`
` 512 6.250
152 6.218
6.258 152.80 1222 -32 1266`
` 768 9.375
230 9.341
9.381 230.85 1847`
`1024 12.500 307
12.498 12.539 307.05 2456 -14
1203`
`1280 15.625 381
15.594 15.636 381.74 3054`
`1536 18.750 455
18.725 18.768 455.58 3645
14 1167`
`1792 21.875 528
21.834 21.876 528.98 4232`
`2048 25.000 603
24.987 25.029 603.31 4826
30 1151`
`2304 28.125 676
28.12 28.162 676.12 5409`
`2560 31.250 750
31.214 31.256 750.86 6007
34 1199`
`2816 34.375 826
34.339 34.380 826.88 6615`
`3072 37.500 904
34.484 37.524 904.99 7240
36 1238`
`3328 40.625 983
40.605 40.645 983.50 7868`
`3584 43.750 1064
43.739 43.777 1064.29 8514 38
1305`
`3840 46.875 1147
46.873 46.910 1147.05 9176`
`4096 50.000 1232
49.998 50.024 1232.08 9857`
`
`` ``
`

**
Table 3. Sine Coefficient Calculation**`
`

`
Degrees Count Sine(c) a
b`
`
-------------------------------------`
`
0.00 0 0
-8 204`
`
11.25 256 100`
`
22.50 512 196
-20 186`
`
33.75 768 284`
`
45.00 1024 362 -34
145`
`
56.25 1280 426`
`
67.50 1536 473 -38
77`
`
78.75 1792 502`
`
90.00 2048 512 -38
-1`
`
101.25 2304 502`
`
112.50 2560 473 -34
-77`
`
123.75 2816 426`
`
135.00 3072 362 -20
-146`
`
146.25 3328 284`
`
157.50 3584 196
-8 -188`
`
168.75 3840 100`
`
180.00 4096 0
8 -204`
`
191.25 4352 -100`
`
202.50 4608 -196
20 -186`
`
213.75 4864 -284`
`
225.00 5120 -362
34 -145`
`
236.25 5376 -426`
`
247.50 5632 -473
38 -77`
`
258.75 5888 -502`
`
270.00 6144 -512
38 1`
`
281.25 6400 -502`
`
292.50 6656 -473
34 77`
`
303.75 6912 -426`
`
315.00 7168 -362
20 146`
`
326.25 7424 -284`
`
337.50 7680 -196
8 188`
`
348.75 7936 -100`
`
360.00 8192 0`

Copyright © 1998 by Jerry Avins