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

Stephan Koehler

 

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