1 function [y]=median(x,DIM)
2 % MEDIAN data elements,
8 % N: median of N-th dimension
9 % default or []: first DIMENSION, with more than 1 element
12 % - can deal with NaN's (missing values)
13 % - accepts dimension argument like in Matlab in Octave, too.
14 % - compatible to Matlab and Octave
16 % see also: SUMSKIPNAN
18 % This program is free software; you can redistribute it and/or modify
19 % it under the terms of the GNU General Public License as published by
20 % the Free Software Foundation; either version 2 of the License, or
21 % (at your option) any later version.
23 % This program is distributed in the hope that it will be useful,
24 % but WITHOUT ANY WARRANTY; without even the implied warranty of
25 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 % GNU General Public License for more details.
28 % You should have received a copy of the GNU General Public License
29 % along with this program; If not, see <http://www.gnu.org/licenses/>.
31 % $Id: median.m 9033 2011-11-08 20:58:07Z schloegl $
32 % Copyright (C) 2000-2003,2009 by Alois Schloegl <alois.schloegl@gmail.com>
33 % This function is part of the NaN-toolbox
34 % http://pub.ist.ac.at/~schloegl/matlab/NaN/
36 global FLAG_NANS_OCCURED;
38 % check dimension of x
41 % find the dimension for median
44 if isempty(DIM), DIM=1; end;
48 sz = [sz,ones(1,DIM-length(sz))];
51 D1 = prod(sz(1:DIM-1));
53 D3 = prod(sz(DIM+1:length(sz)));
54 D0 = [sz(1:DIM-1),1,sz(DIM+1:length(sz))];
56 flag_MexKthElement = exist('kth_element','file')==3;
60 xi = k + l * D1*sz(DIM) + 1 ;
62 t = x(xi+(0:sz(DIM)-1)*D1);
68 elseif flag_MexKthElement,
69 if (D1==1) t = t+0.0; end; % make sure a real copy (not just a reference to x) is used
70 flag_KthE = 0; % fast kth_element can be used, because t does not contain any NaN and there is need to care about in-place sorting
72 y(xo) = sum( kth_element( double(t), n/2 + [0,1], flag_KthE) ) / 2;
74 y(xo) = kth_element(double(t), (n+1)/2, flag_KthE);
79 y(xo) = (t(n/2) + t(n/2+1)) / 2;
86 FLAG_NANS_OCCURED = 1;
92 %!assert(median([1,NaN,3,inf,-inf]),2)
93 %!assert(median([1,NaN,3,inf,4,-inf]),3)