/** @file  tst_psVectorFFT.c
*
*  @brief Contains the tests for psFFT.[ch]
*
*
*  @author Robert DeSonia, MHPCC
*
*  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
*  @date $Date: 2004-08-06 22:34:06 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/

#include <math.h>
#include <float.h>

#include "psTest.h"
#include "pslib.h"

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

static int testVectorFFT( void );
static int testVectorRealImaginary( void );
static int testVectorComplex( void );
static int testVectorConjugate( void );
static int testVectorPowerSpectrum( void );

testDescription tests[] = {
                              {
                                  testVectorFFT, 600, "psVectorFFT", 0, false
                              },
                              {
                                  testVectorRealImaginary, 601, "psVectorRealImaginary", 0, false
                              },
                              {
                                  testVectorComplex, 602, "psVectorComplex", 0, false
                              },
                              {
                                  testVectorConjugate, 603, "psVectorConjugate", 0, false
                              },
                              {
                                  testVectorPowerSpectrum, 604, "psVectorPowerSpectrum", 0, false
                              },
                              {
                                  NULL
                              }
                          };

int main( int argc, char* argv[] )
{
    psLogSetLevel( PS_LOG_INFO );

    return ( ! runTestSuite( stderr, "psFFT", tests, argc, argv ) );
}

int testVectorFFT( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;
    psVector* vec3 = NULL;

    /*
    1. assign a vector to a sinisoid
    2. perform a forward transform
    3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    4. perform a reverse transform
    5. compare to original (should be equal to within a reasonable error)
    */

    // 1. assign a vector to a sinisoid
    vec = psVectorAlloc( 100, PS_TYPE_F32 );
    vec->n = vec->nalloc;
    for ( unsigned int n = 0; n < 100; n++ ) {
        vec->data.F32[ n ] = sinf( ( psF32 ) n / 50.0f * M_PI );
    }

    // 2. perform a forward transform
    vec2 = psVectorFFT( NULL, vec, PS_FFT_FORWARD );
    if ( vec2->type.type != PS_TYPE_C32 ) {
        psError( __func__, "FFT didn't produce complex values?" );
        return 1;
    }


    // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    for ( unsigned int n = 0; n < 100; n++ ) {
        if ( n == 1 || n == 99 ) {
            if ( fabsf( cabsf( vec2->data.C32[ n ] ) - 50.0f ) > 0.1f ) {
                psError( __func__, "FFT didn't work for vector (n=%d)", n );
                return 2;
            }
        } else {
            if ( fabsf( cabsf( vec2->data.C32[ n ] ) ) > 0.1f ) {
                psError( __func__, "FFT didn't work for vector (n=%d)", n );
                return 3;
            }
        }
    }

    // 4. perform a reverse transform
    vec3 = psVectorFFT( NULL, vec2, PS_FFT_REVERSE );
    if ( vec3->type.type != PS_TYPE_C32 ) {
        psError( __func__, "FFT didn't produce complex values?" );
        return 4;
    }
    for ( unsigned int 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 ) {
            psError( __func__, "Reverse FFT didn't give me the original vector back (n=%d) (%.2f vs %.2f)",
                     n, vecVal, val );
            return 5;
        }
    }

    psFree( vec );
    psFree( vec2 );
    psFree( vec3 );

    return 0;
}

int testVectorRealImaginary( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;
    psVector* vec3 = NULL;

    /*
    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
    */

    // 1. create a C32 complex vector with distinctly different real and imaginary parts.
    vec = psVectorAlloc( 100, PS_TYPE_C32 );
    vec->n = vec->nalloc;
    for ( unsigned int n = 0; n < 100; n++ ) {
        vec->data.C32[ n ] = n + I * ( n * 2 );
    }

    // 2. call psVectorReal and psVectorImaginary
    vec2 = psVectorReal( vec2, vec );
    if ( vec2 == NULL ) {
        psError( __func__, "psVectorReal returned a NULL?" );
        return 1;
    }
    if ( vec2->type.type != PS_TYPE_F32 ) {
        psError( __func__, "psVectorReal returned a wrong type (%d)?",
                 vec2->type.type );
        return 2;
    }

    vec3 = psVectorImaginary( vec3, vec );
    if ( vec3 == NULL ) {
        psError( __func__, "psVectorImaginary returned a NULL?" );
        return 3;
    }
    if ( vec3->type.type != PS_TYPE_F32 ) {
        psError( __func__, "psVectorImaginary returned a wrong type (%d)?",
                 vec3->type.type );
        return 4;
    }

    // 3. compare results to the real/imaginary components of input
    for ( unsigned int n = 0; n < 100; n++ ) {
        psF32 r = n;
        psF32 i = ( n * 2 );
        if ( fabsf( vec2->data.F32[ n ] - r ) > FLT_EPSILON ) {
            psError( __func__, "psVectorReal didn't return the real portion at n=%d",
                     n );
            return 5;
        }
        if ( fabsf( vec3->data.F32[ n ] - i ) > FLT_EPSILON ) {
            psError( __func__, "psVectorImaginary didn't return the real portion at n=%d",
                     n );
            return 6;
        }
    }

    psFree( vec );
    psFree( vec2 );
    psFree( vec3 );

    return 0;
}

int testVectorComplex( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;
    psVector* vec3 = NULL;

    /*
    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.
    */

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

    // 2. call psVectorComplex
    vec3 = psVectorComplex( vec3, vec, vec2 );

    // 3. verify that the result is a psC32
    if ( vec3->type.type != PS_TYPE_C32 ) {
        psError( __func__, "Vector Type from psVectorComplex is not complex? (%d)",
                 vec3->type.type );
        return 1;
    }

    // 4. call psVectorReal and psVectorImaginary on step 2 results (not needed, just use crealf/cimagf)
    // 5. compare step 4 results to input.
    for ( unsigned int 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 ) {
            psError( __func__, "psVectorComplex result is invalid (n=%d, %.2f+%.2fi)",
                     n, crealf( vec3->data.C32[ n ] ), cimagf( vec3->data.C32[ n ] ) );
            return 2;
        };
    }


    // 6. create a psF32 and a psF64 vector of the same size
    vec2 = psVectorRecycle( vec2, 100, PS_TYPE_F64 );

    // 7. call psVectorComplex
    psLogMsg( __func__, PS_LOG_INFO, "Following should be an error (type mismatch)." );
    vec3 = psVectorComplex( vec3, vec, vec2 );
    // 8. verify that an appropriate error occurred. (this partially has to be done via inspection)
    if ( vec3 != NULL ) {
        psError( __func__, "psVectorComplex returned a vector though input types mismatched." );
        return 3;
    }

    // 9. create two psf32 vectors of different sizes
    vec2 = psVectorRecycle( vec2, 200, PS_TYPE_F32 );

    // 10. call psVectorComplex
    vec3 = psVectorComplex( vec3, vec, vec2 );

    // 11. verify thet an appropriate error occurred. (actually, it isn't an error...)
    if ( vec3->n != 100 ) {
        psError( __func__, "psVectorComplex returned a larger vector than the input supported (%d).", vec3->n );
        return 4;
    }

    psFree( vec );
    psFree( vec2 );
    psFree( vec3 );

    return 0;
}

int testVectorConjugate( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;

    /*
    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)
    */

    // 1. create a psC32 with unique real and imaginary values.
    vec = psVectorAlloc( 100, PS_TYPE_C32 );
    vec->n = vec->nalloc;
    for ( unsigned int n = 0; n < 100; n++ ) {
        vec->data.C32[ n ] = n + I * ( n * 2 );
    }

    // 2. call psVectorConjugate
    vec2 = psVectorConjugate( vec2, vec );

    // 3. verify result is psC32
    if ( vec2->type.type != PS_TYPE_C32 ) {
        psError( __func__, "the psVectorConjugate didn't return a C32 vector" );
        return 1;
    }

    // 4. verify each value is conjugate of input (a+bi -> a-bi)
    for ( unsigned int 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 ) {
            psError( __func__, "psVectorComplex result is invalid (n=%d, %.2f+%.2fi)",
                     n, crealf( vec2->data.C32[ n ] ), cimagf( vec2->data.C32[ n ] ) );
            return 2;
        };
    }

    psFree( vec );
    psFree( vec2 );

    return 0;
}

int testVectorPowerSpectrum( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;
    psF32 val;

    /*
    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
    */

    // 1. create a psC32 vector with unique real and imaginary components
    vec = psVectorAlloc( 100, PS_TYPE_C32 );
    vec->n = vec->nalloc;
    for ( unsigned int n = 0; n < 100; n++ ) {
        vec->data.C32[ n ] = n + I * sinf( ( ( psF32 ) n ) / 50.f * M_PI );
    }

    // 2. call psVectorPowerSpectrum
    vec2 = psVectorPowerSpectrum( vec2, vec );

    // 3. verify result is psF32
    if ( vec2->type.type != PS_TYPE_F32 ) {
        psError( __func__, "the type was not PS_TYPE_F32." );
        return 1;
    }

    // 4. verify the values are the square of the absolute values of the original
    //   (ADD specifies something else)
    //   P_0 = |C_0|^2/N^2
    //   P_j = (|C_j|^2+|C_N-j|^2)/N^2
    //   P_N/2 = |C_N/2|^2/N^2
    //  where j = 1,2,...,(N/2-1)

    val = cabsf( vec->data.C32[ 0 ] ) * cabsf( vec->data.C32[ 0 ] ) / 100 / 100;
    if ( fabsf( vec2->data.F32[ 0 ] - val ) > FLT_EPSILON ) {
        psError( __func__, "psVectorPowerSpectrum result is invalid (n=0, %.2f %.2f)",
                 vec2->data.F32[ 0 ], val );
        return 2;
    };

    for ( unsigned int 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;

        if ( fabsf( val - vec2->data.F32[ n ] ) > FLT_EPSILON ) {
            psError( __func__, "psVectorPowerSpectrum result is invalid (n=%d, %.2f %.2f)",
                     n, vec2->data.F32[ n ], val );
            return 2;
        };
    }

    val = cabsf( vec->data.C32[ 50 ] ) * cabsf( vec->data.C32[ 50 ] ) / 100 / 100;
    if ( fabsf( vec2->data.F32[ 50 ] - val ) > FLT_EPSILON ) {
        psError( __func__, "psVectorPowerSpectrum result is invalid (n=50, %.2f %.2f)",
                 vec2->data.F32[ 0 ], val );
        return 2;
    };

    psFree( vec );
    psFree( vec2 );

    return 0;
}
