function Y = runmean(X, m, dim, modestr) ; % RUNMEAN - Very fast running mean (aka moving average) filter % For vectors, Y = RUNMEAN(X,M) computes a running mean (also known as % moving average) on the elements of the vector X. It uses a window of % 2*M+1 datapoints. M an positive integer defining (half) the size of the % window. In pseudo code: % Y(i) = sum(X(j)) / (2*M+1), for j = (i-M):(i+M), and i=1:length(X) % % For matrices, Y = RUNMEAN(X,M) or RUNMEAN(X,M,[]) operates on the first % non-singleton dimension of X. RUNMEAN(X,M,DIM) computes the running % mean along the dimension DIM. % % If the total window size (2*M+1) is larger than the size in dimension % DIM, the overall average along dimension DIM is computed. % % As always with filtering, the values of Y can be inaccurate at the % edges. RUNMEAN(..., MODESTR) determines how the edges are treated. MODESTR can be % one of the following strings: % 'edge' : X is padded with first and last values along dimension % DIM (default) % 'zero' : X is padded with zeros % 'mean' : X is padded with the mean along dimension DIM % % X should not contains NaNs, yielding an all NaN result. NaNs can be % replaced by using, e.g., "inpaint_nans" created by John D'Errico. % % Examples % runmean([1:5],1) % % -> 1.33 2 3 4 4.67 % runmean([1:5],1,'mean') % % -> 2 2 3 4 4 % runmean([2:2:10],1,1) % dimension 1 is larger than 2*(M=1)+1 ... % % -> 2 4 6 8 10 % runmean(ones(10,7),3,2,'zero') ; % along columns, using mode 'zero' % runmean(repmat([1 2 4 8 NaN 5 6],5,1),2,2) ; % % -> all NaN result % A = rand(10,10) ; A(2,7) = NaN ; % runmean(A,3,2) ; % % -> column 7 is all NaN % runmean(1:2:10,100) % mean % % -> 5 5 5 5 5 % % This is an incredibly fast implementation of a running mean, since % execution time does not depend on the size of the window. % % See also MEAN, FILTER
% for Matlab R13 % version 3.0 (sep 2006) % Jos van der Geest % email: jos@jasen.nl
% History: % 1.0 (2003) created, after a snippet from Peter Acklam (?) % 1.1 (feb 2006) made suitable for the File Exchange (extended help and % documentation) % 1.2 (feb 2006) added a warning when the window size is too big % 1.3 (feb 2006) improved help section % 2.0 (sep 2006) working across a dimension of a matrix. % 3.0 (sep 2006) several treatments of the edges.
% Acknowledgements: (sep 2006) Thanks to Markus Hahn for the idea of % working in multi-dimensions and the way to treat edges.
error(nargchk(2,4,nargin)) ;
if ~isnumeric(m) || (numel(m) ~= 1) || (m < 0) || fix(m) ~= m, error('The window size (M) should be a positive integer') ; end
if nargin == 2, dim = [] ; modestr = 'edge' ; elseif nargin==3, if ischar(dim), % no dimension given modestr = dim ; dim = [] ; else modestr = 'edge' ; end end
modestr = lower(modestr) ;
% check mode specifier if ~ismember(modestr,{'edge','zero','mean'}), error('Unknown mode') ; end
szX = size(X) ; if isempty(dim), dim = min(find(szX>1)) ; end
if m == 0 || dim > ndims(X), % easy Y = X ; else mm = 2*m+1 ; if mm >= szX(dim), % if the window is larger than X, average all sz2 = ones(size(szX)) ; sz2(dim) = szX(dim) ; Y = repmat(mean(X,dim),sz2) ; else % here starts the real stuff % shift dimensions so that the desired dimensions comes first [X, nshifts] = shiftdim(X, dim-1); szX = size(X) ; % make the rest of the dimensions columns, so we have a 2D matrix % (suggested of Markus Hahn) X = reshape(X,szX(1),[]) ; % select how to pad the matrix switch (modestr), case 'edge' % pad with first and last elements Xfirst = repmat(X(1,:),m,1) ; Xlast = repmat(X(end,:),m,1) ; case 'zero' % pad with zeros Xfirst = zeros(m,1) ; Xlast= zeros(m,1) ; case 'mean', % pad with the average Xfirst = repmat(mean(X,1),m,1) ; Xlast = Xfirst ; end % pad the array Y = [zeros(1,size(X,2)) ; Xfirst ; X ; Xlast] ; % the cumsum trick (by Peter Acklam ?) Y = cumsum(Y,1) ; Y = (Y(mm+1:end,:)-Y(1:end-mm,:)) ./ mm ;
% reshape into original size Y = reshape(Y,szX) ; % and re-shift the dimensions Y = shiftdim(Y,ndims(Y)-nshifts) ; end end
% ===================== % CODE OF VERSION 1.3 % =====================
% function Y = runmean(X,m) ; % % RUNMEAN - Very fast running mean filter for vectors % % Y = RUNMEAN(X,M) computes a running mean on vector X using a window of % % 2*M+1 datapoints. X is a vector, and M an positive integer defining % % (half) the size of the window. In pseudo code: % % Y(i) = sum(X(j)) / (2*M+1), for j = (i-M):(i+M), and i=1:length(X) % % % % If the total window size (2M+1) is larger than the length of the vector, the overall % % average is returned. % % % % Example: % % runmean(1:10,1) % -> % % [1.3333 2 3 4 5 6 7 8 9 9.6667] % % % % This is an incredibly fast implementation of a running average, since % % execution time does not depend on the size of the window. % % % % X should not contains NaNs (a NaN will result in a all NaN result) % % At both ends the values of Y can be inaccurate, as the first and last % % values of X are used multiple times. % % % % See also MEAN % % % for Matlab R13 % % version 1.3 (feb 2006) % % Jos van der Geest % % email: jos@jasen.nl % % % History: % % 1.0 (2003) created, after a snippet from Peter Acklam (?) % % 1.1 (feb 2006) made suitable for the File Exchange (extended help and % % documentation) % % 1.2 (feb 2006) added a warning when the window size is too big % % 1.3 (feb 2006) improved help section % % error(nargchk(2,2,nargin)) ; % % sz = size(X) ; % % if numel(sz) ~= 2 || (min(sz) ~= 1), % error('X should be a vector') ; % end % % if any(isnan(X)), % error('NaNs cannot be dealt with') ; % end % % if ~isnumeric(m) || (numel(m) ~= 1) || (m < 0) || fix(m) ~= m, % error('The window size (M) should be a positive integer') ; % elseif m == 0, % Y = X ; % return ; % end % % mm = 2*m+1 ; % % if mm >= prod(sz), % % if the window is larger than X, average all % warning('Window size is larger than the length of the vector.') % Y = repmat(mean(X),sz) ; % else % % the cumsum trick ... % Y = [repmat(X(1),m,1) ; X(:) ; repmat(X(end),m,1)] ; % Y = [0 ; cumsum(Y)] ; % Y = (Y(mm+1:end)-Y(1:end-mm)) / mm ; % Y = reshape(Y,sz) ; % end