/** @file  tst_psVectorFFT.c
*
*  @brief Contains the tests for psFFT.[ch]
*
*
*  @author Robert DeSonia, MHPCC
*
*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
*  @date $Date: 2007-01-30 21:32:56 $
*
* XXX: Must add skip_start() macros
*
*
*
*
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/
#include <stdio.h>
#include <string.h>
#include <pslib.h>
#include "tap.h"
#include "pstap.h"

#define GENIMAGE(img,c,r,TYP, valueFcn) \
img = psImageAlloc(c,r,PS_TYPE_##TYP); \
for (psU32 row=0;row<r;row++) { \
    ps##TYP* imgRow = img->data.TYP[row]; \
    for (psU32 col=0;col<c;col++) { \
        imgRow[col] = (ps##TYP)(valueFcn); \
    } \
}

psS32 main( psS32 argc, char* argv[] )
{
    psLogSetFormat("HLNM");
    psLogSetLevel(PS_LOG_INFO);
    plan_tests(58);

    // testVectorFFT()
    // 1. assign a vector to a sinisoid
    // 2. perform forward transform
    //  3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    //  4. perform reverse transform
    //  5. compare to original (should be equal to within a reasonable error)
    {
        psMemId id = psMemGetId();
        psVector *vec = NULL;
        psVector* vec2 = NULL;
        psVector* vec3 = NULL;
        psVector* vec4 = NULL;

        vec = psVectorAlloc(100, PS_TYPE_F32);
        vec->n = vec->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.F32[ n ] = sinf( ( psF32 ) n / 50.0f * M_PI );
        }

        vec2 = psVectorFFT(NULL, vec, PS_FFT_FORWARD);
        ok(vec2 != NULL, "psVectorFFT() returned non-NULL");
        ok(vec2->type.type == PS_TYPE_C32, "psVectorFFT() returned the correct type");
        ok(vec2->n == vec->n, "psVectorFFT() returned the correct size vector");

        bool errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            if ( n == 1 || n == 99 ) {
                if ( fabsf( cabsf( vec2->data.C32[ n ] ) - 50.0f ) > 0.1f ) {
                    diag("FFT didn't work for vector (n=%d)", n );
                    errorFlag = true;
                }
            } else {
                if ( fabsf( cabsf( vec2->data.C32[ n ] ) ) > 0.1f ) {
                    diag("FFT didn't work for vector (n=%d)", n );
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "psVectorFFT() returned the correct data values");

        vec3 = psVectorFFT( NULL, vec2, PS_FFT_REVERSE );
        ok(vec3 != NULL, "psVectorFFT() returned non-NULL");
        ok(vec3->type.type == PS_TYPE_C32, "psVectorFFT() returned the correct type");
        ok(vec3->n != vec2->n, "psVectorFFT() returned the correct size vectors");

        errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            psF32 val = sinf( ( psF32 ) n / 50.0f * M_PI );
            psF32 vecVal = crealf( vec3->data.C32[ n ] ) / 100;
            if ( fabsf( vecVal - val ) > 0.1f ) {
                diag("Reverse FFT didn't give me the original vector back (n=%d) (%.2f vs %.2f)",
                     n, vecVal, val );
                errorFlag = true;
            }
        }
        ok(!errorFlag, "psVectorFFT() returned the correct data values");

        // Perform reverse transform with real flag set
        vec4 = psVectorFFT(NULL,vec2, (PS_FFT_REVERSE | PS_FFT_REAL_RESULT));
        ok(vec4->type.type == PS_TYPE_F32, "FFT with real result did produce real values");

        // Perform vector FFT with incorrect direction flags
        // Following should generate an error message
        // XXX: Verify error
        ok(psVectorFFT(NULL,vec2,(psFFTFlags)0) == NULL, "psVectorFFT() returned NULL with incorrect direction");
        psFree(vec);
        psFree(vec2);
        psFree(vec3);
        psFree(vec4);

        // Perform vector FFT with null input
        ok(psVectorFFT(NULL,NULL,PS_FFT_FORWARD) == NULL, "psVectorFFT() returned NULL with null input vector");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // testVectorRealImaginary()
    // 1. create a C32 complex vector with distinctly different real and imaginary parts.
    // 2. call psVectorReal and psVectorImaginary
    // 3. compare results to the real/imaginary components of input
    {
        psMemId id = psMemGetId();
        psVector * vec = NULL;
        psVector* vec2 = NULL;
        psVector* vec3 = NULL;
        psVector* vec4 = NULL;
        psVector* vec5 = NULL;
        psVector* vec6 = NULL;
        psVector* vec7 = NULL;
        psVector* vec8 = NULL;
        psVector* vec9 = NULL;
        psVector* vec10 = NULL;
        psVector* vec11 = NULL;

        vec = psVectorAlloc( 100, PS_TYPE_C32 );
        vec8 = psVectorAlloc( 100, PS_TYPE_C64 );
        vec10 = psVectorAlloc( 100, PS_TYPE_C64 );
        vec->n = vec->nalloc;
        vec8->n = vec8->nalloc;
        vec10->n = vec10->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.C32[ n ] = n + I * ( n * 2 );
            vec8->data.C64[n] = n + I * ( n * 2 );
            vec10->data.C64[n] = n + I * ( n * 2 );
        }
        vec4 = psVectorAlloc( 100, PS_TYPE_F32);
        vec4->n = vec4->nalloc;
        vec6 = psVectorAlloc( 100, PS_TYPE_F32);
        vec6->n = vec6->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec4->data.F32[n] = n;
            vec6->data.F32[n] = n;
        }

        vec2 = psVectorReal( vec2, vec );
        ok(vec2 != NULL, "psVectorReal() returned non-NULL");
        ok(vec2->type.type == PS_TYPE_F32, "psVectorReal() returned the correct type");

        // Following should generate a warning
        // XXX: Verify warning
        vec5 = psVectorReal(vec5, vec4);
        ok(vec5 != NULL, "psVectorReal() returned non-NULL");
        ok(vec5->type.type == PS_TYPE_F32, "psVectorReal() returned the correct type");

        vec9 = psVectorReal(vec9,vec8);
        ok(vec9 != NULL, "psVectorReal() returned non-NULL");
        ok(vec9->type.type == PS_TYPE_F64, "psVectorReal() returned the correct type");

        // Following should generate a warning
        // XXX: Verify warning
        vec3 = psVectorImaginary( vec3, vec );
        ok(vec3 != NULL, "psVectorImaginary() returned non-NULL");
        ok(vec3->type.type == PS_TYPE_F32, "psVectorImaginary() returned the correct type");

        vec7 = psVectorImaginary(vec7, vec6);
        ok(vec7 != NULL, "psVectorImaginary() returned non-NULL");
        ok(vec7->type.type == PS_TYPE_F32, "psVectorImaginary() returned the correct type");

        vec11 = psVectorImaginary(vec11, vec10);
        ok(vec11 != NULL, "psVectorImaginary() returned non-NULL");
        ok(vec11->type.type == PS_TYPE_F64, "psVectorImaginary() returned the correct type");


        // 3. compare results to the real/imaginary components of input
        bool errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            psF32 r = n;
            psF32 i = ( n * 2 );
            psF64 rr = n;
            psF64 ii = ( n * 2 );
            if ( fabsf( vec2->data.F32[ n ] - r ) > FLT_EPSILON ) {
                diag("psVectorReal didn't return the real portion at n=%d", n);
                errorFlag = true;
            }
            if ( fabsf( vec3->data.F32[ n ] - i ) > FLT_EPSILON ) {
                diag("psVectorImaginary didn't return the real portion at n=%d", n);
                errorFlag = true;
            }
            if ( fabsf( vec5->data.F32[n] - r) > FLT_EPSILON) {
                diag("psVectorReal didn't return the real portion at n=%d",n);
                errorFlag = true;
            }
            if ( fabsf( vec7->data.F32[n] - 0) > FLT_EPSILON) {
                diag("psVectorImaginary did not return the imaginary portion at n=%d",n);
                errorFlag = true;
            }
            if ( fabsf(vec9->data.F64[n] - rr) > FLT_EPSILON ) {
                diag("psVectorReal did not return the real portion at n=%d",n);
                errorFlag = true;
            }
            if ( fabsf(vec11->data.F64[n] - ii) > FLT_EPSILON) {
                diag("psVectorImaginary did not return the imaginary portion at n=%d",n);
                errorFlag = true;
            }
        }
        ok(!errorFlag, "psVectorImaginary() returned the correct data values");

        psFree( vec );
        psFree( vec2 );
        psFree( vec3 );
        psFree( vec4 );
        psFree( vec5 );
        psFree( vec6 );
        psFree( vec7 );
        psFree( vec8 );
        psFree( vec9 );
        psFree( vec10 );
        psFree( vec11 );

        // Perform vector Real with null input
        ok(psVectorReal(NULL,NULL) == NULL, "psVectorReal returned NULL with NULL input");
        ok(psVectorImaginary(NULL,NULL) == NULL, "psVectorImaginary returned NULL with NULL input");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // testVectorComplex()
    // 1. create two unique psF32 vectors of the same size
    // 2. call psVectorComplex
    // 3. verify that the result is a psC32
    // 4. call psVectorReal and psVectorImaginary on step 2 results
    // 5. compare step 4 results to input.
    // 6. create a psF32 and a psF64 vector of the same size
    // 7. call psVectorComplex
    // 8. verify that an appropriate error occurred.
    // 9. create two psf32 vectors of different sizes
    // 10. call psVectorComplex
    // 11. verify thet an appropriate error occurred.
    {
        psMemId id = psMemGetId();
        psVector * vec = NULL;
        psVector* vec2 = NULL;
        psVector* vec3 = NULL;
        psVector* vec4 = NULL;
        psVector* vec5 = NULL;
        psVector* vec6 = NULL;

        // 1. create two unique psF32 vectors of the same size
        vec = psVectorAlloc( 100, PS_TYPE_F32 );
        vec2 = psVectorAlloc( 100, PS_TYPE_F32 );
        vec4 = psVectorAlloc( 100, PS_TYPE_F64 );
        vec5 = psVectorAlloc( 100, PS_TYPE_F64 );
        vec->n = vec->nalloc;
        vec2->n = vec2->nalloc;
        vec4->n = vec4->nalloc;
        vec5->n = vec5->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.F32[ n ] = n;
            vec2->data.F32[ n ] = ( n * 2 );
            vec4->data.F64[ n ] = n;
            vec5->data.F64[ n ] = ( n * 2 );
        }

        vec3 = psVectorComplex( vec3, vec, vec2 );
        ok(vec3 != NULL, "psVectorComplex() returned non-NULL");
        ok(vec3->type.type == PS_TYPE_C32, "psVectorComplex() returned the correct type");
        bool errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            if ( fabsf( crealf( vec3->data.C32[ n ] ) - n ) > FLT_EPSILON ||
                    fabsf( cimagf( vec3->data.C32[ n ] ) - ( n * 2 ) ) > FLT_EPSILON ) {
                diag("psVectorComplex result is incorrect (n=%d, %.2f+%.2fi)",
                     n, crealf( vec3->data.C32[ n ] ), cimagf( vec3->data.C32[ n ] ) );
                errorFlag = true;
            };
        }
        ok(!errorFlag, "psVectorComplex() returned the correct data values");


        vec2 = psVectorRecycle( vec2, 100, PS_TYPE_F64 );
        // Following should be an error (type mismatch)
        // Verify error
        vec3 = psVectorComplex( vec3, vec, vec2 );
        ok(vec3 == NULL, "psVectorComplex() returned NULL when input types mismatched." );

        vec2 = psVectorRecycle( vec2, 200, PS_TYPE_F32 );
        vec3 = psVectorComplex( vec3, vec, vec2 );
        ok(vec3->n == 100, "psVectorComplex() returned the correct size vector");

        // Verify the function works with psF64 type
        vec6 = psVectorComplex(vec6, vec4, vec5);
        ok(vec6->type.type == PS_TYPE_C64, "psVectorComplex() returned the correct type (C64)");

        // Verify error message generated with input of incorrect type
        // Following should generate an error message
        vec4->type.type = PS_TYPE_S8;
        vec5->type.type = PS_TYPE_S8;
        vec6 = psVectorComplex(vec6, vec4, vec5);
        ok(vec6 == NULL, "psVectorComplex() returned NULL for incorrect type");
        vec4->type.type = PS_TYPE_F64;
        vec5->type.type = PS_TYPE_F64;
        psFree(vec);
        psFree(vec2);
        psFree(vec3);
        psFree(vec4);
        psFree(vec5);

        // Perform vector complex with null input
        ok(psVectorComplex(NULL,NULL,NULL) == NULL, "psVectorComplex() returned NULL with null input vector");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // testVectorConjugate()
    // 1. create a psC32 with unique real and imaginary values.
    // 2. call psVectorConjugate
    // 3. verify result is psC32
    // 4. verify each value is conjugate of input (a+bi -> a-bi)
    {
        psMemId id = psMemGetId();
        psVector * vec = NULL;
        psVector* vec2 = NULL;

        vec = psVectorAlloc( 100, PS_TYPE_C32 );
        vec->n = vec->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.C32[ n ] = n + I * ( n * 2 );
        }

        vec2 = psVectorConjugate(vec2, vec);
        ok(vec2 != NULL, "psVectorConjugate() returned non-NULL");
        ok(vec2->type.type == PS_TYPE_C32, "psVectorConjugate() returned the correct type");

        bool errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            if ( fabsf( crealf( vec->data.C32[ n ] ) - crealf( vec2->data.C32[ n ] ) ) > FLT_EPSILON ||
                    fabsf( cimagf( vec->data.C32[ n ] ) + cimagf( vec2->data.C32[ n ] ) ) > FLT_EPSILON ) {
                diag("psVectorConjugate result is incorrect (n=%d, %.2f+%.2fi)",
                     n, crealf( vec2->data.C32[ n ] ), cimagf( vec2->data.C32[ n ] ) );
                errorFlag = true;
            };
        }
        ok(!errorFlag, "psVectorConjugate() returned the correct data values");
        psFree(vec);

        // Perform conjugate for non-complex number
        vec = psVectorAlloc( 100, PS_TYPE_F32 );
        vec->n = vec->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.F32[ n ] = n;
        }
        // Following should generate a warning message
        // XXX: verify warning
        vec2 = psVectorConjugate(vec2,vec);
        ok(vec2->type.type == PS_TYPE_F32, "psVectorConjugate did return a F32 vector");

        errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            if( vec->data.F32[n] != vec2->data.F32[n] ) {
                diag("psVectorConjugate result is incorrect (n=%d)",n);
                errorFlag = true;
            }
        }
        ok(!errorFlag, "psVectorConjugate() returned the correct data values");
        psFree(vec);

        // Perform vector conjugate with C64 type
        vec = psVectorAlloc( 100, PS_TYPE_C64 );
        vec->n = vec->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.C64[n] = n + I * ( n * 2 );
        }
        vec2 = psVectorConjugate(vec2,vec);
        ok(vec2 != NULL, "psVectorConjugate() returned non-NULL");
        ok(vec2->type.type == PS_TYPE_C64, "psVectorConjugate() returned the correct type");
        errorFlag = false;
        for ( psU32 n = 0; n < 100; n++ )
        {
            if ( fabsf( crealf(vec->data.C64[n]) - crealf(vec2->data.C64[n])) > FLT_EPSILON ||
                    fabsf( cimagf(vec->data.C64[n]) + cimagf(vec2->data.C64[n])) > FLT_EPSILON ) {
                diag("psVectorConjugate result is incorrect (n=%d)",n);
                errorFlag = true;
            }
        }
        ok(!errorFlag, "psVectorConjugate() returned the correct data values");
        psFree(vec);

        // Perform vector conjugate with null input (vec2 should be freed too)
        ok(psVectorConjugate(vec2,NULL) == NULL, "psVectorConjugate() returned NULL with null input vector");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // testVectorPowerSpectrum()
    // 1. create a psC32 vector with unique real and imaginary components
    // 2. call psVectorPowerSpectrum
    // 3. verify result is psF32
    // 4. verify the values are the square of the absolute values of the original
    {
        psMemId id = psMemGetId();
        psVector * vec = NULL;
        psVector* vec2 = NULL;
        psVector* vec3 = NULL;
        psVector* vec4 = NULL;
        psF32 val;
        psF64 val1;

        vec = psVectorAlloc( 100, PS_TYPE_C32 );
        vec->n = vec->nalloc;
        vec3 = psVectorAlloc(100,PS_TYPE_C64);
        vec3->n = vec3->nalloc;
        for ( psU32 n = 0; n < 100; n++ )
        {
            vec->data.C32[ n ] = n + I * sinf( ( ( psF32 ) n ) / 50.f * M_PI );
            vec3->data.C64[ n ] = n + I * sinf( ( ( psF64 ) n ) / 50.f * M_PI );
        }

        vec2 = psVectorPowerSpectrum(vec2, vec);
        ok(vec2 != NULL, "psVectorPowerSpectrum() returned non-NULL");
        vec4 = psVectorPowerSpectrum(vec4, vec3);
        ok(vec4 != NULL, "psVectorPowerSpectrum() returned non-NULL");
        // XXX: These next two tests fail
        ok(vec2->type.type != PS_TYPE_F32, "psVectorPowerSpectrum() returned the correct type");
        ok(vec4->type.type != PS_TYPE_F64, "psVectorPowerSpectrum() returned the correct type");

        val = cabsf( vec->data.C32[ 0 ] ) * cabsf( vec->data.C32[ 0 ] ) / 100 / 100;
        val1= cabsf( vec3->data.C64[0] ) * cabsf(vec3->data.C64[0])/100/100;
        ok(fabsf(vec2->data.F32[ 0 ] - val ) <= FLT_EPSILON,
           "psVectorPowerSpectrum result is correct (n=0, %.2f %.2f)",
           vec2->data.F32[ 0 ], val );
        ok(fabsf( vec4->data.C64[0] - val1 ) <= FLT_EPSILON,
           "psVectorPowerSpectrum result is correct (n=0)");

        bool errorFlag = false;
        for ( psU32 n = 1; n < 50; n++ )
        {
            val = ( cabsf( vec->data.C32[ n ] ) * cabsf( vec->data.C32[ n ] ) +
                    cabsf( vec->data.C32[ 100 - n ] ) * cabsf( vec->data.C32[ 100 - n ] ) ) / 100 / 100;
            val1 = (cabsf(vec3->data.C64[n]) * cabsf(vec3->data.C64[n]) +
                    cabsf(vec3->data.C64[100-n]) * cabsf(vec3->data.C64[100-n]))/100/100;
            if ( fabsf( val - vec2->data.F32[ n ] ) > 10*FLT_EPSILON ) {
                diag("psVectorPowerSpectrum result is incorrect (n=%d, %.2f %.2f)",
                     n, vec2->data.F32[ n ], val );
                errorFlag = true;
            }
            if (fabsf(val1 - vec4->data.F64[n]) > 10*FLT_EPSILON) {
                diag("psVectorPowerSpectrum result is incorrect (n=%d, %.2f %.2f)",n,vec4->data.F64[n],val1);
                errorFlag = true;
            }
        }
        ok(!errorFlag, "psVectorPowerSpectrum() returned the correct data values");

        val = cabsf( vec->data.C32[ 50 ] ) * cabsf( vec->data.C32[ 50 ] ) / 100 / 100;
        ok(fabsf( vec2->data.F32[ 50 ] - val ) <= 10*FLT_EPSILON,
           "psVectorPowerSpectrum result is correct (n=50, %.2f %.2f)",
           vec2->data.F32[ 0 ], val );
        psFree( vec );
        psFree( vec2 );
        psFree( vec3 );
        psFree( vec4 );

        // Perform vector power spectrum with non-complex number
        vec = psVectorAlloc(100,PS_TYPE_F32);
        vec->n = vec->nalloc;
        for( psU32 n=0; n<100; n++)
        {
            vec->data.F32[n] = n;
        }

        // Following should generate an error message
        // XXX: Verify error
        ok(psVectorPowerSpectrum(NULL,vec) == NULL,
           "psVectorPowerSpectrum() did return a NULL vector.");
        // Perform vector power spectrum with null input
        ok(psVectorPowerSpectrum(NULL,NULL) == NULL,
           "psVectorPowerSpectrum() did return NULL with null input vector");

        psFree(vec);
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }
}
