polyfit函数的C语言实现

void PolyfitCf(int n_poly,int Nwin_length,int Npoly,double * ypoly,double **fitcoef)
{
 int i,j,m;
 int nwin_length=Nwin_length;
 int poly_n=n_poly;
 int npoly=Npoly;

 double *x=NULL;//[nwin_length];
 if (x==NULL)
 {
  x=new double[nwin_length];
 }

 for (int ix=0;ix<nwin_length;ix++)
 {
  x[ix]=ix+1;
 }

 double **y=NULL;
 int iypoly=0;
 y=new double*[nwin_length*npoly];
 for (int ifc=0;ifc<npoly;ifc++)
 {
  y[ifc]=new double[nwin_length];
 }
 for (int ifcx=0;ifcx<npoly;ifcx++)
 {
  for (int ifcy=0;ifcy<nwin_length;ifcy++)
  {
   y[ifcx][ifcy]=ypoly[iypoly];
   iypoly++;
  }
 }

 double apoly[3];
 for (int ixy=0;ixy<npoly;ixy++)
 {
  polyfit(nwin_length,x,y[ixy],poly_n,apoly);
  fitcoef[ixy][0]=apoly[2];
  fitcoef[ixy][1]=apoly[1];
  fitcoef[ixy][2]=apoly[0];
 }
 
 delete [] x;
 x=NULL;

 for (int ifc=0;ifc<npoly;ifc++)
 {
  delete [] y[ifc];
  y[ifc]=NULL;

 }

 delete [] y;
 y=NULL;

}

//*==================polyfit(n,x,y,poly_n,a)===================*/
//*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/
//*=====n是数据个数 xy是数据值 poly_n是多项式的项数======*/
//*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/

void polyfit(int n,double x[],double y[],int poly_n,double a[])
{
 int i,j;
 double *tempx,*tempy,*sumxx,*sumxy,*ata;

 tempx=new double[n];
 sumxx=new double[poly_n*2+1];
 tempy=new double[n];
 sumxy=new double[poly_n+1];
 ata=new double[(poly_n+1)*(poly_n+1)];

 for (i=0;i<n;i++)
 {
  tempx[i]=1;
  tempy[i]=y[i];
 }
 for (i=0;i<2*poly_n+1;i++)
  for (sumxx[i]=0,j=0;j<n;j++)
  {
   sumxx[i]+=tempx[j];
   tempx[j]*=x[j];
  }
  for (i=0;i<poly_n+1;i++)
   for (sumxy[i]=0,j=0;j<n;j++)
   {
    sumxy[i]+=tempy[j];
    tempy[j]*=x[j];
   }
   for (i=0;i<poly_n+1;i++)
    for (j=0;j<poly_n+1;j++)
     ata[i*(poly_n+1)+j]=sumxx[i+j];
   gauss_solve(poly_n+1,ata,a,sumxy);

   delete [] tempx;
   tempx=NULL;
   delete [] sumxx;
   sumxx=NULL;
   delete [] tempy;
   tempy=NULL;
   delete [] sumxy;
   sumxy=NULL;
   delete [] ata;
   ata=NULL;
}

void gauss_solve(int n,double A[],double x[],double b[])
{
 int i,j,k,r;
 double max;
 for (k=0;k<n-1;k++)
 {
  max=fabs(A[k*n+k]); /*find maxmum*/
  r=k;
  for (i=k+1;i<n-1;i++)
   if (max<fabs(A[i*n+i]))
   {
    max=fabs(A[i*n+i]);
    r=i;
   }
   if (r!=k)
    for (i=0;i<n;i++) /*change array:A[k]&A[r] */
    {
     max=A[k*n+i];
     A[k*n+i]=A[r*n+i];
     A[r*n+i]=max;
    }
    max=b[k]; /*change array:b[k]&b[r] */
    b[k]=b[r];
    b[r]=max;
    for (i=k+1;i<n;i++)
    {
     for (j=k+1;j<n;j++)
      A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k];
     b[i]-=A[i*n+k]*b[k]/A[k*n+k];
    }
 }

 for (i=n-1;i>=0;x[i]/=A[i*n+i],i--)
  for (j=i+1,x[i]=b[i];j<n;j++)
   x[i]-=A[i*n+j]*x[j];
}