/**
 * \defgroup libbfin32 32bit algorithm collection of Blackfin DSP
 */




/**
 * \file main.c
 * \brief test main file.
 */

#include <stdio.h>
#include <math.h>
#include "fr32_nco.h"

/**
 * \brief define AVERAGE when calculate the whole average of sine and cosine.
 * \details
 * undef AVERAGE calcs the standard deviation of sine and cosine.
 */
#define AVERAGE

/**
 * \brief test step. 1 tests 4G points, 2 tests 2G points.
 */
#define STEP 1

extern void fract32_sincos( fract32 angle, fract32 *c, fract32 *s );
#define NUM2P31 (LONG_MAX+1LL)



int main(void)
{
    fract32 c, s;
    long long int i,j;
    double theta;
    double co, si;
    double ecsum = 0, essum = 0;
    double maxes = 0, mines = 0, maxec = 0, minec = 0;
    long long int count = 0;
    TNCO32 nco;


    fr32_nco_init(&nco, STEP);
    for( j=LONG_MIN; j<LONG_MAX; j+=0x01000000 )
    {
        for ( i=0; i<0x01000000; i+=STEP)
        {

            theta = (i+j)*M_PI/NUM2P31;
#ifdef AVERAGE
            fr32_nco( &nco, &c, &s, 1 );
            co = cos(theta)*NUM2P31 - c;
            si = sin(theta)*NUM2P31 - s;
#else
            co = (cos(theta)-(float)cos(theta))*NUM2P31;
            si = (sin(theta)-(float)sin(theta))*NUM2P31;
#endif

#ifdef AVERAGE
            ecsum += co;
            essum += si;
#else
            ecsum += co*co;
//            ecsum += (co+1.272776)*(co+1.272776);
            essum += si*si;
#endif
            maxec = ( co > maxec )? co : maxec;
            minec = ( co < minec )? co : minec;
            maxes = ( si > maxes )? si : maxes;
            mines = ( si < mines )? si : mines;
            count++;

        }
#ifdef AVERAGE
        printf("%3d ave cos err = %f, ave sin err = %f, max cos err : %f, min cos err : %f, max sin err : %f, min sin err : %f \n",
                (int)(j>>24),
                ecsum/count,
                essum/count,
                maxec,
                minec,
                maxes,
                mines
                );
#else
        printf("%3d stddev cos err = %f, stddev sin err = %f,  \n",
                (int)(j>>24),
                sqrt(ecsum/count),
                sqrt(essum/count)
                );
#endif
    }


	return 0;
}

