低阶拟合算法的优化：
当阶数比较低时，如直线或者抛物线，解行列式可以优化，不需要做成递推模式，直接强行解就是，而且一般不会出现溢出问题。
还有pow（double x, double y）函数应该自己写替换，原因是在我们的应用中y肯定是整数，而这个函数按照double运算的，效率应该会差很多，所以来个阉割版的pow（double x, int y）.在结束低的时候，直接写两个宏定义。
#define my_pow0(x) (1.0)
#define my_pow1(x) (x)
#define my_pow2(x) (my_pow1(x) * my_pow1(x))
#define my_pow2(x) (my_pow1(x) * my_pow2(x))
#define my_pow4(x) (my_pow2(x) * my_pow2(x))
#define my_pow5(x) (my_pow2(x) * my_pow3(x))
#define my_pow6(x) (my_pow3(x) * my_pow3(x))
直线拟合的简化：
假设最后拟合直线：
y = a * x + b
残差方程：
（a * x + b – y）^2 = x^2 * a^2 + b^2 + y^2 + 2 * (x * a * b – x * y * a – y * b)
残差方程对a偏导：
x^2 * a + x^1 * b – y * x^1
残差方程对a偏导：
x^1 * a + x^0 * b – y * x^0
当残差方程之和取极值时，偏导为零。
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
（1）首先是定义一个2行3列的系数矩阵，并初始化
double Para[2] [3] = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
Para[1][1] = Amount;
for ( ; Amount; Amount--, X++, Y++)
{
Para[0][0] += my_pow2(*X);
Para[0][1] += my_pow1(*X);
//Para[1][1] += my_pow0(*Y);
Para[0][2] += my_pow1(*X) * my_pow1(*Y);
//Para[1][2] += my_pow0(*X) * my_pow1(*Y);
Para[1][2] += my_pow1(*Y);
}
Para[1][0] = Para[0][1];
return 0;
}
（2）解行列式：
由方程组
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
首先可以看到：
A，Para[0][0]和Para[1][1]肯定大于0
B，Para[0][1] = Para[1][0]
C， 可以证明[2][2]行列式的值
Para[0][0] * Para[1][1] - Para[0][1] * Para[1][0]
（∑x^2） * （∑x^0） - （∑x^1） * （∑x^1）> 0
也即是说，方阵正定，偏导为零的点就是最优值的点
可以看出，第1行的系数比第0行的系数要小很多，所以先采用第1行消第0行，再采用第0行消第1行。
int ParaDeal(void)
{
/*先把Para[0][1]消成0，由于有Para[1][1]肯定大于0，所以用乘法（除法)，除法更好一些，数据溢出的问题几乎不存在，适应性会更好*/
Para[0][0] = Para[0][0] - Para[1][0] * (Para[0][1] / Para[1][1];
Para[0][2] = Para[0][2] - Para[1][2] * (Para[0][1] / Para[1][1];
Para[0][1] = 0;
//再把Para[1][0] 消成0
Para[1][2] = Para[1][2] - Para[0][2] * (Para[1][0] / Para[0][0]);
Para[1][0] = 0;
//把Para[0][0]，Para[1][1] 消成1
Para[0][2] /= Para[0][0];//这个就是最后a的值
Para[0][0] = 1.0;
Para[1][2] /= Para[1][1]; //这个就是最后b的值
Para[1][1] = 1.0;
return 0;
}
至此，一阶直线拟合的推导过程和C代码全部完毕，可以看出运算量可以说很小，两个函数均只需调用一次即可得到结果，在源数据的个数和大小不是特别大的时候，完全可以移植到MCU中进行运算，实用价值比较高。
再来个二阶抛物线拟合的优化：
前期推导和一阶直线差不多，略去，系数矩阵方程组：
a * ∑x^4 + b * ∑x^3 + c * ∑x^2 = ∑(y * x^2)
a * ∑x^3 + b * ∑x^2 + c * ∑x^1 = ∑(y * x^1)
a * ∑x^2 + b * ∑x^1 + c * ∑x^0 = ∑(y * x^0)
（1）首先是定义一个3行4列的系数矩阵，并初始化
double Para[3] [4] = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
Para[1][1] = Amount;
for ( ; Amount; Amount--, X++, Y++)
{
Para[0][0] += my_pow4(*X);
Para[0][1] += my_pow3(*X);
Para[0][2] += my_pow2(*X);
Para[1][2] += my_pow1(*X);
//Para[2][2] += my_pow0(*X);
Para[0][3] += my_pow2(*X) * my_pow1(*Y);
Para[1][3] += my_pow1(*X) * my_pow1(*Y);
Para[2][3] += my_pow0(*X) * my_pow1(*Y);
}
Para[1][0] = Para[0][1];
Para[1][1] = Para[0][2];
Para[2][0] = Para[1][1];
Para[2][1] = Para[1][2];
return 0;
}
a * ∑x^4 + b * ∑x^3 + c * ∑x^2 = ∑(y * x^2)
a * ∑x^3 + b * ∑x^2 + c * ∑x^1 = ∑(y * x^1)
a * ∑x^2 + b * ∑x^1 + c * ∑x^0 = ∑(y * x^0)
首先可以看到：
A，Para[0][0],Para[1][1], Para[2][2]肯定大于0
B，左边3*3的方阵是对称矩阵
D， 可以证明[3][3]行列式的值 > 0,也即是说，方阵正定，偏导为零的点就是最优值的点
可以看出，第k + 1行的系数比第k行的系数要小很多，所以先采用第k +1行消第k行，再采用第k行消第k + 1行。
int ParaDeal(void)
{
//用Para[2][2]把Para[0][2]消成0
Para[0][0] = Para[0][0] - Para[2][0] * (Para[0][2] / Para[2][2]);
Para[0][1] = Para[0][1] - Para[2][1] * (Para[0][2] / Para[2][2]);
Para[0][3] = Para[0][3] - Para[2][3] * (Para[0][2] / Para[2][2]);
Para[0][2] = 0;
//用Para[2][2]把Para[1][2] 消成0
Para[1][0] = Para[1][0] - Para[2][0] * (Para[1][2] / Para[2][2]);
Para[1][1] = Para[1][1] - Para[2][1] * (Para[1][2] / Para[2][2]);
Para[1][3] = Para[1][3] - Para[2][3] * (Para[1][2] / Para[2][2]);
Para[1][2] = 0;
//用Para[1][1]把Para[0][1]消成0
Para[0][0] = Para[0][0] - Para[1][0] * (Para[0][1] / Para[1][1]);
Para[0][3] = Para[0][3] - Para[1][3] * (Para[0][1] / Para[1][1] );
Para[0][1] = 0;
//至此，已经成成为三角矩阵
//用Para[0][0]把Para[1][0]消成0
Para[1][3] = Para[1][3] - Para[0][3] * (Para[1][0] / Para[0][0]);
Para[1][0] = 0;
//用Para[0][0]把Para[2][0]消成0
Para[2][3] = Para[2][3] - Para[0][3] * (Para[2][0] / Para[0][0]);
Para[2][0] = 0;
//用Para[1][1]把Para[2][1]消成0
Para[2][3] = Para[2][3] - Para[1][3] * (Para[2][1] / Para[1][1]);
Para[2][1] = 0;
//至此，已经成成为对角矩阵
//把Para[0][0]，Para[1][1]，Para[2][2]消成1
Para[0][3] /= Para[0][0];//这个就是最后a的值
Para[0][0] = 1.0;
Para[1][3] /= Para[1][1]; //这个就是最后b的值
Para[1][1] = 1.0;
Para[2][3] /= Para[2][2]; //这个就是最后c的值
Para[2][2] = 1.0;
return 0;
}
|
| | [tr][/tr]
[tr][/tr]
|
QQ好友和群
QQ空间
腾讯微博
腾讯朋友
收藏
淘帖
顶
踩
