多数据点拟合曲线,最小二乘法,矩阵

最新推荐文章于 2024-07-15 11:02:18 发布

搬砖小码农嘿嘿嘿

最新推荐文章于 2024-07-15 11:02:18 发布

阅读量1.1k

点赞数

分类专栏： C++ 数据处理文章标签： c++

本文链接：https://blog.csdn.net/qq_42985268/article/details/120256911

版权

C++ 同时被 2 个专栏收录

8 篇文章 0 订阅

订阅专栏

数据处理

1 篇文章 0 订阅

订阅专栏

本文介绍了使用C++实现的多数据点拟合曲线算法，通过范德蒙化简方法求解一元四次方程，展示了如何通过多项式拟合和系数计算来拟合给定的数据点。重点讲解了`calmatchline`函数中的范德蒙矩阵操作和费法拉利求解过程。

摘要由CSDN通过智能技术生成

#include <vector>
#include <qdebug.h>
//多数据点拟合曲线
bool calmatchline(const std::vector<int>&xValue ,const vector<float>&yValues)
{
	using namespace std;
	int cloth = 4;//次幂
	int matlength = cloth +1;
	std::vector<float> xvalues;
	//模拟数值
	xvalues.push_back(-65);
	xvalues.push_back(-35);
	xvalues.push_back(-5);
	xvalues.push_back(0);
	xvalues.push_back(5);
	xvalues.push_back(10);
	xvalues.push_back(15);
	xvalues.push_back(20);
	xvalues.push_back(25);
	xvalues.push_back(37);
	vector<float> yvalues;
	yvalues.push_back(0.26749875);
	yvalues.push_back(0.40583125);
	yvalues.push_back(0.85247625);
	yvalues.push_back(1.001101299);
	yvalues.push_back(0.97136125);
	yvalues.push_back(1.1904425);
	yvalues.push_back(1.017414286);
	yvalues.push_back(0.839881818);
	yvalues.push_back(0.759257143);
	yvalues.push_back(0.511851948);
	CvMat  *A=cvCreateMat(matlength,matlength,CV_32FC1);
    CvMat  *X=cvCreateMat(matlength,1,CV_32FC1);
    CvMat  *b=cvCreateMat(matlength,1,CV_32FC1);

	//
	float onevalue = 0;
	for (int i = 0; i < matlength; i++)
	{
		for (int j = 0; j < matlength; j++)
		{
			onevalue = sumxvalue(xvalues,i+j);
			cvmSet(A,i,j,onevalue);
		}
	}

	for (int i = 0; i < matlength; i++)
	{
		onevalue = sumxyvalue(xvalues,yvalues,i);
		cvmSet(b,i,0,onevalue);
	}

	cvSolve(A,b,X,CV_LU);


	vector<float> value;
	for (int i = 0; i < matlength; i++)
	{
		value.push_back(cvmGet(X,i,0));
		qDebug()<<QString::fromLocal8Bit("系数%1:").arg(QString::number(i))<<QString::number(value.at(i))+" ";
	}
	float *y = new float[20];
	float *ptr = y;
	for (int i = 0; i < yvalues.size(); i++)
	{
		float tmp = 0;
		for (int j = 0; j < matlength; j++)
		{
			tmp +=  pow(xvalues.at(i),j)*value.at(j);
		}
		
		*ptr = tmp;
		ptr++;
	}
	ptr = y;
	for (int i = 0; i < xvalues.size(); i++)
	{
		qDebug()<<"x:"<<QString::number(xvalues.at(i))<<"  y:"<<QString::number(*ptr);
		ptr++;
	}
	std::complex<double> x[4];
	x[4] = Ferrari(x,value[1],value[2],value[3],value[4],value[0]-0.5); 

	
	std::cout<<"root1: "<<x[0]<<std::endl
                     <<"root2: "<<x[1]<<std::endl
                     <<"root3: "<<x[2]<<std::endl
                     <<"root4: "<<x[3]<<std::endl;
	return true;

	//范德蒙化简方式
	CvMat  *A1=cvCreateMat(xvalues.size(),matlength,CV_32FC1);
	CvMat  *tA1=cvCreateMat(matlength,xvalues.size(),CV_32FC1);
	CvMat  *sumA=cvCreateMat(xvalues.size(),xvalues.size(),CV_32FC1);
	CvMat  *nA=cvCreateMat(xvalues.size(),xvalues.size(),CV_32FC1);
    CvMat  *X1=cvCreateMat(xvalues.size(),1,CV_32FC1);
    CvMat  *b1=cvCreateMat(xvalues.size(),1,CV_32FC1);
	//
	//Mat A1;
	//Mat tA1;
	//Mat X1;
	//Mat b1;
	for (int i = 0; i < xvalues.size(); i++)
	{
		for (int j = 0; j < matlength; j++)
		{
			onevalue = pow(xvalues.at(i),j);
			cvmSet(A1,i,j,onevalue);
		}
	}
	for (int i = 0; i < xvalues.size(); i++)
	{
		onevalue = yvalues.at(i);
		cvmSet(b1,i,0,onevalue);
	}
	cvTranspose(A1,tA1);
	cvmMul(A1,tA1,sumA);

	cvInvert(sumA,nA);
	
	CvMat  *sum2=cvCreateMat(xvalues.size(),matlength,CV_32FC1);
	cvMul(nA,tA1,sum2);//这里有异常,下次再处理
	//invert(A1*tA1,nA,DECOMP_SVD);
	float tmpvalue;
	for (int i = 0; i < matlength; i++)
	{
		tmpvalue = cvmGet(X,i,0);
		qDebug()<<QString::number(i)+":";
		qDebug()<<QString::number(tmpvalue)+" ";
	}
	return true;
}

float sumxvalue(std::vector<float> value,int cloth)
{
	float sumvalue = 0;
	for (int i = 0; i < value.size(); i++)
	{
		sumvalue += pow(value.at(i),cloth);
	}
	return sumvalue;
}

float sumxyvalue(std::vector<float> xvalue,std::vector<float> yvalue,int cloth)
{
	float sumvalue = 0;
	for (int i = 0; i < xvalue.size(); i++)
	{
		sumvalue += pow(xvalue.at(i),cloth) * yvalue.at(i);
	}
	return sumvalue;
}

std::complex<double>  sqrtn(const std::complex<double>&x,double n)
{
    double r = hypot(x.real(),x.imag()); //模
    if(r > 0.0)
    {
        double a = atan2(x.imag(),x.real()); //辐角
        n = 1.0 / n;
        r = pow(r,n);
        a *= n;
        return std::complex<double>(r * cos(a),r * sin(a));
    }
    return std::complex<double>();
}

//费法拉求解一元四次方程
std::complex<double> Ferrari(std::complex<double> x[4]
,std::complex<double> a
,std::complex<double> b
,std::complex<double> c
,std::complex<double> d
,std::complex<double> e)
{
    a = 1.0 / a;
    b *= a;
    c *= a;
    d *= a;
    e *= a;
    std::complex<double> P = (c * c + 12.0 * e - 3.0 * b * d) / 9.0;
    std::complex<double> Q = (27.0 * d * d + 2.0 * c * c * c + 27.0 * b * b * e - 72.0 * c * e - 9.0 * b * c * d) / 54.0;
    std::complex<double> D = sqrtn(Q * Q - P * P * P,2.0);
    std::complex<double> u = Q + D;
    std::complex<double> v = Q - D;
    if(v.real() * v.real() + v.imag() * v.imag() > u.real() * u.real() + u.imag() * u.imag())
    {
        u = sqrtn(v,3.0);
    }
    else
    {
        u = sqrtn(u,3.0);
    }
    std::complex<double> y;
    if(u.real() * u.real() + u.imag() * u.imag() > 0.0)
    {
        v = P / u;
        std::complex<double> o1(-0.5,+0.86602540378443864676372317075294);
        std::complex<double> o2(-0.5,-0.86602540378443864676372317075294);
        std::complex<double>&yMax = x[0];
        double m2 = 0.0;
        double m2Max = 0.0;
        int iMax = -1;
        for(int i = 0;i < 3;++i)
        {
            y = u + v + c / 3.0;
            u *= o1;
            v *= o2;
            a = b * b + 4.0 * (y - c);
            m2 = a.real() * a.real() + a.imag() * a.imag();
            if(0 == i || m2Max < m2)
            {
                m2Max = m2;
                yMax = y;
                iMax = i;
            }
        }
        y = yMax;
    }
    else
    {//一元三次方程，三重根
        y = c / 3.0;
    }
    std::complex<double> m = sqrtn(b * b + 4.0 * (y - c),2.0);
    if(m.real() * m.real() + m.imag() * m.imag() >= DBL_MIN)
    {
        std::complex<double> n = (b * y - 2.0 * d) / m;
        a = sqrtn((b + m) * (b + m) - 8.0 * (y + n),2.0);
        x[0] = (-(b + m) + a) / 4.0;
        x[1] = (-(b + m) - a) / 4.0;
        a = sqrtn((b - m) * (b - m) - 8.0 * (y - n),2.0);
        x[2] = (-(b - m) + a) / 4.0;
        x[3] = (-(b - m) - a) / 4.0;
    }
    else
    {
        a = sqrtn(b * b - 8.0 * y,2.0);
        x[0] =
        x[1] = (-b + a) / 4.0;
        x[2] =
        x[3] = (-b - a) / 4.0;
    }
return x[4];
}

转载自http://blog.csdn.net/pipisorry/article/details/49804441

https://blog.csdn.net/i_chaoren/article/details/79822574

https://blog.csdn.net/qq_38222947/article/details/118212408?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2_{aggregatepage}first_rank_ecpm_v1~rank_aggregation-1-118212408.pc_agg_rank_aggregation&utm_term=%E5%A4%9A%E9%A1%B9%E5%BC%8F%E6%8B%9F%E5%90%88%E5%8E%9F%E7%90%86&spm=1000.2123.3001.4430