低階擬合算法的優化：
當階數比較低時，如直線或者拋物線，解行列式可以優化，不需要做成遞推模式，直接強行解就是，而且一般不會出現溢出問題。
還有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;
}
|
|