Quadratic Data Fit

AndreyITO · May 5, 2019, 12:20am

Hey guys,

I am quite new at Arduino but already having a lot of fun and advancing a bit.

Anyways... for my project I want to compare several enviroments in terms of light. For that I have a digital light sensor attached to my Arduino which returns serial data in lux.

Well, in order to compare different light enviroments I want to fit the light over time functions to a quadratic curve:

(light in lux) = ax²+bx+c

With x being from 0 to 20.

Why quadratic function? Well I have special light sources and prior experiments without Arduino showed that quadratic fit is quite suitable. Now I want to build a handheld device.

My problem: I found following code here in the forum.

// USAGE:  double g = fit_G( sizeof(x)/sizeof(double), t, x);

// should be declared like that to optimize for the 8bits microcontroller core.
double fit_G(int N_points, double px[], double py[]) __attribute__((__optimize__("O2"))); //Remove this one if you're using DUE.
double fit_G( int N_points, double px[], double py[] ) {

  int i;
  double S00, S10, S20, S30, S40, S01, S11, S21;
  double denom, x, y, a, b, c;
  S00=S10=S20=S30=S40=S01=S11=S21=0;
  
  for (i=0; i<N_points; i++) {  
    x = px[i];
    y = py[i];
    //S00 += 1; // x^0+y^0
    S10 += x;
    S20 += x * x;
    S30 += x * x * x;
    S40 += x * x * x * x;
    S01 += y;
    S11 += x * y;
    S21 += x * x * y;
  }
  S00 = N_points;
  
  denom =   S00*(S20*S40 - S30*S30) - S10*(S10*S40 - S20*S30) + S20*(S10*S30 - S20*S20);

/* c = (S01*(S20*S40-S30*S30)-S11*(S10*S40-S20*S30)+S21*(S10*S30-S20*S20))/denom;
   b = (S00*(S11*S40-S30*S21)-S10*(S01*S40-S21*S20)+S20*(S01*S30-S11*S20))/denom;*/
   a = (  S00*(S20*S21 - S11*S30) - S10*(S10*S21 - S01*S30) + S20*(S10*S11 - S01*S20) )/denom;
            
  double g = a*2;
  return g;
}

And I have no idea how to implement that.

I even understand the mathematics behind it, but I have no clue how I should tell my Arduino to declare x as 0 to 20 and y as my serial data and then fit the curve.

Aim: I want to receive different a,b,c for different enviroments and compare them.

Thank you all already for the answers!

I'm excited to hear them.

jremington · May 5, 2019, 1:21am

To use that function, your main Arduino program needs to fill two float arrays, one with X values and the other with Y values. Then you call the function with those two arrays as arguments.

One way to create the data arrays is by declaring them and initializing them in the program, for example

float x[] = {0.0, 1.0, 2.0, 3.0};
float y[] = {3.1, 4.5, 5.5, 6.7);
...

float g = fit_G( 4, x, y);

Do you realize that the function you have selected does not return the a, b and c coefficients of the quadratic fit, and does not calculate b and c? The calculation of b and c has been changed into a comment using the /* ... */ symbols.

AndreyITO · May 5, 2019, 9:35am

jremington:
To use that function, your main Arduino program needs to fill two float arrays, one with X values and the other with Y values. Then you call the function with those two arrays as arguments.

One way to create the data arrays is by declaring them and initializing them in the program, for example
float x[] = {0.0, 1.0, 2.0, 3.0};

float y[] = {3.1, 4.5, 5.5, 6.7);
...

float g = fit_G( 4, x, y);




Do you realize that the function you have selected [u]does not[/u] return the a, b and c coefficients of the quadratic fit, and does not calculate b and c? The calculation of b and c has been changed into a comment using the /* ... */ symbols.

Thank you for your answer!

Ok, it's a bit clearer now! But there's the problem that whereas my x values I can define for sure, my y values are different in every loop.

Do I define them in the function? Something like

float y [] = digitalRead();

?

And yes, I realize it returns only g=a*2, but that is due to the reason that the original thread wanted to calculate specific acceleration. If I rewrite it to

return a
return b
return c

and uncomment b and c it would work, right?

Also: In which section of the code do I put what?
Which part goes into the setup? And what into loop?

wvmarle · May 5, 2019, 10:05am

You can pass only one value through the return statement (only the first will be executed, after that you're out of the function).

The others have to be done through the function parameters.

AndreyITO · May 5, 2019, 12:46pm

wvmarle:
You can pass only one value through the return statement (only the first will be executed, after that you're out of the function).

The others have to be done through the function parameters.

Ok, I got that... So I'll need to work something out to return all the coefficients.

But! I still do not really have a clue how to implement all this into a function.
However, I tried something out, obviously it's not quite working... but still, here is my unready code...
(TSL2561 is the luminosity sensor.)

// First try of own fit

#include <Wire.h>
#include <TSL2561.h>

/*
   Setup TSL2561
*/
TSL2561 tsl(TSL2561_ADDR_FLOAT);

// Quadratic fit function:
// USAGE:  double g = fit_G( sizeof(x)/sizeof(double), t, x);

float x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 
             6.0, 7.0, 8.0, 9.0, 10.0, 
             11.0, 12.0, 13.0, 14.0, 15.0,
             16.0, 17.0, 18.0, 19.0, 20.0};
int numBytes // is that right?


// should be declared like that to optimize for the 8bits microcontroller core.
double fit_G(int N_points, double px[], double py[]) __attribute__((__optimize__("O2"))); //Remove this one if you're using DUE.

double fit_G( int N_points, double px[], double py[] ) {

  int i;
  double S00, S10, S20, S30, S40, S01, S11, S21;
  double denom, x, y, a, b, c;
  S00=S10=S20=S30=S40=S01=S11=S21=0;
  
  for (i=0; i<N_points; i++) {  
    x = px[i];
    y = py[i];
    //S00 += 1; // x^0+y^0
    S10 += x;
    S20 += x * x;
    S30 += x * x * x;
    S40 += x * x * x * x;
    S01 += y;
    S11 += x * y;
    S21 += x * x * y;
  }
  S00 = N_points;
  
  denom =   S00*(S20*S40 - S30*S30) - S10*(S10*S40 - S20*S30) + S20*(S10*S30 - S20*S20);

   c = (S01*(S20*S40-S30*S30)-S11*(S10*S40-S20*S30)+S21*(S10*S30-S20*S20))/denom;
   b = (S00*(S11*S40-S30*S21)-S10*(S01*S40-S21*S20)+S20*(S01*S30-S11*S20))/denom;
   a = (  S00*(S20*S21 - S11*S30) - S10*(S10*S21 - S01*S30) + S20*(S10*S11 - S01*S20) )/denom;

// Here I'll need to work yomething out to return a, b and c...
  double g = a*2;
  return g;
}

void setup() {

  // Initialize the serial connection for debugging.
  Serial.begin(9600);
  Wire.begin();
  Serial.println(F("Serial initialized"));
  
  // Initialize the TSL2561.
  if (tsl.begin()) {
    Serial.println("TSL2561 Confirmed");
  } else {
    Serial.println("Error initializing TSL2561");
  }
  // Setting Gain and Integration Time
  tsl.setGain(TSL2561_GAIN_16X);
  tsl.setTiming(TSL2561_INTEGRATIONTIME_402MS);
}

// Main function
void loop() {
  for (int i = 0; i <= 20; i++) { // Measure light 20 times
    // TSL2561 Parameters
    uint16_t y = tsl.getLuminosity(TSL2561_VISIBLE);
    Serial.println(y, DEC);
    delay(150);
  }

  numBytes = Serial.available();
  for (n = 0; n < numBytes; n++) {
    mySerialData[n] = Serial.read();
 float y[] = mySerialData[];  // Is that right?

 float g = fit_G( 20, x, y);  
}

Can someone help?

jremington · May 5, 2019, 2:28pm

To return multiple values from a function, you need another mechanism. A common one is to pass to a function the addresses of (pointers to) variables as follows.

Be careful to observe usage of "*" and "&".

float x[] = {0.0, 1.0, 2.0, 3.0};
float y[] = {3.1, 4.5, 5.5, 6.7);
float a=0.0, b=0.0, c=0.0;
...

fit_G( 4, x, y, &a, &b, &c);

...

double fit_G( int N_points, double px[], double py[], float* a, float* b, float* c ) {

...
   *c = (S01*(S20*S40-S30*S30)-S11*(S10*S40-S20*S30)+S21*(S10*S30-S20*S20))/denom;
   *b = (S00*(S11*S40-S30*S21)-S10*(S01*S40-S21*S20)+S20*(S01*S30-S11*S20))/denom;
   *a = (  S00*(S20*S21 - S11*S30) - S10*(S10*S21 - S01*S30) + S20*(S10*S11 - S01*S20) )/denom;

// no return statement needed

AndreyITO · May 6, 2019, 10:14am

jremington:
To return multiple values from a function, you need another mechanism. A common one is to pass to a function the addresses of (pointers to) variables as follows.

Be careful to observe usage of "*" and "&".
float x[] = {0.0, 1.0, 2.0, 3.0};

float y[] = {3.1, 4.5, 5.5, 6.7);
float a=0.0, b=0.0, c=0.0;
...

fit_G( 4, x, y, &a, &b, &c);

...

double fit_G( int N_points, double px[], double py[], float* a, float* b, float* c ) {

...
c = (S01(S20S40-S30S30)-S11*(S10S40-S20S30)+S21*(S10S30-S20S20))/denom;
b = (S00(S11S40-S30S21)-S10*(S01S40-S21S20)+S20*(S01S30-S11S20))/denom;
a = ( S00(S20S21 - S11S30) - S10*(S10S21 - S01S30) + S20*(S10S11 - S01S20) )/denom;

// no return statement needed

Ok, wow! That's already amazing! Thank you a lot!

My last question, I think, would be how I tell my programm to take the y values from the serial data.
Because in your code y is defined as an array with 3.1, 4.5,...

wvmarle · May 6, 2019, 10:53am

Read the values from the serial input, place them in an array, when complete call that function to do the calculation.

AndreyITO · May 7, 2019, 8:14am

Ok, I understood nearly the whole concept. But I think I'm just not advanced enough...

Something in my declarations is not right, so I get the error message

"no matching function for call to 'fit_G(int, float [21], float [1], float*, float*, float*)' "

Why is my second array only 1 bit long?
Why is there an error?

Full code:

// First try of own fit

#include <Wire.h>
#include <TSL2561.h>

/*
   Setup TSL2561
*/
TSL2561 tsl(TSL2561_ADDR_FLOAT);

// Quadratic fit function:
// USAGE:  double g = fit_G( sizeof(x)/sizeof(double), t, x);

float x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 
             6.0, 7.0, 8.0, 9.0, 10.0, 
             11.0, 12.0, 13.0, 14.0, 15.0,
             16.0, 17.0, 18.0, 19.0, 20.0};

float y[] = {0};
float a=0.0, b=0.0, c=0.0;

float array[] = {0};


// should be declared like that to optimize for the 8bits microcontroller core.
double fit_G(int N_points, double px[], double py[]) __attribute__((__optimize__("O2"))); //Remove this one if you're using DUE.

double fit_G( int N_points, double px[], double py[], float* a, float* b, float* c ) {

  int i;
  double S00, S10, S20, S30, S40, S01, S11, S21;
  double denom, x, y;
  S00=S10=S20=S30=S40=S01=S11=S21=0;
  
  for (i=0; i<N_points; i++) {  
    x = px[i];
    y = py[i];
    //S00 += 1; // x^0+y^0
    S10 += x;
    S20 += x * x;
    S30 += x * x * x;
    S40 += x * x * x * x;
    S01 += y;
    S11 += x * y;
    S21 += x * x * y;
  }
  S00 = N_points;
  
  denom =   S00*(S20*S40 - S30*S30) - S10*(S10*S40 - S20*S30) + S20*(S10*S30 - S20*S20);

   *c = (S01*(S20*S40-S30*S30)-S11*(S10*S40-S20*S30)+S21*(S10*S30-S20*S20))/denom;
   *b = (S00*(S11*S40-S30*S21)-S10*(S01*S40-S21*S20)+S20*(S01*S30-S11*S20))/denom;
   *a = (  S00*(S20*S21 - S11*S30) - S10*(S10*S21 - S01*S30) + S20*(S10*S11 - S01*S20) )/denom;

}

void setup() {

  // Initialize the serial connection for debugging.
  Serial.begin(9600);
  Wire.begin();
  Serial.println(F("Serial initialized"));
  
  // Initialize the TSL2561.
  if (tsl.begin()) {
    Serial.println("TSL2561 Confirmed");
  } else {
    Serial.println("Error initializing TSL2561");
  }
  // Setting Gain and Integration Time
  tsl.setGain(TSL2561_GAIN_16X);
  tsl.setTiming(TSL2561_INTEGRATIONTIME_402MS);
}

// Main function
void loop() {
  for (int i = 0; i <= 20; i++) { // Measure light 20 times
    // TSL2561 Parameters
    uint16_t y = tsl.getLuminosity(TSL2561_VISIBLE);
    Serial.println(y, DEC);
    delay(150);
  }
  if (Serial.available() > 0)
{
    for (int i=0; i <= 20 ; i++)
    {
        array[i] = Serial.read();
    }
    Serial.flush();
}  
  fit_G( 20, x, array, &a, &b, &c);
  }

wvmarle · May 7, 2019, 10:57am

For starters you're filling your array[] with a lot of byte values, read from the Serial console. That is almost certainly not what you want.

Take a look at the Serial input basics tutorial.

system · May 7, 2019, 5:14pm

  if (Serial.available() > 0)
{
    for (int i=0; i <= 20 ; i++)
    {
        array[i] = Serial.read();
    }

If there is one byte to read, it is NOT OK to then read all 21 of them.

When you do not define a value in the [] when declaring an array, the compiler counts the initializers for you. You provided one initializer, so you have a one element array. You can't store 21 values in a one element array.

AndreyITO · May 7, 2019, 6:36pm

Apologies for so many questions, this is my first real Arduino Project

Let me quickly summarize and divide the problem:

My aim is something like that:

declare all necessary variables for the least square calculation and data fit (done?)
read the serial data from the sensor

...
float array[20];
...
if (Serial.available() > 0)
{
    for (int i=0; i <= 20 ; i++)
    {
        array[i] = Serial.read();
    }
   
}
...

Yes, I declared the array size wrong, after further research I found this mistake and now you confirmed my suspicion in the answers. But is this approach wrong? It is what I found in a few different sources. Or do I need to declare all the additional functions stated in the Serial Basics Tutorial? (recvWithStartEndMarkers() etc.) Isn't that for parsing?

call that array into my fit function to be the y of my graph

Here too, I dont really have an idea how to accomplish that since I get an error message no matter what I try...

jremington · May 7, 2019, 7:02pm

A recommended approach is to get each part of a project running and understand how each one works, before combining them.

For example, for data input from the serial monitor, start with this tutorial (examples 4 and 5 for numeric input) and once they are working and you understand the programs, modify example 5 to allow entry of an array of numeric values.

It will be fairly simple to add the quadratic fit function to the working, modified serial input example 5.

wvmarle · May 8, 2019, 1:05am

Data over the Serial interface comes one byte at a time, not 21 in one go. A float value is at least 4 bytes, if represented as ASCII (human readable) data it may be 6-10 bytes per value.

You think you may send lots of bytes in one go on the Serial console, but they arrive one by one. At 9600 bps it's about 1 ms per byte, in that time the loop() may run several times. So it sees bytes come in one at a time.

To properly read the Serial input, you first of all have to know how the data is sent exactly. What format. Does it have start/end markers? Is it binary perhaps? Figure that one out and you can try to read the values. Most likely you have to create a buffer (array of char or byte) in which you store the incoming data, then when a value is complete you parse it to recreate the original float value.

So start with that: read your Serial input, recreate the float value, print it out. When you see those values appear correctly printed back on your console you're a good way in.

wvmarle · May 8, 2019, 1:09am

But writing that I more and more realise:

you never told what the source of the values is - what sensor produces them, and how the data looks like (as well be explained in the data sheet of this mystery sensor).
your code (which is incomplete - the suspiciously named TSL2561_ADDR_FLOAT in your code of #8 has not been declared) only mentions a TSL2561 light sensor. That one works over I2C, not Serial. The library will take care of turning the received values into something useful. I forgot if it's a float or int.

AndreyITO · May 10, 2019, 2:02pm

Thank you for all your input!

I think I got it now....

At least if the fit itself works.

Here is the code:

.....

void loop() {
 for (int i = 0; i <= 20; i++) {// Measure light 20 times
  uint16_t lum = tsl.getLuminosity(TSL2561_VISIBLE);
  y[i] = double(lum);

  delay(250);
 }
  for(int j = 0; j <= 20; j++)
  {
  Serial.println(y[j]);
  }
  delay(250);

  fit_G( 20, x, y, &a, &b, &c);

  Serial.println("a: ");
  Serial.println(a);
  Serial.println("b: ");
  Serial.println(b);
  Serial.println("c: ");
  Serial.println(c);

...
  
}

Comments?

jremington · May 10, 2019, 2:24pm

Could you elaborate on the reason for reading the sensor 20 times, as quickly as possible, and fitting just the first four (or all 20) points of the data with a quadratic function?

What "x" values are you using, and what do they represent?

AndreyITO · May 10, 2019, 2:53pm

jremington:
Could you elaborate on the reason for reading the sensor 20 times, as quickly as possible, and fitting just the first four (or all 20) points of the data with a quadratic function?

What "x" values are you using, and what do they represent?

About the fit:

It was a typo. I changed it so that I fit all the 20 values.

About the general idea:

It's a project for my university where I have to analyze different light situations and I wanted to automatize it.

Topic		Replies	Views
[SOLVED] Fitting a quadratic curve with Arduino. Possible? General Guidance	6	4099	May 5, 2021
Fitting a complex curve Programming	19	2598	May 5, 2021
Exponential curve fitting on the Teensy3.6 General Guidance	7	1701	May 5, 2021
Least Squares Approximation Programming	18	17991	May 5, 2021
Not passing array correctly to function Programming	6	4494	May 5, 2021

Quadratic Data Fit

Related topics