Index: trunk/psLib/src/fft/psImageFFT.c =================================================================== --- trunk/psLib/src/fft/psImageFFT.c (revision 11668) +++ trunk/psLib/src/fft/psImageFFT.c (revision 11703) @@ -1,13 +1,14 @@ -/** @file psImageFFT.c - * - * @brief Contains FFT transform related functions for psImage. - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-02-06 21:36:09 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psImageFFT.c +/// +/// @brief Contains FFT transform related functions for psImage. +/// +/// @author Paul Price, IfA +/// @author Robert DeSonia, MHPCC +/// +/// @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// +/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +/// #ifdef HAVE_CONFIG_H @@ -20,442 +21,254 @@ #include -#include "psImageFFT.h" +#include "psAssert.h" #include "psError.h" #include "psMemory.h" #include "psLogMsg.h" +#include "psConstants.h" #include "psImageStructManip.h" - - - -#define PS_FFTW_PLAN_RIGOR FFTW_ESTIMATE - -static bool p_fftwWisdomImported = false; - -psImage* psImageFFT(psImage* out, const psImage* image, psFFTFlags direction) -{ - psU32 numCols; - psU32 numRows; - psElemType type; - fftwf_plan plan; - - /* got good image data? */ - if (image == NULL) { - psFree(out); +#include "psImageFFT.h" + + +#define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is + +static psBool fftwWisdomImported = false; // Has the FFTW wisdom been imported yet? + +bool psImageForwardFFT(psImage **real, psImage **imag, const psImage *in) +{ + PS_ASSERT_IMAGE_NON_NULL(in, false); + PS_ASSERT_IMAGE_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; + } + + int numCols = in->numCols; // Number of columns + int numRows = in->numRows; // Number of rows + + psF32 *input; // Input data, for FFTW + if (!in->parent) { + // No parent --- can just use the data buffer + input = psMemIncrRefCounter(in->p_rawDataBuffer); + } else { + // Need to copy the data + input = psAlloc(numCols * numRows); + for (int y = 0; y < numRows; y++) { + memcpy(&input[y * numRows], in->data.F32[y], numCols); + } + } + + // Do the FFT + fftwf_complex *out = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // Output data + fftwf_plan plan = fftwf_plan_dft_r2c_2d(numRows, numCols, input, out, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); + psFree(input); + + // Pull the real and imaginary parts out + numCols = numCols/2 + 1; // x dimension is halved by FFTW + *real = psImageRecycle(*real, numCols, numRows, PS_TYPE_F32); + *imag = psImageRecycle(*imag, numCols, numRows, PS_TYPE_F32); + for (int y = 0, index = 0; y < numRows; y++) { + for (int x = 0; x < numCols; index++, x++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + (*real)->data.F32[y][x] = creal(out[index]); + (*imag)->data.F32[y][x] = cimag(out[index]); +#else + // FFTW's backup complex support + (*real)->data.F32[y][x] = out[index][0]; + (*imag)->data.F32[y][x] = out[index][1]; +#endif + } + } + fftwf_free(out); + + return true; +} + +bool psImageBackwardFFT(psImage **out, const psImage *real, const psImage *imag, int origCols) +{ + PS_ASSERT_IMAGE_NON_NULL(real, false); + PS_ASSERT_IMAGE_TYPE(real, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(imag, false); + PS_ASSERT_IMAGE_TYPE(imag, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(real, imag, false); + PS_ASSERT_PTR_NON_NULL(out, false); + + int numCols = real->numCols; // Number of columns + int numRows = real->numRows; // Number of rows + + // Because of the way FFT r2c and c2r work, need the number of columns in the target to be: + // 2 * numCols = 2*(numCols/2 + 1) = origCols % 2 ? origCols + 1 : origCols + 2 + if (2 * numCols != ((origCols % 2) ? origCols + 1 : origCols + 2)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "Number of columns in FFT-ed images (%d) should be origCols (%d) / 2 + 1 = %d", + numCols, origCols, origCols/2 + 1); + return false; + } + + if (*out && (*out)->parent) { + // It has a parent, so we can't write directly into the buffer. + // It had better be the correct size and type, because we don't want to resize a child image + PS_ASSERT_IMAGE_SIZE(*out, origCols, numRows, false); + PS_ASSERT_IMAGE_TYPE(*out, PS_TYPE_F32, false); + } + + + // Make sure the system-level wisdom information is imported. + if (!fftwWisdomImported) { + fftwf_import_system_wisdom(); + fftwWisdomImported = true; + } + + // Stuff the real and imaginary parts in + psF32 *target = psAlloc(2 * numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Target for FFTW + fftwf_complex *in = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // Input data + for (int y = 0, index = 0; y < numRows; y++) { + for (int x = 0; x < numCols; index++, x++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + in[index] = real->data.F32[y][x] + imag->data.F32[y][x] * I; +#else + // FFTW's backup complex support + in[index][0] = real->data.F32[y][x]; + in[index][1] = imag->data.F32[y][x]; +#endif + } + } + + // Do the FFT + fftwf_plan plan = fftwf_plan_dft_c2r_2d(numRows, origCols, in, target, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); + fftwf_free(in); + + if (!(*out) || !(*out)->parent) { + *out = psImageRecycle(*out, origCols, numRows, PS_TYPE_F32); + } + + // Copy the target pixels into the output + // There's a slight offset --- target has the additional padding required by FFTW, + // while the output doesn't. + for (int y = 0, index = 0; y < numRows; y++, index += origCols) { + memcpy(&(*out)->data.F32[y][0], &target[index], numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + } + psFree(target); + + return true; +} + +psImage* psImagePowerSpectrum(psImage *out, const psImage *in) +{ + PS_ASSERT_IMAGE_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL); + + psImage *real = NULL; // Real component of FFT + psImage *imag = NULL; // Imaginary component of FFT + + if (!psImageForwardFFT(&real, &imag, in)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT."); return NULL; } - if ( ((direction & PS_FFT_FORWARD) != 0) ) { - if ((direction & PS_FFT_REVERSE) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options.")); - psFree(out); - return NULL; - } - if ((direction & PS_FFT_REAL_RESULT) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported.")); - psFree(out); - return NULL; - } - } else if ((direction & PS_FFT_REVERSE) == 0) { + int numCols = real->numCols; // Number of columns + int numRows = real->numRows; // Number of rows + + float norm = 1.0 / PS_SQR(in->numCols) / PS_SQR(in->numRows); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // Power spectrum is the square of the complex modulus + real->data.F32[y][x] = norm * (PS_SQR(real->data.F32[y][x]) + PS_SQR(imag->data.F32[y][x])); + } + } + psFree(imag); + + return real; +} + + + +bool psImageComplexMultiply(psImage **outReal, psImage **outImag, + const psImage *in1Real, const psImage *in1Imag, + const psImage *in2Real, const psImage *in2Imag) +{ + PS_ASSERT_IMAGE_NON_NULL(in1Real, false); + PS_ASSERT_IMAGE_NON_NULL(in1Imag, false); + PS_ASSERT_IMAGE_NON_NULL(in2Real, false); + PS_ASSERT_IMAGE_NON_NULL(in2Imag, false); + PS_ASSERT_IMAGE_TYPE(in1Real, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(in1Imag, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(in2Real, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(in2Imag, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(in1Imag, in1Real, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(in2Real, in1Real, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(in2Imag, in1Real, false); + PS_ASSERT_PTR_NON_NULL(outReal, false); + PS_ASSERT_PTR_NON_NULL(outImag, false); + if (*outReal) { + PS_ASSERT_IMAGE_NON_NULL(*outReal, false); + PS_ASSERT_IMAGE_NON_NULL(*outImag, false); + } else if (*outImag) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option.")); - psFree(out); - return NULL; - } - - type = image->type.type; - - /* make sure the system-level wisdom information is imported. */ - if (!p_fftwWisdomImported) { - fftwf_import_system_wisdom(); - p_fftwWisdomImported = true; - } - - numRows = image->numRows; - numCols = image->numCols; - - // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. - psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; - - fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex)); - p_psImageCopyToRawBuffer(outBuffer, image, PS_TYPE_C32); - - plan = fftwf_plan_dft_2d(numRows, numCols, - outBuffer, - outBuffer, - sign, - PS_FFTW_PLAN_RIGOR); - - /* check if a plan exists now -- if not, it is a real problem at this point */ - 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) { - // n.b., we do this instead of using fftwf_plan_dft_c2r because that - // plan requires a half-image, which would require a image reordering - // that is not as simple as performing a normal complex transform and - // then taking the real part of the result. If performance here - // becomes an issue, the use of fftwf_plan_dft_r2c should be considered - // as well as fftwf_plan_dft_c2r. - out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); - int index = 0; - for (int row=0; row < numRows; row++) { - psF32* outRow = out->data.F32[row]; - for (int col=0; col < numCols; col++) { - outRow[col] = crealf(outBuffer[index++]); // take just the real part - } - } + "If the output real part is provided, the output imaginary part must also be provided."); + return false; + } + + int numRows = in1Real->numRows; // Number of rows + int numCols = in1Real->numCols; // Number of columns + + // Need to worry if the outputs are children + psImage *targetReal, *targetImag; // Target real and imaginary parts + + if ((*outReal)->parent) { + // It had better be the correct size and type, because we don't want to resize a child image + PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false); + targetReal = psImageAlloc(numCols, numRows, PS_TYPE_F32); } else { - out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); - int index = 0; - for (int row=0; row < numRows; row++) { - psC32* outRow = out->data.C32[row]; - for (int col=0; col < numCols; col++) { - outRow[col] = outBuffer[index++]; // take just the real part - } - } - // memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex)); - } - - fftwf_free(outBuffer); - - return out; - -} - -psImage* psImageReal(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - /* if not a complex number, this is logically just a copy then */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - return psImageCopy(out, in, type); - } - - if (type == PS_TYPE_C32) { - psF32* outRow; - psC32* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F32[row]; - inRow = in->data.C32[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = crealf(inRow[col]); - } - } - } else if (type == PS_TYPE_C64) { - psF64* outRow; - psC64* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = creal(inRow[col]); - } - } + *outReal = psImageRecycle(*outReal, numCols, numRows, PS_TYPE_F32); + targetReal = psMemIncrRefCounter(*outReal); + } + if ((*outImag)->parent) { + // It had better be the correct size and type, because we don't want to resize a child image + if ((*outReal)->numCols != numCols || (*outReal)->numRows != numRows || + (*outImag)->type.type != PS_TYPE_F32) { + // Plug potential memory leak --- need to free targetReal if there's a problem + psFree(targetReal); + } + PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false); + targetImag = psImageAlloc(numCols, numRows, PS_TYPE_F32); } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - - psFree(out); - return NULL; - } - - return out; -} - -psImage* psImageImaginary(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - /* if not a complex image type, this is logically just zeroed image of same size */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - out = psImageRecycle(out, numCols, numRows, type); - memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows); - return out; - } - - if (type == PS_TYPE_C32) { - psF32* outRow; - psC32* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F32[row]; - inRow = in->data.C32[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = cimagf(inRow[col]); - } - } - } else if (type == PS_TYPE_C64) { - psF64* outRow; - psC64* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = cimag(inRow[col]); - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (real == NULL || imag == NULL) { - psFree(out); - return NULL; - } - - type = real->type.type; - numCols = real->numCols; - numRows = real->numRows; - - 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 psImage type (%s) and imaginary psImage type (%s) must be the same."), - typeStrReal,typeStrImag); - psFree(out); - return NULL; - } - - if (imag->numCols != numCols || imag->numRows != numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Real psImage size (%dx%d) and imaginary psImage size (%dx%d) must be the same."), - numCols, numRows, imag->numCols, imag->numRows); - psFree(out); - return NULL; - } - - if (type == PS_TYPE_F32) { - psC32* outRow; - psF32* realRow; - psF32* imagRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); - - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C32[row]; - realRow = real->data.F32[row]; - imagRow = imag->data.F32[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = realRow[col] + I * imagRow[col]; - } - } - } else if (type == PS_TYPE_F64) { - psC64* outRow; - psF64* realRow; - psF64* imagRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C64[row]; - realRow = real->data.F64[row]; - imagRow = imag->data.F64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = realRow[col] + I * imagRow[col]; - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psImage type, %s, is required to be either psF32 or psF64."), - typeStr); - psFree(out); - return NULL; - } - - - return out; -} - -psImage* psImageConjugate(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - /* if not a complex image, this is logically just a image copy */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - return psImageCopy(out, in, type); - } - - if (type == PS_TYPE_C32) { - psC32* outRow; - psC32* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C32[row]; - inRow = in->data.C32[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]); - } - } - } else if (type == PS_TYPE_C64) { - psC64* outRow; - psC64* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]); - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psImage type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psImage* psImagePowerSpectrum(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - if (type == PS_TYPE_C32) { - psF32* outRow; - psC32* inRow; - psF32 real; - psF32 imag; - psF32 fNumCols = numCols; - psF32 fNumRows = numRows; - psF32 numElementsSquared = fNumCols * fNumCols * fNumRows * fNumRows; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F32[row]; - inRow = in->data.C32[row]; - - for (psU32 col = 0; col < numCols; col++) { - real = crealf(inRow[col]); - imag = cimagf(inRow[col]); - outRow[col] = (real * real + imag * imag) / numElementsSquared; - } - } - } else if (type == PS_TYPE_C64) { - psF64* outRow; - psC64* inRow; - psF64 real; - psF64 imag; - psF64 numElementsSquared = numCols * numCols * numRows * numRows; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - real = creal(inRow[col]); - imag = cimag(inRow[col]); - outRow[col] = (real * real + imag * imag) / numElementsSquared; - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psImage type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; - -} + *outImag = psImageRecycle(*outImag, numCols, numRows, PS_TYPE_F32); + targetImag = psMemIncrRefCounter(*outImag); + } + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i + float real = in1Real->data.F32[y][x] * in2Real->data.F32[y][x] - + in1Imag->data.F32[y][x] * in2Imag->data.F32[y][x]; + float imag = in1Imag->data.F32[y][x] * in2Real->data.F32[y][x] + + in1Real->data.F32[y][x] * in2Imag->data.F32[y][x]; + targetReal->data.F32[y][x] = real; + targetImag->data.F32[y][x] = imag; + } + } + + if ((*outReal)->parent) { + *outReal = psImageCopy(*outReal, targetReal, PS_TYPE_F32); + } + if ((*outImag)->parent) { + *outImag = psImageCopy(*outImag, targetImag, PS_TYPE_F32); + } + + psFree(targetReal); + psFree(targetImag); + + return true; +} Index: trunk/psLib/src/fft/psImageFFT.h =================================================================== --- trunk/psLib/src/fft/psImageFFT.h (revision 11668) +++ trunk/psLib/src/fft/psImageFFT.h (revision 11703) @@ -1,11 +1,12 @@ -/* @file psImageFFT.h - * @brief Contains FFT transform related functions for psImage - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-23 22:47:22 $ - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psImageFFT.h +/// @brief Contains FFT transform related functions for psImage +/// +/// @author Paul Price, IfA +/// @author Robert DeSonia, MHPCC +/// +/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +/// #ifndef PS_IMAGE_FFT_H @@ -13,69 +14,50 @@ #include "psImage.h" -#include "psVectorFFT.h" // for psFFTFlags /// @addtogroup MathOps Mathematical Operations /// @{ -/** Forward and reverse FFT calculations. - * - * This takes as input the image of interest (in) and the direction - * (direction), which is specified by an enumerated type psFftDirection. - * The input image may be of type psF32 or psC32, the result is always - * psC32. If the input vector is psF32, the direction must be forward. - * - * @return psImage* the FFT transformation result - */ -psImage* psImageFFT( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* image, ///< the psImage to apply transform to - psFFTFlags direction ///< the direction of the transform -); +/// Forward FFT of an image +/// +/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts. The FFT +/// is not normalised (a forward followed by a reverse is the original scaled by the total number of pixels). +/// The FFT takes advantage of the fact that the input is purely real; hence the output number of columns is +/// numCols/2 + 1 (with division rounding down). Only implemented for F32 input. +bool psImageForwardFFT(psImage **real, ///< Real part of FFT + psImage **imag, ///< Imaginary part of FFT + const psImage *in///< Input image (F32) + ); -/** extract the real portion of a complex image - * - * @return psImage* real portion of the input image. - */ -psImage* psImageReal( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract real portion from -); +/// Backward FFT of an image +/// +/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result +/// returned. The FFT is not normalised (a forward followed by a reverse is the original scaled by the total +/// number of pixels). The FFT takes advantage of the fact that the output will be purely real; hence the +/// input number of columns is numCols/2 + 1 (with division rounding down); to manage the redundancy (is the +/// original number of columns even or odd?), we need the original number of columns (the number of columns of +/// the image that was input to psImageForwardFFT) to be provided. Only implemented for F32 input. +bool psImageBackwardFFT(psImage **out,///< Output image + const psImage *real, ///< Real input (F32) + const psImage *imag, ///< Imaginary input (F32) + int origCols ///< Original number of columns + ); -/** extract the imaginary portion of a complex image - * - * @return psImage* imaginary portion of the input image. - */ -psImage* psImageImaginary( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract imaginary portion from -); -/** creates a complex image from separate real and imaginary plane images - * - * @return psImage* resulting complex image - */ -psImage* psImageComplex( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* real, ///< the real plane image - const psImage* imag ///< the imaginary plane image -); +/// Power spectrum of an image +/// +/// Generates the power spectrum of an image. Only implemented for F32 input. +psImage* psImagePowerSpectrum(psImage *out, const psImage *in); -/** computes the complex conjugate of an image - * - * @return psImage* the complex conjugate of the 'in' image - */ -psImage* psImageConjugate( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to compute conjugate of -); - -/** computes the power spectrum of an image - * - * @return psImage* the power spectrum of the 'in' image - */ -psImage* psImagePowerSpectrum( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to power spectrum of -); +/// Multiply complex images +/// +/// The input images are the real and imaginary parts of each of two images. The real and imaginary parts +/// of the output image are returned. Only implemented for F32 input. +bool psImageComplexMultiply(psImage **outReal, ///< Real part of output + psimage **outImag, ///< Imaginary part of output + const psImage *in1Real, ///< Real part of input 1 + const psImage *in1Imag, ///< Imaginary part of input 1 + const psImage *in2Real, ///< Real part of input 2 + const psImage *in2Imag ///< Imaginary part of input 2 + ); /// @} 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; - -} Index: trunk/psLib/src/fft/psVectorFFT.h =================================================================== --- trunk/psLib/src/fft/psVectorFFT.h (revision 11668) +++ trunk/psLib/src/fft/psVectorFFT.h (revision 11703) @@ -1,11 +1,12 @@ -/* @file psVectorFFT.h - * @brief Contains FFT transform related functions for psVector - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-23 22:47:22 $ - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psVectorFFT.h +/// @brief Contains FFT transform related functions for psVector +/// +/// @author Paul Price, IfA +/// @author Robert DeSonia, MHPCC +/// +/// @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +/// #ifndef PS_VECTOR_FFT_H @@ -17,77 +18,47 @@ /// @{ -/** Specify direction of FFT */ -typedef enum { - /// psImageFFT/psVectorFFT should perform a forward FFT. - PS_FFT_FORWARD = 1, +/// Forward FFT of a vector +/// +/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts. The FFT +/// is not normalised (a forward followed by a reverse is the original scaled by the number of elements). The +/// FFT takes advantage of the fact that the input is purely real; hence the output size is N/2 + 1 (with +/// division rounding down). Only implemented for F32 input. +bool psVectorForwardFFT(psVector **real,///< Real part of FFT + psVector **imag,///< Imaginary part of FFT + const psVector *in ///< Input vector (F32) + ); - /// psImageFFT/psVectorFFT should perform a reverse FFT. - PS_FFT_REVERSE = 2, +/// Backward FFT of a vector +/// +/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result +/// returned. The FFT is not normalised (a forward followed by a reverse is the original scaled by the number +/// of elements). The FFT takes advantage of the fact that the output will be purely real; hence the input +/// size is N/2 + 1 (with division rounding down); to manage the redundancy (is the original size even or +/// odd?), we need the original size (the size of the array that was input to psVectorForwardFFT) to be +/// provided. Only implemented for F32 input. +bool psVectorBackwardFFT(psVector **out,///< Output vector + const psVector *real, ///< Real input (F32) + const psVector *imag, ///< Imaginary input (F32) + int origN ///< Original number of elements + ); - /// psImageFFT/psVectorFFT should perform a reverse FFT with a real result. - PS_FFT_REAL_RESULT = 4 -} psFFTFlags; +/// Power spectrum of a vector +/// +/// Generates the power spectrum of a vector. Only implemented for F32 input. +psVector *psVectorPowerSpectrum(psVector *out, ///< Output power spectrum, or NULL + const psVector* in ///< Input vector (F32) + ); - -/** 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 - psFFTFlags 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 -); +/// Multiply complex vectors +/// +/// The input vectors are the real and imaginary parts of each of two vectors. The real and imaginary parts +/// of the output vector are returned. Only implemented for F32 input. +bool psVectorComplexMultiply(psVector **outReal, ///< Real part of output + psVector **outImag, ///< Imaginary part of output + const psVector *in1Real, ///< Real part of input 1 + const psVector *in1Imag, ///< Imaginary part of input 1 + const psVector *in2Real, ///< Real part of input 2 + const psVector *in2Imag ///< Imaginary part of input 2 + ); /// @}