本方法转自https://github.com/natedomin/polyfit
polyfit.c:
/*------------------------------------------------------------
METHOD: polyfit
INPUTS: dependentValues[0..(countOfElements-1)] //即xData
independentValues[0...(countOfElements-1)] //即yData
countOfElements //即元素个数
order - Order of the polynomial fitting //要拟合的最高阶数
OUTPUTS: coefficients[0..order] - indexed by term
(the (coef*x^3) is coefficients[3])
//拟合的各项系数的结果,x^3 的系数为coefficients[3]
-------------------------------------------------------------*/
int polyfit(const double* const dependentValues ,
const double* const independentValues ,
unsigned int countOfElements ,
unsigned int order ,
double* coefficients)
{
// Declarations...
// ----------------------------------
enum { maxOrder = 5 };
double B[maxOrder + 1] = {0.0f};
double P[((maxOrder + 1) * 2) + 1] = {0.0f};
double A[(maxOrder + 1) * 2 * (maxOrder + 1)] = {0.0f};
double x , y , powx;
unsigned int ii , jj , kk;
// Verify initial conditions....
// ----------------------------------
// This method requires that the countOfElements >
// (order+1)
if (countOfElements <= order)
return -1;
// This method has imposed an arbitrary bound of
// order <= maxOrder. Increase maxOrder if necessary.
if (order > maxOrder)
return -1;
// Begin Code...
// ----------------------------------
// Identify the column vector
for (ii = 0; ii < countOfElements; ii++)
{
x = dependentValues[ii];
y = independentValues[ii];
powx = 1;
for (jj = 0; jj < (order + 1); jj++)
{
B[jj] = B[jj] + (y * powx);
powx = powx * x;
}
}
// Initialize the PowX array
P[0] = countOfElements;
// Compute the sum of the Powers of X
for (ii = 0; ii < countOfElements; ii++)
{
x = dependentValues[ii];
powx = dependentValues[ii];
for (jj = 1; jj < ((2 * (order + 1)) + 1); jj++)
{
P[jj] = P[jj] + powx;
powx = powx * x;
}
}
// Initialize the reduction matrix
//
for (ii = 0; ii < (order + 1); ii++)
{
for (jj = 0; jj < (order + 1); jj++)
{
A[(ii * (2 * (order + 1))) + jj] = P[ii + jj];
}
A[(ii * (2 * (order + 1))) + (ii + (order + 1))] = 1;
}
// Move the Identity matrix portion of the redux matrix
// to the left side (find the inverse of the left side
// of the redux matrix
for (ii = 0; ii < (order + 1); ii++)
{
x = A[(ii * (2 * (order + 1))) + ii];
if (x != 0)
{
for (kk = 0; kk < (2 * (order + 1)); kk++)
{
A[(ii * (2 * (order + 1))) + kk] =
A[(ii * (2 * (order + 1))) + kk] / x;
}
for (jj = 0; jj < (order + 1); jj++)
{
if ((jj - ii) != 0)
{
y = A[(jj * (2 * (order + 1))) + ii];
for (kk = 0; kk < (2 * (order + 1)); kk++)
{
A[(jj * (2 * (order + 1))) + kk] =
A[(jj * (2 * (order + 1))) + kk] -
y * A[(ii * (2 * (order + 1))) + kk];
}
}
}
}
else
{
// Cannot work with singular matrices
return -1;
}
}
// Calculate and Identify the coefficients
for (ii = 0; ii < (order + 1); ii++)
{
for (jj = 0; jj < (order + 1); jj++)
{
x = 0;
for (kk = 0; kk < (order + 1); kk++)
{
x = x + (A[(ii * (2 * (order + 1))) + (kk + (order + 1))] *
B[kk]);
}
coefficients[ii] = x;
}
}
return 0;
}
验证该函数的正确性,代码如下:
main.c
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#define ELE_COUNT 5
int main(void)
{
const unsigned int order = 3;
int result;
// These inputs should result in the following approximate coefficients:
// 0.5 2.5 1.0 3.0
// y = (0.5 * x^3) + (2.5 * x^2) + (1.0 * x) + 3.0
double xData[ELE_COUNT] = {12.0, 77.8, 44.1, 23.6, 108.2};
double yData[ELE_COUNT] = {1239.00, 250668.38, 47792.19, 7991.13, 662740.98};
double coefficients[order + 1]; // resulting array of coefs
// Perform the polyfit
result = polyfit(xData ,
yData ,
ELE_COUNT ,
order ,
coefficients);
printf("\n------------------------------\npolyfit result is:\n");
printf("%lfx^3+%lfx^2+%lfx+%f" , coefficients[3],coefficients[2],coefficients[1],
coefficients[0]);
printf("\n------------------------------\n");
return 0;
}
进行编译运行后,结果如下,可看出结果基本一致,也可自行拿其他函数进行验证。