00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <string>
00037 #include <cstring>
00038 #include "emfft.h"
00039 #include "log.h"
00040
00041 #include <iostream>
00042 using std::cout;
00043 using std::endl;
00044
00045 #include "util.h"
00046
00047 #ifdef EMAN2_USING_CUDA
00048 #include "cuda/cuda_emfft.h"
00049 #endif
00050
00051 #ifdef DJBFFT
00052 extern "C" {
00053 #include <fftr4.h>
00054 }
00055 #endif //DJBFFT
00056
00057 using namespace EMAN;
00058
00059 namespace {
00060 int get_rank(int ny, int nz)
00061 {
00062 int rank = 3;
00063 if (ny == 1) {
00064 rank = 1;
00065 }
00066 else if (nz == 1) {
00067 rank = 2;
00068 }
00069 return rank;
00070 }
00071 }
00072 #ifdef _WIN32
00073 MUTEX fft_mutex;
00074 #else
00075 pthread_mutex_t fft_mutex=PTHREAD_MUTEX_INITIALIZER;
00076 #endif
00077
00078 #ifdef FFTW_PLAN_CACHING
00079
00080 const int EMfft::EMAN2_REAL_2_COMPLEX = 1;
00081 const int EMfft::EMAN2_COMPLEX_2_REAL = 2;
00082
00083
00084 const int EMfft::EMAN2_FFTW2_INPLACE = 1;
00085 const int EMfft::EMAN2_FFTW2_OUT_OF_PLACE=0;
00086
00087
00088 #ifdef FFTW3
00089 EMfft::EMfftw3_cache::EMfftw3_cache() :
00090 num_plans(0)
00091 {
00092 for(int i = 0; i < EMFFTW3_CACHE_SIZE; ++i)
00093 {
00094 rank[i] = 0;
00095 plan_dims[i][0] = 0; plan_dims[i][1] = 0; plan_dims[i][2] = 0;
00096 r2c[i] = -1;
00097 ip[i] = -1;
00098 fftwplans[i] = NULL;
00099 }
00100 }
00101
00102 void EMfft::EMfftw3_cache::debug_plans()
00103 {
00104 for(int i = 0; i < EMFFTW3_CACHE_SIZE; ++i)
00105 {
00106 cout << "Plan " << i << " has dims " << plan_dims[i][0] << " "
00107 << plan_dims[i][1] << " " <<
00108 plan_dims[i][2] << ", rank " <<
00109 rank[i] << ", rc flag "
00110 << r2c[i] << ", ip flag " << ip[i] << endl;
00111 }
00112 }
00113
00114 EMfft::EMfftw3_cache::~EMfftw3_cache()
00115 {
00116 for(int i = 0; i < EMFFTW3_CACHE_SIZE; ++i)
00117 {
00118 if (fftwplans[i] != NULL)
00119 {
00120 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00121 fftwf_destroy_plan(fftwplans[i]);
00122 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00123 fftwplans[i] = NULL;
00124 }
00125 }
00126 }
00127
00128 fftwf_plan EMfft::EMfftw3_cache::get_plan(const int rank_in, const int x, const int y, const int z, const int r2c_flag, const int ip_flag, fftwf_complex* complex_data, float* real_data )
00129 {
00130
00131 if ( rank_in > 3 || rank_in < 1 ) throw InvalidValueException(rank_in, "Error, can not get an FFTW plan using rank out of the range [1,3]");
00132 if ( r2c_flag != EMAN2_REAL_2_COMPLEX && r2c_flag != EMAN2_COMPLEX_2_REAL ) throw InvalidValueException(r2c_flag, "The real two complex flag is not supported");
00133
00134
00135
00136
00137 int dims[3];
00138 dims[0] = z;
00139 dims[1] = y;
00140 dims[2] = x;
00141
00142
00143 int i;
00144 for (i=0; i<num_plans; i++) {
00145 if (plan_dims[i][0]==x && plan_dims[i][1]==y && plan_dims[i][2]==z
00146 && rank[i]==rank_in && r2c[i]==r2c_flag && ip[i]==ip_flag) return fftwplans[i];
00147 }
00148
00149 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00150
00151 fftwf_plan plan;
00152
00153 if ( y == 1 && z == 1 )
00154 {
00155 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00156 plan = fftwf_plan_dft_r2c_1d(x, real_data, complex_data, FFTW_ESTIMATE);
00157 else
00158 plan = fftwf_plan_dft_c2r_1d(x, complex_data, real_data, FFTW_ESTIMATE);
00159 }
00160 else
00161 {
00162 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00163 plan = fftwf_plan_dft_r2c(rank_in, dims + (3 - rank_in), real_data, complex_data, FFTW_ESTIMATE);
00164 else
00165 plan = fftwf_plan_dft_c2r(rank_in, dims + (3 - rank_in), complex_data, real_data, FFTW_ESTIMATE);
00166 }
00167
00168 if (fftwplans[EMFFTW3_CACHE_SIZE-1] != NULL )
00169 {
00170 fftwf_destroy_plan(fftwplans[EMFFTW3_CACHE_SIZE-1]);
00171 fftwplans[EMFFTW3_CACHE_SIZE-1] = NULL;
00172 }
00173
00174 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00175
00176 int upper_limit = num_plans;
00177 if ( upper_limit == EMFFTW3_CACHE_SIZE ) upper_limit -= 1;
00178 for (int i=upper_limit-1; i>0; i--)
00179 {
00180 fftwplans[i]=fftwplans[i-1];
00181 rank[i]=rank[i-1];
00182 r2c[i]=r2c[i-1];
00183 ip[i]=ip[i-1];
00184 plan_dims[i][0]=plan_dims[i-1][0];
00185 plan_dims[i][1]=plan_dims[i-1][1];
00186 plan_dims[i][2]=plan_dims[i-1][2];
00187 }
00188
00189
00190 plan_dims[0][0]=x;
00191 plan_dims[0][1]=y;
00192 plan_dims[0][2]=z;
00193 r2c[0]=r2c_flag;
00194 ip[0]=ip_flag;
00195 fftwplans[0] = plan;
00196 rank[0]=rank_in;
00197 if (num_plans<EMFFTW3_CACHE_SIZE) num_plans++;
00198
00199
00200
00201
00202 return fftwplans[0];
00203
00204 }
00205
00206
00207 EMfft::EMfftw3_cache EMfft::plan_cache;
00208
00209 #endif // FFTW3
00210
00211 #ifdef FFTW2
00212
00213
00214 EMfft::EMfftw2_cache_nd::EMfftw2_cache_nd() :
00215 num_plans(0)
00216 {
00217 for(int i = 0; i < EMFFTW2_ND_CACHE_SIZE; ++i)
00218 {
00219 rank[i] = 0;
00220 plan_dims[i][0] = 0; plan_dims[i][1] = 0; plan_dims[i][2] = 0;
00221 r2c[i] = -1;
00222 ip[i] = -1;
00223 rfftwnd_plans[i] = NULL;
00224 }
00225 }
00226
00227 void EMfft::EMfftw2_cache_nd::debug_plans()
00228 {
00229 for(int i = 0; i < EMFFTW2_ND_CACHE_SIZE; ++i)
00230 {
00231 cout << "Plan " << i << " has dims " << plan_dims[i][0] << " " << plan_dims[i][1] << " " << plan_dims[i][2] << ", rank " <<
00232 rank[i] << ", rc flag " << r2c[i] << ", ip flag " << ip[i] << endl;
00233 }
00234 }
00235
00236 EMfft::EMfftw2_cache_nd::~EMfftw2_cache_nd()
00237 {
00238 for(int i = 0; i < EMFFTW2_ND_CACHE_SIZE; ++i)
00239 {
00240 if (rfftwnd_plans[i] != NULL)
00241 {
00242 rfftwnd_destroy_plan(rfftwnd_plans[i]);
00243 rfftwnd_plans[i] = NULL;
00244 }
00245 }
00246 }
00247
00248 rfftwnd_plan EMfft::EMfftw2_cache_nd::get_plan(const int rank_in, const int x, const int y, const int z, const int r2c_flag, const int ip_flag)
00249 {
00250
00251 if ( rank_in > 3 || rank_in < 1 ) throw InvalidValueException(rank_in, "Error, can not get an FFTW2 plan using rank out of the range [2,3]");
00252 if ( r2c_flag != EMAN2_REAL_2_COMPLEX && r2c_flag != EMAN2_COMPLEX_2_REAL ) throw InvalidValueException(r2c_flag, "The real two complex flag is not supported");
00253
00254 int dims[3];
00255 dims[0] = z;
00256 dims[1] = y;
00257 dims[2] = x;
00258
00259
00260 int i;
00261 for (i=0; i<num_plans; i++) {
00262 if (plan_dims[i][0]==x && plan_dims[i][1]==y && plan_dims[i][2]==z
00263 && rank[i]==rank_in && r2c[i]==r2c_flag && ip[i] == ip_flag) break;
00264 }
00265
00266
00267 if (i!=num_plans) {
00268 return rfftwnd_plans[i];
00269 }
00270 else
00271 {
00272 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00273
00274 rfftwnd_plan plan;
00275
00276 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00277 {
00278 if (ip_flag == EMAN2_FFTW2_INPLACE)
00279 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE|FFTW_IN_PLACE);
00280 else
00281 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00282 }
00283 else
00284 {
00285 if (ip_flag == EMAN2_FFTW2_INPLACE)
00286 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE|FFTW_IN_PLACE);
00287 else
00288 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE);
00289 }
00290
00291
00292 if (rfftwnd_plans[EMFFTW2_ND_CACHE_SIZE-1] != NULL )
00293 {
00294 rfftwnd_destroy_plan(rfftwnd_plans[EMFFTW2_ND_CACHE_SIZE-1]);
00295 rfftwnd_plans[EMFFTW2_ND_CACHE_SIZE-1] = NULL;
00296 }
00297 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00298
00299 int upper_limit = num_plans;
00300 if ( upper_limit == EMFFTW2_ND_CACHE_SIZE ) upper_limit -= 1;
00301 for (int i=upper_limit-1; i>0; i--)
00302 {
00303 rfftwnd_plans[i]=rfftwnd_plans[i-1];
00304 rank[i]=rank[i-1];
00305 r2c[i]=r2c[i-1];
00306 ip[i]=ip[i-1];
00307 plan_dims[i][0]=plan_dims[i-1][0];
00308 plan_dims[i][1]=plan_dims[i-1][1];
00309 plan_dims[i][2]=plan_dims[i-1][2];
00310 }
00311
00312
00313 plan_dims[0][0]=x;
00314 plan_dims[0][1]=y;
00315 plan_dims[0][2]=z;
00316 r2c[0]=r2c_flag;
00317 rfftwnd_plans[0] = plan;
00318 rank[0]=rank_in;
00319 ip[0]=ip_flag;
00320 if (num_plans<EMFFTW2_ND_CACHE_SIZE) num_plans++;
00321
00322 return rfftwnd_plans[0];
00323 }
00324 }
00325
00326 EMfft::EMfftw2_cache_1d::EMfftw2_cache_1d() :
00327 num_plans(0)
00328 {
00329 for(int i = 0; i < EMFFTW2_1D_CACHE_SIZE; ++i)
00330 {
00331 plan_dims[i] = 0;
00332 r2c[i] = -1;
00333 rfftw1d_plans[i] = NULL;
00334 }
00335 }
00336
00337 void EMfft::EMfftw2_cache_1d::debug_plans()
00338 {
00339 for(int i = 0; i < EMFFTW2_1D_CACHE_SIZE; ++i)
00340 {
00341 cout << "Plan " << i << " has length " << plan_dims[i] << ", rc flag " << r2c[i] << endl;
00342 }
00343 }
00344
00345 EMfft::EMfftw2_cache_1d::~EMfftw2_cache_1d()
00346 {
00347 for(int i = 0; i < EMFFTW2_1D_CACHE_SIZE; ++i)
00348 {
00349 if (rfftw1d_plans[i] != NULL)
00350 {
00351 rfftw_destroy_plan(rfftw1d_plans[i]);
00352 rfftw1d_plans[i] = NULL;
00353 }
00354 }
00355 }
00356
00357 rfftw_plan EMfft::EMfftw2_cache_1d::get_plan(const int x, const int r2c_flag )
00358 {
00359
00360 if ( r2c_flag != EMAN2_REAL_2_COMPLEX && r2c_flag != EMAN2_COMPLEX_2_REAL ) throw InvalidValueException(r2c_flag, "The real two complex flag is not supported");
00361
00362
00363 int i;
00364 for (i=0; i<num_plans; i++) {
00365 if (plan_dims[i]==x && r2c[i]==r2c_flag ) return rfftw1d_plans[i];
00366 }
00367
00368 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00369
00370 rfftw_plan plan;
00371
00372 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00373 plan = rfftw_create_plan(x, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00374 else
00375 plan = rfftw_create_plan(x, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00376
00377
00378
00379 if (rfftw1d_plans[EMFFTW2_1D_CACHE_SIZE-1] != NULL )
00380 {
00381 rfftw_destroy_plan(rfftw1d_plans[EMFFTW2_1D_CACHE_SIZE-1]);
00382 rfftw1d_plans[EMFFTW2_1D_CACHE_SIZE-1] = NULL;
00383 }
00384
00385 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00386
00387 int upper_limit = num_plans;
00388 if ( upper_limit == EMFFTW2_1D_CACHE_SIZE ) upper_limit -= 1;
00389 for (int i=upper_limit-1; i>0; i--)
00390 {
00391 rfftw1d_plans[i]=rfftw1d_plans[i-1];
00392 plan_dims[i]=plan_dims[i-1];
00393 r2c[i] = r2c[i-1];
00394 }
00395
00396
00397 plan_dims[0]=x;
00398 r2c[0]=r2c_flag;
00399 rfftw1d_plans[0] = plan;
00400 if (num_plans<EMFFTW2_1D_CACHE_SIZE) num_plans++;
00401
00402 return rfftw1d_plans[0];
00403
00404 }
00405
00406
00407 EMfft::EMfftw2_cache_nd EMfft::plan_nd_cache;
00408 EMfft::EMfftw2_cache_1d EMfft::plan_1d_cache;
00409
00410 #endif // FFTW2
00411 #endif // FFTW_PLAN_CACHING
00412
00413 #ifdef FFTW2
00414
00415 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00416 {
00417 #ifdef DJBFFT
00418 switch(n)
00419 {
00420 if ( n==2 || n==4 || n==8 || n==16 || n==32 || n==64 || n==128
00421 || n==256 || n==512 || n==1024 || n==2048 || n==4096 || n==8192 )
00422 {
00423 memcpy( complex_data, real_data, n * sizeof(float) );
00424 }
00425
00426 case 2:
00427 fftr4_2( (real4 *)complex_data );
00428
00429 case 4:
00430 fftr4_4( (real4 *)complex_data );
00431
00432 case 8:
00433 fftr4_8( (real4 *)complex_data );
00434
00435 case 16:
00436 fftr4_16( (real4 *)complex_data );
00437
00438 case 32:
00439 fftr4_32( (real4 *)complex_data );
00440
00441 case 64:
00442 fftr4_64( (real4 *)complex_data );
00443
00444 case 128:
00445 fftr4_128( (real4 *)complex_data );
00446
00447 case 256:
00448 fftr4_256( (real4 *)complex_data );
00449
00450 case 512:
00451 fftr4_512( (real4 *)complex_data );
00452
00453 case 1024:
00454 fftr4_1024( (real4 *)complex_data );
00455
00456 case 2048:
00457 fftr4_2048( (real4 *)complex_data );
00458
00459 case 4096:
00460 fftr4_4096( (real4 *)complex_data );
00461
00462 case 8192:
00463 fftr4_8192( (real4 *)complex_data );
00464
00465 default:
00466 const int complex_n = n + 2 - n%2;
00467 float * fft_data = new float[n];
00468 rfftw_plan p = rfftw_create_plan(n, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00469 rfftw_one(p, (fftw_real *) real_data, (fftw_real *) complex_data);
00470 rfftw_destroy_plan(p);
00471 for(int i=0; i<complex_n; ++i) {
00472 if(i%2==0) {
00473 complex_data[i] = fft_data[i/2];
00474 }
00475 else {
00476 if(i==1) {
00477 complex_data[i] = 0.0f;
00478 }
00479 else {
00480 if(n%2 == 0 && i == complex_n-1 ) {
00481 complex_data[i] = 0.0f;
00482 }
00483 else {
00484 complex_data[i] = fft_data[n-i/2];
00485 }
00486 }
00487 }
00488 }
00489
00490 delete [] fft_data;
00491 }
00492 #else
00493 const int complex_n = n + 2 - n%2;
00494 float * fft_data = new float[n];
00495 #ifdef FFTW_PLAN_CACHING
00496 rfftw_plan plan_1d = plan_1d_cache.get_plan(n,EMAN2_REAL_2_COMPLEX);
00497 rfftw_one(plan_1d, (fftw_real *) real_data, (fftw_real *) fft_data);
00498 #else
00499 rfftw_plan p = rfftw_create_plan(n, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00500 rfftw_one(p, (fftw_real *) real_data, (fftw_real *) fft_data);
00501 rfftw_destroy_plan(p);
00502 #endif //FFTW_PLAN_CACHING
00503
00504 for(int i=0; i<complex_n; ++i) {
00505 if(i%2==0) {
00506 complex_data[i] = fft_data[i/2];
00507 }
00508 else {
00509 if(i==1) {
00510 complex_data[i] = 0.0f;
00511 }
00512 else {
00513 if(n%2 == 0 && i == complex_n-1 ) {
00514 complex_data[i] = 0.0f;
00515 }
00516 else {
00517 complex_data[i] = fft_data[n-i/2];
00518 }
00519 }
00520 }
00521 }
00522
00523 delete [] fft_data;
00524 #endif //DJBFFT
00525
00526 return 0;
00527 }
00528
00529 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00530 {
00531 #ifdef DJBFFT
00532 switch(n)
00533 {
00534 if ( n==2 || n==4 || n==8 || n==16 || n==32 || n==64 || n==128
00535 || n==256 || n==512 || n==1024 || n==2048 || n==4096 || n==8192 )
00536 {
00537 memcpy( real_data, complex_data, n * sizeof(float) );
00538 }
00539
00540 case 2:
00541 fftr4_scale2( (real4 *)real_data );
00542 fftr4_un2( (real4 *)real_data );
00543
00544 case 4:
00545 fftr4_scale4( (real4 *)real_data );
00546 fftr4_un4( (real4 *)real_data );
00547
00548 case 8:
00549 fftr4_scale8( (real4 *)real_data );
00550 fftr4_un8( (real4 *)real_data );
00551
00552 case 16:
00553 fftr4_scale16( (real4 *)real_data );
00554 fftr4_un16( (real4 *)real_data );
00555
00556 case 32:
00557 fftr4_scale32( (real4 *)real_data );
00558 fftr4_un32( (real4 *)real_data );
00559
00560 case 64:
00561 fftr4_scale64( (real4 *)real_data );
00562 fftr4_un64( (real4 *)real_data );
00563
00564 case 128:
00565 fftr4_scale128( (real4 *)real_data );
00566 fftr4_un128( (real4 *)real_data );
00567
00568 case 256:
00569 fftr4_scale256( (real4 *)real_data );
00570 fftr4_un256( (real4 *)real_data );
00571
00572 case 512:
00573 fftr4_scale512( (real4 *)real_data );
00574 fftr4_un512( (real4 *)real_data );
00575
00576 case 1024:
00577 fftr4_scale1024( (real4 *)real_data );
00578 fftr4_un1024( (real4 *)real_data );
00579
00580 case 2048:
00581 fftr4_scale2048( (real4 *)real_data );
00582 fftr4_un2048( (real4 *)real_data );
00583
00584 case 4096:
00585 fftr4_scale4096( (real4 *)real_data );
00586 fftr4_un4096( (real4 *)real_data );
00587
00588 case 8192:
00589 fftr4_scale8192( (real4 *)real_data );
00590 fftr4_un8192( (real4 *)real_data );
00591
00592 default:
00593 const int complex_n = n + 2 - n%2;
00594 float * fft_data = new float[n];
00595
00596 for(int i=0; i<complex_n; ++i) {
00597 if(i%2 == 0) {
00598 fft_data[i/2] = complex_data[i];
00599 }
00600 else {
00601 if(i==1) {continue;}
00602 if(!(n%2 == 0 && i == complex_n-1)) {
00603 fft_data[n-i/2] = complex_data[i];
00604 }
00605 }
00606 }
00607
00608 rfftw_plan p = rfftw_create_plan(n, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00609 rfftw_one(p, (fftw_real *) fft_data, (fftw_real *) real_data);
00610 rfftw_destroy_plan(p);
00611 }
00612 #else
00613 const int complex_n = n + 2 - n%2;
00614 float * fft_data = new float[n];
00615
00616 for(int i=0; i<complex_n; ++i) {
00617 if(i%2 == 0) {
00618 fft_data[i/2] = complex_data[i];
00619 }
00620 else {
00621 if(i==1) {continue;}
00622 if(!(n%2 == 0 && i == complex_n-1)) {
00623 fft_data[n-i/2] = complex_data[i];
00624 }
00625 }
00626 }
00627 #ifdef FFTW_PLAN_CACHING
00628 rfftw_plan plan_1d = plan_1d_cache.get_plan(n,EMAN2_COMPLEX_2_REAL);
00629 rfftw_one(plan_1d, (fftw_real *) fft_data, (fftw_real *) real_data);
00630 #else
00631 rfftw_plan p = rfftw_create_plan(n, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00632 rfftw_one(p, (fftw_real *) fft_data, (fftw_real *) real_data);
00633 rfftw_destroy_plan(p);
00634 #endif // FFTW_PLAN_CACHING
00635 delete [] fft_data;
00636 #endif //DJBFFT
00637
00638 return 0;
00639
00640 }
00641
00642
00643
00644 int EMfft::complex_to_complex_1d(float *complex_data_in, float *complex_data_out, int n)
00645 {
00646 fftw_complex *in=(fftw_complex *) complex_data_in;
00647 fftw_complex *out=(fftw_complex *) complex_data_out;
00648 fftw_plan p = fftw_create_plan(n/2, FFTW_FORWARD, FFTW_ESTIMATE);
00649 fftw_one(p,in,out);
00650 fftw_destroy_plan(p);
00651 return 0;
00652 }
00653
00654
00655 int EMfft::complex_to_complex_nd(float *complex_data_in, float *complex_data_out, int nx,int ny,int nz)
00656 {
00657 const int rank = get_rank(ny, nz);
00658 int dims[3];
00659 dims[0] = nz;
00660 dims[1] = ny;
00661 dims[2] = nx;
00662
00663 switch(rank) {
00664 case 1:
00665 complex_to_complex_1d(complex_data_in, complex_data_out, nx);
00666 break;
00667
00668 case 2:
00669 case 3:
00670 {
00671 fftw_complex *in=(fftw_complex *) complex_data_in;
00672 fftw_complex *out=(fftw_complex *) complex_data_out;
00673 fftwnd_plan plan_nd;
00674 if(complex_data_out == complex_data_in) {
00675 plan_nd = fftwnd_create_plan(rank, dims+(3-rank),FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
00676 }
00677 else {
00678 plan_nd = fftw3d_create_plan(nx/2,ny,nz,FFTW_FORWARD, FFTW_ESTIMATE);
00679 }
00680 fftwnd_one(plan_nd, in, out);
00681 fftwnd_destroy_plan(plan_nd);
00682 }
00683 }
00684 return 0;
00685 }
00686
00687
00688
00689
00690
00691 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00692 {
00693
00694 const int rank = get_rank(ny, nz);
00695 int dims[3];
00696 dims[0] = nz;
00697 dims[1] = ny;
00698 dims[2] = nx;
00699
00700 switch(rank) {
00701 case 1:
00702 real_to_complex_1d(real_data, complex_data, nx);
00703 break;
00704
00705 case 2:
00706 case 3:
00707 {
00708 #ifdef FFTW_PLAN_CACHING
00709 bool ip = ( complex_data == real_data );
00710 rfftwnd_plan plan_nd = plan_nd_cache.get_plan(rank,nx,ny,nz,EMAN2_REAL_2_COMPLEX,ip);
00711 rfftwnd_one_real_to_complex(plan_nd, (fftw_real *) real_data, (fftw_complex *) complex_data);
00712 #else
00713 rfftwnd_plan plan_nd;
00714 if(real_data == complex_data) {
00715 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00716 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
00717 }
00718 else {
00719 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00720 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00721 }
00722 rfftwnd_one_real_to_complex(plan_nd, (fftw_real *) real_data,
00723 (fftw_complex *) complex_data);
00724 rfftwnd_destroy_plan(plan_nd);
00725 #endif // FFTW_PLAN_CACHING
00726 }
00727 break;
00728
00729 default:
00730 LOGERR("Should NEVER be here!!!");
00731 break;
00732 }
00733
00734 return 0;
00735 }
00736
00737 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00738 {
00739
00740 const int rank = get_rank(ny, nz);
00741 int dims[3];
00742 dims[0] = nz;
00743 dims[1] = ny;
00744 dims[2] = nx;
00745
00746 switch(rank) {
00747 case 1:
00748 complex_to_real_1d(complex_data, real_data, nx);
00749 break;
00750
00751 case 2:
00752 case 3:
00753 {
00754 #ifdef FFTW_PLAN_CACHING
00755 bool ip = ( complex_data == real_data );
00756 rfftwnd_plan plan_nd = plan_nd_cache.get_plan(rank,nx,ny,nz,EMAN2_COMPLEX_2_REAL,ip);
00757 rfftwnd_one_complex_to_real(plan_nd, (fftw_complex *) complex_data, (fftw_real *) real_data);
00758 #else
00759 rfftwnd_plan plan_nd;
00760 if(real_data == complex_data) {
00761 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00762 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
00763 }
00764 else {
00765 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00766 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00767 }
00768 rfftwnd_one_complex_to_real(plan_nd, (fftw_complex *) complex_data,
00769 (fftw_real *) real_data);
00770 rfftwnd_destroy_plan(plan_nd);
00771 #endif // FFTW_PLAN_CACHING
00772 }
00773 break;
00774
00775 default:
00776 LOGERR("Should NEVER be here!!!");
00777 break;
00778 }
00779
00780 return 0;
00781 }
00782 #endif //FFTW2
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807 #ifdef FFTW3
00808
00809 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00810 {
00811 #ifdef FFTW_PLAN_CACHING
00812 bool ip = ( complex_data == real_data );
00813 fftwf_plan plan = plan_cache.get_plan(1,n,1,1,EMAN2_REAL_2_COMPLEX,ip,(fftwf_complex *) complex_data, real_data);
00814
00815 fftwf_execute_dft_r2c(plan, real_data,(fftwf_complex *) complex_data);
00816 #else
00817 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00818 fftwf_plan plan = fftwf_plan_dft_r2c_1d(n, real_data, (fftwf_complex *) complex_data,
00819 FFTW_ESTIMATE);
00820 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00821
00822 fftwf_execute(plan);
00823 mrt = Util::MUTEX_LOCK(&fft_mutex);
00824 fftwf_destroy_plan(plan);
00825 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00826 #endif // FFTW_PLAN_CACHING
00827 return 0;
00828 };
00829
00830 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00831 {
00832 #ifdef FFTW_PLAN_CACHING
00833 bool ip = ( complex_data == real_data );
00834 fftwf_plan plan = plan_cache.get_plan(1,n,1,1,EMAN2_COMPLEX_2_REAL,ip,(fftwf_complex *) complex_data, real_data);
00835
00836 fftwf_execute_dft_c2r(plan, (fftwf_complex *) complex_data, real_data);
00837 #else
00838 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00839 fftwf_plan plan = fftwf_plan_dft_c2r_1d(n, (fftwf_complex *) complex_data, real_data,
00840 FFTW_ESTIMATE);
00841 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00842 fftwf_execute(plan);
00843 mrt = Util::MUTEX_LOCK(&fft_mutex);
00844 fftwf_destroy_plan(plan);
00845 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00846 #endif // FFTW_PLAN_CACHING
00847
00848 return 0;
00849 }
00850
00851
00852
00853 int EMfft::complex_to_complex_1d(float *complex_data_in, float *complex_data_out, int n)
00854 {
00855 fftwf_plan p;
00856 fftwf_complex *in=(fftwf_complex *) complex_data_in;
00857 fftwf_complex *out=(fftwf_complex *) complex_data_out;
00858 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00859 p=fftwf_plan_dft_1d(n/2,in,out, FFTW_FORWARD, FFTW_ESTIMATE);
00860 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00861 fftwf_execute(p);
00862 mrt = Util::MUTEX_LOCK(&fft_mutex);
00863 fftwf_destroy_plan(p);
00864 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00865 return 0;
00866 }
00867
00868
00869 int EMfft::complex_to_complex_nd(float *in, float *out, int nx,int ny,int nz)
00870 {
00871 const int rank = get_rank(ny, nz);
00872 int dims[3];
00873 dims[0] = nz;
00874 dims[1] = ny;
00875 dims[2] = nx;
00876
00877 switch(rank) {
00878 case 1:
00879 complex_to_complex_1d(in, out, nx);
00880 break;
00881
00882 case 2:
00883 case 3:
00884 {
00885 fftwf_plan p;
00886
00887 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00888
00889 if(out == in) {
00890 p=fftwf_plan_dft_3d(nx/2,ny,nz,(fftwf_complex *) in,(fftwf_complex *) out, FFTW_FORWARD, FFTW_ESTIMATE);
00891 }
00892 else {
00893
00894 p=fftwf_plan_dft_3d(nx/2,ny,nz,(fftwf_complex *) in,(fftwf_complex *) out, FFTW_FORWARD, FFTW_ESTIMATE);
00895 }
00896 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00897
00898 fftwf_execute(p);
00899
00900 mrt = Util::MUTEX_LOCK(&fft_mutex);
00901 fftwf_destroy_plan(p);
00902 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00903
00904 }
00905 }
00906 return 0;
00907 }
00908
00909
00910
00911 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00912 {
00913 const int rank = get_rank(ny, nz);
00914 int dims[3];
00915 dims[0] = nz;
00916 dims[1] = ny;
00917 dims[2] = nx;
00918
00919 switch(rank) {
00920 case 1:
00921 real_to_complex_1d(real_data, complex_data, nx);
00922 break;
00923
00924 case 2:
00925 case 3:
00926 {
00927 #ifdef FFTW_PLAN_CACHING
00928 bool ip = ( complex_data == real_data );
00929 fftwf_plan plan = plan_cache.get_plan(rank,nx,ny,nz,EMAN2_REAL_2_COMPLEX,ip,(fftwf_complex *) complex_data, real_data);
00930
00931 fftwf_execute_dft_r2c(plan, real_data,(fftwf_complex *) complex_data );
00932 #else
00933 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00934 fftwf_plan plan = fftwf_plan_dft_r2c(rank, dims + (3 - rank),
00935 real_data, (fftwf_complex *) complex_data, FFTW_ESTIMATE);
00936 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00937
00938 fftwf_execute(plan);
00939
00940 mrt = Util::MUTEX_LOCK(&fft_mutex);
00941 fftwf_destroy_plan(plan);
00942 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00943
00944 #endif // FFTW_PLAN_CACHING
00945 }
00946 break;
00947
00948 default:
00949 LOGERR("Should NEVER be here!!!");
00950 break;
00951 }
00952
00953 return 0;
00954 }
00955
00956 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00957 {
00958 const int rank = get_rank(ny, nz);
00959 int dims[3];
00960 dims[0] = nz;
00961 dims[1] = ny;
00962 dims[2] = nx;
00963
00964 switch(rank) {
00965 case 1:
00966 complex_to_real_1d(complex_data, real_data, nx);
00967 break;
00968
00969 case 2:
00970 case 3:
00971 {
00972 #ifdef FFTW_PLAN_CACHING
00973 bool ip = ( complex_data == real_data );
00974 fftwf_plan plan = plan_cache.get_plan(rank,nx,ny,nz,EMAN2_COMPLEX_2_REAL,ip,(fftwf_complex *) complex_data, real_data);
00975
00976 fftwf_execute_dft_c2r(plan, (fftwf_complex *) complex_data, real_data);
00977 #else
00978 int mrt = Util::MUTEX_LOCK(&fft_mutex);
00979 fftwf_plan plan = fftwf_plan_dft_c2r(rank, dims + (3 - rank),
00980 (fftwf_complex *) complex_data, real_data, FFTW_ESTIMATE);
00981 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00982
00983 fftwf_execute(plan);
00984
00985 mrt = Util::MUTEX_LOCK(&fft_mutex);
00986 fftwf_destroy_plan(plan);
00987 mrt = Util::MUTEX_UNLOCK(&fft_mutex);
00988
00989 #endif // FFTW_PLAN_CACHING
00990
00991
00992 }
00993 break;
00994
00995 default:
00996 LOGERR("Should NEVER be here!!!");
00997 break;
00998 }
00999
01000 return 0;
01001 }
01002
01003 #endif //FFTW3
01004
01005 #ifdef NATIVE_FFT
01006 #include "sparx/native_fft.h"
01007 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
01008 {
01009
01010 float * work = (float*) malloc((2*n+15)*sizeof(float));
01011 if (!work) {
01012 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
01013 LOGERR("real_to_complex_1d: failed to allocate work\n");
01014 }
01015
01016 rffti(n, work);
01017 memcpy(&complex_data[1], real_data, n * sizeof(float));
01018 rfftf(n, &complex_data[1], work);
01019 complex_data[0] = complex_data[1] ;
01020 complex_data[1] = 0.0f ;
01021 if (n%2 == 0) complex_data[n+1] = 0.0f ;
01022
01023 free(work);
01024 return 0;
01025 }
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
01065 {
01066
01067 float * work = (float*) malloc((2*n+15)*sizeof(float));
01068 if (!work) {
01069 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
01070 LOGERR("complex_to_real_1d: failed to allocate work\n");
01071 }
01072
01073 rffti(n, work);
01074
01075 memcpy(&real_data[1], &complex_data[2], (n-1) * sizeof(float));
01076 real_data[0] = complex_data[0] ;
01077 rfftb(n, real_data, work);
01078
01079 float nrm = 1.0f/float(n);
01080 for (int i = 0; i<n; i++) real_data[i] *= nrm;
01081 free(work);
01082 return 0;
01083 }
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
01102 {
01103 const int rank = get_rank(ny, nz);
01104 const int complex_nx = nx + 2 - nx%2;
01105
01106 switch(rank) {
01107 case 1:
01108 real_to_complex_1d(real_data, complex_data, nx);
01109 return 0;
01110 case 2:
01111 {
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132 int status = Nativefft::ftp_2rf(real_data, complex_data, complex_nx, nx, ny);
01133 if (status !=0) {
01134 fprintf(stderr, "real_to_complex_nd(2df): status = %d\n", status);
01135 LOGWARN("real_to_complex_nd(2df): status = %d\n", status);
01136 }
01137 return 0;
01138 }
01139 case 3:
01140 {
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155 int status = Nativefft::ftp_3rf(real_data, complex_data, complex_nx, nx, ny, nz);
01156 if (status !=0) {
01157 fprintf(stderr, "real_to_complex_nd(3df): status = %d\n", status);
01158 LOGWARN("real_to_complex_nd(3df): status = %d\n", status);
01159 }
01160
01161
01162 return 0;
01163 }
01164 default:
01165 LOGERR("Never should be here...\n");
01166 return -1;
01167 }
01168 }
01169
01170 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
01171 {
01172 const int rank = get_rank(ny, nz);
01173 const int complex_nx = nx + 2 - nx%2;
01174
01175 switch(rank) {
01176 case 1:
01177 complex_to_real_1d(complex_data, real_data, nx);
01178 return 0;
01179 case 2:
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201 {
01202 memcpy(real_data, complex_data, complex_nx*ny*sizeof(float));
01203
01204
01205 int status = Nativefft::ftp_2rb(real_data, complex_nx, nx, ny);
01206
01207 if (status !=0) {
01208 fprintf(stderr, "complex_to_real_nd(2db): status = %d\n", status);
01209 LOGWARN("complex_to_real_nd(2db): status = %d\n", status);
01210 }
01211 return 0;
01212 }
01213 case 3:
01214 {
01215 memcpy(real_data, complex_data, complex_nx*ny*nz*sizeof(float));
01216
01217
01218 int status = Nativefft::ftp_3rb(real_data, complex_nx, nx, ny, nz);
01219 if (status !=0) {
01220 fprintf(stderr, "complex_to_real_nd(3db): status = %d\n", status);
01221 LOGWARN("complex_to_real_nd(3db): status = %d\n", status);
01222 }
01223 return 0;
01224 }
01225 default:
01226 LOGERR("Never should be here...\n");
01227 return -1;
01228 }
01229 }
01230 #endif //NATIVE_FFT
01231
01232 #ifdef ACML
01233 #include <iostream>
01234 using std::cout;
01235 using std::endl;
01236
01237 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
01238 {
01239 const int complex_n = n + 2 - n%2;
01240 int info;
01241
01242
01243 float * comm = (float *)malloc((3*n+100)*sizeof(float));
01244
01245 float * fft_data = (float *)malloc(n * sizeof(float));
01246
01247
01248 memcpy(fft_data, real_data, n * sizeof(float));
01249
01250
01251 scfft(0, n, fft_data, comm, &info);
01252 if(info != 0) {
01253 LOGERR("Error happening in Initialize communication work array: %d", info);
01254 }
01255
01256
01257 scfft(1, n, fft_data, comm, &info);
01258 if(info != 0) {
01259 LOGERR("Error happening in Initialize communication work array: %d", info);
01260 }
01261
01268 transform(fft_data, fft_data+n, fft_data, time_sqrt_n(n));
01269
01270 for(int i=0; i<complex_n; ++i) {
01271 if(i%2==0) {
01272 complex_data[i] = fft_data[i/2];
01273 }
01274 else {
01275 if(i==1) {
01276 complex_data[i] = 0.0f;
01277 }
01278 else {
01279 if(n%2 == 0 && i == complex_n-1 ) {
01280 complex_data[i] = 0.0f;
01281 }
01282 else {
01283 complex_data[i] = fft_data[n-i/2];
01284 }
01285 }
01286 }
01287 }
01288
01289 free(fft_data);
01290 free(comm);
01291 return 0;
01292 }
01293
01294 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
01295 {
01296 int complex_n = n + 2 - n%2;
01297 int info;
01298
01299
01300 float * comm = (float *)malloc((3*n+100)*sizeof(float));
01301
01302 for(int i=0; i<complex_n; ++i) {
01303 if(i%2 == 0) {
01304 real_data[i/2] = complex_data[i];
01305 }
01306 else {
01307 if(i==1) {continue;}
01308 if(!(n%2 == 0 && i == complex_n-1)) {
01309 real_data[n-i/2] = complex_data[i];
01310 }
01311 }
01312 }
01313 transform(real_data, real_data+n, real_data, divide_sqrt_n(n));
01314
01315
01316 csfft(0, n, real_data, comm, &info);
01317 if(info != 0) {
01318 LOGERR("Error happening in Initialize communication work array: %d", info);
01319 }
01320
01321
01322 for (int j = n/2+1; j < n; j++) {
01323 real_data[j] = -real_data[j];
01324 }
01325
01326
01327 csfft(1, n, real_data, comm, &info);
01328 if(info != 0) {
01329 LOGERR("Error happening in Initialize communication work array: %d", info);
01330 }
01331
01332 free(comm);
01333 return 0;
01334 }
01335
01336 int EMfft::real_to_complex_2d(float *real_data, float *complex_data, int nx, int ny)
01337 {
01338 const int complex_nx = nx + 2 - nx%2;
01339 int info;
01340
01341 float * comm = (float *)malloc((3*nx+100)*sizeof(float));
01342
01343
01344 float * fft_data = (float *)malloc(complex_nx * ny * sizeof(float));
01345 cout << "fft_data after allocation:" << endl;
01346 for(int j=0; j<ny; ++j) {
01347 for(int i=0; i<complex_nx; ++i) {
01348 cout << "data[" << i << "][" << j << "] = " << fft_data[i+j*complex_nx] << "\t";
01349 }
01350 cout << endl;
01351 }
01352
01353
01354 for(int i=0; i<ny; ++i) {
01355 memcpy(fft_data+i*complex_nx, real_data+i*nx, nx*sizeof(float));
01356 }
01357
01358 cout << endl << "real_data array: " << endl;
01359 for(int j=0; j<ny; ++j) {
01360 for(int i=0; i<nx; ++i) {
01361 cout << real_data[i+j*nx] << "\t";
01362 }
01363 cout << endl;
01364 }
01365 cout << endl << "the fft_data array: " << endl;
01366 for(int j=0; j<ny; ++j) {
01367 for(int i=0; i<complex_nx; ++i) {
01368 cout << fft_data[i+j*complex_nx] << "\t";
01369 }
01370 cout << endl;
01371 }
01372
01373
01374 scfftm(ny, nx, fft_data, comm, &info);
01375
01376 cout << endl << "the fft_data array after x dim transform: " << endl;
01377 for(int j=0; j<ny; ++j) {
01378 for(int i=0; i<complex_nx; ++i) {
01379 cout << fft_data[i+j*complex_nx] << "\t";
01380 }
01381 cout << endl;
01382 }
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394 float * fft_data2 = (float *)malloc(complex_nx * ny * sizeof(float));
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405 if(info != 0) {
01406 LOGERR("Error happening in scfftm: %d", info);
01407 }
01408
01409 return 0;
01410 }
01411
01412 int EMfft::complex_to_real_2d(float *complex_data, float *real_data, int nx, int ny)
01413 {
01414 return 0;
01415 }
01416
01417 int EMfft::real_to_complex_3d(float *real_data, float *complex_data, int nx, int ny, int nz)
01418 {
01419 return 0;
01420 }
01421
01422 int EMfft::complex_to_real_3d(float *complex_data, float *real_data, int nx, int ny, int nz)
01423 {
01424 return 0;
01425 }
01426
01427 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
01428 {
01429 const int rank = get_rank(ny, nz);
01430
01431 switch(rank) {
01432 case 1:
01433 return real_to_complex_1d(real_data, complex_data, nx);
01434 case 2:
01435 return real_to_complex_2d(real_data, complex_data, nx, ny);
01436 case 3:
01437 return real_to_complex_3d(real_data, complex_data, nx, ny, nz);
01438 default:
01439 LOGERR("Never should be here...\n");
01440 return -1;
01441 }
01442 }
01443
01444 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
01445 {
01446 const int rank = get_rank(ny, nz);
01447
01448 switch(rank) {
01449 case 1:
01450 return complex_to_real_1d(complex_data, real_data, nx);
01451 case 2:
01452 return complex_to_real_2d(complex_data, real_data, nx, ny);
01453 case 3:
01454 return complex_to_real_3d(complex_data, real_data, nx, ny, nz);
01455 default:
01456 LOGERR("Never should be here...\n");
01457 return -1;
01458 }
01459 }
01460
01461
01462 #endif //ACML