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.

                                        Listing One

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 ;
 

                          Listing Two

: 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