C语言实现最小二乘法抛物线拟合

archi9966

如题~求大虾们指点~

背负拉天的云

主要是你想拟合成什么样的？用直线拟合曲线的话 是使得误差最小的直线  曲线拟合直线的话 一般二阶就可以了

zouzhichao

matlab 实现，再转C
或者直接推导公式实现吧

zouzhichao

y = a * x^2 + b * x + c
Error = (a * x^2 + b * x + c - y)^2
Error = a^2 * x^4 +b^2 * x^2 + C^2 + y^2 + 2 * (a * b * x^3 + a * c * x^2 - a * y * x^2 + b * c * x - b * y * y - c * y)

Error(a) = 2 * x^4 * a + 2 * b * x^3 + 2 * c * x^2 - 2 * y * x^2
Error(b) = 2 * x^2 * b + 2 * a * x^3 + 2 * c * x - 2 * y * x
Error(c) = 2 * c + 2 * a * x^2 + 2 * b * x - 2 * y

zouzhichao

a * x^4 + b * x^3 + c * x^2 = y * x^2
a * x^3 + b * x^2 + c * x = y * x
a * x^2 + b * x + c = y

LearningASM

算法运算量注定很大，用在电脑上就可以，如果是在MSP上的话。。。。。

zouzhichao

C语言实现，已通过matlab对比验证：
/***********************************************************************************
从txt文件里读取double型XY数据
***********************************************************************************/
int GetXY(const char* FileName, double* X, double* Y, int* Amount)
{
        FILE* File = fopen(FileName, "r");
        if (!File)
                return -1;
        for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
                if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
                        break;
        fclose(File);
        return 0;
}

/***********************************************************************************
从double型XY数据中获取系数矩阵
***********************************************************************************/
int GetK(double* K, const double* X, const double* Y, int Amount)
{
        *(K + 10) = Amount;
        *K = 0;
        *(K + 1) = 0;
        *(K + 2) = 0;
        *(K + 3) = 0;
        *(K + 6) = 0;
        *(K + 7) = 0;
        *(K + 11) = 0;
        for ( ; Amount; Amount--, X++, Y++)
        {
                *K += (*X) * (*X) * (*X) * (*X);
                *(K + 1) += (*X) * (*X) * (*X);
                *(K + 2) += (*X) * (*X);
                *(K + 3) += (*X) * (*X) * (*Y);
                *(K + 6) += (*X);
                *(K + 7) += (*X) * (*Y);
                *(K + 11) += (*Y);
        }
        *(K + 4) = *(K + 1);
        *(K + 5) = *(K + 2);
        *(K + 8) = *(K + 2);
        *(K + 9) = *(K + 6);
        return 0;
}

/***********************************************************************************
系数矩阵变换，解矩阵方程
***********************************************************************************/
int TransK(double* K)
{
        //
        *(K + 5) = (*(K + 5)) * (*K) - (*(K + 4)) * (*(K + 1));
        *(K + 6) = (*(K + 6)) * (*K) - (*(K + 4)) * (*(K + 2));
        *(K + 7) = (*(K + 7)) * (*K) - (*(K + 4)) * (*(K + 3));
        *(K + 4) = 0;
        *(K + 9) = (*(K + 9)) * (*K) - (*(K + 8)) * (*(K + 1));
        *(K + 10) = (*(K + 10)) * (*K) - (*(K + 8)) * (*(K + 2));
        *(K + 11) = (*(K + 11)) * (*K) - (*(K + 8)) * (*(K + 3));
        *(K + 8) = 0;
        //
        //
        *(K + 10) = (*(K + 10)) * (*(K + 5)) - (*(K + 9)) * (*(K + 6));
        *(K + 11) = (*(K + 11)) * (*(K + 5)) - (*(K + 9)) * (*(K + 7));
        *(K + 9) = 0;
        //
        //
        *(K + 5) = (*(K + 5)) * (*(K + 10));
        *(K + 7) = (*(K + 7)) * (*(K + 10)) - (*(K + 6)) * (*(K + 11));
        *(K + 6) = 0;
        //
        //
        *K = (*K) * (*(K + 5));
        *(K + 2) = (*(K + 2)) * (*(K + 5));
        *(K + 3) = (*(K + 3)) * (*(K + 5)) - (*(K + 1)) * (*(K + 7));
        *(K + 1) = 0;
        //
        //
        *K = (*K) * (*(K + 10));
        *(K + 3) = (*(K + 3)) * (*(K + 10)) - (*(K + 2)) * (*(K + 11));
        *(K + 2) = 0;
        //
        //
        if (0 != *(K + 00))
        {
                *(K + 3) /= *(K + 00);
                *K = 1.0;
        }
        if (0 != *(K + 5))
        {
                *(K + 7) /= *(K + 5);
                *(K + 5) = 1.0;
        }
        if (0 != *(K + 10))
        {
                *(K + 11) /= *(K + 10);
                *(K + 10) = 1.0;
        }
        //
        return 0;
}

/***********************************************************************************
***********************************************************************************/
int Cal(const char* FileName, double* ParaA, double* ParaB, double* ParaC)
{
        double ParaK[12], BufferX[1024], BufferY[1024];
        int Amount;
        GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount);
        GetK((double*)ParaK, (const double*)BufferX, (const double*)BufferY, Amount);
        TransK((double*)ParaK);
        *ParaA = ParaK[3];
        *ParaB = ParaK[7];
        *ParaC = ParaK[11];
        return 0;
}

int main(int argc, char* argv[])
{
        double ParaA, ParaB, ParaC;
        Cal((const char*)"test.txt", &ParaA, &ParaB, &ParaC);
        printf("ParaA = %lf\r\nParaB = %lf\r\nParaC = %lf\r\n", ParaA, ParaB, ParaC);
        return 0;
}

zouzhichao

使用了大量浮点乘法运算（也可以换成整型运算，但是精度下降），尽量避免了使用除法运算，运算量还算不大

shufang0520

路过学习一下

archi9966

zouzhichao 发表于 2013-5-19 20:35
C语言实现，已通过matlab对比验证：
/***************************************************************** ...

人才啊， 多谢啦

cumtyap

学习学习啊

ljc2012

早点有就好了 ！当初也是要做个曲线拟合 ，但不知道怎么搞..多谢了！

zpc2013

是呀，求解

emobbs

曲线拟合