文章目录
原文链接:
/weixin_44344462/article/details/88850409
/xsz591541060/article/details/107532611
配置Eigen矩阵运算库
后续计算需要利用矩阵运算来求解拟合系数,用到Eigen矩阵运算库,配置方法可自行搜索或MacOs可以参考Mac配置Eigen库进行配置。
拟合原理
以二次曲线拟合为例,拟合曲线应满足以下方程:
即有:
由上面最后一个等式利用矩阵的转置与求逆,则可以得出其拟合曲线的系数W矩阵。
C++最小二乘拟合代码
以下代码理论上可以完成N阶的曲线拟合,需要对参数N进行设置,并提供正确的数据点。
#include <Eigen/Core>#include <Eigen/Dense>#include <Eigen/Geometry>#include <Eigen/Eigenvalues>#include <iostream>#include <fstream>#include <string>#include <vector>#include <cmath>using namespace std;int main() {// txt点数据文件路径const string fileName = "./points.txt";// 设置是几次拟合const int N = 2;// 创建两个vectorvector<float> x, y;// 读取文件ifstream f(fileName);if (!f) {cout << "数据文件打开失败" << endl;exit(EXIT_FAILURE);}float tempx, tempy;while (f >> tempx >> tempy) {x.push_back(tempx);y.push_back(tempy);}if (x.size() != y.size()) {cout << "数据文件内容有误" << endl;exit(EXIT_FAILURE);}// 创建A矩阵Eigen::MatrixXd A(x.size(), N + 1);for (unsigned int i = 0; i < x.size(); ++i) {// 遍历所有点for (int n = N, dex = 0; n >= 1; --n, ++dex) {// 遍历N到1阶A(i, dex) = pow(x[i], n);}A(i, N) = 1; //}// 创建B矩阵Eigen::MatrixXd B(y.size(), 1);for (unsigned int i = 0; i < y.size(); ++i) {B(i, 0) = y[i];}// 创建矩阵WEigen::MatrixXd W;W = (A.transpose() * A).inverse() * A.transpose() * B;// 打印W结果cout << W << endl;}
与设置的系数相吻合,但这是建立在噪声为0的情况下的结果。
#include <QCoreApplication>#include <QtCore>#include <QObject>#include <QVector>#include <Eigen/Core>#include <Eigen/Dense>#include <Eigen/Geometry>#include <Eigen/Eigenvalues>#include <iostream>#include <fstream>#include <string>#include <vector>#include <cmath>using namespace Eigen;using namespace std;bool polyfit(const QVector<double> &x, const QVector<double> &y, quint8 dim, QVector<double> &z){Q_ASSERT(x.size() == y.size());quint8 N = 3;// 设置是几次拟合if(dim < 0 && dim > 20){return false;}N = dim;// 创建A矩阵Eigen::MatrixXd A(x.size(), N + 1);for (unsigned int i = 0; i < x.size(); ++i) {// 遍历所有点for (int n = N, dex = 0; n >= 1; --n, ++dex) {// 遍历N到1阶A(i, dex) = pow(x[i], n);}A(i, N) = 1; //}// 创建B矩阵Eigen::MatrixXd B(y.size(), 1);for (unsigned int i = 0; i < y.size(); ++i) {B(i, 0) = y[i];}// 创建矩阵WEigen::MatrixXd W;W = (A.transpose() * A).inverse() * A.transpose() * B;// 打印W结果cout << W << endl;cout << "rows:"<<W.rows() << " cols:" << W.cols() << endl;double t;z.clear();for(int i = 0; i< W.rows()*W.cols(); i++){t = W(i,0);z.push_back(t);}return true;}int main(int argc, char *argv[]){QCoreApplication a(argc, argv);// 创建两个vectorQVector<double> x, y, z;x.push_back(1);x.push_back(2);x.push_back(3);x.push_back(4);x.push_back(5);x.push_back(6);y.push_back(10);y.push_back(13);y.push_back(15);y.push_back(21);y.push_back(22);y.push_back(28);bool ret = polyfit(x, y, 5, z);qDebug()<<"QVector z: "<<z;return a.exec();}
如果觉得《C++最小二乘拟合 (高阶最小二乘拟合)》对你有帮助,请点赞、收藏,并留下你的观点哦!