fncmb.m

演示matlab曲线拟和与插直的基本方法

function fn = fncmb(fn1,fnorsc,fn2,sc2)
%FNCMB Arithmetic with function(s).
%
%   FNCMB(function,operation) operation applied to function; specifically:
%     FNCMB(function,scalar)    multiplies function by scalar,
%     FNCMB(function,vector)    translates function value by vector,
%         (if function is scalar-valued, use
%         FNCMB(function,'+',scalar) to translate function value by scalar)
%     FNCMB(function,matrix)    applies matrix to coefficients of function.
%     FNCMB(function,string)    applies the m-file or function specified
%                               by string to coefficients of function.
%
%   FNCMB(function,function)  sum of the  two  functions of same form
%   FNCMB(function,matrix,function)  same as
%                        FNCMB(fncmb(function,matrix),function)
%   FNCMB(function,matrix,function,matrix)  same as
%                        FNCMB(fncmb(function,matrix),fncmb(function,matrix))
%
%   FNCMB(function,op,function)  ppform of the sum (op is '+'), difference
%                             (op is '-'), or pointwise product (op is '*') of
%                             the two functions possibly of different forms.
%                             In particular, in case of addition/subtraction,
%                             the second function may be just a point in the
%                             target of the first (i.e., a constant function).
%
%   At present, all functions must be UNIVARIATE except for the call
%   FNCMB(function,operation).
%
%   Examples:
%
%      fncmb( sp1, '+', sp2 );
%
%   returns the (pointwise) sum of the function in SP1 and that in SP2, while
%
%      fncmb( spmak( augknt(4:9,4), eye(8) ), [1:8] )
%
%   is a complicated way to construct sum_{j=1:8} j*B_j for the eight cubic
%   B-splines B_j for the knot sequence AUGKNT(4:9,4).
%
%   If SP contains a spline in B-form, then
%
%      spa = fncmb( sp, 'abs' );
%
%   changes all its coefficients to their absolute value.
%
%   If FN contains a 3-vector-valued function (i.e., a map into R^3), then
%
%      fncmb( fn, [1 0 0; 0 0 1] )
%
%   provides the projection onto the (x,z)-plane, while
%
%      fncmb( fncmb (fn, [1 0 0; 0 0 -1; 0 1 0] ), [1;2;3] )
%
%   rotates the image of the function in FN 90 degrees around the x-axis and
%   then translates that by the vector (1,2,3).

%   Copyright 1987-2005 C. de Boor and The MathWorks, Inc.
%   $Revision: 1.20.4.2 $

if nargin<2, fn = fn1; return, end

if ~isstruct(fn1)
   try
      fn1 = fn2fm(fn1);
   catch
      error('SPLINES:FNCMB:unknownfrstfn',...
      'First argument is of unknown function type.')
   end
end

try
   fn1form = fnbrk(fn1,'form');
catch
   error('SPLINES:FNCMB:unknownfrstfn',...
         'First argument is of unknown function type.')
end

if nargin==2&&~isstruct(fnorsc) % we are to apply the operation specified
                                % by FNORSC to the coefficients of FN1
   fn1form(3:end) = [];
   [coefs,d] = fnbrk(fn1,'coefs','dim');
   sizeval = d; if length(d)>1,  d = prod(d); end
   if ischar(fnorsc) %the operation is given by some m-file or built-in function
      eval(['coefs = ',fnorsc,'(coefs);'])
   elseif length(fnorsc)==1% we multiply the coefficients by the scalar FNORSC
      switch fn1form
      case {'rp','rB'} % in rational case, only multiply the numerator
         coefs(1:d,:) = fnorsc*coefs(1:d,:);
      otherwise
         coefs = fnorsc*coefs;
      end
   else %FNORSC is either an array with more dimensions than fn1.dim and 
        % with its last dimensions matching exactly fn1.dim;
        % or the first dimensions of FNORSC exactly match fn1.dim and the rest
        % are trivial; anything else is an error.
      sizesc = size(fnorsc);
        % will have to deal with the problem of vanishing unit dimensions
      diffl = length(sizesc)-length(sizeval);
      if diffl>0&&isequal(sizesc(diffl+1:end),sizeval)
            % we are to apply the `matrix' FNORSC to the coefficients
         sizeval = sizesc(1:diffl); r = prod(sizeval);
         fnorsc = reshape(fnorsc,r,d);
         sizec = size(coefs);
         l = 1;
         if fn1form(2)=='p'
            l = fnbrk(fn1,'pieces');
            if length(l)>1, l = 1; end
         end
         coefs = reshape(coefs,sizec(1)/l,l*prod(sizec(2:end)));
         switch fn1form
         case {'rB','rp'}
            coefs = reshape([fnorsc*coefs(1:d,:); coefs(d+1,:)], ...
                            [(r+1)*l,sizec(2:end)]);
         otherwise
            coefs = reshape(fnorsc*coefs,[r*l,sizec(2:end)]);
         end
         d = r;
      elseif prod(sizesc)==d&&...
             ((diffl>=0&&isequal(sizesc(1:length(sizeval)),sizeval))||...
             (diffl<0&&isequal(sizesc,sizeval(1:length(sizec)))))
         % we are to translate the appropriate coefficients by FNORSC
         fnorsc = reshape(fnorsc,d,1);
         switch fn1form
         case {'rB','rp'}
            l = 1;
            if fn1form(2)=='p'
               l = fnbrk(fn1,'pieces');
               if length(l)>1, l = 1; end
            end
            sizec = size(coefs); ncoefs = l*prod(sizec(2:end));
            coefs = reshape(coefs,sizec(1)/l,ncoefs);
            coefs(1:d,:) = repmat(fnorsc,1,ncoefs).* ...
                           repmat(coefs(d+1,:),d,1) + ...
                           coefs(1:d,:);
            coefs = reshape(coefs,sizec);
         case 'pp'         % now only the constant terms get translated.
                           % In this case, the matrix COEFS is, equivalently,
                           % an array, of size [d,l1*k1,...,lm*km], with
                           % m>1 the dimension of the domain of the function.
                           % This means that FNORSC is added to each column of
                           % COEFS(:,(1:l1)*k1,...,(1:lm)*km). Our problem is
                           % that  m  is a variable. We take the coward's way
                           % out: a loop.
                           % Another complication: for m==1, COEFS is of size
                           % [d*l,k].
            k = fnbrk(fn1,'order'); l = fnbrk(fn1,'pieces');
            m = length(k); len(m:-1:1) = cumprod(l(m:-1:1));
            index = (k(m)-1)*l(m)+(1:l(m)).';
            for j=m-1:-1:1
               index = reshape(repmat((l(j)*k(j))*(index-1),1,l(j)) + ...
                       repmat((k(j)-1)*l(j)+(1:l(j)),len(j+1),1),len(j),1);
            end
            prodk = prod(k);
            coefs = reshape(coefs,d,len(1)*prodk);
            coefs(:,index) = coefs(:,index) + repmat(fnorsc,1,len(1));
            if m>1
               coefs = reshape(coefs,[d,l.*k]);
            else
               coefs = reshape(coefs,d*l,k);
            end

         case {'B-','BB'}
            coefs = coefs + repmat(fnorsc,[1,fnbrk(fn1,'number')]);
         case 'st' % in this case, only the coefficient of the constant term
                   % gets translated, i.e., the last coefficient
                   % (this is no good for 'st-tp' ...):
            coefs(:,end) = coefs(:,end) + fnorsc;
         end
      else
         error('SPLINES:FNCMB:unknownsecfn',...
               'Second argument is of unknown function type.')
      end
   end
      % generate fn to be of the same form as fn1
   switch fn1form
   case 'pp'
      breaks = fnbrk(fn1,'breaks');
      if iscell(breaks)
         fn = ppmak(breaks,coefs, ...
               [size(coefs,1),fnbrk(fn1,'pieces').*fnbrk(fn1,'order')]);
      else
         fn = ppmak(breaks,coefs,d);
      end
   case {'B-','BB'}
      fn = spmak(fnbrk(fn1,'knots'),coefs,[size(coefs,1),fnbrk(fn1,'number')]);
      fn.form = fn1.form;
   case 'rB'
      fn = rsmak(fnbrk(fn1,'knots'),coefs,[size(coefs,1),fnbrk(fn1,'number')]);
   case 'rp'
      breaks = fnbrk(fn1,'breaks');
      if iscell(breaks)
         fn = rpmak(breaks,coefs, ...
               [size(coefs,1),fnbrk(fn1,'pieces').*fnbrk(fn1,'order')]);
      else
         fn = rpmak(breaks,coefs,d);
      end
   case 'st'
      fn = stmak(fnbrk(fn1,'centers'),coefs,fn1.form(4:end), ...
                 fnbrk(fn1,'interv'));
   otherwise
      % cannot happen since the earlier fnbrk(fn1,'coefs') would have failed
      error('SPLINES:FNCMB:impossible',...
            'The world is coming to an end; notify [e-mail]')
   end
   if length(sizeval)>1, fn = fnchg(fn,'dz',sizeval); end
   return
end

if fnbrk(fn1,'var')>1&&(nargin>2||(nargin==2&&isstruct(fnorsc)))
   error('SPLINES:FNCMB:onlyuni',...
         'At present, this works only for univariate functions.'), end

if ischar(fnorsc) % we are in the case fn1 'op' fn2, with 'op' +-* and the
                 % two functions of possibly different forms.
   if length(fnorsc)~=1||~any('+-*'==fnorsc)
      error('SPLINES:FNCMB:unknownop',...
           ['The operation ',fnorsc,' is not recognizable.'])
   end
   if nargin~=3
      error('SPLINES:FNCMB:needthree',...
            'Expect three input arguments for this case.')
   end

   if ~isstruct(fn2) 
      try fn2 = fn2fm(fn2);% check whether FN2 is a function written the old way
      catch                % otherwise, FN2 must be a point in FN1's target
       % must ascertain that the constant FN2 is in the target of FN1
       % or else a scalar (in which case it will be taken to be the 
       % corresponding constant element in the target of FN1).
       sizeval = fnbrk(fn1,'dim'); d = prod(sizeval); fn2 = fn2(:);
       if length(fn2)>1 && length(fn2)~=prod(sizeval)
          error('SPLINES:FNCMB:constwrongsize',...
                'The constant, fn2, must lie in the target of fn1.')
       end
       switch fnorsc
       case '*'
          if length(fn2)==1, fn = fncmb(fn1,fn2); clear sizeval
   	  else % as a quickie, handle pointwise multiplication as matrix mult
  	     fn = fncmb(fn1,reshape(diag(fn2), [sizeval,sizeval]));
  	  end
  
       otherwise  % it's addition of a constant and, since we can't be sure
                  % that the relevant B-splines form a partition of unity,