Chris@366: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@366: 
Chris@366: #include "VectorOpsComplex.h"
Chris@366: 
Chris@366: #include <iostream>
Chris@366: #include <cstdlib>
Chris@366: 
Chris@366: #include <time.h>
Chris@366: 
Chris@366: using namespace std;
Chris@366: 
Chris@366: namespace breakfastquay {
Chris@366: 
Chris@366: namespace Test {
Chris@366: 
Chris@366: #ifdef _WIN32
Chris@366: #define drand48() (-1+2*((float)rand())/RAND_MAX)
Chris@366: #endif
Chris@366: 
Chris@366: bool
Chris@366: testMultiply()
Chris@366: {
Chris@366:     cerr << "testVectorOps: testing v_multiply complex" << endl;
Chris@366: 
Chris@366:     const int N = 1024;
Chris@366:     turbot_complex_sample_t target[N];
Chris@366:     turbot_complex_sample_t src1[N];
Chris@366:     turbot_complex_sample_t src2[N];
Chris@366: 
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	src1[i].re = drand48();
Chris@366: 	src1[i].im = drand48();
Chris@366: 	src2[i].re = drand48();
Chris@366: 	src2[i].im = drand48();
Chris@366:     }
Chris@366: 
Chris@366:     turbot_sample_t mean, first, last, total = 0;
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366:         turbot_complex_sample_t result;
Chris@366:         c_multiply(result, src1[i], src2[i]);
Chris@366: 	if (i == 0) first = result.re;
Chris@366: 	if (i == N-1) last = result.im;
Chris@366: 	total += result.re;
Chris@366: 	total += result.im;
Chris@366:     }
Chris@366:     mean = total / (N*2);
Chris@366:     cerr << "Naive method: mean = " << mean << ", first = " << first
Chris@366: 	 << ", last = " << last << endl;
Chris@366: 
Chris@366:     v_multiply(target, src1, src2, N);
Chris@366:     total = 0;
Chris@366: 
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	if (i == 0) first = target[i].re;
Chris@366: 	if (i == N-1) last = target[i].im;
Chris@366: 	total += target[i].re;
Chris@366: 	total += target[i].im;
Chris@366:     }
Chris@366:     mean = total / (N*2);
Chris@366:     cerr << "v_multiply: mean = " << mean << ", first = " << first
Chris@366: 	 << ", last = " << last << endl;
Chris@366: 
Chris@366:     int iterations = 50000;
Chris@366:     cerr << "Iterations: " << iterations << endl;
Chris@366: 	
Chris@366:     cerr << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << endl;
Chris@366:     float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@366: 
Chris@366:     clock_t start = clock();
Chris@366: 
Chris@366:     for (int j = 0; j < iterations; ++j) {
Chris@366: 	for (int i = 0; i < N; ++i) {
Chris@366:             c_multiply(target[i], src1[i], src2[i]);
Chris@366: 	}
Chris@366:     }
Chris@366:     
Chris@366:     clock_t end = clock();
Chris@366: 
Chris@366:     cerr << "Time for naive method: " << float(end - start)/divisor << endl;
Chris@366: 
Chris@366:     start = clock();
Chris@366: 
Chris@366:     for (int j = 0; j < iterations; ++j) {
Chris@366:         v_multiply(target, src1, src2, N);
Chris@366:     }
Chris@366:     
Chris@366:     end = clock();
Chris@366: 
Chris@366:     cerr << "Time for v_multiply: " << float(end - start)/divisor << endl;
Chris@366: 
Chris@366:     return true;
Chris@366: }
Chris@366: 
Chris@366: bool
Chris@366: testPolarToCart()
Chris@366: {
Chris@366:     cerr << "testVectorOps: testing v_polar_to_cartesian" << endl;
Chris@366: 
Chris@366:     const int N = 1024;
Chris@366:     turbot_complex_sample_t target[N];
Chris@366:     turbot_sample_t mag[N];
Chris@366:     turbot_sample_t phase[N];
Chris@366: 
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	mag[i] = drand48();
Chris@366: 	phase[i] = (drand48() * M_PI * 2) - M_PI;
Chris@366:     }
Chris@366: 
Chris@366:     turbot_sample_t mean, first, last, total = 0;
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	turbot_sample_t real = mag[i] * cos(phase[i]);
Chris@366: 	turbot_sample_t imag = mag[i] * sin(phase[i]);
Chris@366: 	if (i == 0) first = real;
Chris@366: 	if (i == N-1) last = imag;
Chris@366: 	total += real;
Chris@366: 	total += imag;
Chris@366:     }
Chris@366:     mean = total / (N*2);
Chris@366:     cerr << "Naive method: mean = " << mean << ", first = " << first
Chris@366: 	 << ", last = " << last << endl;
Chris@366: 
Chris@366:     v_polar_to_cartesian(target, mag, phase, N);
Chris@366: 
Chris@366:     total = 0;
Chris@366: 
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	if (i == 0) first = target[i].re;
Chris@366: 	if (i == N-1) last = target[i].im;
Chris@366: 	total += target[i].re;
Chris@366: 	total += target[i].im;
Chris@366:     }
Chris@366:     mean = total / (N*2);
Chris@366:     cerr << "v_polar_to_cartesian: mean = " << mean << ", first = " << first
Chris@366: 	 << ", last = " << last << endl;
Chris@366: 
Chris@366:     int iterations = 10000;
Chris@366:     cerr << "Iterations: " << iterations << endl;
Chris@366: 	
Chris@366:     cerr << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << endl;
Chris@366:     float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@366: 
Chris@366:     clock_t start = clock();
Chris@366: 
Chris@366:     for (int j = 0; j < iterations; ++j) {
Chris@366: 	for (int i = 0; i < N; ++i) {
Chris@366: 	    target[i].re = mag[i] * cos(phase[i]);
Chris@366: 	    target[i].im = mag[i] * sin(phase[i]);
Chris@366: 	}
Chris@366:     }
Chris@366:     
Chris@366:     clock_t end = clock();
Chris@366: 
Chris@366:     cerr << "Time for naive method: " << float(end - start)/divisor << endl;
Chris@366: 
Chris@366:     start = clock();
Chris@366: 
Chris@366:     for (int j = 0; j < iterations; ++j) {
Chris@366: 	v_polar_to_cartesian(target, mag, phase, N);
Chris@366:     }
Chris@366:     
Chris@366:     end = clock();
Chris@366: 
Chris@366:     cerr << "Time for v_polar_to_cartesian: " << float(end - start)/divisor << endl;
Chris@366: 
Chris@366:     return true;
Chris@366: }
Chris@366: 
Chris@366: bool
Chris@366: testPolarToCartInterleaved()
Chris@366: {
Chris@366:     cerr << "testVectorOps: testing v_polar_interleaved_to_cartesian" << endl;
Chris@366: 
Chris@366:     const int N = 1024;
Chris@366:     turbot_complex_sample_t target[N];
Chris@366:     turbot_sample_t source[N*2];
Chris@366: 
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	source[i*2] = drand48();
Chris@366: 	source[i*2+1] = (drand48() * M_PI * 2) - M_PI;
Chris@366:     }
Chris@366: 
Chris@366:     turbot_sample_t mean, first, last, total = 0;
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	turbot_sample_t real = source[i*2] * cos(source[i*2+1]);
Chris@366: 	turbot_sample_t imag = source[i*2] * sin(source[i*2+1]);
Chris@366: 	if (i == 0) first = real;
Chris@366: 	if (i == N-1) last = imag;
Chris@366: 	total += real;
Chris@366: 	total += imag;
Chris@366:     }
Chris@366:     mean = total / (N*2);
Chris@366:     cerr << "Naive method: mean = " << mean << ", first = " << first
Chris@366: 	 << ", last = " << last << endl;
Chris@366: 
Chris@366:     v_polar_interleaved_to_cartesian(target, source, N);
Chris@366: 
Chris@366:     total = 0;
Chris@366: 
Chris@366:     for (int i = 0; i < N; ++i) {
Chris@366: 	if (i == 0) first = target[i].re;
Chris@366: 	if (i == N-1) last = target[i].im;
Chris@366: 	total += target[i].re;
Chris@366: 	total += target[i].im;
Chris@366:     }
Chris@366:     mean = total / (N*2);
Chris@366:     cerr << "v_polar_interleaved_to_cartesian: mean = " << mean << ", first = " << first
Chris@366: 	 << ", last = " << last << endl;
Chris@366: 
Chris@366:     int iterations = 10000;
Chris@366:     cerr << "Iterations: " << iterations << endl;
Chris@366: 	
Chris@366:     cerr << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << endl;
Chris@366:     float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@366: 
Chris@366:     clock_t start = clock();
Chris@366: 
Chris@366:     for (int j = 0; j < iterations; ++j) {
Chris@366: 	for (int i = 0; i < N; ++i) {
Chris@366: 	    target[i].re = source[i*2] * cos(source[i*2+1]);
Chris@366: 	    target[i].im = source[i*2] * sin(source[i*2+1]);
Chris@366: 	}
Chris@366:     }
Chris@366:     
Chris@366:     clock_t end = clock();
Chris@366: 
Chris@366:     cerr << "Time for naive method: " << float(end - start)/divisor << endl;
Chris@366: 
Chris@366:     start = clock();
Chris@366: 
Chris@366:     for (int j = 0; j < iterations; ++j) {
Chris@366: 	v_polar_interleaved_to_cartesian(target, source, N);
Chris@366:     }
Chris@366:     
Chris@366:     end = clock();
Chris@366: 
Chris@366:     cerr << "Time for v_polar_interleaved_to_cartesian: " << float(end - start)/divisor << endl;
Chris@366: 
Chris@366:     return true;
Chris@366: }
Chris@366: 
Chris@366: bool
Chris@366: testVectorOps()
Chris@366: {
Chris@366:     if (!testMultiply()) return false;
Chris@366:     if (!testPolarToCart()) return false;
Chris@366:     if (!testPolarToCartInterleaved()) return false;
Chris@366:     
Chris@366:     return true;
Chris@366: }
Chris@366: 
Chris@366: }
Chris@366: 
Chris@366: }
Chris@366: