www.pudn.com > LeastSquare.rar > LeastSquare.h

#define max1 10 ///拟合多项式的最高幂数或线性方程组的最大阶数 
#define max2 1000///最小二乘拟合时的最大离散点数 
////最小二乘曲线拟合 
BOOL CDataProcess::LeastSquare(double x[], double y[], int m, int n, double a[]) 
{ 
	//x,y -- X,Y两轴的坐标 
	//m 多项式的最高幂 
	//n 离散点数 
	//a 多项式的系数 
	int i,j,k,r;  
	double f[MaxExp*2]; 
	double A[MaxExp][MaxExp]; 
	double Y[MaxExp]; 
	double s=0; 
	double max3,temp,xt; 
	 
	if(m>MaxExp) 
	{ 
		printf("多项式幂数不能大于%d",MaxExp); 
		exit(0); 
	} 
	if(n>MaxFitPt) 
	{ 
		printf("输入点数不能大于%d",MaxFitPt); 
		exit(0); 
	} 
	if(m>n) 
	{ 
		m=n;//多项式的最高幂不能超过离散点数 
	} 
	for(k=0;k<=2*m;k++)//计算x的次幂 
	{ 
		f[k]=0; 
		for(i=0;iMaxExp) 
			{ 
				max3=fabs(A[i][k]); 
				r=i; 
			  
			 
			} 
		} 
		if(r!=k)//如果r和k不相等就说明需要进行两行交换 
		{ 
			for(int j=k;j<=m;j++) 
			{ 
			   temp=A[k][j];//下面三句是交换A数组 
			   A[k][j]=A[r][j]; 
			   A[r][j]=temp; 
			} 
			temp=Y[k];//下面三句是交换Y数组 
			Y[k]=Y[r]; 
			Y[r]=temp; 
		} 
		for(i=k+1;i<=m;i++) 
		{ 
			xt=A[i][k]/A[k][k];//交换后置零 
			Y[i]=Y[i]-xt*Y[k]; 
			for(j=k+1;j<=m;j++) 
				A[i][j]=A[i][j]-xt*A[k][j]; 
			A[i][k]=0; 
		} 
	} 
 
	a[m]=Y[m]/A[m][m];//求的是最后一个x值 
	for (i=m-1;i>=0;i--)//求前面的几个x值 
	{ 
		temp=0; 
		for (j=i+1;j<=m;j++) 
		{ 
			temp=temp+A[i][j]*a[j]; 
		} 
		a[i]=(Y[i]-temp)/A[i][i]; 
	} 
	return TRUE; 
} 
 
BOOL CDataProcess::Gauss(double A[][MaxExp], double Y[], int m, double x[]) 
{ 
	double max3,temp,xt; 
	int r,i,j; 
	if(m>MaxExp) 
	{ 
		printf("矩阵阶数不能大于%d",MaxExp); 
		exit(0); 
	} 
 
	for(int k=0;k<=m;k++) 
	{ 
		r=k; 
		max3=fabs(A[k][k]); 
		for( i=k;i<=m;i++) 
		{ 
			if(fabs(A[i][k])>MaxExp) 
			{ 
				max3=fabs(A[i][k]); 
				r=i; 
			  
			 
			} 
		} 
		if(r!=k)//如果r和k不相等就说明需要进行两行交换 
		{ 
			for(int j=k;j<=m;j++) 
			{ 
			   temp=A[k][j];//下面三句是交换A数组 
			   A[k][j]=A[r][j]; 
			   A[r][j]=temp; 
			} 
			temp=Y[k];//下面三句是交换Y数组 
			Y[k]=Y[r]; 
			Y[r]=temp; 
		} 
		for(i=k+1;i<=m;i++) 
		{ 
			xt=A[i][k]/A[k][k];//交换后置零 
			Y[i]=Y[i]-xt*Y[k]; 
			for(j=k+1;j<=m;j++) 
				A[i][j]=A[i][j]-xt*A[k][j]; 
			A[i][k]=0; 
		} 
	} 
 
	x[m]=Y[m]/A[m][m];//求的是最后一个x值 
	for (i=m-1;i>=0;i--)//求前面的几个x值 
	{ 
		temp=0; 
		for (j=i+1;j<=m;j++) 
		{ 
			temp=temp+A[i][j]*x[j]; 
		} 
		x[i]=(Y[i]-temp)/A[i][i]; 
	} 
	return TRUE; 
}