Polynomial Interpolation Algorithm from 1 to 9 Degree that Allows to Force the Constant Term to 0

• Download source - 188.5 KB

Introduction

For applications like transducers linearization, it is mandatory to have that when the transducer output is 0 the linearizated output is also 0. This is the case, for example, of a load cell or pressure transducer. When the load or pressure is 0, the linearized output must be 0 too.

Standard polynomial algorithms compute all the polynomial coefficients included the constant term that normal is different from 0. If applied, this leads to an error. The same happens if you just decide to not consider the constant term. This program allows to compute the coefficient of the polynomial forcing the constant term to 0 from a polynomial of degree 1 to 9. Just for convenience, the program also computes the complete polynomial with the constant term.

Background

• Visual C++ using MFC

Using the Code

This article can be useful for those who must linearize a transducer (like load cell or pressure transducer) that are typically not enough linear to use them directly when accurate values are needed.

To correct the output, it frequently uses a polynomial correction. The degree of the polynomial depends on how much the transducer is not linear.

Normally, a second or third degree polynomial is enough if the transducer has just one or two non linear sections, but in some cases, a polynomial of higher degree is necessary. See the picture below where a third degree polynomial has been computed:

Here you can see that the fitting curve does not intersect the X-Y axis in the points (0,0). If we compute the polynomial coefficients by forcing the constant term to 0, we obtain a curve fitting as below:

The first step in the calibration process is to measure the transducer output (x values) and compare it with a reference transducers (y values) that are assumed as the correct values. We obtain a grid X-Y like below: some points are repeated in the calibration process because more cycles have been performed.

The linearization equation is computed as below (I assumed a third degree polynomial).

y = ax^3+bx^2+cx+d

where a, b, c, d are the polynomial coefficients computed by the algorithm.

If d is different from 0 when x = 0, then y=d that is wrong (of course x=0 must involve y=0). A different error arises if you decide just to neglect the term 'd'. The curve is just shifted (below or above) and a measurement error will occur. So the idea is to compute the polynomial forcing the term 'd' to be 0.

This leads to a different algorithm and you will have:

y = a₁x^3+b₁x^2+c₁x

a_1, b_1, c₁ are new coefficients that solve the problem that arises above.

The program has been developed in Visual C++ using MFC in Visual Studio 13.

For data input, I used an old grid control from Microsoft FlexGrid. In any case, it was included in the zip file (see the folder <componenti>) so if you don't have it installed on your computer, you need just to register it and use it.

To register, use the command:

regsvr32 msflxgrd.ocx

The program can be used just as it is because it allows to load a data table (csv file), save it, compute the polynomial coefficient, both with the constant term forced to 0 or not, evaluate graphically the data fitting and save the coefficient in a file. It is possible to choose the degree (from 1 to 9) of the interpolating polynomial.

It is possible to show a graph where the errors among the fitting curve and the measured points are plotted to better evaluate the curve fitting.

The two algorithms used are reported below:

• The first one "compute_terminenoto()" computes all the coefficients (including the constant term).
• The second one 'compute_no_terminenoto()' computes the coefficients by forcing the constant term to 0.

Both the routines use the variable 'ncol' to select the degree of the polynomial and 'nrow' to select the number of points considered. Actually 'ncol' is the number of the coefficients to compute (including the constant term) so the degree of the polynomial is ncol-1.

The matrix Val[][] contains the input values:

Val[][0] are the y values while Val[][1] are the x values.

The resulting coefficients are stored in the w[] array:

/*
This procedure calculates the polynomial only with the coefficient of the term 
known for all degrees of the polynomial
*/

void CMyFitpolDlg::compute_terminenoto()
{
    long double zy[MAX_DEGREE * 2 + 1], zx[MAX_DEGREE * 2 + 1], 
    ZMatr[MAX_DEGREE * 2 + 1][MAX_DEGREE * 2 + 1], ZTemp;
    int i, k, j, Grado;

    Grado = ncol - 1;        
                             
    for (i = 0; i <= Grado * 2; i++)
        zx[i] = 0.0;

    for (i = 0; i <= Grado; i++)
        zy[i] = 0.0;

    for (i = 0; i<nrow; i++)
    {
        for (j = 0; j <= Grado * 2; j++)
            zx[j] += pow(Val[i][1], j);

        for (j = 0; j <= Grado; j++)
            zy[j] += Val[i][0] * pow(Val[i][1], j);

    }

    for (i = 1; i <= Grado + 1; i++)
    {
        k = Grado - 1 + i;
        for (j = 1; j <= Grado + 1; j++)
        {
            ZMatr[i][j] = zx[k];
            k--;
        }
    }

    for (i = 1; i<Grado + 2; i++)
        ZMatr[i][Grado + 2] = zy[i - 1];

    for (i = 1; i <= Grado + 1; i++)
    {
        for (k = (i + 1); k <= Grado + 2; k++)
            ZMatr[i][k] = ZMatr[i][k] / ZMatr[i][i];

        ZMatr[i][i] = 1.0;
        for (j = (i + 1); j <= Grado + 1; j++)
        {
            ZTemp = ZMatr[j][i];
            for (k = 1; k <= Grado + 2; k++)
                ZMatr[j][k] = ZMatr[j][k] - ZMatr[i][k] * ZTemp;
        }
    }

    for (i = Grado + 1; i>0; i--)
    {
        ZMatr[i][Grado + 1] = ZMatr[i][Grado + 1] / ZMatr[i][i];
        ZMatr[i][i] = 1.0;
        for (j = (i - 1); j>0; j--)
        {
            ZMatr[j][Grado + 2] = ZMatr[j][Grado + 2] - ZMatr[i][Grado + 2] * ZMatr[j][i];
            ZMatr[j][i] = 0.0;
        }
    }
    for (i = 0; i < Grado + 1; i++)
        w[Grado - i][0] = ZMatr[i + 1][Grado + 2];
}

/*
This procedure calculates the polynomial without coefficient of the term known 
for all degrees of the polynomial
*/


void CMyFitpolDlg::compute_no_terminenoto()
{
    long double zy[MAX_DEGREE * 2 + 1], zx[MAX_DEGREE * 2 + 1], 
                ZMatr[MAX_DEGREE * 2 + 1][MAX_DEGREE * 2 + 1], ZTemp;
    int i, k, j, Myncol;

    Myncol = ncol - 1;
    for (i = 0; i < Myncol * 2; i++)
        zx[i] = 0.0;

    for (i = 0; i <= Myncol; i++)
        zy[i] = 0.0;

    for (i = 0; i<nrow; i++)
    {
        for (j = 1; j <= Myncol * 2; j++)
            zx[j - 1] += pow(Val[i][1], j);

        for (j = 0; j <= Myncol; j++)
            zy[j] += Val[i][0] * pow(Val[i][1], j);
    }

    for (i = 1; i < Myncol + 1; i++)
    {
        k = Myncol - 1 + i;
        for (j = 1; j < Myncol + 1; j++)
        {
            ZMatr[i][j] = zx[k];
            k--;
        }
    }

    for (i = 1; i<Myncol + 1; i++)
        ZMatr[i][Myncol + 1] = zy[i];

    for (i = 1; i <= Myncol + 1; i++)
    {
        for (k = (i + 1); k <= Myncol + 2; k++)
            ZMatr[i][k] = ZMatr[i][k] / ZMatr[i][i];

        ZMatr[i][i] = 1.0;
        for (j = (i + 1); j <= Myncol + 1; j++)
        {
            ZTemp = ZMatr[j][i];
            for (k = 1; k <= Myncol + 2; k++)
                ZMatr[j][k] = ZMatr[j][k] - ZMatr[i][k] * ZTemp;
        }
    }

    for (i = Myncol; i>0; i--)
    {
        ZMatr[i][Myncol + 1] = ZMatr[i][Myncol + 1] / ZMatr[i][i];
        ZMatr[i][i] = 1.0;
        for (j = (i - 1); j>0; j--)
        {
            ZMatr[j][Myncol + 1] = ZMatr[j][Myncol + 1] - ZMatr[i][Myncol + 1] * ZMatr[j][i];
            ZMatr[j][i] = 0.0;
        }
    }
    for (i = 0; i < Myncol + 1; i++)
        w[Myncol - i][0] = ZMatr[i + 1][Myncol + 1];
    w[0][0] = 0.0;
}

History

• 9^th June, 2020: Initial version