In previous articles I wrote about the advantages of using a moving average over a windowed average, and how a continuous windowed standard deviation can be calculated. Keeping a continuous windowed average can be costly to memory. So the question is, can one keep a standard deviation when running the moving average algorithm?
The answer looks to be yes—mostly. The results of this standard deviation are close to the actual windowed standard deviation, but like the moving average filter itself, it is not an exact. Let us first explore the algorithm I implemented.
an = an-1 + (xn – an-1) / w
bn = bn-1 + (xn – an-1)(xn – an) – bn-1/w
sn= sqrt[ bn / ( w - 1 ) ]
Where a is is the moving average, b reflects the sum of variance, w is the window size, and s the standard deviation at anytime.
I wrote a simple program to generate Gaussian distributed random data with an initial standard deviation and ramping to a final standard deviation. Because the data is random, the actual standard deviation at any point over any number of samples can vary. However with enough samples and a windowed standard deviation there should be a clear trend of increasing standard deviation. I wanted to compare this estimated standard deviation method I created with windowed standard deviation to see how closely they follow one an other. To do this, I generated 1000 points of random data centered at an average of 0, and linearly increasing the standard deviation from 0 to 100.
In this plot we can see the moving standard deviation over a window of 10 samples vs the actual standard deviation over 10 samples. Included are the linearly regression lines for both data sets, showing a steady increase from a standard deviation of 0 to around 100. While the two methods are not exact, they tend to follow one an other. The true windowed standard deviation seems to be a little more erratic than the moving standard deviation. So the estimated seems to be dampened version of windowed standard deviation, and may prove useful to applications that require such things. This damping effect can be shown in a closer look at a set of data.
This plot shows 256 standard deviation periods on a set of data. Notice that the moving standard deviation is on average closer to the overall standard deviation, hovering around 100. The windowed standard deviation for any window along the sample set can vary more wildly. However, the results are similar enough that the moving values should be a decent reflection of actual standard deviation.
This plot shows a rapid change in data, starting at 400, jumping to 900, then back to 400. The blue trace shows the incoming data, and the green the filtered data. In yellow is the true windowed standard deviation for the same time period as the filter. The orange shows the estimated standard deviation. While the standard deviation can ramp up quickly to meet the true deviation, it takes much longer to decay. What this says is the algorithm is quick to respond to increases in standard deviation, but slow for decreases.
Why this algorithm exhibits these characteristics is due to what the algorithm actually is: a squared error accumulator with a moving average drain. Error is accumulated just like normal standard deviation. However, there isn't a way to subtract off the exact amount of error at a latter time like that of a true windowed standard deviation. Instead, old error is removed with a low-pass filter just like that used in a moving average.
Software implementation requires two running values: an average, and a variance accumulator.
uint32_t varianceDrain = ( varianceSum / count );
int32_t lastAverage = average;
int32_t newAverage = average;
newAverage += ( value - lastAverage ) / count;
average = newAverage;
varianceSum += ( value - lastAverage ) * ( value - newAverage );
varianceSum -= varianceDrain;
In this snip of code, the two running variables are average, and varianceSum. They are shown being updated with a new sample called value. The variable windowsSize might be a constant, and for speed can be converted into a power-of-2 so the divides can turn into shifts. The variable average can be used directly. Standard deviation requires one last conversion:
standardDeviation = varianceSum;
standardDeviation /= ( count - 1 );
standardDeviation = sqrt( standardDeviation );
The sqrt function is an integer square root, which takes up to the word bit size shift and addition operations to complete. They are much faster than the floating-point alternative. The square root can be eliminated if one can live with variance rather than standard deviation as a result.
To start the process, using a count variable rather than windowSize may be desired. Or the first sample could seed the average like this:
average = value * windowSize;
For implementation one must consider the worst-case values to avoid issues of overflow. To do this one must know the largest data value (xmax), the highest standard deviation (smax), and the word size (dmax). Note that if doing a saturated summation, only the largest meaningful values need to be considered. If standard deviation is being used as a threshold, then the maximum standard deviation can be the threshold as values above this are not necessary.
To avoid overflowing the average accumulator, insure w·xmax < dmax. For 10-bit data, and a window size of 64 (6-bits), 16-bits are needed for data accumulation which is the largest data type needed. Also be mindful of sine/unsigned. For a 16-bit unsigned integer dmax = 65535, but for a 16-bit signed dmax = 32767.
Variance is a square of error and grows more quickly than does the average. However one does not need to consider the data max (xmax)—only the highest standard deviation (smax). The restriction is that 2·w·smax < word size. So if the worst case standard deviation is 128 (7-bits) and the window size is 4096 (10-bits), the variance sum needs only to be 24-bits. Detecting variance saturation can be useful as well, which will indicate that the standard deviation is at or above it's maximum. If using the standard deviation for a threshold, this might actually be preferred because it means the accumulator only needs to cover up to the threshold but not actually account for anything above. On small microprocessors that allows smaller word sizes to be used, which can greatly improve speed.
//=============================================================================
// Name: movingAverage32.h
// Uses: Interface to compute a moving average and moving standard deviation.
// Date: 2015-08-02
// Author: Andrew Que (http://www.DrQue.net/)
// Revisions:
// 2015-08-02 - QUE - Creation.
//
// Copyright (c) 2015 Andrew Que (http://www.DrQue.net/)
// The MIT License (MIT). See details at end.
//=============================================================================
#ifndef _MOVINGAVERAGE32_H_
#define _MOVINGAVERAGE32_H_
#include <stdint.h>
typedef struct Average32
{
int count; // Number of samples thus far (up to 'windowSize').
int windowSize; // Number of samples to average over.
int32_t average; // Current average (signed or unsigned).
uint32_t varianceSum; // Running variance sum (always unsigned).
} Average32;
//-----------------------------------------------------------------------------
// USES:
// Take the square root of a 32-bit integer.
// INPUT:
// input - Value for which the root is to be taken.
// OUTPUT:
// Integer square root of 'input'.
// NOTE:
// The result of a taking the square root of a 32-bit number always fits
// into a 16-bit number because sqrt( 2^32 ) = 2^(32/2) = 2^16.
//-----------------------------------------------------------------------------
static inline uint16_t sqrt_u32( uint32_t input )
{
uint32_t result = 0;
// The second-to-top bit is set.
uint32_t bit = UINT32_C( 1 ) << ( sizeof( uint32_t) * 8 - 2 );
// "bit" starts at the highest power of four <= the argument.
while ( bit > input )
bit >>= 2;
while ( bit != 0 )
{
if ( input >= result + bit )
{
input -= result + bit;
result = ( result >> 1 ) + bit;
}
else
result >>= 1;
bit >>= 2;
}
return (uint16_t)result;
}
//-----------------------------------------------------------------------------
// USES:
// Setup 32-bit moving average.
// INPUT:
// average - 32-bit moving average instance.
// windowsSize - Number of samples to average over.
// shift - Shift to scale input to add precision. Optional.
// OUTPUT:
// Modified 'average'.
// Nothing is returned.
//-----------------------------------------------------------------------------
static inline void average32_Init
(
Average32 * average,
int32_t windowSize
)
{
average->windowSize = windowSize;
average->count = 0;
average->average = 0;
average->varianceSum = 0;
}
//-----------------------------------------------------------------------------
// USES:
// Add a data point to the moving average.
// INPUT:
// average - 32-bit moving average instance.
// value - New value to incorporate.
// OUTPUT:
// Modified 'average'.
// Nothing is returned.
//-----------------------------------------------------------------------------
static inline void average32_AddValue( Average32 * average, int32_t value )
{
int maxCount = ( 1 << average->windowSize );
int32_t lastAverage = average->average;
int32_t newAverage = average->average;
// Roll-off old data.
newAverage -= ( lastAverage >> average->windowSize );
// Accumulate new data.
newAverage += value;
// Use the very first sample as a representation of all data to come.
// (Optional--only needed if desiring accurate results in start up.)
if ( 0 == average->count )
{
// Make the first value the starting average.
newAverage <<= average->windowSize;
// Set last is the same as first to force variance at 0 for first sample.
lastAverage = newAverage;
}
else
// A change to the most-significant bit means the accumulator has overflowed.
// (Optional if average saturation is needed.)
if ( ( newAverage ^ average->average ) >> ( sizeof( newAverage ) * 8 - 1 ) )
{
// Saturate negative or positive?
if ( average->average < 0 )
newAverage = INT32_MIN;
else
newAverage = INT32_MAX;
}
// Save new average.
average->average = newAverage;
// Calculate variance the new sample adds.
uint32_t variance =
( value - ( newAverage >> average->windowSize ) )
* ( value - ( lastAverage >> average->windowSize ) );
// Roll-off old variance.
// (Optional if. Can be removed if always waiting for at least
// average->windowSize samples before calculating result.)
if ( average->count >= maxCount )
average->varianceSum -= ( average->varianceSum >> average->windowSize );
// Accumulate new variance.
variance += average->varianceSum;
// Unsigned saturate test and correction.
// (Optional if saturation is needed.)
variance |= -( variance < average->varianceSum );
average->varianceSum = variance;
// Keep incrementing count if just getting started.
// (Optional--only needed if desiring accurate results in start up.)
if ( average->count < maxCount )
average->count += 1;
}
//-----------------------------------------------------------------------------
// USES:
// Return moving average of data up to this point.
// INPUT:
// average - 32-bit moving average instance.
// OUTPUT:
// The moving average of all data.
//-----------------------------------------------------------------------------
static inline int32_t average32_Get( Average32 const * average )
{
return average->average >> average->windowSize;
}
//-----------------------------------------------------------------------------
// USES:
// Return approximate standard deviation of data up to this point.
// INPUT:
// average - 32-bit moving average instance.
// OUTPUT:
// An approximation of standard deviation of data.
//-----------------------------------------------------------------------------
static inline uint32_t average32_GetStandardDeviation
(
Average32 const * average
)
{
uint32_t standardDeviation = average->varianceSum;
// Avoid a divide-by-zero.
// (Can be removed if using shifts instead.)
if ( average->count > 1 )
{
// Division can be replaced with shift, but the result will be biased.
// Since the standard deviation is an estimate anyway, that shouldn't be
// a problem.
standardDeviation /= ( average->count - 1 );
//standardDeviation >>= average->windowSize; // Divide by shift.
// The square root can be skipped, which results in variance rather than
// standard deviation.
standardDeviation = sqrt_u32( standardDeviation );
}
return standardDeviation;
}
#endif // _MOVINGAVERAGE32_H_
// The MIT License (MIT)
//
// Copyright (c) 2015 Andrew Que (http://www.DrQue.net/)
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
//
The above is a C implementation for 32-bit signed data with most features written in. It can be cut down depending on the need of the application. For instance, to make this system faster biased variance can be used rather than true standard deviation. This turns the divide into a shift, and removes the square root function. These changes can be made with ease:
-
From signed to unsigned. This allows the average saturation code to be simplified to a single line just like it is for the unsigned saturation of the variance sum.
-
No accumulator saturation. If you don't need it, this will save a couple instructions.
-
No initial condition setting. Requires average to have at least a window size of samples before results are accurate.
-
Biased standard deviation. Change the divide to a shift when calculating standard deviation and increase the speed. Almost always a good idea.
-
Use variance instead of standard deviation. Just remove the square root and save that amount of time.
The code can easily be changed to other bit widths as well. To change to 16-bit code, replace all occurrences of 16 with 8, and all 32 with 16. The switch of 16 to 8 is for the result of the square root function which is moot because it's promoted when used anyway. In the same way, this code can also be scaled up to 64-bits.
