Mathematical functions (i.e. interpolation) ?

For my project I need to interpolate or resample a vector. This is needed to plot an array (physical quantity vs time) on a LCD screen which has less dots than the number of samples I have available. I'm thinking of some spline or running average, but this is not easy (if not impossible at all, at least for me!) to do this through simple algebric functions available in C. Do you know if a mathematical library exists to perform such complex functions? I found math.h, but even this does only quite simple operations (trigonometry, max and min, rounding, logs...). Maybe a library exists oriented to data or image processing... Thanks a lot.

I have not looked at this in detail but it might do what you want http://playground.arduino.cc/Main/Smoothstep

You can do it with simple algebra. Imagine that the y values are areas, h*w, in "bins" between consecutive limits, x0 x1 x2 and so on. Define new x values, r0 r1, r2,... reflecting new sample areas. Now, calculate new y values for each r. Each r "bin" will include 2 of the x bins. However you can easily calculate the areas and simple add them. The summed areas will be the new y values for your vector.

aarg: You can do it with simple algebra. Imagine that the y values are areas, h*w, in "bins" between consecutive limits, x0 x1 x2 and so on. Define new x values, r0 r1, r2,... reflecting new sample areas. Now, calculate new y values for each r. Each r "bin" will include 2 of the x bins. However you can easily calculate the areas and simple add them. The summed areas will be the new y values for your vector.

This is very smart :wink: , actually I didn't think to areas (i.e. integrals) but I thought to average over bins. However, it works (as simple as you said) as far as the number of samples is multiple of the number of bins I need (typically 2 times, but may be done with factor 3 or more). Not so simple if the two are not multiple e.g. 50 samples in 30 bins...

tmmlrd: This is very smart :wink: , actually I didn't think to areas (i.e. integrals) but I thought to average over bins. However, it works (as simple as you said) as far as the number of samples is multiple of the number of bins I need (typically 2 times, but may be done with factor 3 or more). Not so simple if the two are not multiple e.g. 50 samples in 30 bins...

It is not at all necessary for the bin widths to be integer multiples with this method. Actually, x and r could be real numbers and the calculations are the same, and no more complex. I've used it for audio resampling. You just have to calculate and add two areas for each r.

The question is unclear... If you have less output points than input points, this has nothing to do with interpolation.

As a starting point, you could simply try remapping. Remapping is used for windowing and viewporting on all computers since the first GUI. See Map() : http://www.arduino.cc/en/Reference/Map ; this function does exactly what you need : conversion of one coordinate system into another. You could have to convert floats into long integers (floating point into fixed point)

But it also depends on the source array size : it could be more efficient to first resample, and then remap.

Using polynomial functions is not easy : for an accurate representation of n points, you have to calculate a n-1 degree polynom that has all the input points as solutions. (2 points -> degree 1, 3 points -> cubic, and so on).

Splines (particular polynomial functions) are generally used as an approximation, and you have to calculate control points. You also have to to decide a criteria for accuracy : mean square deviation for example.

So :

what are yout data and how much ? Xmin, Xmax, Ymin, Ymax, N what is your viewport : width * height ?

I am pretty sure the answer to your problem is simple.

tmmlrd: For my project I need to interpolate or resample a vector. This is needed to plot an array (physical quantity vs time) on a LCD screen which has less dots than the number of samples I have available.

There are algorithms for reducing the number of points in a polygon curve.

Such like the Ramer–Douglas–Peucker algorithm which is explained by Wikipedia, including some "pseudocode" which perhaps can be translated to C/C++ easily.

Would one of the spline algorithms work?

http://www.codeproject.com/Articles/560163/Csharp-Cubic-Spline-Interpolation

AppCrash: If you have less output points than input points, this has nothing to do with interpolation. [...] I am pretty sure the answer to your problem is simple.

Right. This is not interpolation. KISS principle.

Thanks all !!! Luckily it is holiday here on June 1st and 2nd, I need time to evaluate your suggestions!

P.S. You are right, it is no interpolation. I think the right name is resampling.

It is resampling if you reduce the number of points. But interpolation can be needed :

If the samples count is a multiple of the viewport width, or the amount of data is much larger than the viewport x axis pixels count, it is easy : you drop points, but doing this you lose information, particularly transients.

If the samples count is not a multiple of the viewport width, or the viewport width is not far smaller than the samples count, it is more complex, and then interpolation is probably needed. Again, you lose information. Worse, you create information. It can result in glitches and artefacts.

If the data represent a periodic function, and you have to resample and then interpolate, a commonly used algorithm is "sin(x)/x interpolation". It is used by digital oscilloscopes firmwares and hardwares, but it is very complex, and creates funny glitches in certain circumstances.

You could read this document : http://cdn.teledynelecroy.com/files/whitepapers/wp_interpolation_102203.pdf

Remapping and drawing pixels (or segments) keeps all information, even transients. Even with large amounts of data, it can be faster (depending on your display) as it does not need floating math (or need little). Arduinos are not good at floating math, as math are not hardware math. We forgot this with our computers since mid '90 when microprocessors all cames with an FPU ! Before the 486, we had to buy an expensive 8087, 80287 or 80387 to get a hardware floating point unit.

With Arduinos, interpolation and regression must be avoided when possible.

Another point to be discussed is the amount of memory needed for calculations. You can easily be short in memory with the Atmega328

(sorry for my bad english...)

We're talking about a low-resolution display here, not studio-quality audio. There are many, many very straight-forward pure integer algorithms that will easily do the job, whether you're up-sampling or down-sampling. Google is your friend. When you're generating what is undoubtedly 8-bit output values, you don't need extreme precision. Depending on the data, even simple linear interpolation may well be adequate.

Regards, Ray L.

RayLivingston: We're talking about a low-resolution display here, not studio-quality audio. There are many, many very straight-forward pure integer algorithms that will easily do the job, whether you're up-sampling or down-sampling. Google is your friend. When you're generating what is undoubtedly 8-bit output values, you don't need extreme precision. Depending on the data, even simple linear interpolation may well be adequate.

Regards, Ray L.

I agree (although all what has been said before remains true).

@tmmlrd

Can you give an example of the vector? how many elements does it have? to what number it should be reduced?

I recall from image processing that reducing arrays is easy and fast if length B is a whole divider of length A. e.g and array of 100 ==> array of 25.

However if the lengths are not whole dividers there is "trouble" e.g your original vector is 7 elements and you want to reduce to 3. You always end up with an input bin that needs to be divided over 2 bins. Simplest is to use linear interpolation.

example

10 12 14 18 35 23 7 ==> [10+12+ 1/3*14] [2/3*14 + 18 + 2/3*35] [1/3*35 + 23 + 7] 27 50 42 problem is now that the peaks are all higher than in the original number so you need to multiply * 3/7 12 21 18

If you have memory enough there is a trick by using an intermediate array

make an array with length C = length A * length B (7x3 in this case) and repeat the values length B times.

10 12 14 18 35 23 7 ==> 10 10 10 12 12 12 14 14 14 18 18 18 35 35 35 23 23 23 7 7 7

now reduce it by taking the average of every 7 bins (length A) 80/7= 11 152/7=22 125/7=18

note this can be optimized without the intermediate array. (left as an exercise )

another option is to fill also the "display vector" simultaneously.

robtillaart: @tmmlrd

Can you give an example of the vector? how many elements does it have? to what number it should be reduced?

I recall from image processing that reducing arrays is easy and fast if length B is a whole divider of length A. e.g and array of 100 ==> array of 25.

However if the lengths are not whole dividers there is "trouble" e.g your original vector is 7 elements and you want to reduce to 3. You always end up with an input bin that needs to be divided over 2 bins. Simplest is to use linear interpolation.

example

10 12 14 18 35 23 7 ==> [10+12+ 1/3*14] [2/3*14 + 18 + 2/3*35] [1/3*35 + 23 + 7] 27 50 42 problem is now that the peaks are all higher than in the original number so you need to multiply * 3/7 12 21 18

If you have memory enough there is a trick by using an intermediate array

make an array with length C = length A * length B (7x3 in this case) and repeat the values length B times.

10 12 14 18 35 23 7 ==> 10 10 10 12 12 12 14 14 14 18 18 18 35 35 35 23 23 23 7 7 7

now reduce it by taking the average of every 7 bins (length A) 80/7= 11 152/7=22 125/7=18

note this can be optimized without the intermediate array. (left as an exercise )

OMG: this is applied mathematics!!! :) So many years at university studying analysis and didn't know this trick! :D It sounds very simple and effective.

For my application I'm planning to take one measurement every 30 sec and have to display the data vector collected over, say, 12 hours or so on a 64 pixel display. I know I'll throw most of the information there, but would be nice to have a graphical indication of the history during the past hours beside the instantaneous value of the measured quantity.

P.S. I'm not sure I know how to estimate the amount of memory needed for the long vector, I can deal with integers and could have, say 1440 readings times 64 pixels of the display i.e. 92160 numbers. If an integer takes 16 bit (is this correct for arduino nano?) the needed amount should be 184 kbyte. Is this correct? Is it compatible with the arduino nano internal memory? I'm afraid not... Mmmhhh... more tricky than I thought, mathematics is just the starting point...

For that, simple low-pass filtering, followed by averaging will work perfectly. Anything more elaborate will be a waste of time.

Regards, Ray L.

There is still at least one missing information. What do you want to display ? Peak values (for example min and/or max? mean value ?)

The solutions are simple.

The "x cursor" of your display will be incremented by 1 every 12hour/64 ; (64 = display width)

min/max : putPixel(x, y) : x see above ; y : all the values (scaled to display height) ; every 12hour/64 you increment x by 1 ; you will get a graph with a max cruve, a min curve, and pixels in between (a pixel "cloud") ; you can also retain only min and max, and only display these two curves

mean : you accumulate values during 12hour/64, count them, calculate mean value, and putPixel(x, y) when time is elapsed ; x will be equal to 0, y will be equal to the (scaled) mean value you calculated. Then, you increment y by 1, you reset the accumulator, and you continue

Extremely simple. No math. Just plotting + scaling. You curve will grow by itself.

Think of an oscilloscope in free run mode with a very low basetime !

[EDIT] when x reaches 64, you set it to 0, clear display, and the plotting will restart from the left

putPixel(x,y) or "lineTo(x,y)". putPixel is fine if you collect many samples within 12hour/64. Otherwise, use "lineTo" if only one value is collected.

You can also do nice things : for example, you can can count the occurences of each and every value, and put a pixel more or less bright depending on its number of occurences. You get what is called "Z axis" on digital ocilloscopes. The more a value is displayed, the more it will shine.

Could you post an exemple of your data (a large amount) and the duration ?

Did some testing / coding to rebin an array. Tested with MEGA 2560 + IDE 1.5.8.

The code iterates over all the input elements and sums them. It detects if it should add an input element partially. It does not use an intermediate array, only a sum so it is quite memory efficient. Furthermore the code rounds and scales the resulting elements. The performance of the code depends on the sizes of both arrays. Some sample runs:

600 elements to 12 new bins in 2252 uSec => 3.75 uSec per input element.
600 elements to 43 new bins in 3600 uSec => 6.00 uSec per input element.
600 elements to 500 new bins in 23484 uSec => 39.14 usec per input element.

Note if the ‘from’ array is an whole multiple of the ‘to’ array the code can be simplified and will be 5~10% faster.

//
//    FILE: rebin.ino
//  AUTHOR: Rob Tillaart
// VERSION: 0.1.00
// PURPOSE: rebinning demo
//    DATE: 2015-05-30
//     URL: http://forum.arduino.cc/index.php?topic=326320.msg2254353#msg2254353
//
// Released to the public domain
//

uint32_t start;
uint32_t stop;

#define LENX 600
#define LENY 43 

int x[LENX];
int y[LENY];

void setup()
{
  Serial.begin(115200);
  Serial.println("rebin.ino");

  for (int i = 0; i < LENX; i++) x[i] = random(100);
  //  for (int i = 0; i < LENX; i++)
  //  {
  //    Serial.print(x[i]);
  //    Serial.print(' ');
  //  }
  //  Serial.println();

  Serial.println("REBIN 2");
  start = micros();
  rebin2(x, LENX, y, LENY);
  stop = micros();
  Serial.println(stop - start);
  for (int i = 0; i < LENY; i++)
  {
    Serial.print(y[i]);
    Serial.print(' ');
  }
  Serial.println();
  Serial.println("done...");

}

void loop()
{
}

// preferred - stable & fastest overal
// handle every input only once.
void rebin2(int *A, uint16_t lengthA, int *B, uint16_t lengthB)
{
  long sum = 0;
  uint16_t i = 0;
  uint16_t ii = 0;
  uint16_t j = 0;
  uint16_t t = lengthA * lengthB;
  while (t)
  {
    t -= lengthB;
    ii += lengthB;
    if (ii < lengthA) sum += A[i++] * lengthB;
    else
    {
      ii -= lengthA;
      B[j++] = (sum + A[i] * (lengthB - ii) + lengthA / 2) / lengthA;  // scaling & rounding
      sum = A[i] * ii;
      i++;
    }
  }
}
rebin.ino
REBIN 2
3600
49 43 47 56 37 46 45 49 60 46 54 39 54 48 57 48 57 47 38 42 54 43 53 56 35 42 49 35 51 51 50 40 51 38 45 49 49 39 39 52 41 44 48 
done...

... For my application I'm planning to take one measurement every 30 sec and have to display the data vector collected over, say, 12 hours or so on a 64 pixel display. I know I'll throw most of the information there, but would be nice to have a graphical indication of the history during the past hours beside the instantaneous value of the measured quantity.

P.S. I'm not sure I know how to estimate the amount of memory needed for the long vector, I can deal with integers and could have, say 1440 readings times 64 pixels of the display i.e. 92160 numbers. If an integer takes 16 bit (is this correct for arduino nano?) the needed amount should be 184 kbyte. Is this correct? Is it compatible with the arduino nano internal memory? I'm afraid not... Mmmhhh... more tricky than I thought, mathematics is just the starting point...

Create the display vector while making the measurements and you do not need the intermediate array. the framework of your code could look like:

int displayVector[64];  // used as circular buffer


void loop()
{
  raw = makeMeasurement(system);
  
  cooked = doMath(raw);

  measurement[index] = cooked;
  displayVector[dispIndex] += cooked;    // or min, max, ...

  index++;
  if (index % 1440 == 0) // take next display index every 1440 elements
  {
    dispIndex++;    
    if (dispIndex == 64) dispIndex = 0;  // wrap around to have a circular buffer
    displayVector[dispIndex] = 0;  // reset the next element
  }

  if (keyPressed())
  {
     updateDisplay(displayVector, dispIndex);
  }
}