time   calls  line
                  1 function varargout = interp1(varargin)
                  2 %INTERP1 1-D interpolation (table lookup)
                  3 %   YI = INTERP1(X,Y,XI) interpolates to find YI, the values of the
                  4 %   underlying function Y at the points in the array XI. X must be a
                  5 %   vector of length N.
                  6 %   If Y is a vector, then it must also have length N, and YI is the
                  7 %   same size as XI.  If Y is an array of size [N,D1,D2,...,Dk], then
                  8 %   the interpolation is performed for each D1-by-D2-by-...-Dk value
                  9 %   in Y(i,:,:,...,:).
                 10 %   If XI is a vector of length M, then YI has size [M,D1,D2,...,Dk].
                 11 %   If XI is an array of size [M1,M2,...,Mj], then YI is of size
                 12 %   [M1,M2,...,Mj,D1,D2,...,Dk].
                 13 %
                 14 %   YI = INTERP1(Y,XI) assumes X = 1:N, where N is LENGTH(Y)
                 15 %   for vector Y or SIZE(Y,1) for array Y.
                 16 %
                 17 %   Interpolation is the same operation as "table lookup".  Described in
                 18 %   "table lookup" terms, the "table" is [X,Y] and INTERP1 "looks-up"
                 19 %   the elements of XI in X, and, based upon their location, returns
                 20 %   values YI interpolated within the elements of Y.
                 21 %
                 22 %   YI = INTERP1(X,Y,XI,METHOD) specifies alternate methods.
                 23 %   The default is linear interpolation. Use an empty matrix [] to specify
                 24 %   the default. Available methods are:
                 25 %
                 26 %     'nearest'  - nearest neighbor interpolation
                 27 %     'linear'   - linear interpolation
                 28 %     'spline'   - piecewise cubic spline interpolation (SPLINE)
                 29 %     'pchip'    - shape-preserving piecewise cubic interpolation
                 30 %     'cubic'    - same as 'pchip'
                 31 %     'v5cubic'  - the cubic interpolation from MATLAB 5, which does not
                 32 %                  extrapolate and uses 'spline' if X is not equally
                 33 %                  spaced.
                 34 %
                 35 %   YI = INTERP1(X,Y,XI,METHOD,'extrap') uses the interpolation algorithm
                 36 %   specified by METHOD to perform extrapolation for elements of XI outside
                 37 %   the interval spanned by X. 
                 38 %
                 39 %   YI = INTERP1(X,Y,XI,METHOD,EXTRAPVAL) replaces the values outside of the 
                 40 %   interval spanned by X with EXTRAPVAL.  NaN and 0 are often used for 
                 41 %   EXTRAPVAL.  The default extrapolation behavior with four input arguments 
                 42 %   is 'extrap' for 'spline' and 'pchip' and EXTRAPVAL = NaN for the other 
                 43 %   methods.
                 44 %
                 45 %   PP = INTERP1(X,Y,METHOD,'pp') will use the interpolation algorithm specified
                 46 %   by METHOD to generate the ppform (piecewise polynomial form) of Y. The
                 47 %   method may be any of the above METHOD except for 'v5cubic'. PP may then
                 48 %   be evaluated via PPVAL. PPVAL(PP,XI) is the same as
                 49 %   INTERP1(X,Y,XI,METHOD,'extrap').
                 50 %
                 51 %   For example, generate a coarse sine curve and interpolate over a
                 52 %   finer abscissa:
                 53 %       x = 0:10; y = sin(x); xi = 0:.25:10;
                 54 %       yi = interp1(x,y,xi); plot(x,y,'o',xi,yi)
                 55 %
                 56 %   For a multi-dimensional example, we construct a table of functional
                 57 %   values:
                 58 %       x = [1:10]'; y = [ x.^2, x.^3, x.^4 ];
                 59 %       xi = [ 1.5, 1.75; 7.5, 7.75]; yi = interp1(x,y,xi);
                 60 %
                 61 %   creates 2-by-2 matrices of interpolated function values, one matrix for
                 62 %   each of the 3 functions. yi will be of size 2-by-2-by-3.
                 63 %
                 64 %   Class support for inputs X, Y, XI, EXTRAPVAL:
                 65 %      float: double, single
                 66 %
                 67 %   See also INTERP1Q, INTERPFT, SPLINE, PCHIP, INTERP2, INTERP3, INTERPN, PPVAL.
                 68 
                 69 %   Copyright 1984-2010 The MathWorks, Inc.
                 70 %   $Revision: 5.41.4.17 $  $Date: 2010/11/22 02:46:39 $
                 71 %
                 72 % Determine input arguments.
                 73 
                 74 % The following syntaxes do not have X as the first input:
                 75 % INTERP1(Y,XI)
                 76 % INTERP1(Y,XI,METHOD)
                 77 % INTERP1(Y,XI,METHOD,'extrap')
                 78 % INTERP1(Y,XI,METHOD,EXTRAPVAL)
          1778   79 xOffset = 1; 
  0.01    1778   80 if (nargin==2) || ... 
                 81    (nargin==3 && ischar(varargin{3}))  || ...
                 82    (nargin==4 && ~(ischar(varargin{4}) || isempty(varargin{4})) || ...
                 83    (nargin==4 && strcmp(varargin{4}, 'extrap')));
            12   84     xOffset = 0; 
            12   85 end 
                 86 
          1778   87 ppOutput = false; 
                 88 % PP = INTERP1(X,Y,METHOD,'pp')
          1778   89 if nargin>=4 && ischar(varargin{3}) && isequal('pp',varargin{4}) 
                 90     ppOutput = true;
                 91     if (nargin > 4)
                 92         error(message('MATLAB:interp1:ppOutput'))
                 93     end
                 94 end
                 95 
                 96 % Process Y in INTERP1(Y,...) and INTERP1(X,Y,...)
  0.02    1778   97 y = varargin{1+xOffset}; 
          1778   98 siz_y = size(y); 
                 99 % y may be an ND array, but collapse it down to a 2D yMat. If yMat is
                100 % a vector, it is a column vector.
          1778  101 if isvector(y) 
          1776  102     if isrow(y) 
                103         % Prefer column vectors for y
            15  104         yMat = y.'; 
            15  105         n = siz_y(2); 
  0.01    1761  106     else 
          1761  107         yMat = y; 
          1761  108         n = siz_y(1); 
  0.01    1761  109     end 
          1776  110     ds = 1; 
          1776  111     prodDs = 1; 
             2  112 else 
             2  113     n = siz_y(1); 
             2  114     ds = siz_y(2:end); 
             2  115     prodDs = prod(ds); 
             2  116     yMat = reshape(y,[n prodDs]); 
             2  117 end 
                118 
                119 % Process X in INTERP1(X,Y,...), or supply default for INTERP1(Y,...)
          1778  120 if xOffset 
          1766  121     x = varargin{xOffset}; 
          1766  122     if ~isvector(x) 
                123         error(message('MATLAB:interp1:Xvector'));
                124     end
  0.01    1766  125     if length(x) ~= n 
                126         if isvector(y)
                127             error('MATLAB:interp1:YInvalidNumRows', ...
                128                 'X and Y must be of the same length.')
                129         else
                130             error('MATLAB:interp1:YInvalidNumRows', ...
                131                 'LENGTH(X) and SIZE(Y,1) must be the same.');
                132         end
                133     end
                134     % Prefer column vectors for x
          1766  135     xCol = x(:); 
            12  136 else 
            12  137     xCol = (1:n)'; 
            12  138 end 
                139 
                140 % Process XI in INTERP1(Y,XI,...) and INTERP1(X,Y,XI,...)
                141 % Avoid syntax PP = INTERP1(X,Y,METHOD,'pp')
          1778  142 if ~ppOutput 
          1778  143     xi = varargin{2+xOffset}; 
          1778  144     siz_xi = size(xi); 
                145     % xi may be an ND array, but flatten it to a column vector xiCol
  0.01    1778  146     xiCol = xi(:); 
                147     % The size of the output YI
          1778  148     if isvector(y) 
                149         % Y is a vector so size(YI) == size(XI)
          1776  150         siz_yi = siz_xi; 
             2  151     else 
             2  152         if isvector(xi) 
                153             % Y is not a vector but XI is
             2  154             siz_yi = [length(xi) ds]; 
                155         else
                156             % Both Y and XI are non-vectors
                157             siz_yi = [siz_xi ds];
                158         end
             2  159     end 
          1778  160 end 
                161 
  0.01    1778  162 if xOffset && ~isreal(x) 
                163     error(message('MATLAB:interp1:ComplexX'))
                164 end
                165 
          1778  166 if ~ppOutput && ~isreal(xi) 
                167     error(message('MATLAB:interp1:ComplexInterpPts'))
                168 end
                169 
                170 % Error check for NaN values in X and Y
                171 % check for NaN's
  0.04    1778  172 if xOffset && (any(isnan(xCol))) 
                173     error(message('MATLAB:interp1:NaNinX'));
                174 end
                175 
                176 % NANS are allowed as a value for F(X), since a function may be undefined
                177 % for a given value.
  0.02    1778  178 if any(isnan(yMat(:))) 
                179     warning(message('MATLAB:interp1:NaNinY'));
                180 end
                181 
  0.01    1778  182 if (n < 2) 
                183     if ppOutput || ~isempty(xi)
                184         error(message('MATLAB:interp1:NotEnoughPts'))
                185     else
                186         yi = zeros(siz_yi,superiorfloat(x,y,xi));
                187         varargout{1} = yi;
                188         return
                189     end
                190 end
                191 
                192 % Process METHOD in
                193 % PP = INTERP1(X,Y,METHOD,'pp')
                194 % YI = INTERP1(Y,XI,METHOD,...)
                195 % YI = INTERP1(X,Y,XI,METHOD,...)
                196 % including explicit specification of the default by an empty input.
          1778  197 if ppOutput 
                198     if isempty(varargin{3})
                199         method = 'linear';
                200     else
                201         method = varargin{3};
                202     end
          1778  203 else 
  0.01    1778  204     if nargin >= 3+xOffset && ~isempty(varargin{3+xOffset}) 
  0.01    1727  205         method = varargin{3+xOffset}; 
            51  206     else 
            51  207         method = 'linear'; 
            51  208     end 
  0.01    1778  209 end 
                210     
                211 % The v5 option, '*method', asserts that x is equally spaced.
          1778  212 eqsp = (method(1) == '*'); 
  0.01    1778  213 if eqsp 
                214     method(1) = [];
                215 end
                216 
                217 % INTERP1([X,]Y,XI,METHOD,'extrap') and INTERP1([X,]Y,Xi,METHOD,EXTRAPVAL)
  0.01    1778  218 if ~ppOutput 
  0.01    1778  219     if nargin >= 4+xOffset 
          1700  220         extrapval = varargin{4+xOffset}; 
            78  221     else 
            78  222         switch method(1) 
            78  223             case {'s','p','c'} 
            25  224                 extrapval = 'extrap'; 
            53  225             otherwise 
            53  226                 extrapval = NaN; 
            53  227         end 
            78  228     end 
          1778  229 end 
                230 
                231 % Start the algorithm
                232 % We now have column vector xCol, column vector or 2D matrix yMat and
                233 % column vector xiCol.
          1778  234 if xOffset 
  0.02    1766  235     if ~eqsp 
          1766  236         h = diff(xCol); 
  0.03    1766  237         eqsp = (norm(diff(h),Inf) <= eps(norm(xCol,Inf))); 
  0.01    1766  238         if any(~isfinite(xCol)) 
                239             eqsp = 0; % if an INF in x, x is not equally spaced
                240         end
  0.03    1766  241     end 
          1766  242     if eqsp 
          1509  243         h = (xCol(n)-xCol(1))/(n-1); 
  0.02    1509  244     end 
            12  245 else 
            12  246     h = 1; 
            12  247     eqsp = 1; 
            12  248 end 
  0.01    1778  249 if any(h < 0) 
                250     [xCol,p] = sort(xCol);
                251     yMat = yMat(p,:);
                252     if eqsp
                253         h = -h;
                254     else
                255         h = diff(xCol);
                256     end
                257 end
  0.03    1778  258 if any(h == 0) 
                259     error(message('MATLAB:interp1:RepeatedValuesX'));
                260 end
                261 
                262 % PP = INTERP1(X,Y,METHOD,'pp')
          1778  263 if nargin==4 && ischar(varargin{3}) && isequal('pp',varargin{4})  
                264     % obtain pp form of output
                265     pp = ppinterp;
                266     varargout{1} = pp;
                267     return
                268 end
                269 
                270 % Interpolate
          1778  271 numelXi = length(xiCol); 
  0.03    1778  272 p = []; 
          1778  273 switch method(1) 
          1778  274     case 's'  % 'spline' 
                275         % spline is oriented opposite to interp1
  0.94       7  276         yiMat = spline(xCol.',yMat.',xiCol.').'; 
                277 
  0.01    1771  278     case {'c','p'}  % 'cubic' or 'pchip' 
                279         % pchip is oriented opposite to interp1
  1.28      18  280         yiMat = pchip(xCol.',yMat.',xiCol.').'; 
                281 
          1753  282     otherwise % 'nearest', 'linear', 'v5cubic' 
  0.09    1753  283         yiMat = zeros(numelXi,prodDs,superiorfloat(xCol,yMat,xiCol)); 
  0.03    1753  284         if ~eqsp && any(diff(xiCol) < 0) 
                285             [xiCol,p] = sort(xiCol);
          1753  286         else 
  0.17    1753  287             p = 1:numelXi; 
  0.01    1753  288         end 
                289 
                290         % Find indices of subintervals, x(k) <= u < x(k+1),
                291         % or u < x(1) or u >= x(m-1).
  0.01    1753  292         if isempty(xiCol) 
                293             k = xiCol;
          1753  294         elseif eqsp 
  1.39    1509  295             k = min(max(1+floor((xiCol-xCol(1))/h),1),n-1); 
           244  296         else 
  0.09     244  297             [~,k] = histc(xiCol,xCol); 
  0.01     244  298             k(xiCol<xCol(1) | ~isfinite(xiCol)) = 1; 
  0.01     244  299             k(xiCol>=xCol(n)) = n-1; 
           244  300         end 
                301 
  0.01    1753  302         switch method(1) 
          1753  303             case 'n'  % 'nearest' 
                304                 i = find(xiCol >= (xCol(k)+xCol(k+1))/2);
                305                 k(i) = k(i)+1;
                306                 yiMat(p,:) = yMat(k,:);
                307                 nanidx = find(isnan(xiCol));
                308                 yiMat(p(nanidx),:) = yiMat(p(nanidx),:)*NaN;
          1753  309             case 'l'  % 'linear' 
          1753  310                 if eqsp 
  1.04    1509  311                     s = (xiCol - xCol(k))/h; 
           244  312                 else 
  0.10     244  313                     s = (xiCol - xCol(k))./h(k); 
           244  314                 end 
          1753  315                 for j = 1:prodDs 
  2.74    1761  316                     yiMat(p,j) = yMat(k,j) + s.*(yMat(k+1,j)-yMat(k,j)); 
  0.01    1761  317                 end 
                318 
                319             case 'v'  % 'v5cubic'
                320                 extrapval = NaN;
                321                 if eqsp
                322                     % Data are equally spaced
                323                     s = (xiCol - xCol(k))/h;
                324                     s2 = s.*s;
                325                     s3 = s.*s2;
                326                     % Add extra points for first and last interval
                327                     yMat = [3*yMat(1,:)-3*yMat(2,:)+yMat(3,:); ...
                328                         yMat; ...
                329                         3*yMat(n,:)-3*yMat(n-1,:)+yMat(n-2,:)];
                330                     for j = 1:prodDs
                331                         yiMat(p,j) = (yMat(k,j).*(-s3+2*s2-s) + ...
                332                             yMat(k+1,j).*(3*s3-5*s2+2) + ...
                333                             yMat(k+2,j).*(-3*s3+4*s2+s) + ...
                334                             yMat(k+3,j).*(s3-s2))/2;
                335                     end
                336                 else
                337                     % Data are not equally spaced
                338                     % spline is oriented opposite to interp1
                339                     yiMat = spline(xCol.',yMat.',xiCol.').';
                340                 end
                341             otherwise
                342                 error(message('MATLAB:interp1:InvalidMethod'))
                343         end
  0.01    1753  344 end 
                345 
                346 % Override extrapolation
          1778  347 if ~isequal(extrapval,'extrap') 
            53  348     if ischar(extrapval) 
                349         error(message('MATLAB:interp1:InvalidExtrap'))
            53  350     elseif ~isscalar(extrapval) 
                351         error(message('MATLAB:interp1:NonScalarExtrapValue'))
                352     end
            53  353     if isempty(p) 
                354         p = 1 : numelXi;
                355     end
  0.20      53  356      outOfBounds = xiCol<xCol(1) | xiCol>xCol(n); 
  0.08      53  357      yiMat(p(outOfBounds),:) = extrapval; 
  0.01      53  358 end 
                359 
                360 % Reshape result, possibly to an ND array
  0.01    1778  361 yi = reshape(yiMat,siz_yi); 
  0.03    1778  362 varargout{1} = yi; 
                363 
                364 %-------------------------------------------------------------------------%
                365     function pp = ppinterp
                366         %PPINTERP ppform interpretation.
                367 
                368         switch method(1)
                369             case 'n' % nearest
                370                 breaks = [xCol(1); ...
                371                     (xCol(1:end-1)+xCol(2:end))/2; ...
                372                     xCol(end)].';
                373                 coefs = yMat.';
                374                 pp = mkpp(breaks,coefs,ds);
                375             case 'l' % linear
                376                 breaks = xCol.';
                377                 page1 = (diff(yMat)./repmat(diff(xCol),[1, prodDs])).';
                378                 page2 = (reshape(yMat(1:end-1,:),[n-1, prodDs])).';
                379                 coefs = cat(3,page1,page2);
                380                 pp = mkpp(breaks,coefs,ds);
                381             case {'p', 'c'} % pchip and cubic
                382                 pp = pchip(xCol.',reshape(yMat.',[ds, n]));
                383             case 's' % spline
                384                 pp = spline(xCol.',reshape(yMat.',[ds, n]));
                385             case 'v' % v5cubic
                386                 b = diff(xCol);
                387                 if norm(diff(b),Inf) <= eps(norm(xCol,Inf))
                388                     % data are equally spaced
                389                     a = repmat(b,[1 prodDs]).';
                390                     yReorg = [3*yMat(1,:)-3*yMat(2,:)+yMat(3,:); ...
                391                         yMat; ...
                392                         3*yMat(n,:)-3*yMat(n-1,:)+yMat(n-2,:)];
                393                     y1 = yReorg(1:end-3,:).';
                394                     y2 = yReorg(2:end-2,:).';
                395                     y3 = yReorg(3:end-1,:).';
                396                     y4 = yReorg(4:end,:).';
                397                     breaks = xCol.';
                398                     page1 = (-y1+3*y2-3*y3+y4)./(2*a.^3);
                399                     page2 = (2*y1-5*y2+4*y3-y4)./(2*a.^2);
                400                     page3 = (-y1+y3)./(2*a);
                401                     page4 = y2;
                402                     coefs = cat(3,page1,page2,page3,page4);
                403                     pp = mkpp(breaks,coefs,ds);
                404                 else
                405                     % data are not equally spaced
                406                     pp = spline(xCol.',reshape(yMat.',[ds, n]));
                407                 end
                408             otherwise
                409                 error(message('MATLAB:interp1:ppinterp:UnknownMethod'));
                410         end
                411 
                412         % Even if method is 'spline' or 'pchip', we still need to record that the
                413         % input data Y was oriented according to INTERP1's rules.
                414         % Thus PPVAL will return YI oriented according to INTERP1's rules and
                415         % YI = INTERP1(X,Y,XI,METHOD) will be the same as
                416         % YI = PPVAL(INTERP1(X,Y,METHOD,'pp'),XI)
                417         pp.orient = 'first';
                418     end % PPINTERP
                419 
  0.01    1778  420 end % INTERP1