Changeset 827
- Timestamp:
- Jun 1, 2004, 5:02:48 PM (22 years ago)
- Location:
- trunk/psLib/src
- Files:
-
- 6 edited
-
dataManip/psFFT.c (modified) (11 diffs)
-
dataManip/psFFT.h (modified) (3 diffs)
-
dataManip/psVectorFFT.c (modified) (11 diffs)
-
dataManip/psVectorFFT.h (modified) (3 diffs)
-
fft/psVectorFFT.c (modified) (11 diffs)
-
fft/psVectorFFT.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataManip/psFFT.c
r824 r827 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 6$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-06-0 1 22:42:57$7 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-06-02 03:02:48 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 26 26 static bool p_fftwWisdomImported = false; 27 27 28 psImage* psImageFFT(psImage* out, const psImage* in, int flags) 29 { 30 static fftwf_plan lastForwardPlan = NULL; 31 static int lastForwardPlanType = 0; 32 static int lastForwardPlanRows = 0; 33 static int lastForwardPlanCols = 0; 34 static fftwf_plan lastReversePlan = NULL; 35 static int lastReversePlanResultType = 0; 36 static int lastReversePlanRows = 0; 37 static int lastReversePlanCols = 0; 28 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) 29 { 38 30 unsigned int numCols; 39 31 unsigned int numRows; … … 46 38 } 47 39 48 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {49 psError(__func__,"Can not perform a forward FFT to real result.");50 return NULL;51 }52 53 40 type = in->type.type; 54 41 … … 59 46 } 60 47 61 if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0)) {48 if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) { 62 49 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 63 50 type); … … 74 61 numRows = in->numRows; 75 62 numCols = in->numCols; 76 // n.b., numCols needs adjustment for c->r, see below. 77 78 if (flags & PS_FFT_REVERSE) { // reverse transform 79 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 80 // Assuming image is from a r2c transform. 81 // (FFTW only uses nx/2+1 for r->c transform results) 82 numCols = (numCols-1)*2; 83 84 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 85 86 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 87 || lastReversePlanResultType != PS_TYPE_F32 88 || lastReversePlan == NULL) { 89 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, 90 (fftwf_complex*)in->data.C32[0], 91 out->data.F32[0], 92 P_FFTW_PLAN_RIGOR); 93 } 94 } else { // complex result desired 95 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 96 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 97 || lastReversePlanResultType != PS_TYPE_C32 98 || lastReversePlan == NULL) { 99 lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, 100 (fftwf_complex*)in->data.C32[0], 101 (fftwf_complex*)out->data.C32[0], 102 FFTW_BACKWARD, 103 P_FFTW_PLAN_RIGOR); 104 } 105 } 106 plan = lastReversePlan; 107 } else { // forward transform 108 if (type == PS_TYPE_F32) { // real data input 109 // (FFTW only uses nx/2+1 for r->c transform results) 110 out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); 111 112 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 113 || lastForwardPlanType != type || lastForwardPlan == NULL) { 114 /* need new plan created */ 115 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, 116 in->data.F32[0], 117 (fftwf_complex*)out->data.C32[0], 118 P_FFTW_PLAN_RIGOR); 119 } 120 } else if (type == PS_TYPE_C32) { // complex data input 121 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 122 123 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 124 || lastForwardPlanType != type || lastForwardPlan == NULL) { 125 lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, 126 (fftwf_complex*)in->data.C32[0], 127 (fftwf_complex*)out->data.C32[0], 128 FFTW_FORWARD, 129 P_FFTW_PLAN_RIGOR); 130 } 131 } 132 plan = lastForwardPlan; 133 } 63 64 out = psImageCopy(out,in, PS_TYPE_C32); 65 66 plan = fftwf_plan_dft_2d(numCols, numRows, 67 (fftwf_complex*)out->data.C32[0], 68 (fftwf_complex*)out->data.C32[0], 69 direction, 70 P_FFTW_PLAN_RIGOR); 134 71 135 72 /* check if a plan exists now*/ … … 142 79 /* finally, call FFTW with the plan made above */ 143 80 fftwf_execute(plan); 81 82 fftwf_destroy_plan(plan); 144 83 145 84 return out; … … 475 414 /************************************** Vector Functions ***************************************/ 476 415 477 psVector* psVectorFFT(psVector* out, const psVector* in, int flags) 478 { 479 static fftwf_plan lastForwardPlan = NULL; 480 static int lastForwardPlanType = 0; 481 static int lastForwardPlanElements = 0; 482 static fftwf_plan lastReversePlan = NULL; 483 static int lastReversePlanResultType = 0; 484 static int lastReversePlanElements = 0; 416 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction) 417 { 485 418 unsigned int numElements; 486 419 psElemType type; … … 492 425 } 493 426 494 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {495 psError(__func__,"Can not perform a forward FFT to real result.");496 return NULL;497 }498 499 427 type = in->type.type; 500 428 … … 505 433 } 506 434 507 if ( (type != PS_TYPE_C32) && ( (flags & PS_FFT_REVERSE) != 0) ) {435 if ( (type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE) ) { 508 436 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 509 437 type); … … 519 447 520 448 numElements = in->n; 521 // n.b., numElements needs adjustment for c->r, see below. 522 523 if (flags & PS_FFT_REVERSE) { // reverse transform 524 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 525 // Assuming image is from a r2c transform. 526 // (FFTW only uses nx/2+1 for r->c transform results) 527 numElements = (numElements-1)*2; 528 529 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 530 531 if (lastReversePlanElements != numElements || 532 lastReversePlanResultType != PS_TYPE_F32 || 533 lastReversePlan == NULL) { 534 lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, 535 (fftwf_complex*)in->vec.cf, 536 out->vec.f, 537 P_FFTW_PLAN_RIGOR); 538 } 539 } else { // complex result desired 540 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 541 if (lastReversePlanElements != numElements || 542 lastReversePlanResultType != PS_TYPE_C32 || 543 lastReversePlan == NULL) { 544 lastReversePlan = fftwf_plan_dft_1d(numElements, 545 (fftwf_complex*)in->vec.cf, 546 (fftwf_complex*)out->vec.cf, 547 FFTW_BACKWARD, 548 P_FFTW_PLAN_RIGOR); 549 } 550 } 551 plan = lastReversePlan; 552 } else { // forward transform 553 if (type == PS_TYPE_F32) { // real data input 554 // (FFTW only uses nx/2+1 for r->c transform results) 555 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 556 557 if (lastForwardPlanElements != numElements || 558 lastForwardPlanType != type || 559 lastForwardPlan == NULL) { 560 /* need new plan created */ 561 lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, 562 in->vec.f, 563 (fftwf_complex*)out->vec.cf, 564 P_FFTW_PLAN_RIGOR); 565 } 566 } else if (type == PS_TYPE_C32) { // complex data input 567 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 568 569 if (lastForwardPlanElements != numElements || 570 lastForwardPlanType != type || 571 lastForwardPlan == NULL) { 572 lastForwardPlan = fftwf_plan_dft_1d(numElements, 573 (fftwf_complex*)in->vec.cf, 574 (fftwf_complex*)out->vec.cf, 575 FFTW_FORWARD, 576 P_FFTW_PLAN_RIGOR); 577 } 578 } 579 plan = lastForwardPlan; 580 } 449 450 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 451 452 if (type == PS_TYPE_F32) { 453 // need to convert to complex 454 psC32* outVec = out->vec.cf; 455 psF32* inVec = in->vec.f; 456 for (unsigned int i=0;i<numElements;i++) { 457 outVec[i] = inVec[i]; 458 } 459 } else { 460 psC32* outVec = out->vec.cf; 461 psC32* inVec = in->vec.cf; 462 for (unsigned int i=0;i<numElements;i++) { 463 outVec[i] = inVec[i]; 464 } 465 } 466 467 plan = fftwf_plan_dft_1d(numElements, 468 (fftwf_complex*)out->vec.cf, 469 (fftwf_complex*)out->vec.cf, 470 direction, 471 P_FFTW_PLAN_RIGOR); 581 472 582 473 /* check if a plan exists now*/ … … 590 481 fftwf_execute(plan); 591 482 592 return out; 593 483 fftwf_destroy_plan(plan); 484 485 return out; 594 486 } 595 487 -
trunk/psLib/src/dataManip/psFFT.h
r823 r827 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-06-0 1 22:27:16$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-06-02 03:02:48 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 typedef enum { 24 24 /// psImageFFT/psVectorFFT should perform a forward FFT. 25 PS_FFT_FORWARD = 0,25 PS_FFT_FORWARD = FFTW_FORWARD, 26 26 27 27 ///< psImageFFT/psVectorFFT should perform a reverse FFT. 28 PS_FFT_REVERSE = 1, 28 PS_FFT_REVERSE = FFTW_BACKWARD 29 } psFftDirection; 29 30 30 /** psImageFFT/psVectorFFT should return a real image. This is valid for 31 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter 32 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. 33 */ 34 PS_FFT_REAL_RESULT = 2 35 } psFftFlags; 36 37 psImage* psImageFFT(psImage* out, const psImage* in, int flags); 31 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction); 38 32 psImage* psImageReal(psImage *out, const psImage* in); 39 33 psImage* psImageImaginary(psImage *out, const psImage* in); … … 42 36 psImage* psImagePowerSpectrum(psImage* out, const psImage* in); 43 37 44 psVector* psVectorFFT(psVector* out, const psVector* in, int flags);38 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction); 45 39 psVector* psVectorReal(psVector* out, const psVector* in); 46 40 psVector* psVectorImaginary(psVector* out, const psVector* in); -
trunk/psLib/src/dataManip/psVectorFFT.c
r824 r827 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 6$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-06-0 1 22:42:57$7 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-06-02 03:02:48 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 26 26 static bool p_fftwWisdomImported = false; 27 27 28 psImage* psImageFFT(psImage* out, const psImage* in, int flags) 29 { 30 static fftwf_plan lastForwardPlan = NULL; 31 static int lastForwardPlanType = 0; 32 static int lastForwardPlanRows = 0; 33 static int lastForwardPlanCols = 0; 34 static fftwf_plan lastReversePlan = NULL; 35 static int lastReversePlanResultType = 0; 36 static int lastReversePlanRows = 0; 37 static int lastReversePlanCols = 0; 28 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) 29 { 38 30 unsigned int numCols; 39 31 unsigned int numRows; … … 46 38 } 47 39 48 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {49 psError(__func__,"Can not perform a forward FFT to real result.");50 return NULL;51 }52 53 40 type = in->type.type; 54 41 … … 59 46 } 60 47 61 if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0)) {48 if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) { 62 49 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 63 50 type); … … 74 61 numRows = in->numRows; 75 62 numCols = in->numCols; 76 // n.b., numCols needs adjustment for c->r, see below. 77 78 if (flags & PS_FFT_REVERSE) { // reverse transform 79 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 80 // Assuming image is from a r2c transform. 81 // (FFTW only uses nx/2+1 for r->c transform results) 82 numCols = (numCols-1)*2; 83 84 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 85 86 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 87 || lastReversePlanResultType != PS_TYPE_F32 88 || lastReversePlan == NULL) { 89 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, 90 (fftwf_complex*)in->data.C32[0], 91 out->data.F32[0], 92 P_FFTW_PLAN_RIGOR); 93 } 94 } else { // complex result desired 95 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 96 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 97 || lastReversePlanResultType != PS_TYPE_C32 98 || lastReversePlan == NULL) { 99 lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, 100 (fftwf_complex*)in->data.C32[0], 101 (fftwf_complex*)out->data.C32[0], 102 FFTW_BACKWARD, 103 P_FFTW_PLAN_RIGOR); 104 } 105 } 106 plan = lastReversePlan; 107 } else { // forward transform 108 if (type == PS_TYPE_F32) { // real data input 109 // (FFTW only uses nx/2+1 for r->c transform results) 110 out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); 111 112 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 113 || lastForwardPlanType != type || lastForwardPlan == NULL) { 114 /* need new plan created */ 115 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, 116 in->data.F32[0], 117 (fftwf_complex*)out->data.C32[0], 118 P_FFTW_PLAN_RIGOR); 119 } 120 } else if (type == PS_TYPE_C32) { // complex data input 121 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 122 123 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 124 || lastForwardPlanType != type || lastForwardPlan == NULL) { 125 lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, 126 (fftwf_complex*)in->data.C32[0], 127 (fftwf_complex*)out->data.C32[0], 128 FFTW_FORWARD, 129 P_FFTW_PLAN_RIGOR); 130 } 131 } 132 plan = lastForwardPlan; 133 } 63 64 out = psImageCopy(out,in, PS_TYPE_C32); 65 66 plan = fftwf_plan_dft_2d(numCols, numRows, 67 (fftwf_complex*)out->data.C32[0], 68 (fftwf_complex*)out->data.C32[0], 69 direction, 70 P_FFTW_PLAN_RIGOR); 134 71 135 72 /* check if a plan exists now*/ … … 142 79 /* finally, call FFTW with the plan made above */ 143 80 fftwf_execute(plan); 81 82 fftwf_destroy_plan(plan); 144 83 145 84 return out; … … 475 414 /************************************** Vector Functions ***************************************/ 476 415 477 psVector* psVectorFFT(psVector* out, const psVector* in, int flags) 478 { 479 static fftwf_plan lastForwardPlan = NULL; 480 static int lastForwardPlanType = 0; 481 static int lastForwardPlanElements = 0; 482 static fftwf_plan lastReversePlan = NULL; 483 static int lastReversePlanResultType = 0; 484 static int lastReversePlanElements = 0; 416 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction) 417 { 485 418 unsigned int numElements; 486 419 psElemType type; … … 492 425 } 493 426 494 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {495 psError(__func__,"Can not perform a forward FFT to real result.");496 return NULL;497 }498 499 427 type = in->type.type; 500 428 … … 505 433 } 506 434 507 if ( (type != PS_TYPE_C32) && ( (flags & PS_FFT_REVERSE) != 0) ) {435 if ( (type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE) ) { 508 436 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 509 437 type); … … 519 447 520 448 numElements = in->n; 521 // n.b., numElements needs adjustment for c->r, see below. 522 523 if (flags & PS_FFT_REVERSE) { // reverse transform 524 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 525 // Assuming image is from a r2c transform. 526 // (FFTW only uses nx/2+1 for r->c transform results) 527 numElements = (numElements-1)*2; 528 529 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 530 531 if (lastReversePlanElements != numElements || 532 lastReversePlanResultType != PS_TYPE_F32 || 533 lastReversePlan == NULL) { 534 lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, 535 (fftwf_complex*)in->vec.cf, 536 out->vec.f, 537 P_FFTW_PLAN_RIGOR); 538 } 539 } else { // complex result desired 540 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 541 if (lastReversePlanElements != numElements || 542 lastReversePlanResultType != PS_TYPE_C32 || 543 lastReversePlan == NULL) { 544 lastReversePlan = fftwf_plan_dft_1d(numElements, 545 (fftwf_complex*)in->vec.cf, 546 (fftwf_complex*)out->vec.cf, 547 FFTW_BACKWARD, 548 P_FFTW_PLAN_RIGOR); 549 } 550 } 551 plan = lastReversePlan; 552 } else { // forward transform 553 if (type == PS_TYPE_F32) { // real data input 554 // (FFTW only uses nx/2+1 for r->c transform results) 555 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 556 557 if (lastForwardPlanElements != numElements || 558 lastForwardPlanType != type || 559 lastForwardPlan == NULL) { 560 /* need new plan created */ 561 lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, 562 in->vec.f, 563 (fftwf_complex*)out->vec.cf, 564 P_FFTW_PLAN_RIGOR); 565 } 566 } else if (type == PS_TYPE_C32) { // complex data input 567 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 568 569 if (lastForwardPlanElements != numElements || 570 lastForwardPlanType != type || 571 lastForwardPlan == NULL) { 572 lastForwardPlan = fftwf_plan_dft_1d(numElements, 573 (fftwf_complex*)in->vec.cf, 574 (fftwf_complex*)out->vec.cf, 575 FFTW_FORWARD, 576 P_FFTW_PLAN_RIGOR); 577 } 578 } 579 plan = lastForwardPlan; 580 } 449 450 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 451 452 if (type == PS_TYPE_F32) { 453 // need to convert to complex 454 psC32* outVec = out->vec.cf; 455 psF32* inVec = in->vec.f; 456 for (unsigned int i=0;i<numElements;i++) { 457 outVec[i] = inVec[i]; 458 } 459 } else { 460 psC32* outVec = out->vec.cf; 461 psC32* inVec = in->vec.cf; 462 for (unsigned int i=0;i<numElements;i++) { 463 outVec[i] = inVec[i]; 464 } 465 } 466 467 plan = fftwf_plan_dft_1d(numElements, 468 (fftwf_complex*)out->vec.cf, 469 (fftwf_complex*)out->vec.cf, 470 direction, 471 P_FFTW_PLAN_RIGOR); 581 472 582 473 /* check if a plan exists now*/ … … 590 481 fftwf_execute(plan); 591 482 592 return out; 593 483 fftwf_destroy_plan(plan); 484 485 return out; 594 486 } 595 487 -
trunk/psLib/src/dataManip/psVectorFFT.h
r823 r827 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-06-0 1 22:27:16$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-06-02 03:02:48 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 typedef enum { 24 24 /// psImageFFT/psVectorFFT should perform a forward FFT. 25 PS_FFT_FORWARD = 0,25 PS_FFT_FORWARD = FFTW_FORWARD, 26 26 27 27 ///< psImageFFT/psVectorFFT should perform a reverse FFT. 28 PS_FFT_REVERSE = 1, 28 PS_FFT_REVERSE = FFTW_BACKWARD 29 } psFftDirection; 29 30 30 /** psImageFFT/psVectorFFT should return a real image. This is valid for 31 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter 32 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. 33 */ 34 PS_FFT_REAL_RESULT = 2 35 } psFftFlags; 36 37 psImage* psImageFFT(psImage* out, const psImage* in, int flags); 31 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction); 38 32 psImage* psImageReal(psImage *out, const psImage* in); 39 33 psImage* psImageImaginary(psImage *out, const psImage* in); … … 42 36 psImage* psImagePowerSpectrum(psImage* out, const psImage* in); 43 37 44 psVector* psVectorFFT(psVector* out, const psVector* in, int flags);38 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction); 45 39 psVector* psVectorReal(psVector* out, const psVector* in); 46 40 psVector* psVectorImaginary(psVector* out, const psVector* in); -
trunk/psLib/src/fft/psVectorFFT.c
r824 r827 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 6$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-06-0 1 22:42:57$7 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-06-02 03:02:48 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 26 26 static bool p_fftwWisdomImported = false; 27 27 28 psImage* psImageFFT(psImage* out, const psImage* in, int flags) 29 { 30 static fftwf_plan lastForwardPlan = NULL; 31 static int lastForwardPlanType = 0; 32 static int lastForwardPlanRows = 0; 33 static int lastForwardPlanCols = 0; 34 static fftwf_plan lastReversePlan = NULL; 35 static int lastReversePlanResultType = 0; 36 static int lastReversePlanRows = 0; 37 static int lastReversePlanCols = 0; 28 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) 29 { 38 30 unsigned int numCols; 39 31 unsigned int numRows; … … 46 38 } 47 39 48 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {49 psError(__func__,"Can not perform a forward FFT to real result.");50 return NULL;51 }52 53 40 type = in->type.type; 54 41 … … 59 46 } 60 47 61 if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0)) {48 if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) { 62 49 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 63 50 type); … … 74 61 numRows = in->numRows; 75 62 numCols = in->numCols; 76 // n.b., numCols needs adjustment for c->r, see below. 77 78 if (flags & PS_FFT_REVERSE) { // reverse transform 79 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 80 // Assuming image is from a r2c transform. 81 // (FFTW only uses nx/2+1 for r->c transform results) 82 numCols = (numCols-1)*2; 83 84 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 85 86 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 87 || lastReversePlanResultType != PS_TYPE_F32 88 || lastReversePlan == NULL) { 89 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, 90 (fftwf_complex*)in->data.C32[0], 91 out->data.F32[0], 92 P_FFTW_PLAN_RIGOR); 93 } 94 } else { // complex result desired 95 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 96 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 97 || lastReversePlanResultType != PS_TYPE_C32 98 || lastReversePlan == NULL) { 99 lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, 100 (fftwf_complex*)in->data.C32[0], 101 (fftwf_complex*)out->data.C32[0], 102 FFTW_BACKWARD, 103 P_FFTW_PLAN_RIGOR); 104 } 105 } 106 plan = lastReversePlan; 107 } else { // forward transform 108 if (type == PS_TYPE_F32) { // real data input 109 // (FFTW only uses nx/2+1 for r->c transform results) 110 out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); 111 112 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 113 || lastForwardPlanType != type || lastForwardPlan == NULL) { 114 /* need new plan created */ 115 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, 116 in->data.F32[0], 117 (fftwf_complex*)out->data.C32[0], 118 P_FFTW_PLAN_RIGOR); 119 } 120 } else if (type == PS_TYPE_C32) { // complex data input 121 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 122 123 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 124 || lastForwardPlanType != type || lastForwardPlan == NULL) { 125 lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, 126 (fftwf_complex*)in->data.C32[0], 127 (fftwf_complex*)out->data.C32[0], 128 FFTW_FORWARD, 129 P_FFTW_PLAN_RIGOR); 130 } 131 } 132 plan = lastForwardPlan; 133 } 63 64 out = psImageCopy(out,in, PS_TYPE_C32); 65 66 plan = fftwf_plan_dft_2d(numCols, numRows, 67 (fftwf_complex*)out->data.C32[0], 68 (fftwf_complex*)out->data.C32[0], 69 direction, 70 P_FFTW_PLAN_RIGOR); 134 71 135 72 /* check if a plan exists now*/ … … 142 79 /* finally, call FFTW with the plan made above */ 143 80 fftwf_execute(plan); 81 82 fftwf_destroy_plan(plan); 144 83 145 84 return out; … … 475 414 /************************************** Vector Functions ***************************************/ 476 415 477 psVector* psVectorFFT(psVector* out, const psVector* in, int flags) 478 { 479 static fftwf_plan lastForwardPlan = NULL; 480 static int lastForwardPlanType = 0; 481 static int lastForwardPlanElements = 0; 482 static fftwf_plan lastReversePlan = NULL; 483 static int lastReversePlanResultType = 0; 484 static int lastReversePlanElements = 0; 416 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction) 417 { 485 418 unsigned int numElements; 486 419 psElemType type; … … 492 425 } 493 426 494 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {495 psError(__func__,"Can not perform a forward FFT to real result.");496 return NULL;497 }498 499 427 type = in->type.type; 500 428 … … 505 433 } 506 434 507 if ( (type != PS_TYPE_C32) && ( (flags & PS_FFT_REVERSE) != 0) ) {435 if ( (type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE) ) { 508 436 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 509 437 type); … … 519 447 520 448 numElements = in->n; 521 // n.b., numElements needs adjustment for c->r, see below. 522 523 if (flags & PS_FFT_REVERSE) { // reverse transform 524 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 525 // Assuming image is from a r2c transform. 526 // (FFTW only uses nx/2+1 for r->c transform results) 527 numElements = (numElements-1)*2; 528 529 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 530 531 if (lastReversePlanElements != numElements || 532 lastReversePlanResultType != PS_TYPE_F32 || 533 lastReversePlan == NULL) { 534 lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, 535 (fftwf_complex*)in->vec.cf, 536 out->vec.f, 537 P_FFTW_PLAN_RIGOR); 538 } 539 } else { // complex result desired 540 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 541 if (lastReversePlanElements != numElements || 542 lastReversePlanResultType != PS_TYPE_C32 || 543 lastReversePlan == NULL) { 544 lastReversePlan = fftwf_plan_dft_1d(numElements, 545 (fftwf_complex*)in->vec.cf, 546 (fftwf_complex*)out->vec.cf, 547 FFTW_BACKWARD, 548 P_FFTW_PLAN_RIGOR); 549 } 550 } 551 plan = lastReversePlan; 552 } else { // forward transform 553 if (type == PS_TYPE_F32) { // real data input 554 // (FFTW only uses nx/2+1 for r->c transform results) 555 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 556 557 if (lastForwardPlanElements != numElements || 558 lastForwardPlanType != type || 559 lastForwardPlan == NULL) { 560 /* need new plan created */ 561 lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, 562 in->vec.f, 563 (fftwf_complex*)out->vec.cf, 564 P_FFTW_PLAN_RIGOR); 565 } 566 } else if (type == PS_TYPE_C32) { // complex data input 567 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 568 569 if (lastForwardPlanElements != numElements || 570 lastForwardPlanType != type || 571 lastForwardPlan == NULL) { 572 lastForwardPlan = fftwf_plan_dft_1d(numElements, 573 (fftwf_complex*)in->vec.cf, 574 (fftwf_complex*)out->vec.cf, 575 FFTW_FORWARD, 576 P_FFTW_PLAN_RIGOR); 577 } 578 } 579 plan = lastForwardPlan; 580 } 449 450 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 451 452 if (type == PS_TYPE_F32) { 453 // need to convert to complex 454 psC32* outVec = out->vec.cf; 455 psF32* inVec = in->vec.f; 456 for (unsigned int i=0;i<numElements;i++) { 457 outVec[i] = inVec[i]; 458 } 459 } else { 460 psC32* outVec = out->vec.cf; 461 psC32* inVec = in->vec.cf; 462 for (unsigned int i=0;i<numElements;i++) { 463 outVec[i] = inVec[i]; 464 } 465 } 466 467 plan = fftwf_plan_dft_1d(numElements, 468 (fftwf_complex*)out->vec.cf, 469 (fftwf_complex*)out->vec.cf, 470 direction, 471 P_FFTW_PLAN_RIGOR); 581 472 582 473 /* check if a plan exists now*/ … … 590 481 fftwf_execute(plan); 591 482 592 return out; 593 483 fftwf_destroy_plan(plan); 484 485 return out; 594 486 } 595 487 -
trunk/psLib/src/fft/psVectorFFT.h
r823 r827 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-06-0 1 22:27:16$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-06-02 03:02:48 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 typedef enum { 24 24 /// psImageFFT/psVectorFFT should perform a forward FFT. 25 PS_FFT_FORWARD = 0,25 PS_FFT_FORWARD = FFTW_FORWARD, 26 26 27 27 ///< psImageFFT/psVectorFFT should perform a reverse FFT. 28 PS_FFT_REVERSE = 1, 28 PS_FFT_REVERSE = FFTW_BACKWARD 29 } psFftDirection; 29 30 30 /** psImageFFT/psVectorFFT should return a real image. This is valid for 31 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter 32 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. 33 */ 34 PS_FFT_REAL_RESULT = 2 35 } psFftFlags; 36 37 psImage* psImageFFT(psImage* out, const psImage* in, int flags); 31 psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction); 38 32 psImage* psImageReal(psImage *out, const psImage* in); 39 33 psImage* psImageImaginary(psImage *out, const psImage* in); … … 42 36 psImage* psImagePowerSpectrum(psImage* out, const psImage* in); 43 37 44 psVector* psVectorFFT(psVector* out, const psVector* in, int flags);38 psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction); 45 39 psVector* psVectorReal(psVector* out, const psVector* in); 46 40 psVector* psVectorImaginary(psVector* out, const psVector* in);
Note:
See TracChangeset
for help on using the changeset viewer.
