Linear regression of samples in a circular buffer

Assumes a circular buffer of noisy samples arriving at regular time intervals. Computes the instantaneous slope of the best fit straight line, using linear regression on the rolling buffer of data. Uses only integer arithmetic. Optimized for computation on 16bit microcontroller with hardware multiplier.  The linear regression computation is simplified  considerably by subtracting out the rolling average of the buffer samples.

This computation assumes  the sampling rate is known.

Example application : variometer, a device used by glider pilots to indicate vertical climb rate. Here the Z samples are altitude samples derived from a barometric pressure sensor. The linear regression computes the slope dZ/dt, i.e. vertical climb rate.  This code is in a commercial product, the open-source design Flytepark Microvario, implemented on  a 16bit PIC24F16KA101, computing the linear regression at  26 altitude samples per second with an instruction clock of 250kHz ! (The device runs on a Lithium coincell, so power consumption is an issue).

// Linear regression of samples in a circular sample
// buffer. Uses only integer arithmetic, optimized for
// computation on 16bit microcontroller  with hardware 
// multiplier.  The linear regression computation is
// simplified considerably by subtracting out the rolling
// average of the buffer samples.
// This computation assumes the samples arrive at 
// regular intervals, and this sampling rate is known.
//   Usage : 
//   1. call lr_Init() to initialize gnLRDenominator,
//      gnNumSamples and gnSampleIndex
//   2. get first sample value and initialize gZBuffer 
//      with this value 
//   3. for each subsequent incoming sample 
//   	gZBuffer[gnSampleIndex] = lr_GetNewZSample();
//	gZAverage = lr_CalculateAverage(gZBuffer,gnNumSamples);
//	gSlope = lr_CalculateSlope(gZBuffer, gnNumSamples, gnSampleIndex, gZAverage);
//      gnSampleIndex++;
//	if (gnSampleIndex >= gnNumSamples) gnSampleIndex = 0;
//

typedef signed long    s32;

#define MAX_Z_SAMPLES   80
#define SENSOR_SAMPLES_PER_SEC	26L
#define MAX_SLOPE				2000L

#define CLAMP(x,min,max)       {if ((x) <= (min)) (x) = (min); else if ((x) >= (max)) (x) = (max);}

s32 gnLRDenominator;
int gnSampleIndex, gnNumSamples;
s32 gZBuffer[MAX_Z_SAMPLES];
s32 gZAverage;
s32 gSlope;

void lr_Init(int numSamples) {
	s32 zSample, sumT, sumT2;
	int inx;
	sumT = -(numSamples * (numSamples-1L))/2L;
	sumT2 = (numSamples * (numSamples-1L)*(2L*numSamples-1L))/6L;
	gnLRDenominator = (numSamples*sumT2) - (sumT*sumT);
	gnSampleIndex = 0;
	gnNumSamples = numSamples;
	zSample = lr_GetNewZSample(); // get a sample from the sensor
	inx = gnNumSamples;
	while (inx--) gZBuffer[inx] = zSample;  // fill the ZBuffer with first sample value
	}
 
 
s32 lr_CalculateAverage(s32* pZBuffer, int numSamples ) {
   int inx;
   s32 accumulator, average;
   inx = numSamples;
   accumulator = 0;
   while (inx--)  {
	  accumulator += pZBuffer[inx];
      }
   accumulator = (accumulator >= 0 ? accumulator +numSamples/2 : accumulator - numSamples/2); 
   average = accumulator/numSamples;  // rounded up average
   return average; 
   }

/// Linear regression of samples in buffer to calculate slope.

s32 lr_CalculateSlope(s32* pZBuffer, int numSamples, int currentSampleIndex, int zAverage)   {
   int inx,tRelative;
   s32 z, sumZT,slope;
   
   sumZT = 0;
   inx = numSamples;
   while (inx--)  {
      z = pZBuffer[inx] - zAverage;   // subtract out the average value to simplify the arithmetic
      tRelative = inx - currentSampleIndex; // time origin is the current sample in window
      if (tRelative > 0) {
         tRelative -= numSamples;
         }
      sumZT += ((s32)tRelative*z);
      }

   slope = (sumZT*(s32)(SENSOR_SAMPLES_PER_SEC*numSamples))/gnLRDenominator;
   CLAMP(slope,-MAX_SLOPE,MAX_SLOPE);
   return slope;
   }