trunk/benchmarks/triang.cpp

/*
    Triangle wave generator benchmark

    This is a benchmark for comparison between a integer math, table lookup
    and numeric sine wave harmonics solution.

    Copyright (C) 2005 Christian Schoenebeck <cuse@users.sf.net>
*/

#include <math.h>
#include <time.h>
#include <stdio.h>

// whether we should not show any messages on the console
#ifndef SILENT
# define SILENT 0
#endif

// set to 1 if you want to output the three calculated waves as RAW files
// you can e.g. open it as RAW file in Rezound
// (32 bit SP-FP PCM, mono, little endian, 44100kHz)
#ifndef OUTPUT_AS_RAW_WAVE
# define OUTPUT_AS_RAW_WAVE     0
#endif

// how many sample points should we calculate in one sequence
#ifndef STEPS
# define STEPS                  16384
#endif

// how often should we repeat the benchmark loop of each solution
#ifndef RUNS
# define RUNS                   1000
#endif

// whether the wave should have positive and negative range (signed -> 1) or just positive (unsigned -> 0)
#ifndef SIGNED
# define SIGNED                 1
#endif

// pro forma
#ifndef SAMPLING_RATE
# define SAMPLING_RATE          44100.0f
#endif

// return value of this benchmark
// to indicate the best performing solution
#define INT_MATH_SOLUTION       2  /* we don't start with 1, as this is reserved for unknown errors */
#define DI_HARMONIC_SOLUTION    3
#define TABLE_LOOKUP_SOLUTION   4  /* table lookup solution is currently disabled in this benchmark, see below */
#define INVALID_RESULT          -1

// we use 32 bit single precision floating point as sample point format
typedef float sample_t;

// integer math solution
float int_math(sample_t* pDestinationBuffer, float* pAmp, const int steps, const float frequency) {
    // pro forma
    const float r = frequency / SAMPLING_RATE; // frequency alteration quotient
    unsigned int maxvalue = (unsigned int) -1; // all 0xFFFF...
    #if SIGNED
    const float normalizer = 1.0f / ((float) maxvalue / 4.0f);
    #else // unsigned
    const float normalizer = 1.0f / ((float) maxvalue / 2.0f);
    #endif
    const int c = (int) (maxvalue * r);
    const int signshifts = (sizeof(int) * 8) - 1;

    int iSign;

    clock_t stop_time;
    clock_t start_time = clock();

    for (int run = 0; run < RUNS; run++) {
        int iLevel = 0;
        for (int i = 0; i < steps; ++i) {
            iLevel += c;
            iSign  = (iLevel >> signshifts) | 1;
            #if SIGNED
            pDestinationBuffer[i] = (normalizer * (float) (iSign * iLevel) - 1.0f) * pAmp[i]; // * pAmp[i] just to simulate some memory load
            #else // unsigned
            pDestinationBuffer[i] = normalizer * (float) (iSign * iLevel) * pAmp[i]; // * pAmp[i] just to simulate some memory load
            #endif
        }
    }

    stop_time = clock();
    float elapsed_time = (stop_time - start_time) / (double(CLOCKS_PER_SEC) / 1000.0);
    #if ! SILENT
    printf("int math solution elapsed time: %1.0f ms\n", elapsed_time);
    #endif

    return elapsed_time;
}

// table lookup solution (currently disabled)
//
// This solution is not yet implemented in LinuxSampler yet and probably
// never will, I simply haven't found an architectures / system where this
// turned out to be the best solution and it introduces too many problems
// anyway. If you found an architecture where this seems to be the best
// solution, please let us know!
#if 0
float table_lookup(sample_t* pDestinationBuffer, float* pAmp, const int steps, const float frequency) {
    // pro forma
    const float r = frequency / SAMPLING_RATE; // frequency alteration quotient
    #if SIGNED
    float c = r * 4.0f;
    #else
    float c = r * 2.0f;
    #endif
    const int wl = (int) (SAMPLING_RATE / frequency); // wave length (in sample points)

    // 'volatile' to avoid the cache to fudge the benchmark result
    volatile float* pPrerenderingBuffer = new float[wl];

    // prerendering of the triangular wave
    {
        float level = 0.0f;
        for (int i = 0; i < wl; ++i) {
            level += c;
            #if SIGNED
            if (level >= 1.0f) {
                c = -c;
                level = 1.0f;
            }
            else if (level <= -1.0f) {
                c = -c;
                level = -1.0f;
            }
            #else
            if (level >= 1.0f) {
                c = -c;
                level = 1.0f;
            }
            else if (level <= 0.0f) {
                c = -c;
                level = 0.0f;
            }
            #endif
            pPrerenderingBuffer[i] = level;
        }
    }

    clock_t stop_time;
    clock_t start_time = clock();

    for (int run = 0; run < RUNS; run++) {
        for (int i = 0; i < steps; ++i) {
            pDestinationBuffer[i] = pPrerenderingBuffer[i % wl] * pAmp[i]; // * pAmp[i] just to simulate some memory load
        }
    }

    stop_time = clock();
    float elapsed_time = (stop_time - start_time) / (double(CLOCKS_PER_SEC) / 1000.0);
    #if ! SILENT
    printf("Table lookup solution elapsed time: %1.0f ms\n", elapsed_time);
    #endif

    if (pPrerenderingBuffer) delete[] pPrerenderingBuffer;

    return elapsed_time;
}
#endif

// numeric, di-harmonic solution
float numeric_di_harmonic_solution(sample_t* pDestinationBuffer, float* pAmp, const int steps, const float frequency) {
    // we approximate the triangular wave by adding 2 harmonics
    const float c1   = 2.0f * M_PI * frequency / SAMPLING_RATE;
    const float phi1 = 0.0f; // phase displacement
    const float c2   = 2.0f * M_PI * frequency / SAMPLING_RATE * 3.0f;
    const float phi2 = 0.0f; // phase displacement

    // amplitue for the 2nd harmonic (to approximate the triangular wave)
    const float amp2 = 0.1f;

    // initial values for real and imaginary part
    float real1 = cos(phi1);
    float imag1 = sin(phi1);
    float real2 = cos(phi2);
    float imag2 = sin(phi2);

    clock_t stop_time;
    clock_t start_time = clock();

    for (int run = 0; run < RUNS; run++) {
        for (int i = 0; i < steps; i++) {
            #if SIGNED
            pDestinationBuffer[i] = (real1 + real2 * amp2) * pAmp[i]; // * pAmp[i] just to simulate some memory load
            #else // unsigned
            pDestinationBuffer[i] = ((real1 + real2 * amp2) * 0.5f + 0.5f) * pAmp[i]; // * pAmp[i] just to simulate some memory load
            #endif
            real1 -= c1 * imag1;
            imag1 += c1 * real1;
            real2 -= c2 * imag2;
            imag2 += c2 * real2;
        }
    }

    stop_time = clock();
    float elapsed_time = (stop_time - start_time) / (double(CLOCKS_PER_SEC) / 1000.0);
    #if ! SILENT
    printf("Numeric harmonics solution elapsed time: %1.0f ms\n", elapsed_time);
    #endif

    return elapsed_time;
}

// output calculated values as RAW audio format (32 bit floating point, mono) file
void output_as_raw_file(const char* filename, sample_t* pOutputBuffer, int steps) {
    FILE* file = fopen(filename, "w");
    if (file) {
        fwrite((void*) pOutputBuffer, sizeof(float), steps, file);
        fclose(file);
    } else {
        fprintf(stderr, "Could not open %s\n", filename);
    }
}

int main() {
    const int steps = STEPS;
    const int sinusoidFrequency = 100; // Hz

    #if ! SILENT
    # if SIGNED
    printf("Signed triangular wave benchmark\n");
    printf("--------------------------------\n");
    # else
    printf("Unsigned triangular wave benchmark\n");
    printf("----------------------------------\n");
    # endif
    #endif

    // output buffer for the calculated sinusoid wave
    sample_t* pOutputBuffer = new sample_t[steps];
    // just an arbitrary amplitude envelope to simulate a bit higher memory bandwidth
    float* pAmplitude  = new float[steps];

    // pro forma - an arbitary amplitude envelope
    for (int i = 0; i < steps; ++i)
        pAmplitude[i] = (float) i / (float) steps;

    // how long each solution took (in seconds)
    float int_math_result, /*table_lookup_result,*/ numeric_di_harmonic_result;

    int_math_result = int_math(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
    #if OUTPUT_AS_RAW_WAVE
      output_as_raw_file("bench_int_math.raw", pOutputBuffer, steps);
    #endif
    //table_lookup_result = table_lookup(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
    //#if OUTPUT_AS_RAW_WAVE
    //  output_as_raw_file("bench_table_lookup.raw", pOutputBuffer, steps);
    //#endif
    numeric_di_harmonic_result = numeric_di_harmonic_solution(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
    #if OUTPUT_AS_RAW_WAVE
      output_as_raw_file("bench_numeric_harmonics.raw", pOutputBuffer, steps);
    #endif

    #if ! SILENT
    printf("\nOK, 2nd try\n\n");
    #endif

    int_math_result            += int_math(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
    //table_lookup_result        += table_lookup(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
    numeric_di_harmonic_result += numeric_di_harmonic_solution(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);

    if (pAmplitude)    delete[] pAmplitude;
    if (pOutputBuffer) delete[] pOutputBuffer;

    if (/*int_math_result <= table_lookup_result &&*/ int_math_result <= numeric_di_harmonic_result) return INT_MATH_SOLUTION;
    if (/*numeric_di_harmonic_result <= table_lookup_result &&*/ numeric_di_harmonic_result <= int_math_result) return DI_HARMONIC_SOLUTION;
    //if (table_lookup_result <= int_math_result && table_lookup_result <= numeric_di_harmonic_result) return TABLE_LOOKUP_SOLUTION;

    return INVALID_RESULT; // error
}
1	schoenebeck	714	/*
2			Triangle wave generator benchmark
3
4			This is a benchmark for comparison between a integer math, table lookup
5			and numeric sine wave harmonics solution.
6
7			Copyright (C) 2005 Christian Schoenebeck <cuse@users.sf.net>
8			*/
9
10			#include <math.h>
11			#include <time.h>
12			#include <stdio.h>
13
14			// whether we should not show any messages on the console
15			#ifndef SILENT
16			# define SILENT 0
17			#endif
18
19			// set to 1 if you want to output the three calculated waves as RAW files
20			// you can e.g. open it as RAW file in Rezound
21			// (32 bit SP-FP PCM, mono, little endian, 44100kHz)
22			#ifndef OUTPUT_AS_RAW_WAVE
23			# define OUTPUT_AS_RAW_WAVE 0
24			#endif
25
26			// how many sample points should we calculate in one sequence
27			#ifndef STEPS
28			# define STEPS 16384
29			#endif
30
31			// how often should we repeat the benchmark loop of each solution
32			#ifndef RUNS
33			# define RUNS 1000
34			#endif
35
36			// whether the wave should have positive and negative range (signed -> 1) or just positive (unsigned -> 0)
37			#ifndef SIGNED
38			# define SIGNED 1
39			#endif
40
41			// pro forma
42			#ifndef SAMPLING_RATE
43			# define SAMPLING_RATE 44100.0f
44			#endif
45
46			// return value of this benchmark
47			// to indicate the best performing solution
48			#define INT_MATH_SOLUTION 2 /* we don't start with 1, as this is reserved for unknown errors */
49			#define DI_HARMONIC_SOLUTION 3
50			#define TABLE_LOOKUP_SOLUTION 4 /* table lookup solution is currently disabled in this benchmark, see below */
51			#define INVALID_RESULT -1
52
53			// we use 32 bit single precision floating point as sample point format
54			typedef float sample_t;
55
56			// integer math solution
57			float int_math(sample_t* pDestinationBuffer, float* pAmp, const int steps, const float frequency) {
58			// pro forma
59			const float r = frequency / SAMPLING_RATE; // frequency alteration quotient
60			unsigned int maxvalue = (unsigned int) -1; // all 0xFFFF...
61			#if SIGNED
62			const float normalizer = 1.0f / ((float) maxvalue / 4.0f);
63			#else // unsigned
64			const float normalizer = 1.0f / ((float) maxvalue / 2.0f);
65			#endif
66			const int c = (int) (maxvalue * r);
67			const int signshifts = (sizeof(int) * 8) - 1;
68
69			int iSign;
70
71			clock_t stop_time;
72			clock_t start_time = clock();
73
74			for (int run = 0; run < RUNS; run++) {
75			int iLevel = 0;
76			for (int i = 0; i < steps; ++i) {
77			iLevel += c;
78			iSign = (iLevel >> signshifts) \| 1;
79			#if SIGNED
80			pDestinationBuffer[i] = (normalizer * (float) (iSign * iLevel) - 1.0f) * pAmp[i]; // * pAmp[i] just to simulate some memory load
81			#else // unsigned
82			pDestinationBuffer[i] = normalizer * (float) (iSign * iLevel) * pAmp[i]; // * pAmp[i] just to simulate some memory load
83			#endif
84			}
85			}
86
87			stop_time = clock();
88			float elapsed_time = (stop_time - start_time) / (double(CLOCKS_PER_SEC) / 1000.0);
89			#if ! SILENT
90			printf("int math solution elapsed time: %1.0f ms\n", elapsed_time);
91			#endif
92
93			return elapsed_time;
94			}
95
96			// table lookup solution (currently disabled)
97			//
98			// This solution is not yet implemented in LinuxSampler yet and probably
99			// never will, I simply haven't found an architectures / system where this
100			// turned out to be the best solution and it introduces too many problems
101			// anyway. If you found an architecture where this seems to be the best
102			// solution, please let us know!
103			#if 0
104			float table_lookup(sample_t* pDestinationBuffer, float* pAmp, const int steps, const float frequency) {
105			// pro forma
106			const float r = frequency / SAMPLING_RATE; // frequency alteration quotient
107			#if SIGNED
108			float c = r * 4.0f;
109			#else
110			float c = r * 2.0f;
111			#endif
112			const int wl = (int) (SAMPLING_RATE / frequency); // wave length (in sample points)
113
114			// 'volatile' to avoid the cache to fudge the benchmark result
115			volatile float* pPrerenderingBuffer = new float[wl];
116
117			// prerendering of the triangular wave
118			{
119			float level = 0.0f;
120			for (int i = 0; i < wl; ++i) {
121			level += c;
122			#if SIGNED
123			if (level >= 1.0f) {
124			c = -c;
125			level = 1.0f;
126			}
127			else if (level <= -1.0f) {
128			c = -c;
129			level = -1.0f;
130			}
131			#else
132			if (level >= 1.0f) {
133			c = -c;
134			level = 1.0f;
135			}
136			else if (level <= 0.0f) {
137			c = -c;
138			level = 0.0f;
139			}
140			#endif
141			pPrerenderingBuffer[i] = level;
142			}
143			}
144
145			clock_t stop_time;
146			clock_t start_time = clock();
147
148			for (int run = 0; run < RUNS; run++) {
149			for (int i = 0; i < steps; ++i) {
150			pDestinationBuffer[i] = pPrerenderingBuffer[i % wl] * pAmp[i]; // * pAmp[i] just to simulate some memory load
151			}
152			}
153
154			stop_time = clock();
155			float elapsed_time = (stop_time - start_time) / (double(CLOCKS_PER_SEC) / 1000.0);
156			#if ! SILENT
157			printf("Table lookup solution elapsed time: %1.0f ms\n", elapsed_time);
158			#endif
159
160			if (pPrerenderingBuffer) delete[] pPrerenderingBuffer;
161
162			return elapsed_time;
163			}
164			#endif
165
166			// numeric, di-harmonic solution
167			float numeric_di_harmonic_solution(sample_t* pDestinationBuffer, float* pAmp, const int steps, const float frequency) {
168			// we approximate the triangular wave by adding 2 harmonics
169			const float c1 = 2.0f * M_PI * frequency / SAMPLING_RATE;
170			const float phi1 = 0.0f; // phase displacement
171			const float c2 = 2.0f * M_PI * frequency / SAMPLING_RATE * 3.0f;
172			const float phi2 = 0.0f; // phase displacement
173
174			// amplitue for the 2nd harmonic (to approximate the triangular wave)
175			const float amp2 = 0.1f;
176
177			// initial values for real and imaginary part
178			float real1 = cos(phi1);
179			float imag1 = sin(phi1);
180			float real2 = cos(phi2);
181			float imag2 = sin(phi2);
182
183			clock_t stop_time;
184			clock_t start_time = clock();
185
186			for (int run = 0; run < RUNS; run++) {
187			for (int i = 0; i < steps; i++) {
188			#if SIGNED
189			pDestinationBuffer[i] = (real1 + real2 * amp2) * pAmp[i]; // * pAmp[i] just to simulate some memory load
190			#else // unsigned
191			pDestinationBuffer[i] = ((real1 + real2 * amp2) * 0.5f + 0.5f) * pAmp[i]; // * pAmp[i] just to simulate some memory load
192			#endif
193			real1 -= c1 * imag1;
194			imag1 += c1 * real1;
195			real2 -= c2 * imag2;
196			imag2 += c2 * real2;
197			}
198			}
199
200			stop_time = clock();
201			float elapsed_time = (stop_time - start_time) / (double(CLOCKS_PER_SEC) / 1000.0);
202			#if ! SILENT
203			printf("Numeric harmonics solution elapsed time: %1.0f ms\n", elapsed_time);
204			#endif
205
206			return elapsed_time;
207			}
208
209			// output calculated values as RAW audio format (32 bit floating point, mono) file
210			void output_as_raw_file(const char* filename, sample_t* pOutputBuffer, int steps) {
211			FILE* file = fopen(filename, "w");
212			if (file) {
213			fwrite((void*) pOutputBuffer, sizeof(float), steps, file);
214			fclose(file);
215			} else {
216			fprintf(stderr, "Could not open %s\n", filename);
217			}
218			}
219
220			int main() {
221			const int steps = STEPS;
222			const int sinusoidFrequency = 100; // Hz
223
224			#if ! SILENT
225			# if SIGNED
226			printf("Signed triangular wave benchmark\n");
227			printf("--------------------------------\n");
228			# else
229			printf("Unsigned triangular wave benchmark\n");
230			printf("----------------------------------\n");
231			# endif
232			#endif
233
234			// output buffer for the calculated sinusoid wave
235			sample_t* pOutputBuffer = new sample_t[steps];
236			// just an arbitrary amplitude envelope to simulate a bit higher memory bandwidth
237			float* pAmplitude = new float[steps];
238
239			// pro forma - an arbitary amplitude envelope
240			for (int i = 0; i < steps; ++i)
241			pAmplitude[i] = (float) i / (float) steps;
242
243			// how long each solution took (in seconds)
244			float int_math_result, /table_lookup_result,/ numeric_di_harmonic_result;
245
246			int_math_result = int_math(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
247			#if OUTPUT_AS_RAW_WAVE
248			output_as_raw_file("bench_int_math.raw", pOutputBuffer, steps);
249			#endif
250			//table_lookup_result = table_lookup(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
251			//#if OUTPUT_AS_RAW_WAVE
252			// output_as_raw_file("bench_table_lookup.raw", pOutputBuffer, steps);
253			//#endif
254			numeric_di_harmonic_result = numeric_di_harmonic_solution(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
255			#if OUTPUT_AS_RAW_WAVE
256			output_as_raw_file("bench_numeric_harmonics.raw", pOutputBuffer, steps);
257			#endif
258
259			#if ! SILENT
260			printf("\nOK, 2nd try\n\n");
261			#endif
262
263			int_math_result += int_math(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
264			//table_lookup_result += table_lookup(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
265			numeric_di_harmonic_result += numeric_di_harmonic_solution(pOutputBuffer, pAmplitude, steps, sinusoidFrequency);
266
267			if (pAmplitude) delete[] pAmplitude;
268			if (pOutputBuffer) delete[] pOutputBuffer;
269
270			if (/int_math_result <= table_lookup_result &&/ int_math_result <= numeric_di_harmonic_result) return INT_MATH_SOLUTION;
271			if (/numeric_di_harmonic_result <= table_lookup_result &&/ numeric_di_harmonic_result <= int_math_result) return DI_HARMONIC_SOLUTION;
272			//if (table_lookup_result <= int_math_result && table_lookup_result <= numeric_di_harmonic_result) return TABLE_LOOKUP_SOLUTION;
273
274			return INVALID_RESULT; // error
275			}