最小2乘发曲线拟和.txt

最小二乘法拟和曲线

1.
//最小二乘法直线拟合
BOOL CalculateLineKB(CFoldPointList *m_FoldList,double &k,double &b)
{
 //最小二乘法直线拟合
 //m_FoldList为关键点(x,y)的链表
 //拟合直线方程(Y=kX+b)


 if(m_FoldList==NULL)return FALSE;
 long lCount=m_FoldList->GetCount();
 if(lCount<2)return FALSE;
 CFoldPoint *pFold;
 double mX,mY,mXX,mXY,n;
 mX=mY=mXX=mXY=0;
 n=lCount;
 POSITION pos=m_FoldList->GetHeadPosition();
 while(pos != NULL)
 {
  pFold=m_FoldList->GetNext(pos);
  mX+=pFold->X;
  mY+=pFold->Y;
  mXX+=pFold->X*pFold->X;
  mXY+=pFold->X*pFold->Y;
 }
 if(mX*mX-mXX*n==0)return FALSE;
 k=(mY*mX-mXY*n)/(mX*mX-mXX*n);
 b=(mY-mX*k)/n;
 return TRUE;
}




2.
//最小二乘法曲线拟合
typedef CArray<double,double>CDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A)
{
 //X,Y --  X,Y两轴的坐标
 //M   --  结果变量组数
 //N   --  采样数目
 //A   --  结果参数

 register long i,j,k;
 double Z,D1,D2,C,P,G,Q;
 CDoubleArray B,T,S;
 B.SetSize(N);
 T.SetSize(N);
 S.SetSize(N);
 if(M>N)M=N;
 for(i=0;i<M;i++)
  (*A)[i]=0;
 Z=0;
 B[0]=1;
 D1=N;
 P=0;
 C=0;
 for(i=0;i<N;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;
  T[0]=-P;
  D2=0;
  C=0;
  G=0;
  for(i=0;i<N;i++)
  {
   Q=(*X)[i]-Z-P;
   D2=D2+Q*Q;
   C=(*Y)[i]*Q+C;
   G=((*X)[i]-Z)*Q*Q+G;
  }
  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;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;
  C=0;
  G=0;
  for(i=0;i<N;i++)
  {
   Q=S[j];
   for(k=j-1;k>=0;k--)
    Q=Q*((*X)[i]-Z)+S[k];
   D2=D2+Q*Q;
   C=(*Y)[i]*Q+C;
   G=((*X)[i]-Z)*Q*Q+G;
  }
  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];
  }
 }
 return TRUE;
}