/** @file  tst_psVectorFFT.c
*
*  @brief Contains the tests for psFFT.[ch]
*
*
*  @author Robert DeSonia, MHPCC
*
*  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-01-13 21:25:44 $
*
*  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 (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); \
    } \
}

static psS32 testVectorFFT( void );
static psS32 testVectorRealImaginary( void );
static psS32 testVectorComplex( void );
static psS32 testVectorConjugate( void );
static psS32 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}
                          };

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

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

psS32 testVectorFFT( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;
    psVector* vec3 = NULL;
    psVector* vec4 = 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 ( psU32 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(PS_ERR_UNKNOWN,true, "FFT didn't produce complex values?" );
        return 1;
    }

    // 2a. verify output vector is same size as input
    if ( vec2->n != vec->n ) {
        psError(PS_ERR_UNKNOWN,true, "FFT didn't return vector with same size as input");
        return 10;
    }


    // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    for ( psU32 n = 0; n < 100; n++ ) {
        if ( n == 1 || n == 99 ) {
            if ( fabsf( cabsf( vec2->data.C32[ n ] ) - 50.0f ) > 0.1f ) {
                psError(PS_ERR_UNKNOWN,true, "FFT didn't work for vector (n=%d)", n );
                return 2;
            }
        } else {
            if ( fabsf( cabsf( vec2->data.C32[ n ] ) ) > 0.1f ) {
                psError(PS_ERR_UNKNOWN,true, "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(PS_ERR_UNKNOWN,true, "FFT didn't produce complex values?" );
        return 4;
    }
    if ( vec3->n != vec2->n ) {
        psError(PS_ERR_UNKNOWN,true, "FFT didn't return vector with same size as input");
        return 40;
    }
    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 ) {
            psError(PS_ERR_UNKNOWN,true, "Reverse FFT didn't give me the original vector back (n=%d) (%.2f vs %.2f)",
                    n, vecVal, val );
            return 5;
        }
    }

    // Perform a reverse transform with real flag set
    vec4 = psVectorFFT(NULL,vec2, (PS_FFT_REVERSE | PS_FFT_REAL_RESULT));
    if(vec4->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"FFT with real result did not produce real values");
        return 80;
    }

    // Perform vector FFT with invalid direction flags
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message.");
    if ( psVectorFFT(NULL,vec2,(psFFTFlags)0) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid direction");
        return 70;
    }

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

    // Perform vector FFT with null input
    if ( psVectorFFT(NULL,NULL,PS_FFT_FORWARD) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with null input vector");
        return 64;
    }

    return 0;
}

psS32 testVectorRealImaginary( void )
{
    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;

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

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

    psLogMsg(__func__,PS_LOG_INFO,"Following should generate a warning.");
    vec5 = psVectorReal(vec5, vec4);
    if (vec5 == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psVectorReal returned NULL");
        return 10;
    }
    if ( vec5->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"psVectorReal returned the wrong type");
        return 11;
    }

    vec9 = psVectorReal(vec9,vec8);
    if(vec9 == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psVectorReal returned NULL");
        return 20;
    }
    if(vec9->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN,true,"psVectorReal returned the wrong type");
        return 21;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should generate a warning.");
    vec3 = psVectorImaginary( vec3, vec );
    if ( vec3 == NULL ) {
        psError(PS_ERR_UNKNOWN,true, "psVectorImaginary returned a NULL?" );
        return 3;
    }
    if ( vec3->type.type != PS_TYPE_F32 ) {
        psError(PS_ERR_UNKNOWN,true, "psVectorImaginary returned a wrong type (%d)?",
                vec3->type.type );
        return 4;
    }

    vec7 = psVectorImaginary(vec7, vec6);
    if(vec7 == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psVectorImage returned NULL");
        return 12;
    }
    if(vec7->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"psVectorImaginary returned the wrong type");
        return 13;
    }

    vec11 = psVectorImaginary(vec11, vec10);
    if(vec11 == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psVectorImaginary returned NULL");
        return 14;
    }
    if(vec11->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN,true,"psVectorImaginary returned the wrong type");
        return 15;
    }

    // 3. compare results to the real/imaginary components of input
    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 ) {
            psError(PS_ERR_UNKNOWN,true, "psVectorReal didn't return the real portion at n=%d",
                    n );
            return 5;
        }
        if ( fabsf( vec3->data.F32[ n ] - i ) > FLT_EPSILON ) {
            psError(PS_ERR_UNKNOWN,true, "psVectorImaginary didn't return the real portion at n=%d",
                    n );
            return 6;
        }
        if ( fabsf( vec5->data.F32[n] - r) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN,true,"psVectorReal didn't return the real portion at n=%d",n);
            return 50;
        }
        if ( fabsf( vec7->data.F32[n] - 0) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN,true,"psVectorImaginary did not return the imaginary portion at n=%d",n);
            return 51;
        }
        if ( fabsf(vec9->data.F64[n] - rr) > FLT_EPSILON ) {
            psError(PS_ERR_UNKNOWN,true,"psVectorReal did not return the real portion at n=%d",n);
            return 52;
        }
        if ( fabsf(vec11->data.F64[n] - ii) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN,true,"psVectorImaginary did not return the imaginary portion at n=%d",n);
            return 53;
        }
    }

    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
    if ( psVectorReal(NULL,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with null input vector");
        return 63;
    }
    // Perform vector Imaginary with null input
    if ( psVectorImaginary(NULL,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with null input vector");
        return 73;
    }

    return 0;
}

psS32 testVectorComplex( void )
{
    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
    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 );
    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 );
    }

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

    // 3. verify that the result is a psC32
    if ( vec3->type.type != PS_TYPE_C32 ) {
        psError(PS_ERR_UNKNOWN,true, "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 ( 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 ) {
            psError(PS_ERR_UNKNOWN,true, "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(PS_ERR_UNKNOWN,true, "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(PS_ERR_UNKNOWN,true, "psVectorComplex returned a larger vector than the input supported (%d).", vec3->n );
        return 4;
    }

    // Verify the function works with psF64 type
    vec6 = psVectorComplex(vec6, vec4, vec5);
    if( vec6->type.type != PS_TYPE_C64 ) {
        psError(PS_ERR_UNKNOWN,true,"Vector return type is not complex");
        return 40;
    }

    // Verify error message generated with input of invalid type
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message.");
    vec4->type.type = PS_TYPE_S8;
    vec5->type.type = PS_TYPE_S8;
    vec6 = psVectorComplex(vec6, vec4, vec5);
    if(vec6 != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL for invalid type");
        return 50;
    }
    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
    if ( psVectorComplex(NULL,NULL,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with null input vector");
        return 62;
    }

    return 0;
}

psS32 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 ( psU32 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(PS_ERR_UNKNOWN,true, "the psVectorConjugate didn't return a C32 vector" );
        return 1;
    }

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

    psFree( vec );
    psFree( vec2 );

    // 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;
    }
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate a warning message.");
    vec2 = psVectorConjugate(vec2,vec);
    if(vec2->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"psVectorConjugate did not return a F32 vector");
        return 10;
    }
    for ( psU32 n = 0; n < 100; n++ ) {
        if( vec->data.F32[n] != vec2->data.F32[n] ) {
            psError(PS_ERR_UNKNOWN,true,"psVectorConjugate result is invalid (n=%d)",n);
            return 11;
        }
    }
    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);
    if(vec2->type.type != PS_TYPE_C64) {
        psError(PS_ERR_UNKNOWN,true,"psVectorConjugate did not return a C64 vector");
        return 12;
    }
    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 ) {
            psError(PS_ERR_UNKNOWN,true,"psVectorConjugate result is invalid (n=%d)",n);
            return 13;
        }
    }
    psFree(vec);

    // Perform vector conjugate with null input
    if ( psVectorConjugate(NULL,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with null input vector");
        return 61;
    }

    return 0;
}

psS32 testVectorPowerSpectrum( void )
{
    psVector * vec = NULL;
    psVector* vec2 = NULL;
    psVector* vec3 = NULL;
    psVector* vec4 = NULL;

    psF32 val;
    psF64 val1;

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

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

    // 3. verify result is psF32
    if ( vec2->type.type != PS_TYPE_F32 ) {
        psError(PS_ERR_UNKNOWN,true, "the type was not PS_TYPE_F32." );
        return 1;
    }
    if ( vec4->type.type != PS_TYPE_F64 ) {
        psError(PS_ERR_UNKNOWN,true,"psPowerSpectrum did not return type PS_TYPE_F64 type = %d",vec4->type.type);
        return 20;
    }

    // 3a. verify result is the same size a input
    // Awaiting IfA direction on bug #228
    //    if ( vec2->n != vec->n ) {
    //       psError(PS_ERR_UNKNOWN,true, "Output vector size different(%d) than input vector(%d)",vec2->n, vec->n);
    //       return 10;
    //    }

    // 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;
    val1= cabsf( vec3->data.C64[0] ) * cabsf(vec3->data.C64[0])/100/100;
    if ( fabsf( vec2->data.F32[ 0 ] - val ) > FLT_EPSILON ) {
        psError(PS_ERR_UNKNOWN,true, "psVectorPowerSpectrum result is invalid (n=0, %.2f %.2f)",
                vec2->data.F32[ 0 ], val );
        return 2;
    }
    if ( fabsf( vec4->data.C64[0] - val1 ) > FLT_EPSILON ) {
        psError(PS_ERR_UNKNOWN,true,"psVectorPowerSpectrum result is invalid (n=0)");
        return 21;
    }

    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 ) {
            psError(PS_ERR_UNKNOWN,true, "psVectorPowerSpectrum result is invalid (n=%d, %.2f %.2f)",
                    n, vec2->data.F32[ n ], val );
            return 2;
        }
        if (fabsf(val1 - vec4->data.F64[n]) > 10*FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN,true,"psVectorPowerSpectrum result is invalid (n=%d, %.2f %.2f)",n,vec4->data.F64[n],val1);
            return 22;
        }
    }

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

    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;
    }
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message.");
    if(psVectorPowerSpectrum(NULL,vec) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psVectorPowerSpectrum did not return a null vector.");
        return 10;
    }

    // Perform vector power spectrum with null input
    if ( psVectorPowerSpectrum(NULL,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with null input vector");
        return 60;
    }

    psFree(vec);

    return 0;
}
