Index: /trunk/psLib/test/dataManip/tst_psVectorFFT.c
===================================================================
--- /trunk/psLib/test/dataManip/tst_psVectorFFT.c	(revision 1009)
+++ /trunk/psLib/test/dataManip/tst_psVectorFFT.c	(revision 1009)
@@ -0,0 +1,380 @@
+/** @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: 2004-06-12 01:34:10 $
+ *
+ *  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-testVectorFFT",0},
+                              {testVectorRealImaginary,"601-testVectorRealImaginary",0},
+                              {testVectorComplex,"602-testVectorComplex",0},
+                              {testVectorConjugate,"603-testVectorConjugate",0},
+                              {testVectorPowerSpectrum,"604-testVectorPowerSpectrum",0},
+                              {NULL}
+                          };
+
+int main()
+{
+    psLogSetLevel(PS_LOG_INFO);
+
+    if (! runTestSuite(stderr,"psFFT",tests)) {
+        psAbort(__FILE__,"One or more tests failed");
+    }
+    return 0;
+}
+
+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;
+        }
+    }
+
+    psVectorFree(vec);
+    psVectorFree(vec2);
+    psVectorFree(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;
+        }
+    }
+
+    psVectorFree(vec);
+    psVectorFree(vec2);
+    psVectorFree(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,PS_TYPE_F64, 100);
+
+    // 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,PS_TYPE_F32,200);
+
+    // 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 vector though input sizes mismatched.");
+        return 4;
+    }
+
+    psVectorFree(vec);
+    psVectorFree(vec2);
+    psVectorFree(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;
+        };
+    }
+
+    psVectorFree(vec);
+    psVectorFree(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;
+    };
+
+    psVectorFree(vec);
+    psVectorFree(vec2);
+
+    return 0;
+}