K Thermocouple conversion by polynomial conversion formula

[Vindhyachal.takniki]

Was writing code to convert thermocouple k type conversion formula for best accuracy. Goal is to use polynomial functions and not the table or other method, that I will write in next code. Currently accuracy is concern. Below is my code, any one found any issue in that? 
Currently i have simulated only by putting different millivolt values & cold junction value as 0c.  Any way it can be improved further for accuracy or any other error in code.

Code: [Select]

/*
1. Convert the cold-junction temperature (TCJ) to a voltage (VCJ)
2. Add the cold-junction voltage to the measured thermocouple voltage (VCJ + VTC)
3. Convert the summed cold-junction voltage and thermocouple voltage to the thermocouple temperature
(TTC)
 */


#define float32_t  float


const float32_t ci_i0_below0 = 0.000000000000E+00;
const float32_t ci_i1_below0 = 0.394501280250E-01;
const float32_t ci_i2_below0 = 0.236223735980E-04;
const float32_t ci_i3_below0 = -0.328589067840E-06;
const float32_t ci_i4_below0 = -0.499048287770E-08;
const float32_t ci_i5_below0 = -0.675090591730E-10;
const float32_t ci_i6_below0 = -0.574103274280E-12;
const float32_t ci_i7_below0 = -0.310888728940E-14;
const float32_t ci_i8_below0 = -0.104516093650E-16;
const float32_t ci_i9_below0 = -0.198892668780E-19;
const float32_t ci_i10_below0 = -0.163226974860E-22;

const float32_t ci_i0_above0 = -0.176004136860E-01;
const float32_t ci_i1_above0 = 0.389212049750E-01;
const float32_t ci_i2_above0 = 0.185587700320E-04;
const float32_t ci_i3_above0 = -0.994575928740E-07;
const float32_t ci_i4_above0 = 0.318409457190E-09;
const float32_t ci_i5_above0 = -0.560728448890E-12;
const float32_t ci_i6_above0 = 0.560750590590E-15;
const float32_t ci_i7_above0 = -0.320207200030E-18;
const float32_t ci_i8_above0 = 0.971511471520E-22;
const float32_t ci_i9_above0 = -0.121047212750E-25;

const float32_t a0_above0 = 0.118597600000E+00;
const float32_t a1_above0 = -0.118343200000E-03;
const float32_t a2_above0 = 0.126968600000E+03;



const float32_t d0_i0 = 0.0000000E+00;
const float32_t d0_i1 = 2.5173462E+01;
const float32_t d0_i2 = -1.1662878E+00;
const float32_t d0_i3 = -1.0833638E+00;
const float32_t d0_i4 = -8.9773540E-01;
const float32_t d0_i5 = -3.7342377E-01;
const float32_t d0_i6 = -8.6632643E-02;
const float32_t d0_i7 = -1.0450598E-02;
const float32_t d0_i8 = -5.1920577E-04;
const float32_t d0_i9 = 0.0000000E+00;

const float32_t d1_i0 = 0.000000E+00;
const float32_t d1_i1 = 2.508355E+01;
const float32_t d1_i2 = 7.860106E-02;
const float32_t d1_i3 = -2.503131E-01;
const float32_t d1_i4 = 8.315270E-02;
const float32_t d1_i5 = -1.228034E-02;
const float32_t d1_i6 = 9.804036E-04;
const float32_t d1_i7 = -4.413030E-05;
const float32_t d1_i8 = 1.057734E-06;
const float32_t d1_i9 = -1.052755E-08;

const float32_t d2_i0 = -1.318058E+02;
const float32_t d2_i1 = 4.830222E+01;
const float32_t d2_i2 = -1.646031E+00;
const float32_t d2_i3 = 5.464731E-02;
const float32_t d2_i4 = -9.650715E-04;
const float32_t d2_i5 = 8.802193E-06;
const float32_t d2_i6 = -3.110810E-08;
const float32_t d2_i7 = 0.000000E+00;
const float32_t d2_i8 = 0.000000E+00;
const float32_t d2_i9 = 0.000000E+00;



void setup()
{
  Serial.begin(9600);

}

void loop()
{
    uint32_t err;
    float32_t tempc;

    err = thermocouple_tempc(54.819f, 0.0f , &tempc);

    Serial.print(err);
    Serial.print("     ");
    Serial.println(tempc);

}





uint32_t thermocouple_tempc(float32_t vtc_mV , float32_t cj_tempc , float32_t *tc_tempc_return)
{
    uint32_t err = 0U;
    float32_t vcj_mV = 0.0f;
    float32_t vsum = 0.0f;
    float32_t tempc = 0.0f;
   

//////////////// STEP-A cold-junction temperature (TCJ) to a voltage (VCJ)   
/* convert cj_temp into mv by below equation
*
* E = sum(i=0 to n) c_i t^i.
*
* The equation above 0 °C is of the form
* E = sum(i=0 to n) c_i t^i + a0 exp(a1 (t - a2)^2).
*
*/   
    if(cj_tempc < -270.0f)
    {
        err = 1U;  //cold-junction undertemperatyre
    }
    else if((cj_tempc >= -270.0f) && (cj_tempc < 0.0f))
    {
        vcj_mV =         ( ci_i0_below0  * pow(cj_tempc , 0U) )
                       + ( ci_i1_below0  * pow(cj_tempc , 1U) )
                       + ( ci_i2_below0  * pow(cj_tempc , 2U) )
                       + ( ci_i3_below0  * pow(cj_tempc , 3U) )
                       + ( ci_i4_below0  * pow(cj_tempc , 4U) )
                       + ( ci_i5_below0  * pow(cj_tempc , 5U) )
                       + ( ci_i6_below0  * pow(cj_tempc , 6U) )
                       + ( ci_i7_below0  * pow(cj_tempc , 7U) )
                       + ( ci_i8_below0  * pow(cj_tempc , 8U) )
                       + ( ci_i9_below0  * pow(cj_tempc , 9U) )
                       + ( ci_i10_below0  * pow(cj_tempc , 10U) );
    }   
    else if((cj_tempc >= 0.0f) && (cj_tempc <= 1372.0f))
    {
        vcj_mV =         ( ci_i0_above0  * pow(cj_tempc , 0U) )
                       + ( ci_i1_above0  * pow(cj_tempc , 1U) )
                       + ( ci_i2_above0  * pow(cj_tempc , 2U) )
                       + ( ci_i3_above0  * pow(cj_tempc , 3U) )
                       + ( ci_i4_above0  * pow(cj_tempc , 4U) )
                       + ( ci_i5_above0  * pow(cj_tempc , 5U) )
                       + ( ci_i6_above0  * pow(cj_tempc , 6U) )
                       + ( ci_i7_above0  * pow(cj_tempc , 7U) )
                       + ( ci_i8_above0  * pow(cj_tempc , 8U) )
                       + ( ci_i9_above0  * pow(cj_tempc , 9U) )
                       + ( a0_above0 * ( exp(a1_above0 * pow((cj_tempc - a2_above0), 2U) ) )   );
                       
    } 
    else if(cj_tempc > 1372.0f)
    {
        err = 2U;    //cold-junction overtemperatyre
    }
    else
    {
        err = 3U;    //NAN or unknown error
    }
   
   
//////////////// STEP-B cold-junction voltage to the measured thermocouple voltage (VCJ + VTC)   
    if(0U == err)
    {
        vsum = vcj_mV + vtc_mV;
        if(vsum < -5.891f)
        {
            err = 4U;   //sum voltage falls below convertable range
        }
        else if((vsum >= -5.891f) && (vsum <= 54.886f))
        {
            //err = 0U; this is conevertiable range
        }
        else if(vsum > 54.886f)
        {
            err = 5U;   //sum voltage above convertable range
        }
        else
        {
            err = 6U;  //NAN or unknown error
        }
   
    }
   
   
   
//////////////// STEP-C vsum coversion to temperature
    if(0U == err)
    {
        if(vsum < -5.891f)
        {
            err = 7U;  //under-voltage for final ocnversion
        }
        else if((vsum >= -5.891f) && (vsum < 0.0f))
        {
            tempc =              ( d0_i0  * pow(vsum , 0U) )
                               + ( d0_i1  * pow(vsum , 1U) )
                               + ( d0_i2  * pow(vsum , 2U) )
                               + ( d0_i3  * pow(vsum , 3U) )
                               + ( d0_i4  * pow(vsum , 4U) )
                               + ( d0_i5  * pow(vsum , 5U) )
                               + ( d0_i6  * pow(vsum , 6U) )
                               + ( d0_i7  * pow(vsum , 7U) )
                               + ( d0_i8  * pow(vsum , 8U) )
                               + ( d0_i9  * pow(vsum , 9U) );
        }
        else if((vsum >= 0.0f) && (vsum < 20.644f))
        {
            tempc =              ( d1_i0  * pow(vsum , 0U) )
                               + ( d1_i1  * pow(vsum , 1U) )
                               + ( d1_i2  * pow(vsum , 2U) )
                               + ( d1_i3  * pow(vsum , 3U) )
                               + ( d1_i4  * pow(vsum , 4U) )
                               + ( d1_i5  * pow(vsum , 5U) )
                               + ( d1_i6  * pow(vsum , 6U) )
                               + ( d1_i7  * pow(vsum , 7U) )
                               + ( d1_i8  * pow(vsum , 8U) )
                               + ( d1_i9  * pow(vsum , 9U) );
        }
        else if((vsum >= 20.644f) && (vsum <= 54.886f))
        {
            tempc =              ( d2_i0  * pow(vsum , 0U) )
                               + ( d2_i1  * pow(vsum , 1U) )
                               + ( d2_i2  * pow(vsum , 2U) )
                               + ( d2_i3  * pow(vsum , 3U) )
                               + ( d2_i4  * pow(vsum , 4U) )
                               + ( d2_i5  * pow(vsum , 5U) )
                               + ( d2_i6  * pow(vsum , 6U) )
                               + ( d2_i7  * pow(vsum , 7U) )
                               + ( d2_i8  * pow(vsum , 8U) )
                               + ( d2_i9  * pow(vsum , 9U) );
        }
        else if(vsum > 54.886f)
        {
            err = 8U;   //over-voltage for final ocnversion
        }
        else
        {
            err = 9U;  //NAN or unknown error
        }
    }
   
   
    if(0U == err)
    {
        if(tempc < -200.0f)
        {
            err = 10U;  //final calcualtions result in under temperature
        }
        else if(tempc > 1372.0f)
        {
            err = 11U;  //final calcualtions result in over temperature
        }
        else
        {
            *tc_tempc_return = tempc;
            //err = 0U; //final calcualtions ok
        }
 
    }
   
    return err;
   
}




[radiolistener]

polynomial of order higher than 3 works bad for real world curve interpolation, especially for temperature or pressure sensors.

For industrial platinum and copper sensors, there are special standard approximation formulas.

If you don't have such standard for your sensor, it will be much better to use table with spline interpolation between table points. More large table leads to more precise interpolation

But since you cannot get ideal sensor, for most of the cases 2-nd order polynomial is good enough for temperature sensor.

Also, if you're deal with 10-th order polynomial float32 may have not enough precision. So, there is a sense to try float64.

I have some experience with industrial platinum and copper temperature sensors, but unfortunately I didn't worked with K-type, so I cannot suggest you something regarding K-type sensors...

[beanflying]

Seems you have followed the nist data (without checking to closely) https://srdata.nist.gov/its90/type_k/kcoefficients.html so you are about as close as you can be without using a table lookup.

Given the general accuracy or lack of it with Type K thermocouples your code is most likely better than the device providing or reading the data.

What are you trying to achieve with your code?

[mikerj]

Currently accuracy is concern.

What are the limitations of your cold junction measurement and the signal processing for the thermocouple?  It's pointless trying to use a 10th order polynomial if other errors dominate.  I'd also agree with radiolistener, the rounding errors you will get from single precision floats will make this pretty useless.

[Rerouter]

That is a lot of floating point math... May we know why your opposed to the table approach, over gaps like 2C the difference from straight line interpolation is below a Type K's error budget (you will be struggling to pull out more than 2.5 decimal places with most setups)

Usual table approach is to take the NIST table, re-scale it to make the most of a 16 bit integers size, and then with the extra resolution of the table do a linear interpolation between the 2 values,

E.g. -270C = -6.458mV to 1373C = 54.886mV
so you have to represent 61.344mV, and guess what, a 16 bit integer has 65536 steps, so with a simple offset you can store any part of the NIST table to the nanovolt,

Now a 3.2KB table for every degree C step would be a bit overkill, so you would probably trim down the table to your working range, as well as use 2C steps or larger to really cut down on the data size if that is the issue.

While it is ugly code (could flip it so the mV is that table lookup for almost no processing cost) you then find the nearest 2, get the relative position from the 2 and use that for your fractional temperature,

Edit: In the past I only cared about a 327C range, so using the mV table I mentioned storing instead temperatures in the 16 bit values, letting me store to the nearest 0.005C, after offsetting and scaling the raw ADC value, I divided and modulo'd the value, the division was used as the table lookup, and the modulo was used for the interpolation between the steps, leaving me with a similarly scaled 16 bit value that represented 0.005C resolution

[Vindhyachal.takniki]

The code is only for getting best mathematical accuracy & demonstration by using polynomial function.

Later i will shift code to look up table. But I want to try code on polynomial also, just to use & experiment



So wanted to know if the method used is correctly used or somewhere i have done any mistake like converting cold junction temperature to voltage or any other.



This is just for fun/hobby. Made a loop on arduinio Leonardo to see what is max time it take to convert with fixed cold junction temperature 0C & varying  thermocouple voltage from -5.891 to 54.886 with 0.01 step, all done in code only.

Max time to convert I get is 3776us even with this float intensive operation

Code is O3 optimized.

[mikerj]

Are you using the usual thermocouple arrangement i.e. a single welded junction connected to a terminal block that forms the reference/cold junction?  If so then the cold junction temperature does not determine which set of coefficients you should be using, and the cold junction temperature wouldn't cover such a huge range (it will usually be close to ambient). 

The cold junction temperature will need to be measured using some kind of absolute temperature sensor (e.g. RTD, thermistor, bandgap) which is why I asked about the accuracy you expect from this.

[ejeffrey]

I agree with a 10th order polynomial being useless and also suggest using float64 for any high order polynomial interpolation.

However there is one more thing about polynomial fitting.  You should rescale you x axis variable to be small, typically between -1 and +1 to make the most of your numerical precision.  The higher order fit the more important this is.  Remember, Tc^10 is Enormous when Tc =50 and you are adding that to the linear term.

[Kleinstein]

The POW() function in the pure form is very inefficient and may add some numerical error. So it may be better to rewrite the way the polynomnals are calculated. An efficient way is calculating a polynominal is like
y = a_0 + x*(a_1 + x* (a_2 + ..... )).

The compiler may optimize this, but not sure.

There is rarely need for such high degree polynomials. TC are usually not that accurate anyway.

[Dave]

Flash is cheap, you can just put the whole ITS-90 table in there and use it as a look-up table. The whole thing (with single-precision values) will take about 6,4kB flash. I've done exactly this in a multichannel thermometer that I've developed at work.

Doing high order polynomial approximation will be slow as wet wick on a microcontroller and you'll need to do conversions in both directions to calculate the hot junction temperature.