Index: trunk/psLib/src/fft/psVectorFFT.c =================================================================== --- trunk/psLib/src/fft/psVectorFFT.c (revision 824) +++ trunk/psLib/src/fft/psVectorFFT.c (revision 827) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-06-01 22:42:57 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-02 03:02:48 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,14 +26,6 @@ static bool p_fftwWisdomImported = false; -psImage* psImageFFT(psImage* out, const psImage* in, int flags) -{ - static fftwf_plan lastForwardPlan = NULL; - static int lastForwardPlanType = 0; - static int lastForwardPlanRows = 0; - static int lastForwardPlanCols = 0; - static fftwf_plan lastReversePlan = NULL; - static int lastReversePlanResultType = 0; - static int lastReversePlanRows = 0; - static int lastReversePlanCols = 0; +psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) +{ unsigned int numCols; unsigned int numRows; @@ -46,9 +38,4 @@ } - if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) { - psError(__func__,"Can not perform a forward FFT to real result."); - return NULL; - } - type = in->type.type; @@ -59,5 +46,5 @@ } - if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) { + if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) { psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", type); @@ -74,62 +61,12 @@ numRows = in->numRows; numCols = in->numCols; - // n.b., numCols needs adjustment for c->r, see below. - - if (flags & PS_FFT_REVERSE) { // reverse transform - if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) - // Assuming image is from a r2c transform. - // (FFTW only uses nx/2+1 for r->c transform results) - numCols = (numCols-1)*2; - - out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); - - if (lastReversePlanRows != numRows || lastReversePlanCols != numCols - || lastReversePlanResultType != PS_TYPE_F32 - || lastReversePlan == NULL) { - lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, - (fftwf_complex*)in->data.C32[0], - out->data.F32[0], - P_FFTW_PLAN_RIGOR); - } - } else { // complex result desired - out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); - if (lastReversePlanRows != numRows || lastReversePlanCols != numCols - || lastReversePlanResultType != PS_TYPE_C32 - || lastReversePlan == NULL) { - lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, - (fftwf_complex*)in->data.C32[0], - (fftwf_complex*)out->data.C32[0], - FFTW_BACKWARD, - P_FFTW_PLAN_RIGOR); - } - } - plan = lastReversePlan; - } else { // forward transform - if (type == PS_TYPE_F32) { // real data input - // (FFTW only uses nx/2+1 for r->c transform results) - out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); - - if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols - || lastForwardPlanType != type || lastForwardPlan == NULL) { - /* need new plan created */ - lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, - in->data.F32[0], - (fftwf_complex*)out->data.C32[0], - P_FFTW_PLAN_RIGOR); - } - } else if (type == PS_TYPE_C32) { // complex data input - out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); - - if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols - || lastForwardPlanType != type || lastForwardPlan == NULL) { - lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, - (fftwf_complex*)in->data.C32[0], - (fftwf_complex*)out->data.C32[0], - FFTW_FORWARD, - P_FFTW_PLAN_RIGOR); - } - } - plan = lastForwardPlan; - } + + out = psImageCopy(out,in, PS_TYPE_C32); + + plan = fftwf_plan_dft_2d(numCols, numRows, + (fftwf_complex*)out->data.C32[0], + (fftwf_complex*)out->data.C32[0], + direction, + P_FFTW_PLAN_RIGOR); /* check if a plan exists now*/ @@ -142,4 +79,6 @@ /* finally, call FFTW with the plan made above */ fftwf_execute(plan); + + fftwf_destroy_plan(plan); return out; @@ -475,12 +414,6 @@ /************************************** Vector Functions ***************************************/ -psVector* psVectorFFT(psVector* out, const psVector* in, int flags) -{ - static fftwf_plan lastForwardPlan = NULL; - static int lastForwardPlanType = 0; - static int lastForwardPlanElements = 0; - static fftwf_plan lastReversePlan = NULL; - static int lastReversePlanResultType = 0; - static int lastReversePlanElements = 0; +psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction) +{ unsigned int numElements; psElemType type; @@ -492,9 +425,4 @@ } - if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) { - psError(__func__,"Can not perform a forward FFT to real result."); - return NULL; - } - type = in->type.type; @@ -505,5 +433,5 @@ } - if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) { + if ( (type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE) ) { psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", type); @@ -519,64 +447,27 @@ numElements = in->n; - // n.b., numElements needs adjustment for c->r, see below. - - if (flags & PS_FFT_REVERSE) { // reverse transform - if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) - // Assuming image is from a r2c transform. - // (FFTW only uses nx/2+1 for r->c transform results) - numElements = (numElements-1)*2; - - out = psVectorRecycle(out,PS_TYPE_F32,numElements); - - if (lastReversePlanElements != numElements || - lastReversePlanResultType != PS_TYPE_F32 || - lastReversePlan == NULL) { - lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, - (fftwf_complex*)in->vec.cf, - out->vec.f, - P_FFTW_PLAN_RIGOR); - } - } else { // complex result desired - out = psVectorRecycle(out,PS_TYPE_C32,numElements); - if (lastReversePlanElements != numElements || - lastReversePlanResultType != PS_TYPE_C32 || - lastReversePlan == NULL) { - lastReversePlan = fftwf_plan_dft_1d(numElements, - (fftwf_complex*)in->vec.cf, - (fftwf_complex*)out->vec.cf, - FFTW_BACKWARD, - P_FFTW_PLAN_RIGOR); - } - } - plan = lastReversePlan; - } else { // forward transform - if (type == PS_TYPE_F32) { // real data input - // (FFTW only uses nx/2+1 for r->c transform results) - out = psVectorRecycle(out,PS_TYPE_C32,numElements); - - if (lastForwardPlanElements != numElements || - lastForwardPlanType != type || - lastForwardPlan == NULL) { - /* need new plan created */ - lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, - in->vec.f, - (fftwf_complex*)out->vec.cf, - P_FFTW_PLAN_RIGOR); - } - } else if (type == PS_TYPE_C32) { // complex data input - out = psVectorRecycle(out,PS_TYPE_C32,numElements); - - if (lastForwardPlanElements != numElements || - lastForwardPlanType != type || - lastForwardPlan == NULL) { - lastForwardPlan = fftwf_plan_dft_1d(numElements, - (fftwf_complex*)in->vec.cf, - (fftwf_complex*)out->vec.cf, - FFTW_FORWARD, - P_FFTW_PLAN_RIGOR); - } - } - plan = lastForwardPlan; - } + + out = psVectorRecycle(out,PS_TYPE_C32,numElements); + + if (type == PS_TYPE_F32) { + // need to convert to complex + psC32* outVec = out->vec.cf; + psF32* inVec = in->vec.f; + for (unsigned int i=0;ivec.cf; + psC32* inVec = in->vec.cf; + for (unsigned int i=0;ivec.cf, + (fftwf_complex*)out->vec.cf, + direction, + P_FFTW_PLAN_RIGOR); /* check if a plan exists now*/ @@ -590,6 +481,7 @@ fftwf_execute(plan); - return out; - + fftwf_destroy_plan(plan); + + return out; }