Description of smooth

0001 function [c,ww] = smooth(varargin)
0002 %SMOOTH  Smooth data.
0003 %   Z = SMOOTH(Y) smooths data Y using a 5-point moving average.
0004 %
0005 %   Z = SMOOTH(Y,SPAN) smooths data Y using SPAN as the number of points used
0006 %   to compute each element of Z.
0007 %
0008 %   Z = SMOOTH(Y,SPAN,METHOD) smooths data Y with specified METHOD. The
0009 %   available methods are:
0010 %
0011 %           'moving'   - Moving average (default)
0012 %           'lowess'   - Lowess (linear fit)
0013 %           'loess'    - Loess (quadratic fit)
0014 %           'sgolay'   - Savitzky-Golay
0015 %           'rlowess'  - Robust Lowess (linear fit)
0016 %           'rloess'   - Robust Loess (quadratic fit)
0017 %
0018 %   Z = SMOOTH(Y,METHOD) uses the default SPAN 5.
0019 %
0020 %   Z = SMOOTH(Y,SPAN,'sgolay',DEGREE) and Z = SMOOTH(Y,'sgolay',DEGREE)
0021 %   additionally specify the degree of the polynomial to be used in the
0022 %   Savitzky-Golay method. The default DEGREE is 2. DEGREE must be smaller
0023 %   than SPAN.
0024 %
0025 %   Z = SMOOTH(X,Y,...) additionally specifies the X coordinates.  If X is
0026 %   not provided, methods that require X coordinates assume X = 1:N, where
0027 %   N is the length of Y.
0028 %
0029 %   Notes:
0030 %   1. When X is given and X is not uniformly distributed, the default method
0031 %   is 'lowess'.  The 'moving' method is not recommended.
0032 %
0033 %   2. For the 'moving' and 'sgolay' methods, SPAN must be odd.
0034 %   If an even SPAN is specified, it is reduced by 1.
0035 %
0036 %   3. If SPAN is greater than the length of Y, it is reduced to the
0037 %   length of Y.
0038 %
0039 %   4. In the case of (robust) lowess and (robust) loess, it is also
0040 %   possible to specify the SPAN as a percentage of the total number
0041 %   of data points. When SPAN is less than or equal to 1, it is
0042 %   treated as a percentage.
0043 %
0044 %   For example:
0045 %
0046 %   Z = SMOOTH(Y) uses the moving average method with span 5 and
0047 %   X=1:length(Y).
0048 %
0049 %   Z = SMOOTH(Y,7) uses the moving average method with span 7 and
0050 %   X=1:length(Y).
0051 %
0052 %   Z = SMOOTH(Y,'sgolay') uses the Savitzky-Golay method with DEGREE=2,
0053 %   SPAN = 5, X = 1:length(Y).
0054 %
0055 %   Z = SMOOTH(X,Y,'lowess') uses the lowess method with SPAN=5.
0056 %
0057 %   Z = SMOOTH(X,Y,SPAN,'rloess') uses the robust loess method.
0058 %
0059 %   Z = SMOOTH(X,Y) where X is unevenly distributed uses the
0060 %   'lowess' method with span 5.
0061 %
0062 %   Z = SMOOTH(X,Y,8,'sgolay') uses the Savitzky-Golay method with
0063 %   span 7 (8 is reduced by 1 to make it odd).
0064 %
0065 %   Z = SMOOTH(X,Y,0.3,'loess') uses the loess method where span is
0066 %   30% of the data, i.e. span = ceil(0.3*length(Y)).
0067 %
0068 %   See also SPLINE.
0069 
0070 %   Copyright 2001-2009 The MathWorks, Inc.
0071 %   $Revision: 1.1 $  $Date: 2010/09/15 10:04:55 $
0072 
0073 if nargin < 1
0074     error('curvefit:smooth:needMoreArgs', ...
0075         'SMOOTH needs at least one argument.');
0076 end
0077 
0078 if nargout > 1 % Called from the GUI cftool
0079     ws = warning('off', 'all'); % turn warning off and record the previous warning state.
0080     [lw,lwid] = lastwarn;
0081     lastwarn('');
0082 else
0083     ws = warning('query','all'); % Leave warning state alone but save it so resets are no-ops.
0084 end
0085 
0086 % is x given as the first argument?
0087 if nargin==1 || ( nargin > 1 && (length(varargin{2})==1 || ischar(varargin{2})) )
0088     % smooth(Y) | smooth(Y,span,...) | smooth(Y,method,...)
0089     is_x = 0; % x is not given
0090     y = varargin{1};
0091     y = y(:);
0092     x = (1:length(y))';
0093 else % smooth(X,Y,...)
0094     is_x = 1;
0095     y = varargin{2};
0096     x = varargin{1};
0097     y = y(:);
0098     x = x(:);
0099 end
0100 
0101 % is span given?
0102 span = [];
0103 if nargin == 1+is_x || ischar(varargin{2+is_x})
0104     % smooth(Y), smooth(X,Y) || smooth(X,Y,method,..), smooth(Y,method)
0105     is_span = 0;
0106 else
0107     % smooth(...,SPAN,...)
0108     is_span = 1;
0109     span = varargin{2+is_x};
0110 end
0111 
0112 % is method given?
0113 method = [];
0114 if nargin >= 2+is_x+is_span
0115     % smooth(...,Y,method,...) | smooth(...,Y,span,method,...)
0116     method = varargin{2+is_x+is_span};
0117 end
0118 
0119 t = length(y);
0120 if t == 0
0121     c = y;
0122     ww = '';
0123     if nargout > 1
0124         ww = lastwarn;
0125         lastwarn(lw,lwid);
0126         warning(ws);  % turn warning back to the previous state.
0127     end
0128     return
0129 elseif length(x) ~= t
0130     warning(ws); % reset warn state before erroring
0131     error('curvefit:smooth:XYmustBeSameLength',...
0132         'X and Y must be the same length.');
0133 end
0134 
0135 if isempty(method)
0136     diffx = diff(x);
0137     if uniformx(diffx,x,y)
0138         method = 'moving'; % uniformly distributed X.
0139     else
0140         method = 'lowess';
0141     end
0142 end
0143 
0144 % realize span
0145 if span <= 0
0146     warning(ws); % reset warn state before erroring
0147     error('curvefit:smooth:spanMustBePositive', ...
0148         'SPAN must be positive.');
0149 end
0150 if span < 1, span = ceil(span*t); end % percent convention
0151 if isempty(span), span = 5; end % smooth(Y,[],method)
0152 
0153 idx = 1:t;
0154 
0155 sortx = any(diff(isnan(x))<0);   % if NaNs not all at end
0156 if sortx || any(diff(x)<0) % sort x
0157     [x,idx] = sort(x);
0158     y = y(idx);
0159 end
0160 
0161 c = NaN(size(y));
0162 ok = ~isnan(x);
0163 switch method
0164     case 'moving'
0165         c(ok) = moving(x(ok),y(ok),span);
0166     case {'lowess','loess','rlowess','rloess'}
0167         robust = 0;
0168         iter = 5;
0169         if method(1)=='r'
0170             robust = 1;
0171             method = method(2:end);
0172         end
0173         c(ok) = lowess(x(ok),y(ok),span, method,robust,iter);
0174     case 'sgolay'
0175         if nargin >= 3+is_x+is_span
0176             degree = varargin{3+is_x+is_span};
0177         else
0178             degree = 2;
0179         end
0180         if degree < 0 || degree ~= floor(degree) || degree >= span
0181             warning(ws); % reset warn state before erroring
0182             error('curvefit:smooth:invalidDegree', ...
0183                 'Degree must be an integer between 0 and span-1.');
0184         end
0185         c(ok) = sgolay(x(ok),y(ok),span,degree);
0186     otherwise
0187         warning(ws); % reset warn state before erroring
0188         error('curvefit:smooth:unrecognizedMethod', ...
0189             'SMOOTH: Unrecognized method.');
0190 end
0191 
0192 c(idx) = c;
0193 
0194 if nargout > 1
0195     ww = lastwarn;
0196     lastwarn(lw,lwid);
0197     warning(ws);  % turn warning back to the previous state.
0198 end
0199 
0200 %--------------------------------------------------------------------
0201 function c = moving(x,y, span)
0202 % moving average of the data.
0203 
0204 ynan = isnan(y);
0205 span = floor(span);
0206 n = length(y);
0207 span = min(span,n);
0208 width = span-1+mod(span,2); % force it to be odd
0209 xreps = any(diff(x)==0);
0210 if width==1 && ~xreps && ~any(ynan), c = y; return; end
0211 if ~xreps && ~any(ynan)
0212     % simplest method for most common case
0213     c = filter(ones(width,1)/width,1,y);
0214     cbegin = cumsum(y(1:width-2));
0215     cbegin = cbegin(1:2:end)./(1:2:(width-2))';
0216     cend = cumsum(y(n:-1:n-width+3));
0217     cend = cend(end:-2:1)./(width-2:-2:1)';
0218     c = [cbegin;c(width:end);cend];
0219 elseif ~xreps
0220     % with no x repeats, can take ratio of two smoothed sequences
0221     yy = y;
0222     yy(ynan) = 0;
0223     nn = double(~ynan);
0224     ynum = moving(x,yy,span);
0225     yden = moving(x,nn,span);
0226     c = ynum ./ yden;
0227 else
0228     % with some x repeats, loop
0229     notnan = ~ynan;
0230     yy = y;
0231     yy(ynan) = 0;
0232     c = zeros(n,1);
0233     for i=1:n
0234         if i>1 && x(i)==x(i-1)
0235             c(i) = c(i-1);
0236             continue;
0237         end
0238         R = i;                                 % find rightmost value with same x
0239         while(R<n && x(R+1)==x(R))
0240             R = R+1;
0241         end
0242         hf = ceil(max(0,(span - (R-i+1))/2));  % need this many more on each side
0243         hf = min(min(hf,(i-1)), (n-R));
0244         L = i-hf;                              % find leftmost point needed
0245         while(L>1 && x(L)==x(L-1))
0246             L = L-1;
0247         end
0248         R = R+hf;                              % find rightmost point needed
0249         while(R<n && x(R)==x(R+1))
0250             R = R+1;
0251         end
0252         c(i) = sum(yy(L:R)) / sum(notnan(L:R));
0253     end
0254 end
0255 
0256 %--------------------------------------------------------------------
0257 function c = lowess(x,y, span, method, robust, iter)
0258 % LOWESS  Smooth data using Lowess or Loess method.
0259 %
0260 % The difference between LOWESS and LOESS is that LOWESS uses a
0261 % linear model to do the local fitting whereas LOESS uses a
0262 % quadratic model to do the local fitting. Some other software
0263 % may not have LOWESS, instead, they use LOESS with order 1 or 2 to
0264 % represent these two smoothing methods.
0265 %
0266 % Reference:
0267 % [C79] W.S.Cleveland, "Robust Locally Weighted Regression and Smoothing
0268 %    Scatterplots", _J. of the American Statistical Ass._, Vol 74, No. 368
0269 %    (Dec.,1979), pp. 829-836.
0270 %    http://www.math.tau.ac.il/~yekutiel/MA%20seminar/Cleveland%201979.pdf
0271 
0272 
0273 n = length(y);
0274 span = floor(span);
0275 span = min(span,n);
0276 c = y;
0277 if span == 1
0278     return;
0279 end
0280 
0281 useLoess = false;
0282 if isequal(method,'loess')
0283     useLoess = true;
0284 end
0285 
0286 diffx = diff(x);
0287 
0288 % For problems where x is uniform, there's a faster way
0289 isuniform = uniformx(diffx,x,y);
0290 if isuniform
0291     % For uniform data, an even span actually covers an odd number of
0292     % points.  For example, the four closest points to 5 in the
0293     % sequence 1:10 are {3,4,5,6}, but 7 is as close as 3.
0294     % Therefore force an odd span.
0295     span = 2*floor(span/2) + 1;
0296 
0297     c = unifloess(y,span,useLoess);
0298     if ~robust || span<=2
0299         return;
0300     end
0301 end
0302 
0303 % Turn off warnings when called from command line (already off if called from
0304 % cftool).
0305 ws = warning( 'off', 'MATLAB:rankDeficientMatrix' );
0306 cleanup = onCleanup( @() warning( ws ) );
0307 
0308 
0309 ynan = isnan(y);
0310 anyNans = any(ynan(:));
0311 seps = sqrt(eps);
0312 theDiffs = [1; diffx; 1];
0313 
0314 if isuniform
0315     % We've already computed the non-robust smooth, so in preparation for
0316     % the robust smooth, compute the following arrays directly
0317     halfw = floor(span/2);
0318     
0319     % Each local interval is from |halfw| below the current index to |halfw|
0320     % above
0321     lbound = (1:n)-halfw;
0322     rbound = (1:n)+halfw;
0323     % However, there always has to be at least |span| points to the right of the
0324     % left bound
0325     lbound = min( n+1-span, lbound );
0326     % ... and at least |span| points to the left of the right bound
0327     rbound = max( span, rbound );
0328     % Furthermore, because these bounds index into vectors of length n, they
0329     % must contain valid indices
0330     lbound = max( 1, lbound );
0331     rbound = min( n, rbound );
0332     
0333     % Since the input is uniform we can use natural numbers for the input when
0334     % we need them.
0335     x = (1:numel(x))';
0336 else
0337     if robust
0338         % pre-allocate space for lower and upper indices for each fit,
0339         % to avoid re-computing this information in robust iterations
0340         lbound = zeros(n,1);
0341         rbound = zeros(n,1);
0342     end
0343 
0344     % Compute the non-robust smooth for non-uniform x
0345     for i=1:n
0346         % if x(i) and x(i-1) are equal we just use the old value.
0347         if theDiffs(i) == 0
0348             c(i) = c(i-1);
0349             if robust
0350                 lbound(i) = lbound(i-1);
0351                 rbound(i) = rbound(i-1);
0352             end
0353             continue;
0354         end
0355         
0356         % Find nearest neighbours
0357         idx = iKNearestNeighbours( span, i, x, ~ynan );
0358         if robust
0359             % Need to store neighborhoods for robust loop
0360             lbound(i) = min(idx);
0361             rbound(i) = max(idx);
0362         end
0363         
0364         if isempty(idx)
0365             c(i) = NaN;
0366             continue
0367         end
0368 
0369         x1 = x(idx)-x(i); % center around current point to improve conditioning
0370         d1 = abs(x1);
0371         y1 = y(idx);
0372 
0373         weight = iTricubeWeights( d1 );
0374         if all(weight<seps)
0375             weight(:) = 1;    % if all weights are 0, just skip weighting
0376         end
0377 
0378         v = [ones(size(x1)) x1];
0379         if useLoess
0380             v = [v x1.*x1]; %#ok<AGROW> There is no significant growth here
0381         end
0382         
0383         v = weight(:,ones(1,size(v,2))).*v;
0384         y1 = weight.*y1;
0385         if size(v,1)==size(v,2)
0386             % Square v may give infs in the \ solution, so force least squares
0387             b = [v;zeros(1,size(v,2))]\[y1;0];
0388         else
0389             b = v\y1;
0390         end
0391         c(i) = b(1);
0392     end
0393 end
0394 
0395 % now that we have a non-robust fit, we can compute the residual and do
0396 % the robust fit if required
0397 maxabsyXeps = max(abs(y))*eps;
0398 if robust
0399     for k = 1:iter
0400         r = y-c;
0401         
0402         % Compute robust weights
0403         rweight = iBisquareWeights( r, maxabsyXeps ); 
0404         
0405         % Find new value for each point.
0406         for i=1:n
0407             if i>1 && x(i)==x(i-1)
0408                 c(i) = c(i-1);
0409                 continue;
0410             end
0411             if isnan(c(i)), 
0412                 continue; 
0413             end
0414             
0415             idx = lbound(i):rbound(i);
0416             if anyNans
0417                 idx = idx(~ynan(idx));
0418             end
0419             % check robust weights for removed points
0420             if any( rweight(idx) <= 0 )
0421                 idx = iKNearestNeighbours( span, i, x, (rweight > 0) );
0422             end
0423             
0424             x1 = x(idx) - x(i);
0425             d1 = abs(x1);
0426             y1 = y(idx);
0427 
0428             weight = iTricubeWeights( d1 );
0429             if all(weight<seps)
0430                 weight(:) = 1;    % if all weights 0, just skip weighting
0431             end
0432 
0433             v = [ones(size(x1)) x1];
0434             if useLoess
0435                 v = [v x1.*x1]; %#ok<AGROW> There is no significant growth here
0436             end
0437             
0438             % Modify the weights based on x values by multiplying them by
0439             % robust weights.
0440             weight = weight.*rweight(idx);
0441             
0442             v = weight(:,ones(1,size(v,2))).*v;
0443             y1 = weight.*y1;
0444             if size(v,1)==size(v,2)
0445                 % Square v may give infs in the \ solution, so force least squares
0446                 b = [v;zeros(1,size(v,2))]\[y1;0];
0447             else
0448                 b = v\y1;
0449             end
0450             c(i) = b(1);
0451         end
0452     end
0453 end
0454 
0455 
0456 %--------------------------------------------------------------------
0457 function c=sgolay(x,y,f,k)
0458 % savitziki-golay smooth
0459 % (x,y) are given data. f is the frame length to be taken, should
0460 % be an odd number. k is the degree of polynomial filter. It should
0461 % be less than f.
0462 
0463 % Reference: Orfanidis, S.J., Introduction to Signal Processing,
0464 % Prentice-Hall, Englewood Cliffs, NJ, 1996.
0465 
0466 n = length(x);
0467 f = floor(f);
0468 f = min(f,n);
0469 f = f-mod(f-1,2); % will subtract 1 if frame is even.
0470 diffx = diff(x);
0471 notnan = ~isnan(y);
0472 nomissing = all(notnan);
0473 if f <= k && all(diffx>0) && nomissing, c = y; return; end
0474 hf = (f-1)/2; % half frame length
0475 
0476 idx = 1:n;
0477 if any(diffx<0) % make sure x is monotonically increasing
0478     [x,idx]=sort(x);
0479     y = y(idx);
0480     notnan = notnan(idx);
0481     diffx = diff(x);
0482 end
0483 % note that x is sorted so max(abs(x)) must be abs(x(1)) or abs(x(end));
0484 % already calculated diffx for monotonic case, so use it again. Only
0485 % recalculate if we sort x.
0486 if nomissing && uniformx(diffx,x,y)
0487     v = ones(f,k+1);
0488     t=(-hf:hf)';
0489     for i=1:k
0490         v(:,i+1)=t.^i;
0491     end
0492     [q,ignore]=qr(v,0);
0493     ymid = filter(q*q(hf+1,:)',1,y);
0494     ybegin = q(1:hf,:)*q'*y(1:f);
0495     yend = q((hf+2):end,:)*q'*y(n-f+1:n);
0496     c = [ybegin;ymid(f:end);yend];
0497     return;
0498 end
0499 
0500 % non-uniformly distributed data
0501 c = y;
0502 
0503 % Turn off warnings when called from command line (already off if called from
0504 % cftool).
0505 ws = warning('off', 'all');
0506 [lastwarnmsg,lastwarnid]=lastwarn;
0507 
0508 for i = 1:n
0509     if i>1 && x(i)==x(i-1)
0510         c(i) = c(i-1);
0511         continue
0512     end
0513     L = i; R = i;                          % find leftmost and rightmost values
0514     while(R<n && x(R+1)==x(i))
0515         R = R+1;
0516     end
0517     while(L>1 && x(L-1)==x(i))
0518         L = L-1;
0519     end
0520     HF = ceil(max(0,(f - (R-L+1))/2));     % need this many more on each side
0521 
0522     L = min(n-f+1,max(1,L-HF));            % find leftmost point needed
0523     while(L>1 && x(L)==x(L-1))
0524         L = L-1;
0525     end
0526     R = min(n,max(R+HF,L+f-1));            % find rightmost point needed
0527     while(R<n && x(R)==x(R+1))
0528         R = R+1;
0529     end
0530 
0531     tidx = L:R;
0532     tidx = tidx(notnan(tidx));
0533     if isempty(tidx)
0534         c(i) = NaN;
0535         continue;
0536     end
0537     q = x(tidx) - x(i);   % center to improve conditioning
0538     vrank = 1 + sum(diff(q)>0);
0539     ncols = min(k+1, vrank);
0540     v = ones(length(q),ncols);
0541     for j = 1:ncols-1
0542         v(:,j+1) = q.^j;
0543     end
0544     if size(v,1)==size(v,2)
0545         % Square v may give infs in the \ solution, so force least squares
0546         d = [v;zeros(1,size(v,2))]\[y(tidx);0];
0547     else
0548         d = v\y(tidx);
0549     end
0550     c(i) = d(1);
0551 end
0552 c(idx) = c;
0553 
0554 lastwarn(lastwarnmsg,lastwarnid);
0555 warning(ws);
0556 
0557 % --------------------------------------------
0558 function ys = unifloess(y,span,useLoess)
0559 %UNIFLOESS Apply loess on uniformly spaced X values
0560 
0561 y = y(:);
0562 
0563 % Omit points at the extremes, which have zero weight
0564 halfw = (span-1)/2;              % halfwidth of entire span
0565 d = abs((1-halfw:halfw-1));      % distances to pts with nonzero weight
0566 dmax = halfw;                    % max distance for tri-cubic weight
0567 
0568 % Set up weighted Vandermonde matrix using equally spaced X values
0569 x1 = (2:span-1)-(halfw+1);
0570 weight = (1 - (d/dmax).^3).^1.5; % tri-cubic weight
0571 v = [ones(length(x1),1) x1(:)];
0572 if useLoess
0573     v = [v x1(:).^2];
0574 end
0575 V = v .* repmat(weight',1,size(v,2));
0576 
0577 % Do QR decomposition
0578 [Q,ignore] = qr(V,0);
0579 
0580 % The projection matrix is Q*Q'.  We want to project onto the middle
0581 % point, so we can take just one row of the first factor.
0582 alpha = Q(halfw,:)*Q';
0583 
0584 % This alpha defines the linear combination of the weighted y values that
0585 % yields the desired smooth values.  Incorporate the weights into the
0586 % coefficients of the linear combination, then apply filter.
0587 alpha = alpha .* weight;
0588 ys = filter(alpha,1,y);
0589 
0590 % We need to slide the values into the center of the array.
0591 ys(halfw+1:end-halfw) = ys(span-1:end-1);
0592 
0593 % Now we have taken care of everything except the end effects.  Loop over
0594 % the points where we don't have a complete span.  Now the Vandermonde
0595 % matrix has span-1 points, because only 1 has zero weight.
0596 x1 = 1:span-1;
0597 v = [ones(length(x1),1) x1(:)];
0598 if useLoess
0599     v = [v x1(:).^2];
0600 end
0601 for j=1:halfw
0602     % Compute weights based on deviations from the jth point,
0603     % then compute weights and apply them as above.
0604     d = abs((1:span-1) - j);
0605     weight = (1 - (d/(span-j)).^3).^1.5;
0606     V = v .* repmat(weight(:),1,size(v,2));
0607     [Q,ignore] = qr(V,0);
0608     alpha = Q(j,:)*Q';
0609     alpha = alpha .* weight;
0610     ys(j) = alpha * y(1:span-1);
0611 
0612     % These coefficients can be applied to the other end as well
0613     ys(end+1-j) = alpha * y(end:-1:end-span+2);
0614 end
0615 
0616 % -----------------------------------------
0617 function isuniform = uniformx(diffx,x,y)
0618 %ISUNIFORM True if x is of the form a:b:c
0619 
0620 if any(isnan(y)) || any(isnan(x))
0621     isuniform = false;
0622 else
0623     isuniform = all(abs(diff(diffx)) <= eps*max(abs([x(1),x(end)])));
0624 end
0625 
0626 
0627 %------------------------
0628 function idx = iKNearestNeighbours( k, i, x, in )
0629 % Find the k points from x(in) closest to x(i)
0630 
0631 if nnz( in ) <= k
0632     % If we have k points or fewer, then return them all
0633     idx = find( in );
0634 else
0635     % Find the distance to the k closest point
0636     d = abs( x - x(i) );
0637     ds = sort( d(in) );
0638     dk = ds(k);
0639     
0640     % Find all points that are as close as or closer than the k closest point
0641     close = (d <= dk);
0642     
0643     % The required indices are those points that are both close and "in"
0644     idx = find( close & in );
0645 end
0646 
0647 % -----------------------------------------
0648 % Bi-square (robust) weight function
0649 function delta = iBisquareWeights( r, myeps )
0650 % Convert residuals to weights using the bi-square weight function.
0651 % NOTE that this function returns the square root of the weights
0652 
0653 % Only use non-NaN residuals to compute median
0654 idx = ~isnan( r );
0655 % And bound the median away from zero
0656 s = max( 1e8 * myeps, median( abs( r(idx) ) ) );
0657 % Covert the residuals to weights
0658 delta = iBisquare( r/(6*s) );
0659 % Everything with NaN residual should have zero weight
0660 delta(~idx) = 0;
0661 
0662 function b = iBisquare( x )
0663 % This is this bi-square function defined at the top of the left hand
0664 % column of page 831 in [C79]
0665 % NOTE that this function returns the square root of the weights
0666 b = zeros( size( x ) );
0667 idx = abs( x ) < 1;
0668 b(idx) = abs( 1 - x(idx).^2 );
0669 
0670 
0671 %------------------------
0672 % Tri-cubic weight function
0673 function w = iTricubeWeights( d )
0674 % Convert distances into weights using tri-cubic weight function.
0675 % NOTE that this function returns the square-root of the weights.
0676 %
0677 % Protect against divide-by-zero. This can happen if more points than the span
0678 % are coincident.
0679 maxD = max( d );
0680 if maxD > 0
0681     d = d/max( d );
0682 end
0683 w = (1 - d.^3).^1.5;
smooth

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
