Code covered by the BSD License

# running average and standard deviation of vectors

### Stephan Koehler (view profile)

09 Aug 2013 (Updated )

computes the running average and standard deviation of vectors or matrices for a given window size.

running_mean_std( varargin )
```function varargout = running_mean_std( varargin )
% [running_avg, running_std] = running_mean_std( data, filter_width )
% returns the running mean and standard deviations of either a vecotor or matrix.
%
%for the case of data = vector of length N, it will return a running average and
%standard deviation for a window of width filter_width. The output vectors have length N; however,
% to avoid edge effects NAN's are placed in the front and end: running_avg(
% [1:N, end+1-N:end]) = NAN.
%
%for the case of data = vector of 2XN or NX2, it will return a running average and
%standard deviations for the first row (column) and the second row (colum) for a window of width filter_width.
%The output vectors have length N; however, to avoid edge effects NAN's are placed in the front and end: running_avg([1:N, end+1-N:end]) = NAN.
%
%for the case of data = matrix MXN (i.e. a b/w picture), it will return a running average and
%standard deviations for a filter, which is a filled circle with diameter filter_width.
%
%
%
%to see 3 examples, just run running_mean_std without arguments
if nargin == 0
filter_width = 5;
x = -10:20;
x(x<0) = x(x<0) + rand(1, sum(x<0) );
y = x;
y(x<0) = x(x<0) + rand(1, sum(x<0) );
y( x>=0) = x(x>=0).^2 + rand(1, sum(x>=0) );
figure(1);
disp( 'example 1' );
data = x
filter_width
running_mean_std( y, filter_width );
title( 'example 1' );
figure(2);
disp( 'example 2' );
data = [x; y]
filter_width
running_mean_std( [ x; y], filter_width );
title( 'example 2' );

filter_width = 20;
figure(3);
disp( 'example 3' );
filter_width
imdata = imdata(:,:,1);
running_mean_std( imdata, filter_width );
return
end
if size( varargin{1}, 1 ) <= 2
data = varargin{1}';
else
data = double( varargin{1} );
end
filter_width = varargin{2};

if size( data, 2) == 1
target = ones(filter_width, 1)/filter_width;
running_avg = conv2( data, target, 'same' );
running_std = sqrt( conv2( (data - running_avg ).^2, target, 'same' ) );
% cropping invalid region
running_avg([1:filter_width, end-filter_width+1:end]) = nan;
running_std([1:filter_width, end-filter_width+1:end]) = nan;
elseif size( data, 2) == 2
target = ones(filter_width, 1)/filter_width;
running_avg = [conv2( data(:,1), target, 'same' ), conv2( data(:,2), target, 'same' )];
running_std = [sqrt( conv2( (data(:,1) - running_avg(:,1) ).^2, target, 'same' ) ), sqrt( conv2( (data(:,2) - running_avg(:,2) ).^2, target, 'same' ) ) ];
% cropping invalid region
running_avg([1:filter_width, end-filter_width+1:end], :) = nan;
running_std([1:filter_width, end-filter_width+1:end], :) = nan;
else
[r, c] = meshgrid( -ceil(filter_width/2):ceil(filter_width/2) );
target = double( r.^2+c.^2 <= ceil(filter_width/2)^2 );
target = target/sum(target(:));
running_avg = conv2( data, target, 'same' );
running_std = sqrt( conv2( (data - running_avg ).^2, target, 'same' ) );
% cropping invalid region
running_avg([1:filter_width, end-filter_width+1:end], :) = nan;
running_std([1:filter_width, end-filter_width+1:end], :) = nan;
running_avg( :, [1:filter_width, end-filter_width+1:end]) = nan;
running_std( :,[1:filter_width, end-filter_width+1:end]) = nan;
end

varargout = { running_avg, running_std };
varargout = varargout(1:nargout);
if nargout == 0
clf;
if size( data, 2) == 1
errorbar( 1:numel(data), running_avg, running_std, '.r' );
xlabel( 'pt #' );
ylabel( 'y' );
elseif size( data, 2) == 2
%%         ax = gca;
pt = plot( running_avg(:,1), running_avg(:,2), '.r' );
hold on;
tmp = plot( running_avg(:,[1 1])', [running_avg(:,2) - running_std(:,2), running_avg(:,2) + running_std(:,2)]', '-k' );
pt(end+1) = tmp(1);
tmp = plot( [running_avg(:,1)-running_std(:,1), running_avg(:,1)+running_std(:,1)]', running_avg(:,[2 2])', '-b' );
pt(end+1) = tmp(1);
axis tight;
lg = legend( pt, {'avg', 'std along y', 'std along x'  } );
set( lg, 'location', 'eastoutside' );
else
subplot( 1, 3, 1);
imagesc( data );
title( 'raw' );
subplot( 1,3,2);
imagesc( running_avg )
title( 'running average' );
subplot( 1,3,3);
imagesc( running_std );
title( 'running standard deviation' );
end
end
```