最小二乘法 线性 拟合 程序

#include "stdafx.h"
#include "process.h"

/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x) ((x) * (x))

/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[2][3] = {0.0};

/*******************************************************************************
系数矩阵的初始化，没有撒谎，真的只有6行 + 一个循环
*******************************************************************************/
static int ParaInit(double* X, double* Y, int Amount)
{
        Para[1][1] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
        {
                Para[0][0] += my_pow2(*X);
                Para[0][1] += (*X);
                Para[0][2] += (*X) * (*Y);
                Para[1][2] += (*Y);
        }
        Para[1][0] = Para[0][1];
        return 0;
}

/*******************************************************************************
系数矩阵的运算，没有撒谎，真的只有7行
*******************************************************************************/
static int ParaDeal(void)
{
        Para[0][0] -= Para[1][0] * (Para[0][1] / Para[1][1]);
        Para[0][2] -= Para[1][2] * (Para[0][1] / Para[1][1]);
        Para[0][1] = 0;
        Para[1][2] -= Para[0][2] * (Para[1][0] / Para[0][0]);
        Para[1][0] = 0;
        Para[0][2] /= Para[0][0];
        //Para[0][0] = 1.0;
        Para[1][2] /= Para[1][1];
        //Para[1][1] = 1.0;
        return 0;
}

/***********************************************************************************
从txt文件里读取double型的X，Y数据
***********************************************************************************/
static 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;
}

/*******************************************************************************
*******************************************************************************/
int Test(const char* FileName)
{
        int Amount;
        double BufferX[1024], BufferY[1024];
        if (GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount))
        {
                printf("读取数据文件失败！\r\n");
                return -1;
        }
        printf("读取数据文件成功，共有%d个点！\r\n", Amount);
        ParaInit((double*)BufferX, (double*)BufferY, Amount);
        ParaDeal();
        printf("拟合数据成功，拟合直线为：\r\ny = (%lf) * x + (%lf);\r\n", Para[0][2], Para[1][2]);
        system("pause");
        return 0;
}

int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}