logspeed.c

// test signal *power* to dBFS conversion

// gcc -O3 -o logspeed logspeed.c -Wall -lm
// -ffast-math -fno-finite-math-only

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>

// this tool allocates 12 * DATASIZE bytes
#define DATASIZE (10 * 1024*1024) // 120MB

static inline float log_reference (const float v)
{
	return 10.f * log10f (v);
}

static inline float log_fast (float v)
{
	union {float f; int i;} t;
	t.f = v;
	int* const    exp_ptr =  &t.i;
	int           x = *exp_ptr;
	const int     log_2 = ((x >> 23) & 255) - 128;
	x &= ~(255 << 23);
	x += 127 << 23;
	*exp_ptr = x;

	v = ((-1.0f/3) * t.f + 2) * t.f - 2.0f/3.f;
	return (v + log_2) * 3.0188867925f;
}

static inline float log_quick (const float v)
{
	union {float f; int i;} t;
	t.f = v * v * v;
	return -127.f + 1.192092896e-7 * (float)t.i;
}



int main (int argc, char** argv)
{
	unsigned int i;
	struct timeval start;
	struct timeval end;

	float elapsed_ref;
	float elapsed_fast, err_fast;
	float elapsed_quick, err_quick;
	unsigned int inf_fast;
	unsigned int inf_quick;

	float *input = (float*) malloc (DATASIZE * sizeof(float));
	float *out_ref = (float*) malloc (DATASIZE * sizeof(float));
	float *out_sut = (float*) malloc (DATASIZE * sizeof(float));

	srand (time (NULL));
	for (i = 0; i < DATASIZE; ++i) {
		input [i] = 1.f - rand() / (RAND_MAX / 2.f);
	}

	// reference log
	gettimeofday(&start, NULL);
	for (i = 0; i < DATASIZE; ++i) {
		out_ref[i] = log_reference (input[i] * input[i]);
	}
	gettimeofday(&end, NULL);
	elapsed_ref = (end.tv_sec - start.tv_sec) * 1000.f + (end.tv_usec - start.tv_usec) / 1000.f;


	// fast log
	gettimeofday(&start, NULL);
	for (i = 0; i < DATASIZE; ++i) {
		out_sut[i] = log_fast (input[i] * input[i]);
	}
	gettimeofday(&end, NULL);
	elapsed_fast = (end.tv_sec - start.tv_sec) * 1000.f + (end.tv_usec - start.tv_usec) / 1000.f;

	err_fast = 0;
	inf_fast = 0;
	for (i = 0; i < DATASIZE; ++i) {
		if (!isfinite (out_ref[i]) && isfinite (out_sut[i])) {
			++inf_fast;
			continue;
		}
		err_fast += fabs (out_sut[i] - out_ref[i]);
	}
	err_fast /= (float)DATASIZE;

	// quick log
	gettimeofday(&start, NULL);
	for (i = 0; i < DATASIZE; ++i) {
		out_sut[i] = log_quick (input[i] * input[i]);
	}
	gettimeofday(&end, NULL);
	elapsed_quick = (end.tv_sec - start.tv_sec) * 1000.f + (end.tv_usec - start.tv_usec) / 1000.f;

	err_quick = 0;
	inf_quick = 0;
	for (i = 0; i < DATASIZE; ++i) {
		if (!isfinite (out_ref[i]) && isfinite (out_sut[i])) {
			++inf_quick;
			continue;
		}
		err_quick += fabs (out_sut[i] - out_ref[i]);
	}
	err_quick /= (float)DATASIZE;

	printf("Reference: %6.2f ms (100.0 %%)\n", elapsed_ref);
	printf("Fast:      %6.2f ms (%5.2f %%) err: %f inf: %d\n",
			elapsed_fast, 100.f * elapsed_fast / elapsed_ref,
			err_fast, inf_fast);
	printf("Quick:     %6.2f ms (%5.2f %%) err: %f inf: %d\n",
			elapsed_quick, 100.f * elapsed_quick / elapsed_ref,
			err_quick, inf_quick);

	free (input);
	free (out_ref);
	free (out_sut);

#if 0 // -60dBFS .. 0dBFS
	// 2> log.dat  --  gnuplot
	// plot[.001:1] 'log.dat' u 1:3, 'log.dat' u 1:4, 10 * log(x)/log(10)
	// plot 'log.dat' u 1:($3 - $2), 'log.dat' u 1:($4 -$2)
	for (i = 1; i < 1000; ++i) {
		float v = (i / 1000.f) * (i / 1000.f);
		fprintf (stderr, "%f %f %f %f\n", v, log_reference(v), log_fast(v), log_quick(v));
	}
#endif
	return 0;
}