多项式拟合算法

具体效果参考excel

// Function to perform polynomial fit
// degree      hightest m
// numPoints   num of data
// x           ram data
// y           target data
//coeff        the output data
void polynomialFit(int degree, int numPoints, double x[], double y[], double coeff[]) {
    int i, j, k;
    double X[2 * degree + 1]; // Array to store the values of sigma(x[i]^k)
    for (i = 0; i < 2 * degree + 1; i++) {
        X[i] = 0;
        for (j = 0; j < numPoints; j++) {
            X[i] = X[i] + pow(x[j], i);
        }
    }

    double B[degree + 1][degree + 2]; // Augmented matrix
    double a[degree + 1];             // Output array
    for (i = 0; i <= degree; i++) {
        for (j = 0; j <= degree; j++) {
            B[i][j] = X[i + j]; // Build the augmented matrix
        }
    }
    double Y[degree + 1]; // Array to store the values of sigma(x[i]^k * y[i])
    for (i = 0; i < degree + 1; i++) {
        Y[i] = 0;
        for (j = 0; j < numPoints; j++) {
            Y[i] = Y[i] + pow(x[j], i) * y[j];
        }
    }
    for (i = 0; i <= degree; i++) {
        B[i][degree + 1] = Y[i];
    }
    degree = degree + 1;
    for (i = 0; i < degree; i++) {
        for (k = i + 1; k < degree; k++) {
            if (B[i][i] < B[k][i]) {
                for (j = 0; j <= degree; j++) {
                    double temp = B[i][j];
                    B[i][j] = B[k][j];
                    B[k][j] = temp;
                }
            }
        }
    }
    for (i = 0; i < degree - 1; i++) {
        for (k = i + 1; k < degree; k++) {
            double t = B[k][i] / B[i][i];
            for (j = 0; j <= degree; j++) {
                B[k][j] = B[k][j] - t * B[i][j];
            }
        }
    }
    for (i = degree - 1; i >= 0; i--) {
        a[i] = B[i][degree];
        for (j = 0; j < degree; j++) {
            if (j != i) {
                a[i] = a[i] - B[i][j] * a[j];
            }
        }
        a[i] = a[i] / B[i][i];
    }
    for (i = 0; i < degree; i++) {
        coeff[i] = a[i];
    }
}