最小二乘法三次多项式曲线拟合 算法 C++ 实现

LP2010 2010-11-01 10:42:27
求最小二乘法三次多项式曲线拟合(就像matlab里面的三次多项式拟合那样) 算法 C++实现的,只需要求出三次多项式的系数即可
...全文
6235 15 打赏 收藏 转发到动态 举报
写回复
用AI写文章
15 条回复
切换为时间正序
请发表友善的回复…
发表回复
sxzero120 2012-12-12
  • 打赏
  • 举报
回复
真心希望楼主分享一下;可惜已经两年了,估计人去“楼”空了
wangyunfu1118 2012-05-07
  • 打赏
  • 举报
回复
lZ 可否分享一下
黎明露珠 2011-09-04
  • 打赏
  • 举报
回复
想学习下
wcongyuan 2011-08-02
  • 打赏
  • 举报
回复
那你最后是怎么实现的,可否分享一下
LP2010 2010-12-08
  • 打赏
  • 举报
回复
谢谢各位的意见,我好久没有登录了,问题我自己已经解决,用着还是可以的。
job82824 2010-11-16
  • 打赏
  • 举报
回复
以前弄的一个曲线拟合用的函数,是分类来分析的--有最小二乘法的拟合算法的,没啥子难点可言的,就是简单的矩阵对角化求解线性方程组的函数,你看着公式花点儿时间的也可以写出来的。

这个分类求解是不是太繁琐了?我想应该有更简练的方式,或者可以合并一下的?贴出来大家交流下

//.h
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//多项式拟合用函数

//数值计算
int GainParam4(double a[][4], double b[]); //求解4X4矩阵,返回解存在b中
int GainParam3(double a[][3], double b[]); //求解3X3矩阵, 返回解存在b
int GainParam2(double a[][2], double b[]); //求解2X2矩阵, 返回解存在b

//a、b分别是X Y坐标点集合,len是录入的坐标点数量.t是个中间过程用变量数组--长度由拟合阶数决定
BOOL LineFit3(double *a, double *b, double t[], int len); //使用3次二项式拟合曲线
BOOL LineFit2(double *a, double *b, double t[], int len); //使用2次二项式拟合曲线
BOOL LineFit(double *a, double *b, double t[], int len); //使用线性拟合曲线

//.cpp
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//多项式拟合用函数

int GainParam2(double a[][2], double b[])
{
double x, y;
int i;

y=a[0][0]*a[1][1]-a[0][1]*a[1][0];
if(!y)
{
return 1;
}

x=a[0][0];
if(!x)
{
if(!a[1][0])
{
return 2;
}
for(i=0;i<2;i++)
{
x=a[0][i];
a[0][i]=a[1][i];
a[1][i]=x;
}
y=b[0];
b[0]=b[1];
b[1]=y;
}
x=a[1][0]/a[0][0];
a[1][1]-=x*a[0][1];
b[1]-=b[0]*x;
x=a[0][1]/a[1][1];
b[0]-=b[1]*x;

b[0]/=a[0][0];
b[1]/=a[1][1];

return 0;

}

int GainParam3(double a[][3], double b[])
{
double x,y;
int i, j;
y=a[0][0]*(a[1][1]*a[2][2]-a[1][2]*a[2][1])-a[0][1]*(a[1][0]*a[2][2]-a[1][2]*a[2][0])+
a[0][2]*(a[1][0]*a[2][1]-a[1][1]*a[2][0]);

if(!y)
{
return 1;
}

//行列式对角化

//第一列挑选非零项
for(i=0;i<3;i++)
{
x=a[i][0];
if(x)break;
}
if(i==3)
{
return 2;
}
if(i>0)
{
j=i;
for(i=0;i<3;i++)
{
x=a[0][i];
a[0][i]=a[j][i];
a[j][i]=x;
}
y=b[0];
b[0]=b[j];
b[j]=y;
}
//销零
x=a[0][0];
y=a[1][0]/x;
b[1]-=b[0]*y;
for(i=1;i<3;i++)
a[1][i]-=a[0][i]*y;
y=a[2][0]/x;
b[2]-=b[0]*y;
for(i=1;i<3;i++)
a[2][i]-=a[0][i]*y;

//第二列挑选非零项
for(i=1;i<3;i++)
{
x=a[i][1];
if(x)break;
}
if(i==3)
{
return 3;
}
if(i>1)
{
j=i;
for(i=0;i<3;i++)
{
x=a[1][i];
a[1][i]=a[j][i];
a[j][i]=x;
}
y=b[1];
b[1]=b[j];
b[j]=y;
}

//销零
x=a[1][1];
y=a[0][1]/x;
b[0]-=b[1]*y;
a[0][2]-=a[1][2]*y;
y=a[2][1]/x;
b[2]-=b[1]*y;
a[2][2]-=a[1][2]*y;

if(!a[2][2])
{
return 4;
}
x=a[2][2];
for(i=0;i<2;i++)
{
y=a[i][2]/x;
b[i]-=b[2]*y;
}

//求解得到系数:
for(i=0;i<3;i++)
{
b[i]/=a[i][i];
}

return 0;
}


int GainParam4(double a[][4], double b[])
{
double x, y;
int i, k1;

x=0;
y=0;
x+=a[1][0]*(a[2][1]*a[3][2]-a[2][2]*a[3][1])+a[1][1]*(a[2][2]*a[3][1]-a[2][0]*a[3][2])+
a[1][2]*(a[2][0]*a[3][1]-a[2][1]*a[3][0]);
y+=-x*a[0][3];
x=a[1][1]*(a[2][2]*a[3][3]-a[2][3]*a[3][2])+a[1][2]*(a[2][3]*a[3][1]-a[2][1]*a[3][3])+
a[1][3]*(a[2][1]*a[3][2]-a[2][2]*a[3][1]);
y+=a[0][0]*x;
x=a[1][0]*(a[2][2]*a[3][3]-a[2][3]*a[3][2])-a[1][2]*(a[2][0]*a[3][3]-a[2][3]*a[3][0])+
a[1][3]*(a[2][0]*a[3][2]-a[2][2]*a[3][0]);
y+=-x*a[0][1];
x=a[1][0]*(a[2][2]*a[3][3]-a[2][3]*a[3][1])+a[1][1]*(a[2][3]*a[3][0]-a[2][0]*a[3][3])+
a[1][3]*(a[2][0]*a[3][1]-a[2][1]*a[3][0]);
y+=x*a[0][2];

if(!y)
{
return 1;
}

//行列式对角化

//第一列挑选非零项
for(i=0;i<4;i++)
{
x=a[i][0];
if(x)break;
}
if(i==4)
{
return 2;
}
if(i>0)
{
k1=i;
for(i=0;i<4;i++)
{
x=a[0][i];
a[0][i]=a[k1][i];
a[k1][i]=x;
}
y=b[0];
b[0]=b[k1];
b[k1]=y;
}


//销零
x=a[0][0];
y=a[1][0]/x;
b[1]-=b[0]*y;
for(i=1;i<4;i++)
a[1][i]-=a[0][i]*y;
y=a[2][0]/x;
b[2]-=b[0]*y;
for(i=1;i<4;i++)
a[2][i]-=a[0][i]*y;
y=a[3][0]/x;
b[3]-=b[0]*y;
for(i=1;i<4;i++)
a[3][i]-=a[0][i]*y;


//第二列挑选非零项
for(i=1;i<4;i++)
{
x=a[i][1];
if(x)break;
}
if(i==4)
{
return 3;
}
if(i>1)
{
k1=i;
for(i=0;i<4;i++)
{
x=a[1][i];
a[1][i]=a[k1][i];
a[k1][i]=x;
}
y=b[1];
b[1]=b[k1];
b[k1]=y;
}

//销零
x=a[1][1];
y=a[0][1]/x;
b[0]-=b[1]*y;
for(i=2;i<4;i++)
a[0][i]-=a[1][i]*y;
y=a[2][1]/x;
b[2]-=b[1]*y;
for(i=2;i<4;i++)
a[2][i]-=a[0][i]*y;
y=a[3][1]/x;
b[3]-=b[1]*y;
for(i=2;i<4;i++)
a[3][i]-=a[0][i]*y;

//第三列
for(i=2;i<4;i++)
{
x=a[i][2];
if(x)break;
}
if(i==4)
{
return 4;
}
if(i>2)
{
k1=i;
for(i=0;i<4;i++)
{
x=a[2][i];
a[2][i]=a[k1][i];
a[k1][i]=x;
}
y=b[2];
b[2]=b[k1];
b[k1]=y;
}
//销零
x=a[2][2];
y=a[0][2]/x;
b[0]-=b[2]*y;
a[0][3]-=a[2][3]*y;
y=a[1][2]/x;
b[1]-=b[2]*y;
a[1][3]-=a[2][3]*y;
y=a[3][2]/x;
b[3]-=b[2]*y;
a[3][3]-=a[2][3]*y;

//第四列
if(a[3][3]==0)
{
return 5;
}
x=a[3][3];
for(i=0;i<3;i++)
{
y=a[i][3]/a[3][3];
b[i]-=b[3]*y;
}

//求解得到系数:
for(i=0;i<4;i++)
{
b[i]/=a[i][i];
}

return 0;
}


//累加矩阵

BOOL LineFit(double *a, double *b, double t[], int len)
{
int i;
double x,y;
double m[2][2]={0};
double n[2]={0};

for(i=0;i<len;i++)
{
x=a[i];
y=b[i];
m[0][0]+=x;
m[0][1]+=1;
m[1][0]+=x*x;
n[0]+=y;
n[1]+=x*y;
}
m[1][1]=m[0][0];
t[0]=n[0];
t[1]=n[1];

i=GainParam2(m,t);
if(!i)
return TRUE;
else
{
if(i==1)
AfxMessageBox(_T("线性拟合失败--2维数组行列式为零"));
else
AfxMessageBox(_T("线性拟合失败--2维数组行列式为零"));
return FALSE;
}

}
BOOL LineFit2(double *a, double *b, double t[], int len)
{
int i;
double m[3][3]={0};
double n[3]={0};
double x,y,z;


for(i=0;i<len;i++)
{
x=a[i];
y=b[i];
z=x*x;
m[0][0]+=z;
m[0][1]+=x;
m[0][2]+=1;
m[1][0]+=z*x;
m[2][0]+=z*z;
n[0]+=y;
n[1]+=x*y;
n[2]+=z*y;
}
m[1][1]=m[0][0];
m[1][2]=m[0][1];
m[2][1]=m[1][0];
m[2][2]=m[0][0];
for(i=0;i<3;i++)
{
t[i]=n[i];
}

i=GainParam3(m, t);
if(!i)
return TRUE;
else
{
if(i==1)
AfxMessageBox(_T("系数行列式值为零,二次拟合失败"));
else if(i==2)
AfxMessageBox(_T("矩阵第一列全为零"));
else if(i==3)
AfxMessageBox(_T("矩阵第二列迹为零"));
else
AfxMessageBox(_T("矩阵第三列迹为零"));
}
return FALSE;
}

BOOL LineFit3(double *a, double *b, double t[], int len)
{
int i;
double m[4][4],n[4], x, y, z;

for(i=0;i<4;i++)
{
m[i][0]=m[i][1]=0.0;
m[i][2]=m[i][3]=0.0;
n[i]=0.0;
}
for(i=0;i<len;i++)
{
x=a[i];
y=b[i];
z=x*x*x;
m[0][0]+=z;
m[0][1]+=x*x;
m[0][2]+=x;
m[0][3]+=1;
n[0]+=y;
n[1]+=x*y;
n[2]+=x*x*y;
n[3]+=z*y;
m[1][0]+=z*x;
m[2][0]+=z*x*x;
m[3][0]+=z*z;
}
m[1][1]=m[0][0];
m[1][2]=m[0][1];
m[1][3]=m[0][2];
m[2][1]=m[1][0];
m[2][2]=m[1][1];
m[2][3]=m[0][1];
m[3][1]=m[2][0];
m[3][2]=m[1][0];
m[3][3]=m[0][0];


for(i=0;i<4;i++)
{
t[i]=n[i];
}

i=GainParam4(m,t);
if(!i)
return TRUE;
else
{
if(i==1)
AfxMessageBox(_T("系数行列式为零, 无解或无穷解"));
else if(i==2)
AfxMessageBox(_T("矩阵第一列全为零"));
else if(i==3)
AfxMessageBox(_T("矩阵第二列迹为零"));
else if(i==4)
AfxMessageBox(_T("矩阵第三列迹为零"));
else
AfxMessageBox(_T("无定解,第四列迹为零"));
}
return FALSE;
}
hastings 2010-11-16
  • 打赏
  • 举报
回复
上面的sztCiShu说明有误,应该是:
sztCiShu代表(sztCiShu-1)次多项式,按链接里的曲线拟合例子,sztCiShu为4
hastings 2010-11-16
  • 打赏
  • 举报
回复
http://www.matlabfan.com/thread-850-1-1.html
	typedef std::vector<std::pair<double,double> >  DoubleVec ;
DoubleVec fltq;//fltq储存采样数据
//略去采样数据存入fltq的代码
int n;
Matrix Y(j,1);//j行,1列的矩阵(j为采样数据对数,按链接里的曲线拟合例子,j为5)
Matrix X(j,sztCiShu);//j行,sztCiShu列的矩阵(sztCiShu代表sztCiShu次多项式,按链接里的曲线拟合例子,sztCiShu为3)
for(i=0;i<j;++i)
{
Y(i,0)=fltq[i].second;//first储存采样数据x坐标值,second储存采样数据y坐标值
for(n=0;n<sztCiShu;++n)
X(i,n)=pow(fltq[i].first,sztCiShu-n-1);
}
Matrix Xtrans(X.transpose());//矩阵X转置,X本身不变,返回值为X转置后产生的新矩阵,存入Xtrans
Matrix T(Xtrans.multiply(X));//Xtrans乘与X,Xtrans本身不变,返回值为相乘的结果,存入T
T.invert();//T求逆
Matrix A(T.multiply(Xtrans).multiply(Y));//A = T * Xtrans * Y (A为sztCiShu行,1列的矩阵)
//最后的A第一列(也是唯一的一列)即为拟合结果多项式的系数,正确的话,为 -0.1917 31.5821 -60.3262 35.3400
fandh 2010-11-16
  • 打赏
  • 举报
回复
比较成熟的算法,一般别人都不公布的!
hellosunking 2010-11-16
  • 打赏
  • 举报
回复
楼上的算法好像不稳定啊,同样的x,y,在m不同的时候差距会很大,即使m足够大,很高次项都小到可以忽略但是低次的系数差距很大,这个是什么原因呢?
zgl7903 2010-11-02
  • 打赏
  • 举报
回复

/******************************************
//参考 《常用算法程序集 (C语言描述 第三版)》

//最小二乘法
//x[n] y[n] 已知输入
//n输入点个数
//a[m] 返回m-1次拟合多项式的m个系数
//m 拟合多项式的项数,即拟合多项式的最高次为m-1
//dt[3] dt[0]返回拟合多项式与各数据点误差的平方和,
dt[1]返回拟合多项式与各数据点误差的绝对值之和
dt[2]返回拟合多项式与各数据点误差的绝对值的最大值

//
//拟合多项式的输出
//Y(x) = a0 + a1(x-X) + a2(x-X)^2 + …… am(x-X)^m
// 其中X为已知点x的平均值
******************************************/

#include "math.h"
void pir1(x,y,n,a,m,dt)
int n,m;
double x[],y[],a[],dt[];
{
int i,j,k;
double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
for (i=0; i<=m-1; i++) a[i]=0.0;
if (m>n) m=n;
if (m>20) m=20;
z=0.0;
for (i=0; i<=n-1; i++) z=z+x[i]/(1.0*n);
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
for (i=0; i<=n-1; i++)
{ p=p+(x[i]-z); c=c+y[i];}
c=c/d1; p=p/d1;
a[0]=c*b[0];
if (m>1)
{ t[1]=1.0; t[0]=-p;
d2=0.0; c=0.0; g=0.0;
for (i=0; i<=n-1; i++)
{ q=x[i]-z-p; d2=d2+q*q;
c=c+y[i]*q;
g=g+(x[i]-z)*q*q;
}
c=c/d2; p=g/d2; q=d2/d1;
d1=d2;
a[1]=c*t[1]; a[0]=c*t[0]+a[0];
}
for (j=2; j<=m-1; j++)
{ s[j]=t[j-1];
s[j-1]=-p*t[j-1]+t[j-2];
if (j>=3)
for (k=j-2; k>=1; k--)
s[k]=-p*t[k]+t[k-1]-q*b[k];
s[0]=-p*t[0]-q*b[0];
d2=0.0; c=0.0; g=0.0;
for (i=0; i<=n-1; i++)
{ q=s[j];
for (k=j-1; k>=0; k--)
q=q*(x[i]-z)+s[k];
d2=d2+q*q; c=c+y[i]*q;
g=g+(x[i]-z)*q*q;
}
c=c/d2; p=g/d2; q=d2/d1;
d1=d2;
a[j]=c*s[j]; t[j]=s[j];
for (k=j-1; k>=0; k--)
{ a[k]=c*s[k]+a[k];
b[k]=t[k]; t[k]=s[k];
}
}
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
for (i=0; i<=n-1; i++)
{ q=a[m-1];
for (k=m-2; k>=0; k--)
q=a[k]+q*(x[i]-z);
p=q-y[i];
if (fabs(p)>dt[2]) dt[2]=fabs(p);
dt[0]=dt[0]+p*p;
dt[1]=dt[1]+fabs(p);
}
return;
}
wuhuwy 2010-11-02
  • 打赏
  • 举报
回复
到pudn看看,我记得有源代码
LP2010 2010-11-01
  • 打赏
  • 举报
回复
自己再顶一个,麻烦高手们说详细点,做好能有实现的源码
hastings 2010-11-01
  • 打赏
  • 举报
回复
关键是弄个矩阵类~~

19,468

社区成员

发帖
与我相关
我的任务
社区描述
VC/MFC 图形处理/算法
社区管理员
  • 图形处理/算法社区
加入社区
  • 近7日
  • 近30日
  • 至今
社区公告
暂无公告

试试用AI创作助手写篇文章吧