Code covered by the BSD License  

Highlights from
running average and standard deviation of vectors

image thumbnail

running average and standard deviation of vectors

by

 

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' );
    disp( 'loading stock photo ngc6543a.jpg' );
    filter_width
    imdata = imread('ngc6543a.jpg');
    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

Contact us