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 CUDA_FFT
00048 #include "cuda/cuda_emfft.h"
00049 #endif //CUDA_FFT
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 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00625 {
00626
00627 const int rank = get_rank(ny, nz);
00628 int dims[3];
00629 dims[0] = nz;
00630 dims[1] = ny;
00631 dims[2] = nx;
00632
00633 switch(rank) {
00634 case 1:
00635 real_to_complex_1d(real_data, complex_data, nx);
00636 break;
00637
00638 case 2:
00639 case 3:
00640 {
00641 #ifdef FFTW_PLAN_CACHING
00642 bool ip = ( complex_data == real_data );
00643 rfftwnd_plan plan_nd = plan_nd_cache.get_plan(rank,nx,ny,nz,EMAN2_REAL_2_COMPLEX,ip);
00644 rfftwnd_one_real_to_complex(plan_nd, (fftw_real *) real_data, (fftw_complex *) complex_data);
00645 #else
00646 rfftwnd_plan plan_nd;
00647 if(real_data == complex_data) {
00648 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00649 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
00650 }
00651 else {
00652 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00653 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00654 }
00655 rfftwnd_one_real_to_complex(plan_nd, (fftw_real *) real_data,
00656 (fftw_complex *) complex_data);
00657 rfftwnd_destroy_plan(plan_nd);
00658 #endif // FFTW_PLAN_CACHING
00659 }
00660 break;
00661
00662 default:
00663 LOGERR("Should NEVER be here!!!");
00664 break;
00665 }
00666
00667 return 0;
00668 }
00669
00670 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00671 {
00672
00673 const int rank = get_rank(ny, nz);
00674 int dims[3];
00675 dims[0] = nz;
00676 dims[1] = ny;
00677 dims[2] = nx;
00678
00679 switch(rank) {
00680 case 1:
00681 complex_to_real_1d(complex_data, real_data, nx);
00682 break;
00683
00684 case 2:
00685 case 3:
00686 {
00687 #ifdef FFTW_PLAN_CACHING
00688 bool ip = ( complex_data == real_data );
00689 rfftwnd_plan plan_nd = plan_nd_cache.get_plan(rank,nx,ny,nz,EMAN2_COMPLEX_2_REAL,ip);
00690 rfftwnd_one_complex_to_real(plan_nd, (fftw_complex *) complex_data, (fftw_real *) real_data);
00691 #else
00692 rfftwnd_plan plan_nd;
00693 if(real_data == complex_data) {
00694 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00695 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
00696 }
00697 else {
00698 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank),
00699 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00700 }
00701 rfftwnd_one_complex_to_real(plan_nd, (fftw_complex *) complex_data,
00702 (fftw_real *) real_data);
00703 rfftwnd_destroy_plan(plan_nd);
00704 #endif // FFTW_PLAN_CACHING
00705 }
00706 break;
00707
00708 default:
00709 LOGERR("Should NEVER be here!!!");
00710 break;
00711 }
00712
00713 return 0;
00714 }
00715 #endif //FFTW2
00716
00717 #ifdef CUDA_FFT
00718 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00719 {
00720 return cuda_fft_real_to_complex_1d(real_data,complex_data,n);
00721 }
00722
00723 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00724 {
00725 return cuda_fft_complex_to_real_1d(complex_data,real_data,n);
00726 }
00727
00728 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00729 {
00730 return cuda_fft_real_to_complex_nd(real_data,complex_data,nx,ny,nz);
00731 }
00732
00733 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00734 {
00735 return cuda_fft_complex_to_real_nd(complex_data,real_data,nx,ny,nz);
00736 }
00737
00738 #endif
00739
00740 #ifdef FFTW3
00741
00742 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00743 {
00744 #ifdef FFTW_PLAN_CACHING
00745 bool ip = ( complex_data == real_data );
00746 fftwf_plan plan = plan_cache.get_plan(1,n,1,1,EMAN2_REAL_2_COMPLEX,ip,(fftwf_complex *) complex_data, real_data);
00747
00748 fftwf_execute_dft_r2c(plan, real_data,(fftwf_complex *) complex_data);
00749 #else
00750 fftwf_plan plan = fftwf_plan_dft_r2c_1d(n, real_data, (fftwf_complex *) complex_data,
00751 FFTW_ESTIMATE);
00752 fftwf_execute(plan);
00753 fftwf_destroy_plan(plan);
00754 #endif // FFTW_PLAN_CACHING
00755 return 0;
00756 };
00757
00758 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00759 {
00760 #ifdef FFTW_PLAN_CACHING
00761 bool ip = ( complex_data == real_data );
00762 fftwf_plan plan = plan_cache.get_plan(1,n,1,1,EMAN2_COMPLEX_2_REAL,ip,(fftwf_complex *) complex_data, real_data);
00763
00764 fftwf_execute_dft_c2r(plan, (fftwf_complex *) complex_data, real_data);
00765 #else
00766 fftwf_plan plan = fftwf_plan_dft_c2r_1d(n, (fftwf_complex *) complex_data, real_data,
00767 FFTW_ESTIMATE);
00768 fftwf_execute(plan);
00769 fftwf_destroy_plan(plan);
00770 #endif // FFTW_PLAN_CACHING
00771
00772 return 0;
00773 }
00774
00775 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00776 {
00777 const int rank = get_rank(ny, nz);
00778 int dims[3];
00779 dims[0] = nz;
00780 dims[1] = ny;
00781 dims[2] = nx;
00782
00783 switch(rank) {
00784 case 1:
00785 real_to_complex_1d(real_data, complex_data, nx);
00786 break;
00787
00788 case 2:
00789 case 3:
00790 {
00791 #ifdef FFTW_PLAN_CACHING
00792 bool ip = ( complex_data == real_data );
00793 fftwf_plan plan = plan_cache.get_plan(rank,nx,ny,nz,EMAN2_REAL_2_COMPLEX,ip,(fftwf_complex *) complex_data, real_data);
00794
00795 fftwf_execute_dft_r2c(plan, real_data,(fftwf_complex *) complex_data );
00796 #else
00797 fftwf_plan plan = fftwf_plan_dft_r2c(rank, dims + (3 - rank),
00798 real_data, (fftwf_complex *) complex_data, FFTW_ESTIMATE);
00799 fftwf_execute(plan);
00800 fftwf_destroy_plan(plan);
00801 #endif // FFTW_PLAN_CACHING
00802 }
00803 break;
00804
00805 default:
00806 LOGERR("Should NEVER be here!!!");
00807 break;
00808 }
00809
00810 return 0;
00811 }
00812
00813 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00814 {
00815 const int rank = get_rank(ny, nz);
00816 int dims[3];
00817 dims[0] = nz;
00818 dims[1] = ny;
00819 dims[2] = nx;
00820
00821 switch(rank) {
00822 case 1:
00823 complex_to_real_1d(complex_data, real_data, nx);
00824 break;
00825
00826 case 2:
00827 case 3:
00828 {
00829 #ifdef FFTW_PLAN_CACHING
00830 bool ip = ( complex_data == real_data );
00831 fftwf_plan plan = plan_cache.get_plan(rank,nx,ny,nz,EMAN2_COMPLEX_2_REAL,ip,(fftwf_complex *) complex_data, real_data);
00832
00833 fftwf_execute_dft_c2r(plan, (fftwf_complex *) complex_data, real_data);
00834 #else
00835 fftwf_plan plan = fftwf_plan_dft_c2r(rank, dims + (3 - rank),
00836 (fftwf_complex *) complex_data, real_data, FFTW_ESTIMATE);
00837 fftwf_execute(plan);
00838 fftwf_destroy_plan(plan);
00839 #endif // FFTW_PLAN_CACHING
00840
00841
00842 }
00843 break;
00844
00845 default:
00846 LOGERR("Should NEVER be here!!!");
00847 break;
00848 }
00849
00850 return 0;
00851 }
00852
00853 #endif //FFTW3
00854
00855 #ifdef NATIVE_FFT
00856 #include "sparx/native_fft.h"
00857 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00858 {
00859
00860 float * work = (float*) malloc((2*n+15)*sizeof(float));
00861 if (!work) {
00862 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
00863 LOGERR("real_to_complex_1d: failed to allocate work\n");
00864 }
00865
00866 rffti(n, work);
00867 memcpy(&complex_data[1], real_data, n * sizeof(float));
00868 rfftf(n, &complex_data[1], work);
00869 complex_data[0] = complex_data[1] ;
00870 complex_data[1] = 0.0f ;
00871 if (n%2 == 0) complex_data[n+1] = 0.0f ;
00872
00873 free(work);
00874 return 0;
00875 }
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00915 {
00916
00917 float * work = (float*) malloc((2*n+15)*sizeof(float));
00918 if (!work) {
00919 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
00920 LOGERR("complex_to_real_1d: failed to allocate work\n");
00921 }
00922
00923 rffti(n, work);
00924
00925 memcpy(&real_data[1], &complex_data[2], (n-1) * sizeof(float));
00926 real_data[0] = complex_data[0] ;
00927 rfftb(n, real_data, work);
00928
00929 float nrm = 1.0f/float(n);
00930 for (int i = 0; i<n; i++) real_data[i] *= nrm;
00931 free(work);
00932 return 0;
00933 }
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00952 {
00953 const int rank = get_rank(ny, nz);
00954 const int complex_nx = nx + 2 - nx%2;
00955
00956 switch(rank) {
00957 case 1:
00958 real_to_complex_1d(real_data, complex_data, nx);
00959 return 0;
00960 case 2:
00961 {
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982 int status = Nativefft::ftp_2rf(real_data, complex_data, complex_nx, nx, ny);
00983 if (status !=0) {
00984 fprintf(stderr, "real_to_complex_nd(2df): status = %d\n", status);
00985 LOGWARN("real_to_complex_nd(2df): status = %d\n", status);
00986 }
00987 return 0;
00988 }
00989 case 3:
00990 {
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005 int status = Nativefft::ftp_3rf(real_data, complex_data, complex_nx, nx, ny, nz);
01006 if (status !=0) {
01007 fprintf(stderr, "real_to_complex_nd(3df): status = %d\n", status);
01008 LOGWARN("real_to_complex_nd(3df): status = %d\n", status);
01009 }
01010
01011
01012 return 0;
01013 }
01014 default:
01015 LOGERR("Never should be here...\n");
01016 return -1;
01017 }
01018 }
01019
01020 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
01021 {
01022 const int rank = get_rank(ny, nz);
01023 const int complex_nx = nx + 2 - nx%2;
01024
01025 switch(rank) {
01026 case 1:
01027 complex_to_real_1d(complex_data, real_data, nx);
01028 return 0;
01029 case 2:
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051 {
01052 memcpy(real_data, complex_data, complex_nx*ny*sizeof(float));
01053
01054
01055 int status = Nativefft::ftp_2rb(real_data, complex_nx, nx, ny);
01056
01057 if (status !=0) {
01058 fprintf(stderr, "complex_to_real_nd(2db): status = %d\n", status);
01059 LOGWARN("complex_to_real_nd(2db): status = %d\n", status);
01060 }
01061 return 0;
01062 }
01063 case 3:
01064 {
01065 memcpy(real_data, complex_data, complex_nx*ny*nz*sizeof(float));
01066
01067
01068 int status = Nativefft::ftp_3rb(real_data, complex_nx, nx, ny, nz);
01069 if (status !=0) {
01070 fprintf(stderr, "complex_to_real_nd(3db): status = %d\n", status);
01071 LOGWARN("complex_to_real_nd(3db): status = %d\n", status);
01072 }
01073 return 0;
01074 }
01075 default:
01076 LOGERR("Never should be here...\n");
01077 return -1;
01078 }
01079 }
01080 #endif //NATIVE_FFT
01081
01082 #ifdef ACML
01083 #include <iostream>
01084 using std::cout;
01085 using std::endl;
01086
01087 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
01088 {
01089 const int complex_n = n + 2 - n%2;
01090 int info;
01091
01092
01093 float * comm = (float *)malloc((3*n+100)*sizeof(float));
01094
01095 float * fft_data = (float *)malloc(n * sizeof(float));
01096
01097
01098 memcpy(fft_data, real_data, n * sizeof(float));
01099
01100
01101 scfft(0, n, fft_data, comm, &info);
01102 if(info != 0) {
01103 LOGERR("Error happening in Initialize communication work array: %d", info);
01104 }
01105
01106
01107 scfft(1, n, fft_data, comm, &info);
01108 if(info != 0) {
01109 LOGERR("Error happening in Initialize communication work array: %d", info);
01110 }
01111
01118 transform(fft_data, fft_data+n, fft_data, time_sqrt_n(n));
01119
01120 for(int i=0; i<complex_n; ++i) {
01121 if(i%2==0) {
01122 complex_data[i] = fft_data[i/2];
01123 }
01124 else {
01125 if(i==1) {
01126 complex_data[i] = 0.0f;
01127 }
01128 else {
01129 if(n%2 == 0 && i == complex_n-1 ) {
01130 complex_data[i] = 0.0f;
01131 }
01132 else {
01133 complex_data[i] = fft_data[n-i/2];
01134 }
01135 }
01136 }
01137 }
01138
01139 free(fft_data);
01140 free(comm);
01141 return 0;
01142 }
01143
01144 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
01145 {
01146 int complex_n = n + 2 - n%2;
01147 int info;
01148
01149
01150 float * comm = (float *)malloc((3*n+100)*sizeof(float));
01151
01152 for(int i=0; i<complex_n; ++i) {
01153 if(i%2 == 0) {
01154 real_data[i/2] = complex_data[i];
01155 }
01156 else {
01157 if(i==1) {continue;}
01158 if(!(n%2 == 0 && i == complex_n-1)) {
01159 real_data[n-i/2] = complex_data[i];
01160 }
01161 }
01162 }
01163 transform(real_data, real_data+n, real_data, divide_sqrt_n(n));
01164
01165
01166 csfft(0, n, real_data, comm, &info);
01167 if(info != 0) {
01168 LOGERR("Error happening in Initialize communication work array: %d", info);
01169 }
01170
01171
01172 for (int j = n/2+1; j < n; j++) {
01173 real_data[j] = -real_data[j];
01174 }
01175
01176
01177 csfft(1, n, real_data, comm, &info);
01178 if(info != 0) {
01179 LOGERR("Error happening in Initialize communication work array: %d", info);
01180 }
01181
01182 free(comm);
01183 return 0;
01184 }
01185
01186 int EMfft::real_to_complex_2d(float *real_data, float *complex_data, int nx, int ny)
01187 {
01188 const int complex_nx = nx + 2 - nx%2;
01189 int info;
01190
01191 float * comm = (float *)malloc((3*nx+100)*sizeof(float));
01192
01193
01194 float * fft_data = (float *)malloc(complex_nx * ny * sizeof(float));
01195 cout << "fft_data after allocation:" << endl;
01196 for(int j=0; j<ny; ++j) {
01197 for(int i=0; i<complex_nx; ++i) {
01198 cout << "data[" << i << "][" << j << "] = " << fft_data[i+j*complex_nx] << "\t";
01199 }
01200 cout << endl;
01201 }
01202
01203
01204 for(int i=0; i<ny; ++i) {
01205 memcpy(fft_data+i*complex_nx, real_data+i*nx, nx*sizeof(float));
01206 }
01207
01208 cout << endl << "real_data array: " << endl;
01209 for(int j=0; j<ny; ++j) {
01210 for(int i=0; i<nx; ++i) {
01211 cout << real_data[i+j*nx] << "\t";
01212 }
01213 cout << endl;
01214 }
01215 cout << endl << "the fft_data array: " << endl;
01216 for(int j=0; j<ny; ++j) {
01217 for(int i=0; i<complex_nx; ++i) {
01218 cout << fft_data[i+j*complex_nx] << "\t";
01219 }
01220 cout << endl;
01221 }
01222
01223
01224 scfftm(ny, nx, fft_data, comm, &info);
01225
01226 cout << endl << "the fft_data array after x dim transform: " << endl;
01227 for(int j=0; j<ny; ++j) {
01228 for(int i=0; i<complex_nx; ++i) {
01229 cout << fft_data[i+j*complex_nx] << "\t";
01230 }
01231 cout << endl;
01232 }
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244 float * fft_data2 = (float *)malloc(complex_nx * ny * sizeof(float));
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255 if(info != 0) {
01256 LOGERR("Error happening in scfftm: %d", info);
01257 }
01258
01259 return 0;
01260 }
01261
01262 int EMfft::complex_to_real_2d(float *complex_data, float *real_data, int nx, int ny)
01263 {
01264 return 0;
01265 }
01266
01267 int EMfft::real_to_complex_3d(float *real_data, float *complex_data, int nx, int ny, int nz)
01268 {
01269 return 0;
01270 }
01271
01272 int EMfft::complex_to_real_3d(float *complex_data, float *real_data, int nx, int ny, int nz)
01273 {
01274 return 0;
01275 }
01276
01277 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
01278 {
01279 const int rank = get_rank(ny, nz);
01280
01281 switch(rank) {
01282 case 1:
01283 return real_to_complex_1d(real_data, complex_data, nx);
01284 case 2:
01285 return real_to_complex_2d(real_data, complex_data, nx, ny);
01286 case 3:
01287 return real_to_complex_3d(real_data, complex_data, nx, ny, nz);
01288 default:
01289 LOGERR("Never should be here...\n");
01290 return -1;
01291 }
01292 }
01293
01294 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
01295 {
01296 const int rank = get_rank(ny, nz);
01297
01298 switch(rank) {
01299 case 1:
01300 return complex_to_real_1d(complex_data, real_data, nx);
01301 case 2:
01302 return complex_to_real_2d(complex_data, real_data, nx, ny);
01303 case 3:
01304 return complex_to_real_3d(complex_data, real_data, nx, ny, nz);
01305 default:
01306 LOGERR("Never should be here...\n");
01307 return -1;
01308 }
01309 }
01310
01311
01312 #endif //ACML