增加拉格朗日插值和多项式拟合算法。
Signed-off-by: lion.chan <cy187lion@sina.com>
This commit is contained in:
parent
efdd665289
commit
d74958ebeb
51
Algorithm/DSP/LagrangeInterpolation.c
Normal file
51
Algorithm/DSP/LagrangeInterpolation.c
Normal file
@ -0,0 +1,51 @@
|
||||
static double domega(const double* xn, const unsigned int n, const unsigned int k, const double x);
|
||||
|
||||
/**
|
||||
* @brief 拉格朗日插值算法。
|
||||
*
|
||||
* @param xn 输入：已知采样点的 x 坐标。
|
||||
* @param yn 输入：已知采样点的 y 坐标。
|
||||
* @param n 输入：采样点数。
|
||||
* @param x 输入：插值点 x 坐标。
|
||||
* @return double 输出：插值点的 y 坐标值。
|
||||
*/
|
||||
double Lagrange(const double* xn, const double* yn, const unsigned int n, const double x) {
|
||||
unsigned int i;
|
||||
double l = 0.0;
|
||||
|
||||
for (i=0; i<n; i++)
|
||||
l += yn[i]*domega(xn, n, i, x);
|
||||
|
||||
return l;
|
||||
}
|
||||
|
||||
static double domega(const double* xn, const unsigned int n, const unsigned int k, const double x) {
|
||||
unsigned int i;
|
||||
double omega=1.0;
|
||||
double omega_deri = 1.0;
|
||||
|
||||
for (i=0; i<n; i++) {
|
||||
if (i==k)
|
||||
continue;
|
||||
|
||||
omega *= (x-xn[i]);
|
||||
omega_deri *= (xn[k] - xn[i]);
|
||||
}
|
||||
return omega/omega_deri;
|
||||
}
|
||||
|
||||
|
||||
/** Demo */
|
||||
static double Xn[] = {0, 0.3000, 0.6000, 0.9000, 1.2000, 1.5000, 1.8000, 2.1000, 2.4000, 2.7000, 3.0000};
|
||||
static double Yn[] = {2.0000, 2.3780, 3.9440, 7.3460, 13.2320, 22.2500, 35.0480, 52.2740, 74.5760, 102.6020, 137.0000};
|
||||
|
||||
/**
|
||||
* @brief 拉格朗日插值示例。
|
||||
*
|
||||
* @return int
|
||||
*/
|
||||
int main(void) {
|
||||
double out = 0.0;
|
||||
out = Lagrange(Xn, Yn, 11, 1.5);
|
||||
return 0;
|
||||
}
|
||||
164
Algorithm/DSP/Polyfit.c
Normal file
164
Algorithm/DSP/Polyfit.c
Normal file
@ -0,0 +1,164 @@
|
||||
/**
|
||||
* @brief 多项式拟合，求解拟合系数。
|
||||
*
|
||||
* @param dependentValues 输入：自变量。
|
||||
* @param independentValues 输入：因变量。
|
||||
* @param countOfElements 输入：采样点数量。
|
||||
* @param order 输入：序数。
|
||||
* @param coefficients 输出：拟合系数。
|
||||
* @return int
|
||||
*/
|
||||
int Polyfit(const double* const dependentValues,
|
||||
const double* const independentValues,
|
||||
unsigned int countOfElements,
|
||||
unsigned int order,
|
||||
double* coefficients)
|
||||
{
|
||||
// Declarations...
|
||||
// ----------------------------------
|
||||
enum {maxOrder = 5};
|
||||
|
||||
double B[maxOrder+1] = {0.0};
|
||||
double P[((maxOrder+1) * 2)+1] = {0.0};
|
||||
double A[(maxOrder + 1)*2*(maxOrder + 1)] = {0.0};
|
||||
|
||||
double x, y, powx;
|
||||
|
||||
unsigned int ii, jj, kk;
|
||||
|
||||
// Verify initial conditions....
|
||||
// ----------------------------------
|
||||
|
||||
// This method requires that the countOfElements >
|
||||
// (order+1)
|
||||
if (countOfElements <= order)
|
||||
return -1;
|
||||
|
||||
// This method has imposed an arbitrary bound of
|
||||
// order <= maxOrder. Increase maxOrder if necessary.
|
||||
if (order > maxOrder)
|
||||
return -1;
|
||||
|
||||
// Begin Code...
|
||||
// ----------------------------------
|
||||
|
||||
// Identify the column vector
|
||||
for (ii = 0; ii < countOfElements; ii++)
|
||||
{
|
||||
x = dependentValues[ii];
|
||||
y = independentValues[ii];
|
||||
powx = 1;
|
||||
|
||||
for (jj = 0; jj < (order + 1); jj++)
|
||||
{
|
||||
B[jj] = B[jj] + (y * powx);
|
||||
powx = powx * x;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize the PowX array
|
||||
P[0] = countOfElements;
|
||||
|
||||
// Compute the sum of the Powers of X
|
||||
for (ii = 0; ii < countOfElements; ii++)
|
||||
{
|
||||
x = dependentValues[ii];
|
||||
powx = dependentValues[ii];
|
||||
|
||||
for (jj = 1; jj < ((2 * (order + 1)) + 1); jj++)
|
||||
{
|
||||
P[jj] = P[jj] + powx;
|
||||
powx = powx * x;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize the reduction matrix
|
||||
//
|
||||
for (ii = 0; ii < (order + 1); ii++)
|
||||
{
|
||||
for (jj = 0; jj < (order + 1); jj++)
|
||||
{
|
||||
A[(ii * (2 * (order + 1))) + jj] = P[ii+jj];
|
||||
}
|
||||
|
||||
A[(ii*(2 * (order + 1))) + (ii + (order + 1))] = 1;
|
||||
}
|
||||
|
||||
// Move the Identity matrix portion of the redux matrix
|
||||
// to the left side (find the inverse of the left side
|
||||
// of the redux matrix
|
||||
for (ii = 0; ii < (order + 1); ii++)
|
||||
{
|
||||
x = A[(ii * (2 * (order + 1))) + ii];
|
||||
if (x != 0.0)
|
||||
{
|
||||
for (kk = 0; kk < (2 * (order + 1)); kk++)
|
||||
{
|
||||
A[(ii * (2 * (order + 1))) + kk] =
|
||||
A[(ii * (2 * (order + 1))) + kk] / x;
|
||||
}
|
||||
|
||||
for (jj = 0; jj < (order + 1); jj++)
|
||||
{
|
||||
if ((jj - ii) != 0)
|
||||
{
|
||||
y = A[(jj * (2 * (order + 1))) + ii];
|
||||
for (kk = 0; kk < (2 * (order + 1)); kk++)
|
||||
{
|
||||
A[(jj * (2 * (order + 1))) + kk] =
|
||||
A[(jj * (2 * (order + 1))) + kk] -
|
||||
y * A[(ii * (2 * (order + 1))) + kk];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Cannot work with singular matrices
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate and Identify the coefficients
|
||||
for (ii = 0; ii < (order + 1); ii++)
|
||||
{
|
||||
for (jj = 0; jj < (order + 1); jj++)
|
||||
{
|
||||
x = 0;
|
||||
for (kk = 0; kk < (order + 1); kk++)
|
||||
{
|
||||
x = x + (A[(ii * (2 * (order + 1))) + (kk + (order + 1))] *
|
||||
B[kk]);
|
||||
}
|
||||
coefficients[ii] = x;
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
/** Demo */
|
||||
static double Xn[] = {0, 0.3000, 0.6000, 0.9000, 1.2000, 1.5000, 1.8000, 2.1000, 2.4000, 2.7000, 3.0000};
|
||||
static double Yn[] = {2.0000, 2.3780, 3.9440, 7.3460, 13.2320, 22.2500, 35.0480, 52.2740, 74.5760, 102.6020, 137.0000};
|
||||
|
||||
/**
|
||||
* @brief 求拟合系数及多项式拟合示例。
|
||||
*
|
||||
* @return int
|
||||
*/
|
||||
int main(void)
|
||||
{
|
||||
double out = 0.0;
|
||||
unsigned int i;
|
||||
double coff[3];
|
||||
|
||||
polyfit(Xn, Yn, 11, 2, coff);
|
||||
|
||||
for (i=0; i<3; i++)
|
||||
{
|
||||
out += coff[i]*pow(Xn[4], i);
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user