[music-dsp] floating point running sum optimization

[koen vos]

Hi,

I implemented Ross' original idea of calculating and adjusting for the
rounding error after each sum_update.

Just combine the rounding errors that occur in adding/subtracting the
newest/oldest sample to/from the sum, with the rounding error that was left
after the previous call to sum_update. Then subtract this combined rounding
error from the running sum, and retain the rounding error that occurs during
this subtraction itself for the next call to sum_update. After each
sum_update, the running sum is accurate to the precision of the format (in
this case, float).

For me this works even when the samples in the window are negative, but
perhaps I'm just lucky that my CPU happens to be in the right rounding
mode..

Make sure your compiler doesn't inline the add(..) and sub(..) functions.
An assembler implementation would obviously be more efficient for a
practical application.

koen.



#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>

/* state and cyclic buffer */
float *buf;
float sum = 0.0f;
float sum_error = 0.0f;
int   index = 0;


float add( float a, float b )
{
	return a + b;
}

float sub( float a, float b )
{
	return a - b;
}

/* add x to sum, and return rounding error */
float update_add( float *sum_ptr, float x )
{
	float sum_old = *sum_ptr;

	/* update sum */
	*sum_ptr= add( *sum_ptr, x );

	/* return rounding error */
	return sub( sub( *sum_ptr, sum_old ), x );
}

/* subtract x from sum, and return rounding error */
float update_sub( float *sum_ptr, float x )
{
	float sum_old = *sum_ptr;

	/* update sum */
	*sum_ptr= sub( *sum_ptr, x );

	/* return rounding error */
	return add( sub( *sum_ptr, sum_old ), x );
}

void update_sum( float x )
{
	/* subtract oldest sample */
	sum_error += update_sub( &sum, buf[index] );

	/* add new sample */
	sum_error += update_add( &sum, x );

	/* try to subtract rounding error */
	sum_error =  update_sub( &sum, sum_error );

	/* store new sample */
	buf[index] = x;

	/* update cyclic buffer index */
	index = (index + 1) & 1023;
}

int main( int argc, char* argv[] )
{
	int k;

	/* allocate and reset memory */
	buf = calloc( 1024, sizeof(float) );

	/* random seed */
	srand( (unsigned)time( NULL ) );

	for( k = 0; k < 1000000; k++ )
		update_sum( (float)( rand() - RAND_MAX/2 ) );

	/* clear sum (accumulated rounding error remains) */
	for( k = 0; k < 1024; k++ )
		update_sum( 0.0f );

	printf( "\n sum: %f \n\n", sum );
}