I wrote previously about how a windowed average can be used to do filtering, and made mention that a windowed average also allows for the calculation of standard deviation. In this article, I will address some acceleration for calculating standard deviation.

The first question one might ask is what having standard deviation is good for when filtering data. Standard deviation in a filter is a reflection of how much noise there is on the signal being filtered. The higher the standard deviation, the less signal to noise ratio. Knowing standard deviation helps select a filter window size so the signal can be extracted from the noisy data. That can often be calculated without having standard deviation in the finial code. But an adaptive filter might use standard deviation to determine the filter window.

Standard deviation might be used be used to indicate the incoming data is becoming unusable. It might also be used to detect sudden changes in the value being filtered. If the average signal level changes abruptly, standard deviation will rise as the filter transitions from the old value to the new.

Here we have a windowed average of 5 samples. The signal change abruptly to 100, and then back to 0. The windowed average takes 5 samples to fully respond to this change. During that time, the standard deviation rises, and falls. In this way it can act as an indicator of such a change before the average reacts.

There are other uses for knowing standard deviation. By considering what standard deviation means about the data being acquired, addition uses can be conceived.

The classic equation for calculating standard deviation is as follows:

Here *x* is an array of data points, *N* the number of data points, *s* is the corrected sample standard deviation—what most people call simply call standard deviation, and *x* is the average of *x*. There is something one will notice about this calculation: the average must be known before standard deviation can be computed.

I've shown in previous articles a simple windowed average template class that can keep the average without having to do a complete calculation over the entire array. For standard deviation I didn't think this was possible because the average needed to be known before calculating standard deviation. Turns out I was mistaken. There is an iterative method. First, the iterative version of an average:

Then standard deviation sum becomes:

And standard deviation at any point *n* is:

That is nice for doing standard deviation in a single pass, but doesn't help us much with a continuous windowed standard deviation. In order to do this windowed, we need to be able to somehow subtract off older values from the running sum. This isn't possible here because of the divide by *n*-1 at the end.

Failing to find a solution to this, I went in search of an other method for doing windowed standard deviation and came across this article. The article presents this form for standard deviation:

This can be decomposed into two running sums:

Now that we have two running sums, we can make those circular.

Where *w* is the size of the window. We assume that *a*_{n} = 0 for all *n* < 0. That is, the circular buffers start filled with zero. This has an other advantage in that the average is kept along side the standard deviation.

Here is a template class that does both a windowed average and standard deviation. It is derived from my earlier work that just did the average. As before, the template can do any type of input as long as it supports basic arithmetic functions, and has an overloaded version of the square root function *sqrt*.

//==============================================================================

// Template class for keeping a windowed average and standard deviation.

// Assumes math operations on 'TYPE' are computationally expensive, so used

// sparingly.

// For standard deviation to work, sqrt must be specialized to handle the

// template type.

//==============================================================================

template< class TYPE >

class AverageStandardDeviation

{

protected:

// Storage of all the data points that make the average.

TYPE * array;

// A representation of zero in the type of this template.

TYPE const zero;

// Index into 'array' to save next data point.

unsigned currentIndex;

// How many data points are currently in array.

unsigned top;

// Length of 'array'.

unsigned length;

// Running sum used to compute average.

TYPE sum;

// Running sum of squares used for standard deviation.

TYPE squareSum;

public:

//------------------------------------------------------------------------

// Create moving average with 'lengthParameter' number of samples to be

// averaged together.

//------------------------------------------------------------------------

AverageStandardDeviation( unsigned lengthParameter )

:

zero( static_cast< TYPE >( 0 ) )

{

length = lengthParameter;

// Allocate storage for

array = new TYPE[ length ];

// We need a definition of zero, so make a cast to get this.

// NOTE: This is important if TYPE is truly a class.

reset();

}

//------------------------------------------------------------------------

// Reset all data in average.

// NOTE: Average result will be the equivalent of zero.

//------------------------------------------------------------------------

void reset()

{

// Empty array.

for ( currentIndex = 0; currentIndex < length; ++currentIndex )

array[ currentIndex ] = zero;

currentIndex = 0;

top = 0;

sum = zero;

squareSum = zero;

}

//------------------------------------------------------------------------

// Resize the number of samples used in average.

// If the size becomes larger, then the average will not change until

// more data is added.

// If the size becomes smaller, then the last 'newLength' samples will be

// copied to new average.

//------------------------------------------------------------------------

void resize( unsigned newLength )

{

// Allocate memory for new array.

TYPE * newArray = new TYPE[ newLength ];

// Reset sums.

sum = zero;

squareSum = zero;

for ( unsigned count = newLength; count > 0; --count )

{

unsigned index = count - 1;

// Do we have data to copy from the old array at this index?

if ( index < top )

{

// Get previous data index.

if ( currentIndex == 0 )

currentIndex = length - 1;

else

currentIndex -= 1;

// Save data in new array.

newArray[ index ] = array[ currentIndex ];

// Sum this data.

sum += newArray[ index ];

squareSum += newArray[ index ] * newArray[ index ];

}

else

// Set this position to zero.

newArray[ index ] = zero;

}

// Delete old array holder and switch to new array.

delete array;

array = newArray;

// Move top if the new length is lower than previous top.

if ( top > newLength )

top = newLength;

// New length.

length = newLength;

// Set new index location.

currentIndex = top;

if ( currentIndex >= length )

currentIndex = 0;

}

//------------------------------------------------------------------------

// Return the computed average.

// NOTE: Uses the data in the average thus far. If average isn't full,

// then only those data points in the array are used.

//------------------------------------------------------------------------

TYPE getAverage() const

{

TYPE result;

if ( top > 0 )

result = sum / static_cast< TYPE >( top );

else

result = zero;

return result;

}

//------------------------------------------------------------------------

// Return the computed standard deviation.

//------------------------------------------------------------------------

TYPE getStandardDeviation() const

{

TYPE result;

if ( top > 1 )

{

using std::sqrt;

result = sum * sum;

result /= static_cast< TYPE >( top );

result = squareSum - result;

result /= static_cast< TYPE >( top - 1 );

result = sqrt( result );

}

else

result = zero;

return result;

}

//------------------------------------------------------------------------

// Return the number of data points in the moving average so far.

//------------------------------------------------------------------------

unsigned getNumberOfDataPoints() const

{

return top;

}

//------------------------------------------------------------------------

// Add new data to average.

//------------------------------------------------------------------------

void push( TYPE newData )

{

TYPE oldPoint = array[ currentIndex ];

// Remove old data from average sum.

sum -= oldPoint;

// Remove old data from square sum.

squareSum -= oldPoint * oldPoint;

// Add new data to average sum.

sum += newData;

// Add new data to average square sum.

squareSum += newData * newData;

// Save new data.

array[ currentIndex ] = newData;

// Advance index.

++currentIndex;

// If this marks the top...

if ( currentIndex > top )

// Update the top.

top = currentIndex;

// Wrap-around index.

if ( currentIndex >= length )

currentIndex = 0;

}

};

What's nice about this code is that adding a new sample to the average only requires a few operations. An array of the input data does need to be kept to it can be subtracted off, but both the average and standard deviation are kept as sums that can be quickly called upon for an exact result at any point. Standard deviation requires a few operations to translate the sums into a result, but this system is much faster than having to run through all the data points every time in order to get a result.