一、涉及的工具和库
我编译用的工具是Qt,不过本文中未设计任何Qt相关的类,用的是C++和eigen库。
二、头文件
#ifndef FITCURVE_H
#define FITCURVE_H
#include "FitCurve_global.h"
#include "vector"
#include "numeric"
#include "Eigen/Dense"
using namespace std;
class FitCurve
{
public:
FitCurve();
double quadraticFitCurve(vector<float> &X,vector<float>
};
#endif // FITCURVE_H
三、源文件
函数的X,Y为散点坐标值列表,orders代表多项式最高次数,coef为存放多项式系数的数组,coef数组长度为orders+1
返回值是方差R²
#include "fitcurve.h"
FitCurve::FitCurve()
{
}
double quadraticFitCurve(vector<float> &X, vector<float> &Y, uint8_t orders, double *coef)
{
if (X.size() < 2 || Y.size() < 2 || X.size() != Y.size() || orders < 1)
exit(EXIT_FAILURE);
// map sample data from STL vector to eigen vector
Eigen::Map<Eigen::VectorXf> sampleX(X.data(), X.size());
Eigen::Map<Eigen::VectorXf> sampleY(Y.data(), Y.size());
Eigen::MatrixXf mtxVandermonde(X.size(), orders + 1); // Vandermonde matrix of X-axis coordinate vector of sample data
Eigen::VectorXf colVandermonde = sampleX; // Vandermonde column
// construct Vandermonde matrix column by column
for (size_t i = 0; i < orders + 1; ++i)
{
if (0 == i)
{
mtxVandermonde.col(0) = Eigen::VectorXf::Constant(X.size(), 1, 1);
continue;
}
if (1 == i)
{
mtxVandermonde.col(1) = colVandermonde;
continue;
}
colVandermonde = colVandermonde.array()*sampleX.array();
mtxVandermonde.col(i) = colVandermonde;
}
// calculate coefficients vector of fitted polynomial
Eigen::VectorXf result = (mtxVandermonde.transpose()*mtxVandermonde).inverse()*(mtxVandermonde.transpose())*sampleY;
for(uint8_t j = 0; j < orders + 1;j++){
coef[j]=result[j];
}
float sum = std::accumulate(std::begin(Y),std::end(Y),0.0);
float mean = sum/Y.size();//求Y的平均值
double SSR=0;
double SST=0;
for(uint8_t j=0;j<Y.size();j++){
float yi = 0;
for(uint8_t k=0;k<orders+1;k++){
yi += result[k]*pow(X[j],k);
}
SSR += pow(yi-Y[j],2);
SST += pow(Y[j]-mean,2);
}
double R = 1-SSR/SST;
return R;
}
四、调用
#include <vector>
vector<float> x,y;
x.push_back(192);
x.push_back(190);
x.push_back(187);
x.push_back(184);
x.push_back(178);
x.push_back(170);
x.push_back(163);
x.push_back(153);
x.push_back(138);
x.push_back(120);
y.push_back(191);
y.push_back(189);
y.push_back(187);
y.push_back(185);
y.push_back(181);
y.push_back(175);
y.push_back(170);
y.push_back(163);
y.push_back(153);
y.push_back(139);
double xs[3];
double R=quadraticFitCurve(x,y,2,xs);