www.pudn.com > pcir.rar > pcir.cpp

//最小二乘法 
void pcir(double X[],double Y[],double A[],int N,int M,double &DT1,double &DT2,double &DT3) 
{ 
	// X[] Y[] 使需要拟合的点，A[]是拟合曲线的系数， 
	double S[21],T[21],B[21]; 
	double Z,D1,P,C,D2,G,Q,DT; 
	int i,j,jk; 
    for(i=1;i<=M;i++)//	DO 5 i=1,M 
	{ 
		A[i]=0.0; 
	} 
 
	if (M>=N) M=N; 
	if (M>=20) M=20; 
	Z=0.0; 
/*	for(i=1;i<=N;i++)//DO 10 i=1,N 
	{ 
		Z=Z+X[i]/N; 
	}//10	Z=Z+X(i)/N*/ 
	B[1]=1.0; 
	D1=N; 
	P=0.0; 
	C=0.0; 
	for(i=1;i<=N;i++) 
	{ 
		P=P+(X[i]-Z); 
		C=C+Y[i]; 
	} 
	C=C/D1; 
	P=P/D1; 
	A[1]=C*B[1]; 
	if(M>1) 
	{ 
	  T[2]=1.0; 
	  T[1]=-P; 
	  D2=0.0; 
	  C=0.0; 
	  G=0.0; 
	  for(i=1;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[2]=C*T[2]; 
	  A[1]=C*T[1]+A[1]; 
	} 
	for(j=3;j<=M;j++) 
	{ 
	  S[j]=T[j-1]; 
	  S[j-1]=-P*T[j-1]+T[j-2]; 
	  if(j>=4) 
		  for(jk=j-2;jk>=2;jk--) 
            S[jk]=-P*T[jk]+T[jk-1]-Q*B[jk]; 
	  S[1]=-P*T[1]-Q*B[1]; 
	  D2=0.0; 
	  C=0.0; 
	  G=0.0; 
	  for(i=1;i<=N;i++) 
	  { 
	    Q=S[j]; 
		for(jk=j-1;jk>=1;jk--) 
		{     
			Q=Q*(X[i]-Z)+S[jk]; 
		} 
	    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(jk=j-1;jk>=1;jk--) 
	  { 
	    A[jk]=C*S[jk]+A[jk]; 
	    B[jk]=T[jk]; 
	    T[jk]=S[jk]; 
	  } 
	} 
	DT1=0.0; 
	DT2=0.0; 
	DT3=0.0; 
	 
	for(i=1;i<=N;i++) 
	{ 
	  Q=A[M]; 
	  for(jk=M-1;jk>=1;jk--) 
	  { 
		  Q=Q*(X[i]-Z)+A[jk]; 
	  } 
	  DT=Q-Y[i]; 
	  if (fabs(DT)>DT3) DT3=fabs(DT); 
	  DT1=DT1+DT*DT; 
	  DT2=DT2+fabs(DT); 
	} 
 
//AfxMessageBox("hello"); 
 
}