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