Index: trunk/psLib/src/fft/psVectorFFT.c =================================================================== --- trunk/psLib/src/fft/psVectorFFT.c (revision 11668) +++ trunk/psLib/src/fft/psVectorFFT.c (revision 11703) @@ -3,8 +3,9 @@ * @brief Contains FFT transform related functions for psVector * + * @author Paul Price, IfA * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-02-06 21:36:09 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-02-08 04:17:58 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -12,5 +13,5 @@ #ifdef HAVE_CONFIG_H -# include "config.h" +#include "config.h" #endif @@ -21,401 +22,161 @@ #include -#include "psVectorFFT.h" +#include "psAssert.h" #include "psError.h" #include "psMemory.h" #include "psLogMsg.h" +#include "psConstants.h" +#include "psVectorFFT.h" + +#define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is + +static psBool fftwWisdomImported = false; // Has the system wisdom been imported? + +bool psVectorForwardFFT(psVector **real, psVector **imag, const psVector *in) +{ + PS_ASSERT_VECTOR_NON_NULL(in, false); + PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, false); + PS_ASSERT_PTR_NON_NULL(real, false); + PS_ASSERT_PTR_NON_NULL(imag, false); + + // Make sure the system-level wisdom information is imported. + if (!fftwWisdomImported) { + fftwf_import_system_wisdom(); + fftwWisdomImported = true; + } + + long num = in->n; // Number of elements + + // Do the FFT + fftwf_complex *out = fftwf_malloc(num * sizeof(fftwf_complex)); // Output data + fftwf_plan plan = fftwf_plan_dft_r2c_1d(num, in->data.F32, out, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); + + // Pull the real and imaginary parts out + *real = psVectorRecycle(*real, num/2 + 1, PS_TYPE_F32); + *imag = psVectorRecycle(*imag, num/2 + 1, PS_TYPE_F32); + for (int i = 0; i < num; i++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + (*real)->data.F32[i] = creal(out[i]); + (*imag)->data.F32[i] = cimag(out[i]); +#else + // FFTW's backup complex support + (*real)->data.F32[i] = out[i][0]; + (*imag)->data.F32[i] = out[i][1]; +#endif + } + fftwf_free(out); + + return true; +} + +bool psVectorBackwardFFT(psVector **out, const psVector *real, const psVector *imag) +{ + PS_ASSERT_VECTOR_NON_NULL(real, false); + PS_ASSERT_VECTOR_TYPE(real, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_NON_NULL(imag, false); + PS_ASSERT_VECTOR_TYPE(imag, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(real, imag, false); + PS_ASSERT_PTR_NON_NULL(out, false); + + // Make sure the system-level wisdom information is imported. + if (!fftwWisdomImported) { + fftwf_import_system_wisdom(); + fftwWisdomImported = true; + } + + long num = real->n; // Number of elements + + // Stuff the real and imaginary parts in + fftwf_complex *in = fftwf_malloc(num * sizeof(fftwf_complex)); // Input data + for (int i = 0; i < num; i++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + in[i] = real->data.F32[i] + imag->data.F32[i] * I; +#else + // FFTW's backup complex support + in[i][0] = real->data.F32[i]; + in[i][1] = imag->data.F32[i]; +#endif + } + // Do the FFT + *out = psVectorRecycle(*out, 2 * num, PS_TYPE_F32); + fftwf_plan plan = fftwf_plan_dft_c2r_1d(2 * num, in, (*out)->data.F32, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); -#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE + fftwf_free(in); -static bool p_fftwWisdomImported = false; + return true; +} -psVector* psVectorFFT(psVector* out, const psVector* in, psFFTFlags direction) +psVector *psVectorPowerSpectrum(psVector* out, const psVector* in) { - psU32 numElements; - psElemType type; - fftwf_plan plan; + PS_ASSERT_VECTOR_NON_NULL(in, NULL); + PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, NULL); - /* got good image data? */ - if (in == NULL) { - psFree(out); + psVector *real = NULL; // Real component of FFT + psVector *imag = NULL; // Imaginary component of FFT + + if (!psVectorForwardFFT(&real, &imag, in)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT."); return NULL; } - type = in->type.type; + int num = real->n; // Number of elements - /* make sure the system-level wisdom information is imported. */ - if (!p_fftwWisdomImported) { - fftwf_import_system_wisdom(); - p_fftwWisdomImported = true; + float norm = 1.0 / PS_SQR(in->n); + for (int i = 0; i < num; i++) { + // Power spectrum is the square of the complex modulus + real->data.F32[i] = norm * (PS_SQR(real->data.F32[i]) + PS_SQR(imag->data.F32[i])); + } + psFree(imag); + + return real; +} + +bool psVectorComplexMultiply(psVector **outReal, psVector **outImag, + const psVector *in1Real, const psVector *in1Imag, + const psVector *in2Real, const psVector *in2Imag) +{ + PS_ASSERT_VECTOR_NON_NULL(in1Real, false); + PS_ASSERT_VECTOR_NON_NULL(in1Imag, false); + PS_ASSERT_VECTOR_NON_NULL(in2Real, false); + PS_ASSERT_VECTOR_NON_NULL(in2Imag, false); + PS_ASSERT_VECTOR_TYPE(in1Real, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(in1Imag, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(in2Real, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(in2Imag, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(in1Imag, in1Real, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(in2Real, in1Real, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(in2Imag, in1Real, false); + PS_ASSERT_PTR_NON_NULL(outReal, false); + PS_ASSERT_PTR_NON_NULL(outImag, false); + if (*outReal) { + PS_ASSERT_VECTOR_NON_NULL(*outReal, false); + PS_ASSERT_VECTOR_NON_NULL(*outImag, false); + } else if (*outImag) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "If the output real part is provided, the output imaginary part must also be provided."); + return false; } - numElements = in->n; + int num = in1Real->n; // Number of elements - out = psVectorCopy(out, in, PS_TYPE_C32); - out->n = numElements; + *outReal = psVectorRecycle(*outReal, num, PS_TYPE_F32); + *outImag = psVectorRecycle(*outImag, num, PS_TYPE_F32); - if ((direction & PS_FFT_FORWARD) != 0) { - plan = fftwf_plan_dft_1d(numElements, - (fftwf_complex *) out->data.C32, - (fftwf_complex *) out->data.C32, FFTW_FORWARD, P_FFTW_PLAN_RIGOR); - } else if ((direction & PS_FFT_REVERSE) != 0) { - plan = fftwf_plan_dft_1d(numElements, - (fftwf_complex *) out->data.C32, - (fftwf_complex *) out->data.C32, FFTW_BACKWARD, P_FFTW_PLAN_RIGOR); - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Must specify the direction as either PS_FFT_FORWARD or PS_FFT_REVERSE.")); - psFree(out); - return NULL; + for (int i = 0; i < num; i++) { + // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i + (*outReal)->data.F32[i] = in1Real->data.F32[i] * in2Real->data.F32[i] - + in1Imag->data.F32[i] * in2Imag->data.F32[i]; + (*outImag)->data.F32[i] = in1Imag->data.F32[i] * in2Real->data.F32[i] + + in1Real->data.F32[i] * in2Imag->data.F32[i]; } - /* check if a plan exists now */ - if (plan == NULL) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Could not create a valid FFT plan to perform the transform.")); - psFree(out); - return NULL; - } - - /* finally, call FFTW with the plan made above */ - fftwf_execute(plan); - - fftwf_destroy_plan(plan); - - if ((direction & PS_FFT_REAL_RESULT) != 0) { - for (psS32 i = 0; i < numElements; i++) { - out->data.F32[i] = out->data.C32[i]; - } - out->type.type = PS_TYPE_F32; - out->data.U8 = psRealloc(out->data.U8,PSELEMTYPE_SIZEOF(PS_TYPE_F32)*out->nalloc); - } - - return out; + return true; } - -psVector* psVectorReal(psVector* out, const psVector* in) -{ - psElemType type; - psU32 numElements; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numElements = in->n; - - /* if not a complex number, this is logically just a copy */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - // Warn user, as this is probably not expected - psLogMsg(__func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. " - "Just a vector copy was performed."); - out = psVectorRecycle(out, numElements, type); - out->n = numElements; - memcpy(out->data.U8, in->data.U8, numElements * PSELEMTYPE_SIZEOF(type)); - return out; - } - - if (type == PS_TYPE_C32) { - psF32* outVec; - psC32* inVec = in->data.C32; - - out = psVectorRecycle(out, numElements, PS_TYPE_F32); - out->n = numElements; - outVec = out->data.F32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = crealf(inVec[i]); - } - } else if (type == PS_TYPE_C64) { - psF64* outVec; - psC64* inVec = in->data.C64; - - out = psVectorRecycle(out, numElements, PS_TYPE_F64); - out->n = numElements; - outVec = out->data.F64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = creal(inVec[i]); - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psVector type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorImaginary(psVector* out, const psVector* in) -{ - psElemType type; - psU32 numElements; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numElements = in->n; - - /* if not a complex number, this is logically just zeroed image of same size */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - // Warn user, as this is probably not expected - psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. " - "A zeroed vector was returned."); - out = psVectorRecycle(out, numElements, type); - out->n = numElements; - memset(out->data.U8, 0, PSELEMTYPE_SIZEOF(type) * numElements); - return out; - } - - if (type == PS_TYPE_C32) { - psF32* outVec; - psC32* inVec = in->data.C32; - - out = psVectorRecycle(out, numElements, PS_TYPE_F32); - out->n = numElements; - outVec = out->data.F32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = cimagf(inVec[i]); - } - } else if (type == PS_TYPE_C64) { - psF64* outVec; - psC64* inVec = in->data.C64; - - out = psVectorRecycle(out, numElements, PS_TYPE_F64); - out->n = numElements; - outVec = out->data.F64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = cimag(inVec[i]); - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psVector type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorComplex(psVector* out, const psVector* real, const psVector* imag) -{ - psElemType type; - psU32 numElements; - - if (real == NULL || imag == NULL) { - psFree(out); - return NULL; - } - - type = real->type.type; - if (real->n < imag->n) { - numElements = real->n; - } else { - numElements = imag->n; - } - - if (imag->type.type != type) { - char* typeStrReal; - char* typeStrImag; - PS_TYPE_NAME(typeStrReal,type); - PS_TYPE_NAME(typeStrImag,imag->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Real psVector type, %s, and imaginary psVector type, %s, must be the same."), - typeStrReal,typeStrImag); - psFree(out); - return NULL; - } - - if (type == PS_TYPE_F32) { - psC32* outVec; - psF32* realVec = real->data.F32; - psF32* imagVec = imag->data.F32; - - out = psVectorRecycle(out, numElements, PS_TYPE_C32); - out->n = numElements; - outVec = out->data.C32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = realVec[i] + I * imagVec[i]; - } - } else if (type == PS_TYPE_F64) { - psC64* outVec; - psF64* realVec = real->data.F64; - psF64* imagVec = imag->data.F64; - - out = psVectorRecycle(out, numElements, PS_TYPE_C64); - out->n = numElements; - outVec = out->data.C64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = realVec[i] + I * imagVec[i]; - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psVector type, %s, is required to be either psF32 or psF64."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorConjugate(psVector* out, const psVector* in) -{ - psElemType type; - psU32 numElements; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numElements = in->n; - - /* if not a complex number, this is logically just a image copy */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - // Warn user, as this is probably not expected - psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. " - "Vector copy was performed instead."); - - out = psVectorRecycle(out, numElements, type); - out->n = numElements; - memcpy(out->data.U8, in->data.U8, PSELEMTYPE_SIZEOF(type) * numElements); - return out; - } - - if (type == PS_TYPE_C32) { - psC32* outVec; - psC32* inVec = in->data.C32; - - out = psVectorRecycle(out, numElements, PS_TYPE_C32); - out->n = numElements; - outVec = out->data.C32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]); - } - } else if (type == PS_TYPE_C64) { - psC64* outVec; - psC64* inVec = in->data.C64; - - out = psVectorRecycle(out, numElements, PS_TYPE_C64); - out->n = numElements; - outVec = out->data.C64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = creal(inVec[i]) - I * cimag(inVec[i]); - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psVector type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorPowerSpectrum(psVector* out, const psVector* in) -{ - psElemType type; - psU32 outNumElements; - psU32 inNumElements; - psU32 inHalfNumElements; - psU32 inNumElementsSquared; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - inNumElements = in->n; - inNumElementsSquared = inNumElements * inNumElements; - inHalfNumElements = inNumElements / 2; - outNumElements = inHalfNumElements + 1; - - if (type == PS_TYPE_C32) { - psF32* outVec; - psC32* inVec = in->data.C32; - psF32 inAbs1; - psF32 inAbs2; - - out = psVectorRecycle(out, outNumElements, PS_TYPE_F32); - out->n = outNumElements; - outVec = out->data.F32; - - // from ADD: P_0 = |C_0|^2/N^2 - inAbs1 = cabsf(inVec[0]); - outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared; - - // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1) - for (psU32 i = 1; i < inHalfNumElements; i++) { - inAbs1 = cabsf(inVec[i]); - inAbs2 = cabsf(inVec[inNumElements - i]); - outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared; - } - - // from ADD: P_N/2 = |C_N/2|^2/N^2 - inAbs1 = cabsf(inVec[inHalfNumElements]); - outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared; - } else if (type == PS_TYPE_C64) { - psF64* outVec; - psC64* inVec = in->data.C64; - psF64 inAbs1; - psF64 inAbs2; - - out = psVectorRecycle(out, outNumElements, PS_TYPE_F64); - out->n = outNumElements; - outVec = out->data.F64; - - // from ADD: P_0 = |C_0|^2/N^2 - inAbs1 = cabs(inVec[0]); - outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared; - - // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1) - for (psU32 i = 1; i < inHalfNumElements; i++) { - inAbs1 = cabs(inVec[i]); - inAbs2 = cabs(inVec[inNumElements - i]); - outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared; - } - - // from ADD: P_N/2 = |C_N/2|^2/N^2 - inAbs1 = cabs(inVec[inHalfNumElements]); - outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared; - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psVector type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; - -}