对 n 个三维坐标（x,y,z) ， 拟合成这样一个函数
Z(x,y)= aa(i,j)*x^i*y^j, i 从0-->p ,j 从0-->q 求和。
具体的看里面的readme，说的很详细- function A=leftmatrix(x,p,y,q)
- % A*a=B a 即为系数列矩阵
- % A为左边(p-1)(q-1) 乘 (p-1)(q-1) 的矩阵
- % x,y 为长度一样的列矩阵 也就是给定离散点的x,y坐标
- % p,q为拟合的函数中x,y的指数
- m=length(x);
- if (nargin~=4) & (m~=length(y)), error('error check check!'); end
- A_length=p*q; % A 为p*q阶的方阵
- A=zeros(A_length,A_length); % 赋值0
- for i=1 : p*q
- for j= 1 : p*q
- x_z=quotient(j-1,q)+quotient(i-1,q); % x 的指数 quotient为求商
- y_z=mod(j-1,q)+mod(i-1,q); % y 的指数
- A(i,j)=qiuhe(x,x_z,y,y_z);
- end
- end
- function he=qiuhe(x,p,y,q,z)
- % he x^p*y^q 从1－>m的和
- % x,y 行向量 长度相同; p,q 为x,y系数
- %x=[2 3]; y=[ 3 4]; p=2; q= 2;
- m=length(x);
- if (nargin<4 )&(m~=length(y)) %输入量至少为四，x,y行向量长度必需一样
- error('error check check!');
- end
- if nargin==4, z=ones(m,1); end %没有 z , 默认为单位行向量
- he=0;
- for i=1:m
- he=he+x(i)^p * y(i)^q*z(i); % 1-->m 求和
- end
- function sh=quotient(x,y)
- % sh 为 x/y 的商
- sh=(x-mod(x,y))/y;
- function B=rightmatrix(x,p,y,q,z)
- % A*a=B
- % B 为一个列向量 长为p*q
- % x y z 为点的坐标 , p q 为x y指数
- if nargin~=5, error('error check check! rightmatrix'); end
- B=zeros(p*q,1);
- for i=1 : p*q
- x_z=quotient(i-1,q); y_z=mod(i-1,q); B(i,1)=qiuhe(x,x_z,y,y_z,z);
- end
- function ff=main(x,p,y,q,z,xx,yy)
- % x y z 坐标向量 长度要一样
- % p ,q 为拟合函数中x,y 的系数
- % xx yy 为 需要拟合的数据 给出(x,y) 坐标 求z
- A=leftmatrix (x,p,y,q); % A*a_n=B
- B=rightmatrix(x,p,y,q,z);
- %a_n=inv(A)*B;
- a_n=A\B; % 求a_n inv(A)*B 效果不好 （存疑）
- for i=1 : length(a_n) % 把长为p*q 的 a_n 列向量 转还 成p x q 的 aa 矩阵
- ii=quotient(i-1,q)+1; % quotient求商
- jj=mod(i-1,q)+1; aa(ii,jj)=a_n(i,1);
- end
-
- ff=0; % ff 是 xx,yy 带入所拟合的函数 求出 z
- for i=1 : p % 函数为 aa(i,j)*x^i*y^j (i=0...p ,j=0...q)
- for j=1 : q % aa 为系数 p x q 的矩阵
- ff=ff+aa(i,j) * xx^(i-1) * yy^(j-1);
- end
- end
[ 本帖最后由 ChaChing 于 2009-4-7 13:12 编辑 ]
 readme.doc |
