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];
}
