MATLAB二维多项式拟合

itmwk

函数如下：
function z = polyval2d(c,x,y)
% z = POLYVAL2D(V,sx,sy) Two dimensional polynomial evaluation.
%
% If V is a matrix whose elements are the coefficients of a
% polynomial function of 2 variables, then POLYVAL2D(V,sx,sy)
% is the value of the polynomial evaluated at [sx,sy].  Row
% numbers in V correspond to powers of x, while column numbers
% in V correspond to powers of y. If sx and sy are matrices or
% vectors,the polynomial function is evaluated at all points
% in [sx, sy].
%
% If V is one dimensional, POLYVAL2D returns the same result as
% POLYVAL.
%
% Use POLYFIT2D to generate appropriate polynomial matrices from
% f(x,y) data using a least squares method.


% Perry W. Stout  June 28, 1995
% 4829 Rockland Way
% Fair Oaks, CA  95628
% (916) 966-0236
% Based on the Matlab function POLYVAL.

%        Polynomial evaluation c(x,y) is implemented using Horner's slick
% method.  Note use of the filter function to speed evaluation when
% the ordered pair [sx,sy] is single valued.

if nargin==2
  y=ones(x); end

if any(size(x) ~= size(y))
error('x and y must have the same dimensions.')
end

[m,n] = size(x);
[rown,coln]= size(c);

if (m+n) == 2
        % Use the built-in filter function when [sx,sy] is single valued
% to implement Hoerner's method.

  z= 0;
for i1= 1:coln
         ccol= c(:,coln+1-i1);
  w = filter(1,[1 -x],ccol);
         w = w(rown)*(y^(i1-1));
  z= z+w;
end
return
end % Of the scalar computation

% Do general case where X and Y are arrays
z = zeros(m,n);
for i1=1:coln
   ccol= c(:,coln+1-i1);
   w= zeros(m,n);
  for j1=1:rown
           w = x.*w + ccol(j1) * ones(m,n);
  end
z= z+w.*(y.^(i1-1));
end



function p = polyfit2d(x,y,z,n,m)
% P= POLYFIT2D(x,y,z,n,m) finds the coefficients of a
%  polynomial function of 2 variables formed from the
%  data in vectors x and y of degrees n and m, respectively,
%  that fit        the data in vector z in a least-squares sense.  
%
% The regression problem is formulated in matrix format as:
%
%        z = A*P   So that if the polynomial is cubic in  
%             x and linear in y, the problem becomes:
%
%        z = [y.*x.^3  y.*x.^2  y.*x  y x.^3  x.^2  x  ones(length(x),1)]*
%            [p31 p21 p11 p01 p30 p20 p10 p00]'                     
%               
%  Note that the various xy products are column vectors of length(x).
%
%  The coefficents of the output p   
%  matrix are arranged as shown:
%
%      p31 p30
%      p21 p20
%      p11 p10
%      p01 p00
%
% The indices on the elements of p correspond to the
% order of x and y associated with that element.
%
% For a solution to exist, the number of ordered
% triples [x,y,z] must equal or exceed (n+1)*(m+1).
% Note that m or n may be zero.
%
% To evaluate the resulting polynominal function,
% use POLYVAL2D.

% Perry W. Stout  June 29, 1995
% 4829 Rockland Way
% Fair Oaks, CA  95628
% (916) 966-0236
% Based on the Matlab function polyfit.

if any((size(x) ~= size(y)) | (size(z) ~= size(y)))
        error('X, Y,and Z vectors must be the same size')
end

x = x(; y = y(; z= z(;  % Switches vectors to columns--matrices, too

if length(x) < (n+1)*(m+1)
error('Number of points must equal or exceed order of polynomial function.')
end

n = n + 1;
m = m + 1; % Increments n and m to equal row, col numbers of p.

a = zeros(max(size(x)),n*m);

% Construct the extended Vandermonde matrix, containing all xy products.

for i1= 1:m
   for j1=1:n
             a(:,j1+(i1-1)*n) = (x.^(n-j1)).*(y.^(m-i1));
   end
end
p1 = (a\z);
% Reform p as a matrix.

p=[];
for i1=1:m
p=[p, p1((n*(i1-1)+1)n*i1))];
end

itmwk

% M-file "polyfit2dtest"
disp('A Test of ''polyfit2d'' and ''polyval2d.''')
disp(' ')
disp('POLYFIT2D and POLYVAL2D are used to form polynomial data fits to')
disp('functions of two variables, and compute the outputs of the resulting')
disp('polynomial forms.')
disp(' ')
disp('We will use the PEAKS function as the source of our data. Since this')
disp('function is composed of both polynomial and exponential terms, we are')
disp('assured this a ''reasonable'' test--That is, the polynomial can never')
disp('match the output exactly, and consequently there will always be some')
disp('mismatch between PEAKS and the fitting polynomial.')
disp(' ')
disp('Get enough points from PEAKS to define a smooth surface.')
disp('Concentrate on the area near the origin, before PEAKS dies out.')
disp(' ')
echo on
[xs,ys]=meshgrid(-2:.2:2,-2:.2:2);
z=peaks(xs,ys);
echo off
disp(' ')
disp('Plot this section of peaks:')
disp(' ')
figure(1)
clf
echo on
surf(xs,ys,z)
echo off
disp(' ')
disp('Now use POLYFIT2D to form a 5th order polynomial in both x and y,')
disp('using the ordered triples xs, ys, z. Note that POLYFIT2D automatically')
disp('''unwraps'' the matrices xs, ys, and z and forms appropriate vectors')
disp('for the least squares computation--you don''t have to do it yourself.')
disp(' ')
echo on
p=polyfit2d(xs,ys,z,5,5)
echo off
disp(' ')
disp('Now check the fit, using a grid twice as dense as the one used to sample PEAKS:')
disp(' ')
figure(2)
clf
echo on
[xp,yp]=meshgrid(-2:.1:2,-2:.1:2);
zp=polyval2d(p,xp,yp);
surf(xp,yp,zp)
echo off

chinamaker

你的第二个帖子的倒数12行 p=polyfit2d(xs,ys,z,5,5)
polyfit2d好像matlab不存在，请给说明一下。

chinamaker

谢谢，我知道了，对不起。刚才没有看清。

itmwk

没什么

TeresaHer

我是一个初学者，刚好这两天急用二维拟合，只是没太看懂上面的程序，还请大人仔细讲解一下相关函数和具体算法，请速回，谢谢了！！谢谢！

favorapple

太厉害了！[em06]

guanweihaha

十分感谢，但还想请教一下在哪儿可以查到这种matlab的原代码，有专门的网站？谢谢

深蓝世界

找了好久的曲面拟合,先顶一个