Fast lookup plus interpolation computation of non-linear functions

This method can be used to simplify the computation of any numeric value where the actual equation uses special functions (log, pow etc.) that require floating point computations and slow library functions. It is optimized for microcontrollers without floating point units or integer hardware divide.

It basically fits a piecewise linear curve to the original highly non-linear function, and does a linear interpolation on the appropriate linear segment to approximate the function value.  Accuracy is limited only by the number of piecewise linear segments in the region of interest, i.e. by the amount of code memory you are willing to allocate for the LUT. The computation time is independent of the number of piecewise segments.

The example code is for computing altitude above sea level from a barometric pressure reading. Pressure P(z) = Psealevel * e ^(-constant*z)  where z is the altitude. 

With a 16bit PIC24 microcontroller with hardware multiplier, observed speed up from ~150mSecs using floating point and pow() calculations, to ~3mSecs.

// Fast integer arithmetic lookup+interpolation method for converting 
// barometric pressure readings to altitude readings.

// Each Lookup Table (LUT) entry is the altitude in centimeters above
// sea level, corresponding to an implicit pressure value, 
// calculated as [PA_INIT - 1024*LUTindex] in Pascals.
// The region of interest is approximately 460metres below sea level,
// to 10000 metres above sea level.

typedef signed long s32;
 
#define PZLUT_ENTRIES   80
#define PA_INIT         106956L
#define PA_DELTA        1024L

#define Z_LIMIT_LO   -99999L
#define Z_LIMIT_HI   99999L

const s32 gPZTbl[PZLUT_ENTRIES] = {
-45853,
-37662,
-29407,
-21087,
-12700,
-4245,
4279,
12874,
21540,
30279,
...  // values removed for brevity
959708,
984147,
1009345
};

// Calculates the altitude in centimeters above sea level, given the barometric
// sensor pressure reading in pascals. The nearest lower LUT index is computed.
// The altitude is then linearly interpolated from the corresponding altitude
// values at the lower and next higher LUT index. Computation is optimized by
// ensuring the difference between LUT entries are spaced by a power of 2, in
// this case 2^10 (1024), so no integer division is required.
// Returns the error values Z_LIMIT_LO or Z_LIMIT_HI if
// the pressure data exceeds the LUT index limits.

s32 sns_Pa2Cm(s32 pa)  {
   	s32 inx,pa1,z1,z2,z;
    
   	if (pa > PA_INIT) {  
      	z = Z_LIMIT_LO;  
      	}
   	else {
      	inx = (PA_INIT - pa)>>10;      
      	if (inx >= PZLUT_ENTRIES-1) {
         	z = Z_LIMIT_HI;
         	}
      	else {
         	pa1 = PA_INIT - (inx<<10);
         	z1 = gPZTbl[inx];
         	z2 = gPZTbl[inx+1];
         	z = z1 + (((pa1-pa)*(z2-z1))>>10);
         	}
      	}
   	return z;
   	}