time   calls  line
                  1 function y = median(x,dim)
                  2 %MEDIAN Median value.
                  3 %   For vectors, MEDIAN(a) is the median value of the elements in a.
                  4 %   For matrices, MEDIAN(A) is a row vector containing the median
                  5 %   value of each column.  For N-D arrays, MEDIAN(A) is the median
                  6 %   value of the elements along the first non-singleton dimension
                  7 %   of A.
                  8 %
                  9 %   MEDIAN(A,DIM) takes the median along the dimension DIM of A.
                 10 %
                 11 %   Example: If A = [1 2 4 4; 3 4 6 6; 5 6 8 8; 5 6 8 8];
                 12 %
                 13 %   then median(A) is [4 5 7 7] and median(A,2) 
                 14 %   is [3 5 7 7].'
                 15 %
                 16 %   Class support for input A:
                 17 %      float: double, single
                 18 %
                 19 %   See also MEAN, STD, MIN, MAX, VAR, COV, MODE.
                 20 
                 21 %   Copyright 1984-2009 The MathWorks, Inc.
                 22 %   $Revision: 5.15.4.7 $  $Date: 2009/09/03 05:19:03 $
                 23 
                 24 %   Calculation method for even lists: b - (b-a)/2.
                 25 %   This method reduces the likelihood of rounding errors.
                 26 
           254   27 nInputs = nargin; 
                 28 
           254   29 if isempty(x) 
                 30   if nInputs == 1
                 31 
                 32     % The output size for [] is a special case when DIM is not given.
                 33     if isequal(x,[]), y = nan(1,class(x)); return; end
                 34 
                 35     % Determine first nonsingleton dimension
                 36     dim = find(size(x)~=1,1);
                 37 
                 38   end
                 39 
                 40   s = size(x);
                 41   if dim <= length(s)
                 42      s(dim) = 1;                  % Set size to 1 along dimension
                 43   end
                 44   y = nan(s,class(x));
                 45 
           254   46 elseif nInputs == 1 && isvector(x) 
                 47   % If input is a vector, calculate single value of output.
  0.45     206   48   x = sort(x); 
  0.01     206   49   nCompare = numel(x); 
           206   50   half = floor(nCompare/2); 
           206   51   y = x(half+1); 
           206   52   if 2*half == nCompare        % Average if even number of elements 
  0.02     131   53     y = meanof(x(half),y); 
           131   54   end 
  0.01     206   55   if isnan(x(nCompare))        % Check last index for NaN 
                 56     y = nan(class(x));
                 57   end
            48   58 else 
            48   59   if nInputs == 1              % Determine first nonsingleton dimension 
            48   60     dim = find(size(x)~=1,1); 
            48   61   end 
                 62 
            48   63   s = size(x); 
                 64 
  0.01      48   65   if dim > length(s)           % If dimension is too high, just return input. 
                 66     y = x;
                 67     return
                 68   end
                 69 
                 70   % Sort along given dimension
  0.04      48   71   x = sort(x,dim); 
                 72 
            48   73   nCompare = s(dim);           % Number of elements used to generate a median 
            48   74   half = floor(nCompare/2);    % Midway point, used for median calculation 
                 75 
            48   76   if dim == 1 
                 77     % If calculating along columns, use vectorized method with column
                 78     % indexing.  Reshape at end to appropriate dimension.
            48   79     y = x(half+1,:); 
            48   80     if 2*half == nCompare 
            39   81       y = meanof(x(half,:),y); 
            39   82     end 
                 83 
            48   84     y(isnan(x(nCompare,:))) = NaN;   % Check last index for NaN 
                 85 
                 86   elseif dim == 2 && length(s) == 2
                 87     % If calculating along rows, use vectorized method only when possible.
                 88     % This requires the input to be 2-dimensional.   Reshape at end.
                 89     y = x(:,half+1);
                 90     if 2*half == nCompare
                 91       y = meanof(x(:,half),y);
                 92     end
                 93 
                 94     y(isnan(x(:,nCompare))) = NaN;   % Check last index for NaN
                 95 
                 96   else
                 97     % In all other cases, use linear indexing to determine exact location
                 98     % of medians.  Use linear indices to extract medians, then reshape at
                 99     % end to appropriate size.
                100     cumSize = cumprod(s);
                101     total = cumSize(end);            % Equivalent to NUMEL(x)
                102     numMedians = total / nCompare;
                103 
                104     numConseq = cumSize(dim - 1);    % Number of consecutive indices
                105     increment = cumSize(dim);        % Gap between runs of indices
                106     ixMedians = 1;
                107 
                108     y = repmat(x(1),numMedians,1);   % Preallocate appropriate type
                109 
                110     % Nested FOR loop tracks down medians by their indices.
                111     for seqIndex = 1:increment:total
                112       for consIndex = half*numConseq:(half+1)*numConseq-1
                113         absIndex = seqIndex + consIndex;
                114         y(ixMedians) = x(absIndex);
                115         ixMedians = ixMedians + 1;
                116       end
                117     end
                118 
                119     % Average in second value if n is even
                120     if 2*half == nCompare
                121       ixMedians = 1;
                122       for seqIndex = 1:increment:total
                123         for consIndex = (half-1)*numConseq:half*numConseq-1
                124           absIndex = seqIndex + consIndex;
                125           y(ixMedians) = meanof(x(absIndex),y(ixMedians));
                126           ixMedians = ixMedians + 1;
                127         end
                128       end
                129     end
                130 
                131     % Check last indices for NaN
                132     ixMedians = 1;
                133     for seqIndex = 1:increment:total
                134       for consIndex = (nCompare-1)*numConseq:nCompare*numConseq-1
                135         absIndex = seqIndex + consIndex;
                136         if isnan(x(absIndex))
                137           y(ixMedians) = NaN;
                138         end
                139         ixMedians = ixMedians + 1;
                140       end
                141     end
                142   end
                143 
                144   % Now reshape output.
            48  145   s(dim) = 1; 
            48  146   y = reshape(y,s); 
            48  147 end 

Other subfunctions in this file are not included in this listing.