Click here to Skip to main content
14,642,548 members
Articles » General Programming » Algorithms & Recipes » Math
Tip/Trick
Posted 9 Jun 2020

Stats

5.6K views
176 downloads
6 bookmarked

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

Rate this:
4.83 (3 votes)
Please Sign up or sign in to vote.
4.83 (3 votes)
10 Jun 2020CPOL
Polynomial interpolation algorithm from 1 to 9 degree that allow forcing constant term to 0
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.

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:

Image 1

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:

Image 2

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.

Image 3

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 = a1x^3+b1x^2+c1x

a1, b1, c1 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

  • 9th June, 2020: Initial version

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

fernando amatulli
Technical Lead
Italy Italy
No Biography provided
Group type: Organisation (No members)



Comments and Discussions

 
Questionconstant term... Pin
Stefan_Lang14-Jul-20 5:55
mveStefan_Lang14-Jul-20 5:55 
Questionseems good as i look at it but... Pin
Member 156562714-Jun-20 3:47
MemberMember 156562714-Jun-20 3:47 
AnswerRe: seems good as i look at it but... Pin
fernando amatulli23-Jun-20 2:03
professionalfernando amatulli23-Jun-20 2:03 
AnswerRe: seems good as i look at it but... Pin
fernando amatulli27-Jun-20 2:53
professionalfernando amatulli27-Jun-20 2:53 
GeneralRe: seems good as i look at it but... Pin
Member 15656271-Jul-20 13:30
MemberMember 15656271-Jul-20 13:30 
QuestionSolution question Pin
avisal11-Jun-20 7:21
professionalavisal11-Jun-20 7:21 
QuestionA Very Good Start Pin
Rick York9-Jun-20 9:49
mveRick York9-Jun-20 9:49 
This is a very good start but, in my opinion, it is lacking a bit and for that reason I have withheld voting on it.

I think you should present some comparative analysis. For example, define a random polynomial function of arbitrary order and then have it compute a set of a few dozen points or however many you want to give it. Then run the points through your function and let it calculate the coefficients. If they match exactly, great! If not then evaluate the error percentages.

In my line of work I have had to do a LOT of filtering, interpolation, and curve fitting of data so I have evaluated quite a number of these algorithms and that is how I do it. Sometimes I give it a random set of points and then evaluate the error and compare with alternatives.

What you have now is a good start but I think you should take it a step or two farther to give readers a bigger picture.
"They have a consciousness, they have a life, they have a soul! Damn you! Let the rabbits wear glasses! Save our brothers! Can I get an amen?"

AnswerRe: A Very Good Start Pin
Nelek9-Jun-20 19:49
protectorNelek9-Jun-20 19:49 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.