%Calculates the centered moving average of an n-dimensional matrix in any direction.
% 1) data = The matrix to be averaged.
% 2) window = The window size. This works best as an odd number. If an even
% number is entered then it is rounded down to the next odd number.
% 3) dim = The dimension in which you would like do the moving average.
% This is an optional input. To leave blank use  place holder in
% function call. Defaults to 1.
% 4) option = which solution algorithm to use. The default option works
% best in most situations, but option 2 works better for wide
% matrices (i.e. 1000000 x 1 or 10 x 1000 x 1000) when solved in the
% shorter dimension. Data size where option 2 is more efficient will
% vary from computer to computre.This is an optional input. To leave
% blank use  place holder in function call. Defaults to 1.
% Calculate column moving average of 10000 x 10 matrix with a window size
% of 5 in the 1st dimension using algorithm option 1.
% Moving mean for each element uses data centered on that element and
% incorporates (window-1)/2 elements before and after the element.
% Function is broken into two parts. The 1d-2d solution, and the
% n-dimensional solution. The 1d-2d solution is the fastest that I have
% been able to come up with, whereas the n-dimensional solution trades
% speed for versatility.
% Function includes some code at the end so that the user can do their
% own speed testing using the TIMEIT function.
% Has been heavily tested in 1d-2d case, and lightly tested in 3d case.
% Should work in n-dimensional space, but has not been tested in more
% than 3 dimensions other than to make sure it does not return an error.
% Has not been tested for complex inputs.
%errors for leaving out data or window inputs
error('no input data')
error('no window size specified')
%defaults the dimension and option to 1 if it is not specified
if nargin<3 || isempty(dim)
if nargin<4 || isempty(option)
%rounds even window sizes down to next lowest odd number
%Calculates the number of elements in before and after the central element
%to incorporate in the moving mean. Round command is just present to deal
%with the potential problem of division leaving a very small decimal, ie.
%calculates the size of the input data set
%returns error messages for incorrect inputs
error('dimension number must be integer')
if ((ndims(data)<=2 && (option<1 || option>3)) || (ndims(data)>=3 && (option<1 ||option>2))) ...
error('invalid algorithm option')
error('window size must be integer')
error('dimension number exceeds number of dimensions in input matrix')
error('dimension number must be >=1')
error('window is too large')
%Solution for 1d-2d situation. Uses vector operations to optimize
%To simplify algorithm the problem is always solved with dim=1.
%If user input is dim=2 then the data is transposed to calculate the
%solution in dim=1
%The three best solutions I came up to for the 1d-2d problem. I kept
%them in here to preserve the code incase I want to use some of it
%again. I have found that option 1 is consistenetly the fastest.
%option 1, works best for most data sets
%Computes the beginning and ending column for each moving
%average compuation. Divide is the number of elements that are
%incorporated in each moving average.
%Calculates the moving average by calculating the sum of elements
%from the start row to the stop row for each central element,
%and then dividing by the number of elements used in that sum
%to get the average for that central element.
%Implemented by calculating the moving sum of the full data
%set. Cumulative sum for each central element is calculated by
%subtracting the cumulative sum for the row before the start row
%from the cumulative sum for the stop row. Row references are
%adusted to take into account the fact that you can now
%reference a row<1. Divides the series of cumulative sums for
%by the number of elements in each sum to get the moving
%option 2, Can be faster than option 1 in very wide matrices
%(i.e. 100 x 1000 or l000 x 10000) when solving in the shorter
%dimension, but in general it is slower than option 1.
%Uses a for loop to calculate the data from the start row to
%the stop row, and then divides by the number of rows
%option 3, Creates a sparse matrix to indicate which rows to include
%data from for the moving average for each central element.
%Uses matrix multiplication to calculte the cumulative sum for
%each central element. Then divides by the number of elements
%used in each cumulative sum to get the average for each
%Elegant but slow. I think this is because bsxfun has to do
%a lot of matrix replication. Also, it is easy for the use matrix to
%exceed the maximum matrix dimensions for MATLAB.
use=(bsxfun(@ge,baseline,start) & bsxfun(@le,baseline,stop));
%uses matrix multiplication instead of repmat or bsxfun, slower than
%use=(baseline_mat>=start_mat & baseline_mat<=stop_mat);
%Transposes the matrix if the problem was solved in dimension=2 so that
%the ouput is correctly oriented. Undoes the initial transposition.
%Solution of the n-dimensional problem. This solution is
%slower than the 1d-2d solution in a 1d-2d problem, but depending on the
%dimensions of the input array can be faster for the same number of
%elements in a 3d+ problem (ie. 1000000 x 1 vs 100 x 100 x 100).
%Either way it is much more versatile. The algorithm can solve for a
%moving average in any dimension of an n-dimension input matrix.
%Two solutions of the problem. In general option 1 is faster than
%option 2, but I think that for some matrix sizes option 2 will be
%Option 1 is based on option 1 for the 1d-2d solution. In
%general this is the faster solution. For a more extensive
%explanation of the code see option 1 in the 1d-2d section.
%Calculates the start and stop positions, and creates a vector
%of how many elements are in each moving sum.
%Reorients the start, stop, and divide vectors in the direction
%of the specified dimension.
%Creates a vector to use to reshape the start, stop, and divide
%vectors. Puts one in for the each dimension excpet for the
%specified dimension, which is set to the size of the input
%matrix in that direction. ex. [1 1 100 1] for a
%100 x 100 x 100 x 100 input matrix, solving in the 3rd
%Uses the vector created above to reorient the start, stop, and
%divide vectors in the correct direction.
%Calculates the cumulative sum of the data in the specified
%I have found the most versatile way of writing an algorithm
%that will solve in any dimension of an unknown dimensional
%input matrix is to built some functions with text that is
%based on the input parameters.
%This text builds the function that subtracts the cumulative
%sum at the stop position from the cumulative sum at one
%element before the start position for each central element.
%Compensates for not have an index less than 1.
function_text=['@(CumulativeSum,start,stop) CumulativeSum(' repmat(':,',1,dim-1) 'stop,' repmat(':,',1,ndims(data)-dim)];
function_text=[function_text ')-CumulativeSum(' repmat(':,',1,dim-1) 'max(1,start-1),' repmat(':,',1,ndims(data)-dim)];
%Converts the string above into an anonymous function with the
%Runs the function built above to determine the cumulative sum
%between start and stop.
%Text to make a function that corrects for the fact that the
%cumulative sum for moving sums that start at position 1 have
%the first data position subtracted out. This is becasue you
%can not have a matrix index less than 1. Uses bsxfun to
%multiply start==1 by the first plane of data (in the
%appropriate dimension) and then adding that to temp_sum. Need
%to do it this way vs using something based on reshaping
%data(repmat(start,inv_temp_dims)==1) because the correction
%needs to be the same size as temp_sum. This is needed b/c I
%do not know how to use text to set the variable elements that
%will accept the addition (i.e.
correction_text=['@(temp_sum,start,data) temp_sum+bsxfun(@times,start==1,data(' repmat(':,',1,dim-1) '1,' repmat(':,',1,ndims(data)-dim)];
%Converts the above text into a function.
%Applies the correction describe previously to the temp_sum
%Divides the temp_sum matrix by the divide vector to
%calculate the moving average for each central elemant.
%An older solution to the n-dimensional moving average problem.
%Based on the approach in option 2 in the 1d-2d solution.
%Depending on matrix dimensions this might be faster than
%option 1 in certain situations (ex input matrix size is 10 x
%1000 x 1000 and the average is done in dimension 1).
%Builds a string for the function used in the for loop.
%Basiscally it builds something like sum(data(:,:,start:stop,:),3) but
%is able to work in any dimensions and have start:stop be in
%any of those dimensions.
%Converts the text above into a function
%For loop determines the start position, stop position, and
%number of elements for the moving average of each central
%element. Calculates the moving sum for each central element
%in that position, and divides by the number of elements to get
%the moving average.
%This was the most versatile way of reassembling all of the
%moving average matrices back into the final output matrix,
%although it is slow b/c the result matrix keeps resizing.
%code for speed testing
%generates a series of random vectors with 10 to 10000000 cells and then
%times the moving average at each size. Time for each input size is
%1d-2d timer test
%data_power=[1 2 3 4 5 6 7];
% f=@() movingmean(data,window,dim);
%3d timer test
%data_power=[1 2 3];
% f=@() movingmean(data,window,dim);