Index: unk/psLib/src/dataManip/psFFT.c =================================================================== --- /trunk/psLib/src/dataManip/psFFT.c (revision 1627) +++ (revision ) @@ -1,350 +1,0 @@ -/** @file psVectorFFT.c - * - * @brief Contains FFT transform related functions for psVector - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-25 21:10:08 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#include -#include -#include -#include -#include - -#include "psVectorFFT.h" -#include "psError.h" -#include "psMemory.h" -#include "psLogMsg.h" -#include "psImageExtraction.h" - -#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE - -static bool p_fftwWisdomImported = false; - -psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction) -{ - unsigned int numElements; - psElemType type; - fftwf_plan plan; - - /* got good image data? */ - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - - if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) { - psError(__func__, "Input image must be a 32-bit float or complex image (type=%d)", type); - psFree(out); - return NULL; - } - - if ((type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE)) { - psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type); - psFree(out); - return NULL; - - } - - if ((type != PS_TYPE_F32) && (direction == PS_FFT_FORWARD)) { - psError(__func__, "Input image must be real image for forward FFT (type=%d).", type); - psFree(out); - return NULL; - } - - /* make sure the system-level wisdom information is imported. */ - if (!p_fftwWisdomImported) { - fftwf_import_system_wisdom(); - p_fftwWisdomImported = true; - } - - numElements = in->n; - - out = psVectorRecycle(out, numElements, PS_TYPE_C32); - out->n = numElements; - - if (type == PS_TYPE_F32) { - // need to convert to complex - psC32* outVec = out->data.C32; - psF32* inVec = in->data.F32; - - for (unsigned int i = 0; i < numElements; i++) { - outVec[i] = inVec[i]; - } - } else { - psC32* outVec = out->data.C32; - psC32* inVec = in->data.C32; - - for (unsigned int i = 0; i < numElements; i++) { - outVec[i] = inVec[i]; - } - } - - plan = fftwf_plan_dft_1d(numElements, - (fftwf_complex *) out->data.C32, - (fftwf_complex *) out->data.C32, direction, P_FFTW_PLAN_RIGOR); - - /* check if a plan exists now */ - if (plan == NULL) { - psError(__func__, "Failed to create FFTW plan."); - psFree(out); - return NULL; - } - - /* finally, call FFTW with the plan made above */ - fftwf_execute(plan); - - fftwf_destroy_plan(plan); - - return out; -} - -psVector* psVectorReal(psVector* out, const psVector* in) -{ - psElemType type; - unsigned int 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.V, in->data.V, 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 (unsigned int i = 0; i < numElements; i++) { - outVec[i] = crealf(inVec[i]); - } - } else { - psError(__func__, "Can not extract real component from given vector type (%d).", type); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorImaginary(psVector* out, const psVector* in) -{ - psElemType type; - unsigned int 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.V, 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 (unsigned int i = 0; i < numElements; i++) { - outVec[i] = cimagf(inVec[i]); - } - } else { - psError(__func__, "Can not extract imaginary component from given vector type (%d).", type); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorComplex(psVector* out, const psVector* real, const psVector* imag) -{ - psElemType type; - unsigned int 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) { - psError(__func__, "The inputs to psVectorComplex must be the same type."); - psFree(out); - return NULL; - } - - if (PS_IS_PSELEMTYPE_COMPLEX(type)) { - psError(__func__, "The inputs to psVectorComplex can not be complex."); - 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 (unsigned int i = 0; i < numElements; i++) { - outVec[i] = realVec[i] + I * imagVec[i]; - } - } else { - psError(__func__, "Can not merge real and imaginary portions for given vector type (%d).", type); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorConjugate(psVector* out, const psVector* in) -{ - psElemType type; - unsigned int 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.V, in->data.V, 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 (unsigned int i = 0; i < numElements; i++) { - outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]); - } - } else { - psError(__func__, "Can not compute complex conjugate for given vector type (%d).", type); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorPowerSpectrum(psVector* out, const psVector* in) -{ - psElemType type; - unsigned int outNumElements; - unsigned int inNumElements; - unsigned int inHalfNumElements; - unsigned int 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 not a complex number, this is not implemented */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - psError(__func__, "Power Spectrum for non-complex inputs is not implemented."); - psFree(out); - return NULL; - } - - 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 (unsigned int 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 { - psError(__func__, "Can not power spectrum for given vector type (%d).", type); - psFree(out); - return NULL; - } - - return out; - -} Index: unk/psLib/src/dataManip/psFFT.h =================================================================== --- /trunk/psLib/src/dataManip/psFFT.h (revision 1627) +++ (revision ) @@ -1,96 +1,0 @@ -/** @file psVectorFFT.h - * - * @brief Contains FFT transform related functions for psVector - * - * @ingroup Transform - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-25 21:10:08 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PS_VECTOR_FFT_H -#define PS_VECTOR_FFT_H - -#include "psVector.h" - -/// @addtogroup Transform -/// @{ - -/** Specify direction of FFT */ -typedef enum { - /// psImageFFT/psVectorFFT should perform a forward FFT. - PS_FFT_FORWARD = (-1), - - /// psImageFFT/psVectorFFT should perform a reverse FFT. - PS_FFT_REVERSE = (+1) -} psFftDirection; - - -/** Forward and reverse FFT calculations. - * - * This takes as input the vector of interest (in) and the direction - * (direction), which is specified by an enumerated type psFftDirection. - * The input vector may be of type psF32 or psC32, the result is always - * psC32. If the input vector is psF32, the direction must be forward. - * - * @return psVector* the FFT transformation result - */ -psVector* psVectorFFT( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in, ///< the vector to apply transform to - psFftDirection direction ///< the direction of the transform -); - -/** extract the real portion of a complex vector - * - * @return psVector* real portion of the input vector. - */ -psVector* psVectorReal( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to extract real portion from -); - -/** extract the imaginary portion of a complex vector - * - * @return psVector* imaginary portion of the input vector. - */ -psVector* psVectorImaginary( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to extract imaginary portion from -); - -/** creates a complex vector from separate real and imaginary vectors - * - * @return psVector* resulting complex vector - */ -psVector* psVectorComplex( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* real, ///< the real vector - const psVector* imag ///< the imaginary vector -); - -/** computes the complex conjugate of a vector - * - * @return psVector* the complex conjugate of the 'in' vector - */ -psVector* psVectorConjugate( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to compute conjugate of -); - -/** computes the power spectrum of a vector - * - * @return psVector* the power spectrum of the 'in' vector - */ -psVector* psVectorPowerSpectrum( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to power spectrum of -); - -/// @} - -#endif