fncmb.m

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

                  % it's easiest just to convert to ppform and then operate
                  % on the constant coefficients.
          if strcmp(fn1form(1:2),'st')
              error('SPLINES:FNCMB:notforst',...
                    'This does not work for st functions (yet).')
          end
          if d>1&&length(fn2)==1, fn2 = repmat(fn2,d,1); end
          if fn1form(1)=='r'
             if fn1form(2)=='B', fn1 = fn2fm(fn1,'rp'); end
             [breaks,c,l,k,d] = ppbrk(fn2fm(fn1,'pp')); d = d-1;
             % s/w + c = (c*w+s)/w
             temp = reshape(c,[d+1,l,k]);
             eval(['temp(1:d,:) = temp(1:d,:)',fnorsc,'fn2*temp(d+1,:);'])
             fn = rpmak(breaks,reshape(temp, (d+1)*l,k),d);
          else
             if fn1form(1)=='B', fn1 = sp2pp(fn1); end
             [breaks,c,l,k] = ppbrk(fn1);
             %tmp = reshape(1:d,d,1); tmp = reshape(tmp(:,ones(1,l)),d*l,1);
             tmp = reshape(repmat(reshape(1:d,d,1),1,l),d*l,1);
             eval(['c(:,k) = c(:,k)',fnorsc,'fn2(tmp);'])
             fn = ppmak(breaks,c,d);
          end
       end
       if exist('sizeval','var')&&length(sizeval)>1
  	   fn = fnchg(fn,'dz',sizeval);
       end
       return
      end
   end

   %  else,  convert both FN1 and FN2 to ppform, if need be.

      % know from earlier test that fn2 is struct; now check that it is
      % actually a form we know: 
   try
      fn2form = fnbrk(fn2,'form');
   catch
      error('SPLINES:FNCMB:unknownsecfun',...
            'Second argument is of unknown function type.')
   end

      % at present, this does not work for the stform
   if isequal(fn1form(1:2),'st')||isequal(fn2form(1:2),'st')
      error('SPLINES:FNCMB:notforst',...
            'This does not work for functions in stform (yet).')
   end

      % next, establish the basic interval of the result as the smallest
      % interval containing both basic intervals and extend both fns to it.
   intervs = [fnbrk(fn1,'interv'); fnbrk(fn2,'interv')];
   interv = [min(intervs(:,1)), max(intervs(:,2))];
   fn1 = fnbrk(fn1,interv); fn2 = fnbrk(fn2,interv);

      %  next, convert rationals to their equivalent pp-version
   if fn1form(1)=='r', fn1 = fn2fm(fn1,'pp'); end
   if fn2form(1)=='r', fn2 = fn2fm(fn2,'pp'); end

   if fn1form(1)=='B', fn1 = sp2pp(fn1); end
   if fn2form(1)=='B', fn2 = sp2pp(fn2); end

      % establish the common refinement of the two break sequences
   tmp = sort([ppbrk(fn1,'b') ppbrk(fn2,'b')]);
   breaks = tmp([find(diff(tmp)>0),end]);
     % refine breaks in both FN1 and FN2
   fn1 = pprfn(fn1,breaks); fn2 = pprfn(fn2,breaks);
   [ignored,c1,l1,k1,d1] = ppbrk(fn1); [ignored,c2,l2,k2,d2] = ppbrk(fn2);
   if (fnorsc=='+'||fnorsc=='-')
      if fn1form(1)=='r'||fn2form(1)=='r'
         % c1 = reshape(c1,[d1,l1,k1]); c2 = reshape(c2,[d2,l2,k2]);
         k = k1+k2-1;
         if fn1form(1)~='r' % only the second is rational
            if d1~=d2-1
               error('SPLINES:FNCMB:targetsdontmatch',...
                     'The two functions must have the same target.'), end
            % s1 op s2/w2 = (s1*w2 op s2)/w2
            c = reshape([zeros(d2*l2,k1-1),c2],[d2,l2,k]);
            temp =matconv(c1,reshape(repmat(c(d2,:,k1:end),[d1,1,1]),d1*l1,k2));
            eval(['c(1:d1,:,:) = reshape(temp,[d1,l1,k])',fnorsc, ...
                                       'c(1:d1,:,:);'])
            d1 = d1+1;
         elseif fn2form(1)~='r' % only the first is rational
            if d1-1~=d2
               error('SPLINES:FNCMB:targetsdontmatch',...
                     'The two functions must have the same target.'), end
            % s1/w1 op s2 = (s1 op s2*w1)/w1
            c = reshape([zeros(d1*l1,k2-1),c1],[d1,l1,k]);
            temp =matconv(c2,reshape(repmat(c(d1,:,k2:end),[d2,1,1]),d2*l2,k1));
            eval(['c(1:d2,:,:) = c(1:d2,:,:)',fnorsc, ...
                                       'reshape(temp,[d2,l2,k]);'])
         else                   % both are rational
            c1 = reshape(c1,[d1,l1,k1]); c2 = reshape(c2,[d2,l2,k2]);
            if d1~=d2
               error('SPLINES:FNCMB:targetsdontmatch',...
                     'The two functions must have the same target.'), end
            % s1/w1 op s2/w2 = (s1*w2 op s2*w1)/(w1*w2)
            c = zeros([d1,l1,k]); d = d1-1;
            eval(['c(1:d,:,:)=reshape(matconv(reshape(c1(1:d,:,:),d*l1,k1)',...
                   ',reshape(repmat(c2(d1,:,:),[d,1,1]),d*l2,k2))', fnorsc, ...
                                'matconv(reshape(c2(1:d,:,:),d*l2,k2)',...
                   ',reshape(repmat(c1(d1,:,:),[d,1,1]),d*l1,k1)),[d,l1,k]);'])
            c(d1,:,:) = reshape(matconv(reshape(c1(d1,:,:),l1,k1), ...
                                        reshape(c2(d1,:,:),l2,k2)),[1,l1,k]);
         end
      else
         if d1~=d2
            error('SPLINES:FNCMB:targetsdontmatch',...
                  'The two functions must have the same target.'), end
         k = k1;
         if k1<k2, k = k2; c1 = [zeros(d1*l1,k2-k1) c1];
         elseif k1>k2, c2 = [zeros(d2*l2,k1-k2) c2]; end
         eval(['c = c1',fnorsc,'c2;'])
      end
      c = reshape(c,[d1*l1,k]);
   else  % we pointwise multiply:
      k = k1+k2-1;
      if (fn1form(1)=='r'&&fn2form(1)=='r')||(fn1form(1)~='r'&&fn2form(1)~='r')
         c = matconv(c1,c2);
      else
         if fn1form(1)=='r'
            if d1-1~=d2
               error('SPLINES:FNCMB:targetsdontmatch',...
                     'The two functions must have the same target.'), end
            c = reshape([zeros(d1*l1,k2-1),c1],[d1,l1,k]);
            c(1:d2,:,:) = reshape(...
                         matconv(reshape(c(1:d2,:,k2:end),d2*l1,k1),c2), ...
                         [d2,l1,k]);
            c = reshape(c,[d1*l1,k]);
         else
            if d1~=d2-1
               error('SPLINES:FNCMB:targetsdontmatch',...
                     'The two functions must have the same target.'), end
            c = reshape([zeros(d2*l2,k1-1),c2],[d2,l2,k]);
            c(1:d1,:,:) = reshape(...
                         matconv(reshape(c(1:d1,:,k1:end),d1*l2,k2),c1), ...
                         [d1,l2,k]);
            c = reshape(c,[d2*l2,k]);
            d1 = d1+1;
         end
      end
   end
   if fn1form(1)=='r'||fn2form(1)=='r'
      fn = rpmak(breaks,c,d1-1);
   else
      fn = ppmak(breaks,c,d1);
   end
   return
end

%%%%%%%  at this point, we know that  FNORSC  is a matrix (iff nargin>2)
      %  or a function (iff nargin==2)

% reduce it to the case of adding fn1 and fn2
if nargin==2
   fn2 = fnorsc;
else
   fn1 = fncmb(fn1,fnorsc);
   if nargin>3
      fn2 = fncmb(fn2,sc2);
   end
end

try
   coefs2 = fnbrk(fn2,'c');
catch
   error('SPLINES:FNCMB:unknownsecfn',...
         'Second argument is of unknown function type.')
end

if ~isequal(fn1form,fnbrk(fn2,'form'))
  error('SPLINES:FNCMB:formsdontmatch',...
       ['The two functions being added are of different forms.\n',...
	'Had you meant to use fncmb(fn1,''+'',fn2) ?'],0)
end

switch fn1form(1)
case 'B'
   [knots,coefs] = spbrk(fn1);
   fn = spmak(knots,sumcoefs(coefs,coefs2));
case 'p'
   [breaks,coefs,l,k,d] = ppbrk(fn1);
   coefs = sumcoefs(coefs,coefs2);
   if length(k)>1
      fn = ppmak(breaks,coefs);
   else
      fn = ppmak(breaks,coefs,[d,l,k]);
   end

case 'r'

   switch fn1form(2)
   case 'B'
      knots = fnbrk(fn1,'knots'); c1 = fnbrk(fn1,'coefs');
      if any(size(c1)-size(coefs2))
         fn = fncmb(fn1,'+',fn2);
      else
         dw = c1(end,:)-coefs2(end,:);
         if max(abs(dw))>1e-12*max(abs(c1(end,:)))
            fn = fncmb(fn1,fn2);
         else
            c1(1:end-1,:) = c1(1:end-1,:)+coefs2(1:end-1,:);
            fn = rsmak(knots, c1);
         end
      end
   case 'p'
      fn = fncmb(fn1,'+',fn2);
   otherwise
      error('SPLINES:FNCMB:unknownfrstfn',...
            'First argument is of unknown function type.')
   end
otherwise
   error('SPLINES:FNCMB:unknownfrstfn',...
         'First argument is of unknown function type.')
end

function asb = matconv(a,b)
%MATCONV row-by-row convolution of the two matrices a and b
%
% asb(j,:) = conv(a(j,:),b(j,:)), all j
%
% An error will occur if a and b fail to have the same number of rows.

[ra,ca] = size(a); cb = size(b,2);

   % make sure that a is the one with fewer columns
if ca>cb
   temp = b; b = a; a = temp; temp = ca; ca = cb; cb = temp;
end

asb = zeros(ra,ca+cb-1);
for j=1:ca
   asb(:,j-1+(1:cb)) = asb(:,j-1+(1:cb))+repmat(a(:,j),1,cb).*b;
end

function coefs = sumcoefs(coef1,coef2)
%SUMCOEFS if sizes of coef1,2 match, return sum; else, print error message

if ~isequal(size(coef1),size(coef2))
  error('SPLINES:FNCMB:incompatiblecoefs', ...
        ['The coefficient arrays of the two functions\nbeing added have',...
	' incompatible sizes.\nHad you meant to use fncmb(fn1,''+'',fn2) ?'])
end
coefs = coef1 + coef2;