Floating Point Inaccuracy in Polynomial Regression on Atmega328

Go To Last Post
15 posts / 0 new
Author
Message
#1
  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Hello!! I'm trying to implement a polynomic regression on an Atmega to calculate coefficients for a line, a quadratic and a qubic polynom; all this from a set of values obtained from the ADC. I am using the Least Squares method to calculate the coefficients and I have already found the formulas for the three equation, I have also written the code to calculate the coefficients and it works fine on the computer (I have already checked doing the regression on Excel). However once I program the Microcontroller and get the coefficients using the same set of data and code, they are very different from the ones I get using the computer.

 

I've read about the error that one gets using floating point math and I'm trying to avoid it. But it's very difficult since the values tend to get too big too fast and go above the limit of 32 bit signed integers. So I was wondering if there's a way to avoid using floating point math or a way of using it so that it does't affect so much the acuraccy of the results.

 

To show what I'm referring to, I add an snippet of the code I'm using. The x array would be the values I'd get from the ADC, and the values on y array won't go higher than 200. As you see the variables s1 to s7, are the sums of the products of the arrays. Therefore they tend to surpass the limit of the 32bit integer. And the floating point math is only  really needed when the coefficients a1 to a3 are calculated

 

float x[16]={448.0,407.0,374.0,342.0,321.0,310.0,289.0,267.0,258.0,241.0,226.0,212.0,198.0,179.0,163.0,134.0};
float y[16]={10.0,15.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0,55.0,60.0,65.0,70.0,75.0,80.0,85.0};

uint_16 n=16;

s1 = s2 = s3 = s4 = s5 = s6 = s7 = s = stot = 0;
   for (i=0; i<n; i++)
   {  x = px[i];
      y = py[i];
      s1 += x;
      s2 += x * x;
      s3 += x * x * x;
      s4 += x * x * x * x;
      s5 += y;
      s6 += x * y;
      s7 += x * x * y;
   }
   if (denom = n  * (s2 * s4 - s3 * s3) -
               s1 * (s1 * s4 - s2 * s3) +
               s2 * (s1 * s3 - s2 * s2) )
   {  a1 = (s5 * (s2 * s4 - s3 * s3) -
            s6 * (s1 * s4 - s2 * s3) +
            s7 * (s1 * s3 - s2 * s2)) / denom;
      a2 = (n  * (s6 * s4 - s3 * s7) -
            s1 * (s5 * s4 - s7 * s2) +
            s2 * (s5 * s3 - s6 * s2)) / denom;
      a3 = (n  * (s2 * s7 - s6 * s3) -
            s1 * (s1 * s7 - s5 * s3) +
            s2 * (s1 * s6 - s5 * s2)) / denom;
            

 

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Your simplest solution might be 64-bit integers.

IIRC avr-gcc handles them rather slowly.

Another possibility is to do the math and rearrange your formulas to avoid things

like the difference of two nearly equal quantities when one or both is not exact.

Yet another possibility might be 48-bit signed bitfields.

Don't know how avr-gcc would handle them.

IIRC avr-gcc has 24-bit integers, but no 48-bit types.

You might emulate 48-bit integers "by hand", say as 3 16-bit integers.

If you got that route, try to make the compiler do as much of the work as you can.

Note, for example, that adding a sequence of hand-made integers can be done

by adding all the components separately and combining them at the end.

"Demons after money.
Whatever happened to the still beating heart of a virgin?
No one has any standards anymore." -- Giles

Last Edited: Fri. Nov 11, 2016 - 01:49 AM
  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

I add an snippet of the code I'm using.

Can you post a big enough snippet that will compile and produce different results on a PC and an AVR?

 

Obviously going from a 32bit integer number (that uses most of the bits) to a 32bit floating point number loses a lot of precision (a float only has 24 bits of mantissa.)  Probably with an IEEE 64bit float, yu don't lose any of those bits (though you can lose some elsewhere.)

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0
if (denom = n  * (s2 * s4 - s3 * s3) - ...

Don't you want '==' instead of '='?

 

 

Edit: Please post real code, that can be compiled.

Last Edited: Thu. Nov 10, 2016 - 09:00 AM
  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0
if (denom = n  * (s2 * s4 - s3 * s3) - ...

Don't you want '==' instead of '='?

 

I think it's checking for zero divisor.

I'd much prefer to make it explicit, at least:

if (0 != (denom = n  * (s2 * s4 - s3 * s3) - ... )
  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Hello! Thank you all for your answers.

 

skeeve: I'm still chewing on some of the options that you gave me. I'll try to do some of them and will reply you as soon as I have.

 

 westfw and snigelen:

 

I'm attaching a c code that I used to get the quadratic coefficients. It is a Least Squares Analysis, which takes a function and from the data that it produces, tries to obtain the linear, quadratic, n exponential, logarithmic, and power coefficients to compare them and determine which adjust better. You can use it with the functions it comes with, or you can add a set of points. I didn't write this code. In fact, it seems to have a Copyright, or so the author claims on the header. I got it from another forum and prove that the math on the quadratic and linear part was ok. You can compile it and run it in a computer and if you want on a microcontroller too (with some adjustments). It's possible that the quadratic part works fine on a uC for some set of points, but it differs from the result I want with the data I have used. In the case of the cubic polynom the difference tends to grow bigger. I understand that as you said (westfw) a floating point 32bits loses a lot of precision. But for avr-gcc the double and float type of variables are both of 32bits.

 

Thank you again, I'll be posting my advances as soon as I can

Attachment(s): 

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

I am wondering what the resolution of your ADC might be...

There is no point in using math that has more precision than the resolution of the input values (except for some rare instances).

After all, if the ADC value changes by ±1 bit, how much does that affect the final calculation?

Only if the difference is less than the ~7 digit precision of a float should the accuracy of the math come into play.

David (aka frog_jr)

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

bananafish wrote:
However once I program the Microcontroller and get the coefficients using the same set of data and code, they are very different from the ones I get using the computer.

So, what values >>do<< you get on the computer?  On the computer, are you using "float"?  How many bits are in a "float" (single) on your platform?  What values do you get from Excel?  What values do you get from the "microcontroller"?  What toolchain are you using?

 

 

You can put lipstick on a pig, but it is still a pig.

I've never met a pig I didn't like, as long as you have some salt and pepper.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

The posted code is a bit puzzling.  (The function does all the analyses?)  Anyway, I see no declaration for s1, s2, ... ?

 

Perhaps post your complete PC and AVR test program.

You can put lipstick on a pig, but it is still a pig.

I've never met a pig I didn't like, as long as you have some salt and pepper.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

If the OP said "why" he is going through all this math gymnastics, a better way to solve the actual problem could be suggested! 

 

 

Jim

 

Click Link: Get Free Stock: Retire early! PM for strategy

share.robinhood.com/jamesc3274
get $5 free gold/silver https://www.onegold.com/join/713...

 

 

 

 

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

It is still interesting to me that (in the provided code):

   double c1 = exp(1), c2 = 3.1415927, c3 = sqrt(2);

Although declared a double, pi is only defined with the precision of a float. Therefore, any of the equations using that value (c2) are only good to (at best) ~7 digits.

The results on a PC may provide more digits in the result, however, they are still only valid to ~7 decimal places.

 

Further, in the data provided in post #1, there appear to be only 3 to 4 digits available, so answers may seem to be extremely precise, but they may not be nearly so accurate.

 

So I ask again, what is the resolution of the input data?

David (aka frog_jr)

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

How much of a difference are you seeing?  (and on what data?)  I modified your fabls.c (attached) so it will run on my AVR Arduinos, and I see pretty good agreement (not perfect) with the sample data that's in there when compiled on the Mac...  If you're complaining about the same level of "error", that's one thing, but if you're somehow seeing much larger errors, that could be something else!

 

Mac: Linear:      y = (2.718282) x + 3.141593; s = 0.000000
AVR: Linear:      y = (2.7183) x + 3.1416; s = 0.0000
Mac: Quadratic:   y = (0.000000) x^2 + (2.718282) x + 3.141593; s = 0.000000
AVR: Quadratic:   y = (0.00) x^2 + (2.7182) x + 3.1416; s = 0.0001
Mac: Exponential: y = exp(0.258940 x + 1.583046); s = 0.063365
AVR: Exponential: y = exp(0.2589x +1.5830); s = 0.0634
Mac: Logarithmic: y = (6.582515) ln(x) + 4.993691; s = 0.986435
AVR: Logarithmic: y = (6.5825) ln(x) +4.9937; s = 0.9864
Mac: Power:       y = (5.683863) x ^ (0.649854); s = 0.371944
AVR: Power:       y = (5.6839) x ^ (0.6499); s = 0.3719
The best fit is a linear function

Mac: Linear:      y = (19.451284) x - 17.613759; s = 5.085440
AVR: Linear:      y = (19.4513) x - 17.6138; s = 5.0854
Mac: Quadratic:   y = (2.718282) x^2 + (3.141593) x + 1.414214; s = 0.000000
AVR: Quadratic:   y = (2.72) x^2 + (3.1417) x + 1.4139; s = 0.0017
Mac: Exponential: y = exp(0.604827 x + 1.578474); s = 0.169417
AVR: Exponential: y = exp(0.6048x +1.5785); s = 0.1694
Mac: Logarithmic: y = (45.322808) ln(x) - 2.656422; s = 11.921653
AVR: Logarithmic: y = (45.3228) ln(x) -2.6564; s = 11.9217
Mac: Power:       y = (6.903959) x ^ (1.525705); s = 2.544455
AVR: Power:       y = (6.9040) x ^ (1.5257); s = 2.5445
The best fit is a quadratic function

Mac: Linear:      y = (3820408.378552) x - 7500931.583125; s = 5458477.959685
AVR: Linear:      y = (3820404.5000) x - 7500923.0000; s = 5458473.5000
Mac: Quadratic:   y = (2543110.201169) x^2 - (11438252.828462) x + 10300839.825058; s = 2675639.811868
AVR: Quadratic:   y = (2543104.00) x^2 - (11438220.0000) x + 10300831.0000; s = 2675637.2500
Mac: Exponential: y = exp(2.718282 x + 3.141593); s = 0.000000
AVR: Exponential: y = exp(2.7183x +3.1416); s = 0.0000
Mac: Logarithmic: y = (7793587.435647) ln(x) - 3502053.546431; s = 6461623.615646
AVR: Logarithmic: y = (7793580.5000) ln(x) -3502049.5000; s = 6461618.0000
Mac: Power:       y = (147.479711) x ^ (6.582515); s = 6305801.674243
AVR: Power:       y = (147.4798) x ^ (6.5825); s = 6305786.0000
The best fit is an exponential function

Mac: Linear:      y = (1.063398) x + 2.554149; s = 0.396479
AVR: Linear:      y = (1.0634) x + 2.5541; s = 0.3965
Mac: Quadratic:   y = (-0.205384) x^2 + (2.295700) x + 1.116463; s = 0.097761
AVR: Quadratic:   y = (-0.21) x^2 + (2.2957) x + 1.1164; s = 0.0978
Mac: Exponential: y = exp(0.206312 x + 1.085508); s = 0.119828
AVR: Exponential: y = exp(0.2063x +1.0855); s = 0.1198
Mac: Logarithmic: y = (2.718282) ln(x) + 3.141593; s = 0.000000
AVR: Logarithmic: y = (2.7183) ln(x) +3.1416; s = 0.0000
Mac: Power:       y = (3.269923) x ^ (0.542746); s = 0.236027
AVR: Power:       y = (3.2699) x ^ (0.5427); s = 0.2360
The best fit is a logarithmic function

Mac: Linear:      y = (103.577513) x - 160.551786; s = 60.855602
AVR: Linear:      y = (103.5775) x - 160.5518; s = 60.8556
Mac: Quadratic:   y = (32.268206) x^2 - (90.031722) x + 65.325654; s = 7.685993
AVR: Quadratic:   y = (32.27) x^2 - (90.0305) x + 65.3257; s = 7.6860
Mac: Exponential: y = exp(1.228998 x + 0.321075); s = 0.458222
AVR: Exponential: y = exp(1.2290x +0.3211); s = 0.4582
Mac: Logarithmic: y = (230.362017) ln(x) - 70.390499; s = 95.352013
AVR: Logarithmic: y = (230.3619) ln(x) -70.3904; s = 95.3520
Mac: Power:       y = (2.718282) x ^ (3.141593); s = 0.000000
AVR: Power:       y = (2.7183) x ^ (3.1416); s = 0.0000
The best fit is a power function

Attachment(s): 

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

westfw wrote:
If you're complaining about the same level of "error", that's one thing, but if you're somehow seeing much larger errors, that could be something else!

 

Op said:

bananafish wrote:
However once I program the Microcontroller and get the coefficients using the same set of data and code, they are very different from the ones I get using the computer.

So I [indirectly] asked what "very" means:

theusch wrote:
So, what values >>do<< you get on the computer? On the computer, are you using "float"? How many bits are in a "float" (single) on your platform? What values do you get from Excel? What values do you get from the "microcontroller"? What toolchain are you using?

 

And I don't think any of my questions were answered.

You can put lipstick on a pig, but it is still a pig.

I've never met a pig I didn't like, as long as you have some salt and pepper.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Hello again!!

 

westfw:

How much of a difference are you seeing?  (and on what data?)

Try adding the cubic expression, that would be:

s1 = s2 = s3 = s4 = s5 = s6 = s7 = s8 = s9 = s10 = s = stot = 0;
   for (i=0; i<n; i++)
   {  x = px[i];
      y = py[i];
      s1 += x;
      s2 += x * x;
      s3 += x * x * x;
      s4 += x * x * x * x;
      s5 += x * x * x * x * x;
      s6 += x * x * x * x * x * x;
      s7 += y;
      s8 += x * y;
      s9 += x * x * y;
      s10 += x * x * x * y;
   }

if (denom = (s6 * (n * (s2 * s4 - s3 * s3) + s1 * (s3 * s2 - s1 * s4) + s2 * (s1 * s3 - s2 * s2)) -
                s5 * (n * (s5 * s2 - s4 * s3) + s1 * (s3 * s3 - s1 * s5) + s2 * (s4 * s1 - s2 * s3)) +
                s4 * (n * (s5 * s3 - s4 * s4) + s1 * (s3 * s4 - s2 * s5) + s2 * (s2 * s4 - s3 * s3)) -
			    s3 * (s1 * (s5 * s3 - s4 * s4) + s2 * (s3 * s4 - s2 * s5) + s3 * (s2 * s4 - s3 * s3))))

   {  a3 = (s10 * (n * (s2 * s4 - s3 * s3) + s1 * (s3 * s2 - s1 * s4) + s2 * (s1 * s3 - s2 * s2)) -
            s9 *  (n * (s5 * s2 - s4 * s3) + s3 * (s3 * s1 - s2 * s2) + s1 * (s4 * s2 - s1 * s5)) +
            s8 *  (n * (s5 * s3 - s4 * s4) + s1 * (s3 * s4 - s2 * s5) + s2 * (s2 * s4 - s3 * s3)) -
		    s7 *  (s1 * (s5 * s3 - s4 * s4) + s2 * (s3 * s4 - s2 * s5) + s3 * (s2 * s4 - s3 * s3))) / denom;

      a2 = (s6 * (n * (s2 * s9 - s3 * s8) + s1 * (s2 * s8 - s1 * s9) + s7 * (s1 * s3 - s2 * s2)) -
            s5 * (n * (s2 * s10 - s4 * s8) + s1 * (s3 * s8 - s1 * s10) + s7 * (s1 * s4 - s2 * s3)) +
            s4 * (n * (s3 * s10 - s4 * s9) + s1 * (s3 * s9 - s2 * s10) + s7 * (s2 * s4 - s3 * s3)) -
		    s3 * (s1 * (s10 * s3 - s4 * s9) + s2 * (s3 * s9 - s2 * s10) + s8 * (s2 * s4 - s3 * s3))) / denom;

      a1 = (s6 * (n * (s4 * s8 - s3 * s9) + s7 * (s2 * s3 - s1 * s4) + s2 * (s1 * s9 - s2 * s8)) -
            s5 * (n * (s5 * s8 - s3 * s10) + s7 * (s3 * s3 - s1 * s5) + s2 * (s1 * s10 - s3 * s8)) +
            s4 * (n * (s5 * s9 - s4 * s10) + s7 * (s3 * s4 - s2 * s5) + s2 * (s2 * s10 - s3 * s9)) -
		    s3 * (s1 * (s5 * s9 - s4 * s10) + s8 * (s3 * s4 - s2 * s5) + s3 * (s2 * s10 - s3 * s9))) / denom;

	  a0 = (s6 * (s7 * (s2 * s4 - s3 * s3) + s1 * (s9 * s3 - s8 * s4) + s2 * (s3 * s8 - s2 * s9)) -
            s5 * (s7 * (s2 * s5 - s3 * s4) + s1 * (s3 * s10 - s8 * s5) + s2 * (s4 * s8 - s2 * s10)) +
            s4 * (s7 * (s3 * s5 - s4 * s4) + s1 * (s10 * s4 - s9 * s5) + s2 * (s4 * s9 - s3 * s10)) -
		    s3 * (s8 * (s5 * s3 - s4 * s4) + s2 * (s10 * s4 - s9 * s5) + s3 * (s4 * s9 - s3 * s10))) / denom;

      for (i=0; i<n; i++)
      {  x = px[i];
         dy = py[i] - (a3 * x * x * x  + a2 * x * x + a1 * x + a0);
		 dy2 = py[i] - s7/n;
         s += dy * dy;
		 stot += dy2 * dy2;
      }
      // s = sqrt(s / r);
	  s = 1 - (s/stot);
      sign  = (a1 < 0) ? '-' : '+';
      sign2 = (a2 < 0) ? '-' : '+';
	  sign3 = (a0 < 0) ? '-' : '+';
      printf("3rd Degree:   y = (%f) x^3 %c (%f) x^2 %c (%f) x %c %f; r^2 = %f\n",
             a3,sign2,fabs(a2),sign,fabs(a1),sign3,fabs(a0),s);

I'm using this set of data: 

 

double x[16]={448.0,407.0,374.0,342.0,321.0,310.0,289.0,267.0,258.0,241.0,226.0,212.0,198.0,179.0,163.0,134.0};
double y[16]={10.0,15.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0,55.0,60.0,65.0,70.0,75.0,80.0,85.0};

And getting the next results on pc and Atmega 328p:

 

PC:

Quadratic:   y = (0.000357) x^2 - (0.468406) x + 146.065785; r^2 = 0.992688
cubic:   y = (0.000002) x^3 - (0.001781) x^2 + (0.111254) x + 97.518867; r^2 = 0.998062

uC:

Quadratic: y = 0.000357 x^2 -0.468266 x +146.042083; r^2 = 0.992687
Cubic: y = 0.000002 x^3 -0.001130 x^2 +0.117809 x +70.221809; r^2 = 0.233806

For the Quadratic polynom the difference is not too big, but for the Cubic it is enough to not to be of use at all.

 

(I also replaced the s, for the r^2) 

s = 1 - (s/stot);

 

frogjr:

 

I am wondering what the resolution of your ADC might be

It is a 10 bit ADC. But I'm planning to implement it on a xmega with a 12 bit resoultion, but I'd use it with a 11 bit because I'd use a differential input with gain. So the input x would be in a range of

[0 - 2048]

 

Jim:

If the OP said "why" he is going through all this math gymnastics, a better way to solve the actual problem could be suggested! 

I'm trying to characterize (I'm not sure if thats the correct term in english) a sensor. I mean, to obtain the characteristic curve of response of a sensor, which I know from observation and data, that is not linear. I've already done a lookup table and linear interpolation between points, which I've found is the common way of doing that on a uC. I just wanted to obtain an equation which I thought would be more practical, and perhaps accurate.

 

theush:

So, what values >>do<< you get on the computer? On the computer, are you using "float"? How many bits are in a "float" (single) on your platform? What values do you get from Excel? What values do you get from the "microcontroller"? What toolchain are you using?

 

Maybe your questions have been already answered in the code and results I put above. Except for the toolchain that I'm using which is Codevision. And the results that Excel gives, which I'm not putting, because I'm lazy ;) But trust me they are very much alike to what I'm getting from the pc, and different from what the uC gives me. I also tried the same data on wathever web polynomial regression tool, google gives when asked, with the same results.

Last Edited: Fri. Nov 11, 2016 - 02:37 PM
  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Sorry but did you just say you are comparing the mathematical accuracy of the 32bit IEEE754 code in CodeVision to that of Excel? What exactly were you expecting? 32bit IEEE is not some kind of magic with infinite precision. A 23 bit (well really 24 bit because of the normalisation) mantissa is accurate for about 7.5 decimal places.

 

Observe, this is using gcc on a powerful PC:

$ cat fptest.c
#include <stdio.h>

int main(void) {
    volatile float f = 0.0f;
    int n;

    for (n = 0; n < 10000; n++) {
        f = f + 0.1f;
    }
    printf("result = %f\n", (double)f);
    return 0;
}
$ gcc -Os fptest.c -o fptest
$ ./fptest 
result = 999.902893

It adds 0.1 10,000 times. In a perfect world the result would have been 1000.0. 32bit IEEE754 does not make a perfect world.

 

Now I could switch to 64 bit IEE754...

$ cat fptest.c
#include <stdio.h>

int main(void) {
    volatile double d = 0.0;
    int n;

    for (n = 0; n < 10000; n++) {
        d = d + 0.1;
    }
    printf("result = %f\n", d);
    return 0;
}
$ gcc -Os fptest.c -o fptest
$ ./fptest 
result = 1000.000000

But on an AVR you don't get that luxury unless you pay $3,000 for IAR or $800 (is it?) for the "Pro" version of Imagecraft.

 

You need to find a different numeric representation system with more accuracy. "float" on an AVR (or a powerful PC or anything for that matter where it is 32 bit) is not going to cut this.