最小二乘法 线性 拟合 程序

1200阅读　0评论2015-11-19　number007cool
分类：C/C++


点击(此处)折叠或打开

  1. #include "stdafx.h"
  2. #include "process.h"

  4. /*******************************************************************************
  5. 阶乘函数宏定义
  6. *******************************************************************************/
  7. #define my_pow2(x) ((x) * (x))

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

  14. /*******************************************************************************
  15. 系数矩阵的初始化，没有撒谎，真的只有6行 + 一个循环
  16. *******************************************************************************/
  17. static int ParaInit(double* X, double* Y, int Amount)
  18. {
  19.         Para[1][1] = Amount;
  20.         for ( ; Amount; Amount--, X++, Y++)
  21.         {
  22.                 Para[0][0] += my_pow2(*X);
  23.                 Para[0][1] += (*X);
  24.                 Para[0][2] += (*X) * (*Y);
  25.                 Para[1][2] += (*Y);
  26.         }
  27.         Para[1][0] = Para[0][1];
  28.         return 0;
  29. }

  31. /*******************************************************************************
  32. 系数矩阵的运算，没有撒谎，真的只有7行
  33. *******************************************************************************/
  34. static int ParaDeal(void)
  35. {
  36.         Para[0][0] -= Para[1][0] * (Para[0][1] / Para[1][1]);
  37.         Para[0][2] -= Para[1][2] * (Para[0][1] / Para[1][1]);
  38.         Para[0][1] = 0;
  39.         Para[1][2] -= Para[0][2] * (Para[1][0] / Para[0][0]);
  40.         Para[1][0] = 0;
  41.         Para[0][2] /= Para[0][0];
  42.         //Para[0][0] = 1.0;
  43.         Para[1][2] /= Para[1][1];
  44.         //Para[1][1] = 1.0;
  45.         return 0;
  46. }

  48. /***********************************************************************************
  49. 从txt文件里读取double型的X，Y数据
  50. ***********************************************************************************/
  51. static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
  52. {
  53.         FILE* File = fopen(FileName, "r");
  54.         if (!File)
  55.                 return -1;
  56.         for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
  57.                 if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
  58.                         break;
  59.         fclose(File);
  60.         return 0;
  61. }

  63. /*******************************************************************************
  64. *******************************************************************************/
  65. int Test(const char* FileName)
  66. {
  67.         int Amount;
  68.         double BufferX[1024], BufferY[1024];
  69.         if (GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount))
  70.         {
  71.                 printf("读取数据文件失败！\r\n");
  72.                 return -1;
  73.         }
  74.         printf("读取数据文件成功，共有%d个点！\r\n", Amount);
  75.         ParaInit((double*)BufferX, (double*)BufferY, Amount);
  76.         ParaDeal();
  77.         printf("拟合数据成功，拟合直线为：\r\ny = (%lf) * x + (%lf);\r\n", Para[0][2], Para[1][2]);
  78.         system("pause");
  79.         return 0;
  80. }

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

上一篇：xml配置文件
下一篇：VC 模拟量标定