00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00027
00028 #include <math.h>
00029 #include <stdlib.h>
00030 #include <string.h>
00031 #include <complex.h>
00032
00033
00034 #include "nfft3util.h"
00035
00036
00037 #include "nfft3.h"
00038
00039 #include "infft.h"
00040
00041
00042 #include "legendre.h"
00043
00044
00045 #include "api.h"
00046
00047
00057 #define NFSFT_DEFAULT_NFFT_CUTOFF 6
00058
00064 #define NFSFT_DEFAULT_THRESHOLD 1000
00065
00071 #define NFSFT_BREAK_EVEN 5
00072
00079 static struct nfsft_wisdom wisdom = {false,0U,-1,-1,0,0,0,0,0};
00080
00103 static inline void c2e(nfsft_plan *plan)
00104 {
00105 int k;
00106 int n;
00107 double _Complex last;
00108 double _Complex act;
00109 double _Complex *xp;
00110 double _Complex *xm;
00111 int low;
00112 int up;
00113 int lowe;
00114 int upe;
00116
00117 memset(plan->f_hat_intern,0U,(2*plan->N+2)*sizeof(double _Complex));
00118
00119
00120 lowe = -plan->N + (plan->N%2);
00121 upe = -lowe;
00122
00123
00124 for (n = lowe; n <= upe; n += 2)
00125 {
00126
00127
00128 xm = &(plan->f_hat_intern[NFSFT_INDEX(-1,n,plan)]);
00129 xp = &(plan->f_hat_intern[NFSFT_INDEX(+1,n,plan)]);
00130 for(k = 1; k <= plan->N; k++)
00131 {
00132 *xp *= 0.5;
00133 *xm-- = *xp++;
00134 }
00135
00136
00137 *xm = 0.0;
00138 }
00139
00140
00141 low = -plan->N + (1-plan->N%2);
00142 up = -low;
00143
00144
00145
00146 for (n = low; n <= up; n += 2)
00147 {
00148
00149
00150
00151 plan->f_hat_intern[NFSFT_INDEX(0,n,plan)] *= 2.0;
00152 xp = &(plan->f_hat_intern[NFSFT_INDEX(-plan->N-1,n,plan)]);
00153
00154
00155 *xp++ = 0.0;
00156 xm = &(plan->f_hat_intern[NFSFT_INDEX(plan->N,n,plan)]);
00157 last = *xm;
00158 *xm = 0.5 * _Complex_I * (0.5*xm[-1]);
00159 *xp++ = -(*xm--);
00160 for (k = plan->N-1; k > 0; k--)
00161 {
00162 act = *xm;
00163 *xm = 0.5 * _Complex_I * (0.5*(xm[-1] - last));
00164 *xp++ = -(*xm--);
00165 last = act;
00166 }
00167 *xm = 0.0;
00168 }
00169 }
00170
00181 static inline void c2e_transposed(nfsft_plan *plan)
00182 {
00183 int k;
00184 int n;
00185 double _Complex last;
00186 double _Complex act;
00187 double _Complex *xp;
00188 double _Complex *xm;
00189 int low;
00190 int up;
00191 int lowe;
00192 int upe;
00194
00195 lowe = -plan->N + (plan->N%2);
00196 upe = -lowe;
00197
00198
00199 for (n = lowe; n <= upe; n += 2)
00200 {
00201
00202
00203 xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]);
00204 xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]);
00205 for(k = 1; k <= plan->N; k++)
00206 {
00207 *xp += *xm--;
00208 *xp++ *= 0.5;;
00209 }
00210 }
00211
00212
00213 low = -plan->N + (1-plan->N%2);
00214 up = -low;
00215
00216
00217 for (n = low; n <= up; n += 2)
00218 {
00219
00220
00221 xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]);
00222 xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]);
00223 for(k = 1; k <= plan->N; k++)
00224 {
00225 *xp++ -= *xm--;
00226 }
00227
00228 plan->f_hat[NFSFT_INDEX(0,n,plan)] =
00229 -0.25*_Complex_I*plan->f_hat[NFSFT_INDEX(1,n,plan)];
00230 last = plan->f_hat[NFSFT_INDEX(1,n,plan)];
00231 plan->f_hat[NFSFT_INDEX(1,n,plan)] =
00232 -0.25*_Complex_I*plan->f_hat[NFSFT_INDEX(2,n,plan)];
00233
00234 xp = &(plan->f_hat[NFSFT_INDEX(2,n,plan)]);
00235 for (k = 2; k < plan->N; k++)
00236 {
00237 act = *xp;
00238 *xp = -0.25 * _Complex_I * (xp[1] - last);
00239 xp++;
00240 last = act;
00241 }
00242 *xp = 0.25 * _Complex_I * last;
00243
00244 plan->f_hat[NFSFT_INDEX(0,n,plan)] *= 2.0;
00245 }
00246 }
00247
00248 void nfsft_init(nfsft_plan *plan, int N, int M)
00249 {
00250
00251 nfsft_init_advanced(plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
00252 NFSFT_MALLOC_F_HAT);
00253 }
00254
00255 void nfsft_init_advanced(nfsft_plan* plan, int N, int M,
00256 unsigned int flags)
00257 {
00258
00259 nfsft_init_guru(plan, N, M, flags, PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00260 FFT_OUT_OF_PLACE, NFSFT_DEFAULT_NFFT_CUTOFF);
00261 }
00262
00263 void nfsft_init_guru(nfsft_plan *plan, int N, int M, unsigned int flags,
00264 unsigned int nfft_flags, int nfft_cutoff)
00265 {
00266 int *nfft_size;
00267 int *fftw_size;
00268
00269
00270 plan->flags = flags;
00271
00272
00273 plan->N = N;
00274 plan->M_total = M;
00275
00276
00277
00278
00279
00280
00281
00282 plan->N_total = (2*plan->N+2)*(2*plan->N+2);
00283
00284
00285
00286 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00287 {
00288 plan->f_hat_intern = (double _Complex*) nfft_malloc(plan->N_total*
00289 sizeof(double _Complex));
00290 }
00291
00292
00293 if (plan->flags & NFSFT_MALLOC_F_HAT)
00294 {
00295 plan->f_hat = (double _Complex*) nfft_malloc(plan->N_total*
00296 sizeof(double _Complex));
00297 }
00298
00299
00300 if (plan->flags & NFSFT_MALLOC_F)
00301 {
00302 plan->f = (double _Complex*) nfft_malloc(plan->M_total*sizeof(double _Complex));
00303 }
00304
00305
00306 if (plan->flags & NFSFT_MALLOC_X)
00307 {
00308 plan->x = (double*) nfft_malloc(plan->M_total*2*sizeof(double));
00309 }
00310
00311
00312 if (plan->flags & NFSFT_NO_FAST_ALGORITHM)
00313 {
00314 }
00315 else
00316 {
00317 nfft_size = (int*)nfft_malloc(2*sizeof(int));
00318 fftw_size = (int*)nfft_malloc(2*sizeof(int));
00319
00321 nfft_size[0] = 2*plan->N+2;
00322 nfft_size[1] = 2*plan->N+2;
00323 fftw_size[0] = 4*plan->N;
00324 fftw_size[1] = 4*plan->N;
00325
00327 nfft_init_guru(&plan->plan_nfft, 2, nfft_size, plan->M_total, fftw_size,
00328 nfft_cutoff, nfft_flags,
00329 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
00330
00331
00332 plan->plan_nfft.x = plan->x;
00333
00334 plan->plan_nfft.f = plan->f;
00335
00336 plan->plan_nfft.f_hat = plan->f_hat;
00337
00340
00341
00342
00343
00344 nfft_free(nfft_size);
00345 nfft_free(fftw_size);
00346 }
00347
00348 plan->mv_trafo = (void (*) (void* ))nfsft_trafo;
00349 plan->mv_adjoint = (void (*) (void* ))nfsft_adjoint;
00350 }
00351
00352 void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags,
00353 unsigned int fpt_flags)
00354 {
00355 int n;
00356
00357
00358 if (wisdom.initialized == true)
00359 {
00360 return;
00361 }
00362
00363
00364 wisdom.flags = nfsft_flags;
00365
00366
00367
00368 nfft_next_power_of_2_exp(N,&wisdom.N_MAX,&wisdom.T_MAX);
00369
00370
00371 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00372 {
00373 wisdom.alpha = NULL;
00374 wisdom.beta = NULL;
00375 wisdom.gamma = NULL;
00376 }
00377 else
00378 {
00379
00380 wisdom.alpha = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
00381 sizeof(double));
00382 wisdom.beta = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
00383 sizeof(double));
00384 wisdom.gamma = (double*) nfft_malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
00385 sizeof(double));
00387
00388
00389 alpha_al_all(wisdom.alpha,wisdom.N_MAX);
00390 beta_al_all(wisdom.beta,wisdom.N_MAX);
00391 gamma_al_all(wisdom.gamma,wisdom.N_MAX);
00392 }
00393
00394
00395 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00396 {
00397 }
00398 else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN)
00399 {
00400
00401
00402
00403 if (wisdom.alpha != NULL)
00404 {
00405
00406
00407 wisdom.set = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX,
00408 fpt_flags | FPT_AL_SYMMETRY | FPT_PERSISTENT_DATA);
00409 for (n = 0; n <= wisdom.N_MAX; n++)
00410 {
00411
00412
00413
00414 fpt_precompute(wisdom.set,n,&wisdom.alpha[ROW(n)],&wisdom.beta[ROW(n)],
00415 &wisdom.gamma[ROW(n)],n,kappa);
00416 }
00417 }
00418 else
00419 {
00420
00421 wisdom.alpha = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
00422 wisdom.beta = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
00423 wisdom.gamma = (double*) nfft_malloc((wisdom.N_MAX+2)*sizeof(double));
00424 wisdom.set = fpt_init(wisdom.N_MAX+1, wisdom.T_MAX,
00425 fpt_flags | FPT_AL_SYMMETRY);
00426 for (n = 0; n <= wisdom.N_MAX; n++)
00427 {
00428
00429
00430
00431
00432 alpha_al_row(wisdom.alpha,wisdom.N_MAX,n);
00433 beta_al_row(wisdom.beta,wisdom.N_MAX,n);
00434 gamma_al_row(wisdom.gamma,wisdom.N_MAX,n);
00435
00436
00437 fpt_precompute(wisdom.set,n,wisdom.alpha,wisdom.beta,wisdom.gamma,n,
00438 kappa);
00439 }
00440
00441 nfft_free(wisdom.alpha);
00442 nfft_free(wisdom.beta);
00443 nfft_free(wisdom.gamma);
00444 wisdom.alpha = NULL;
00445 wisdom.beta = NULL;
00446 wisdom.gamma = NULL;
00447 }
00448 }
00449
00450
00451 wisdom.initialized = true;
00452 }
00453
00454 void nfsft_forget(void)
00455 {
00456
00457 if (wisdom.initialized == false)
00458 {
00459
00460 return;
00461 }
00462
00463
00464 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00465 {
00466 }
00467 else
00468 {
00469
00470 nfft_free(wisdom.alpha);
00471 nfft_free(wisdom.beta);
00472 nfft_free(wisdom.gamma);
00473 wisdom.alpha = NULL;
00474 wisdom.beta = NULL;
00475 wisdom.gamma = NULL;
00476 }
00477
00478
00479 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00480 {
00481 }
00482 else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN)
00483 {
00484
00485 fpt_finalize(wisdom.set);
00486 }
00487
00488
00489 wisdom.initialized = false;
00490 }
00491
00492
00493 void nfsft_finalize(nfsft_plan *plan)
00494 {
00495 if (!plan)
00496 return;
00497
00498
00499 nfft_finalize(&plan->plan_nfft);
00500
00501
00502
00503 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00504 {
00505 nfft_free(plan->f_hat_intern);
00506 }
00507
00508
00509 if (plan->flags & NFSFT_MALLOC_F_HAT)
00510 {
00511
00512 nfft_free(plan->f_hat);
00513 }
00514
00515
00516 if (plan->flags & NFSFT_MALLOC_F)
00517 {
00518
00519 nfft_free(plan->f);
00520 }
00521
00522
00523 if (plan->flags & NFSFT_MALLOC_X)
00524 {
00525
00526 nfft_free(plan->x);
00527 }
00528 }
00529
00530 void ndsft_trafo(nfsft_plan *plan)
00531 {
00532 int m;
00533 int k;
00534 int n;
00535 int n_abs;
00536 double *alpha;
00537
00538
00539 double *gamma;
00540
00541
00542 double _Complex *a;
00543 double _Complex it1;
00544 double _Complex it2;
00545 double _Complex temp;
00546 double _Complex f_m;
00547
00548 double stheta;
00549 double sphi;
00550
00551 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00552 {
00553 return;
00554 }
00555
00556
00557 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00558 {
00559 memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total*
00560 sizeof(double _Complex));
00561 }
00562 else
00563 {
00564 plan->f_hat_intern = plan->f_hat;
00565 }
00566
00567
00568
00569
00570 if (plan->flags & NFSFT_NORMALIZED)
00571 {
00572
00573 for (k = 0; k <= plan->N; k++)
00574 {
00575 for (n = -k; n <= k; n++)
00576 {
00577
00578 plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *=
00579 sqrt((2*k+1)/(4.0*PI));
00580 }
00581 }
00582 }
00583
00584
00585 if (plan->N == 0)
00586 {
00587
00588
00589
00590 for (m = 0; m < plan->M_total; m++)
00591 {
00592 plan->f[m] = plan->f_hat_intern[NFSFT_INDEX(0,0,plan)];
00593 }
00594 }
00595 else
00596 {
00597
00598
00599
00600
00601
00602
00603
00604 for (m = 0; m < plan->M_total; m++)
00605 {
00606
00607 stheta = cos(2.0*PI*plan->x[2*m+1]);
00608
00609 sphi = 2.0*PI*plan->x[2*m];
00610
00611
00612 f_m = 0.0;
00613
00614
00615
00616
00617
00618 for (n = -plan->N; n <= plan->N; n++)
00619 {
00620
00621 a = &(plan->f_hat_intern[NFSFT_INDEX(0,n,plan)]);
00622
00623
00624 n_abs = abs(n);
00625
00626
00627 alpha = &(wisdom.alpha[ROW(n_abs)]);
00628 gamma = &(wisdom.gamma[ROW(n_abs)]);
00629
00630
00631 it2 = a[plan->N];
00632 it1 = a[plan->N-1];
00633 for (k = plan->N; k > n_abs + 1; k--)
00634 {
00635 temp = a[k-2] + it2 * gamma[k];
00636 it2 = it1 + it2 * alpha[k] * stheta;
00637 it1 = temp;
00638 }
00639
00640
00641 if (n_abs < plan->N)
00642 {
00643 it2 = it1 + it2 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
00644 }
00645
00646
00647
00648
00649
00650
00651
00652 f_m += it2 * wisdom.gamma[ROW(n_abs)] *
00653 pow(1- stheta * stheta, 0.5*n_abs) * cexp(_Complex_I*n*sphi);
00654 }
00655
00656
00657 plan->f[m] = f_m;
00658 }
00659 }
00660 }
00661
00662 void ndsft_adjoint(nfsft_plan *plan)
00663 {
00664 int m;
00665 int k;
00666 int n;
00667 int n_abs;
00668 double *alpha;
00669
00670
00671 double *gamma;
00672
00673
00674 double _Complex it1;
00675 double _Complex it2;
00676 double _Complex temp;
00677 double stheta;
00678 double sphi;
00679
00680 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00681 {
00682 return;
00683 }
00684
00685
00686 memset(plan->f_hat,0U,plan->N_total*sizeof(double _Complex));
00687
00688
00689 if (plan->N == 0)
00690 {
00691
00692
00693
00694 for (m = 0; m < plan->M_total; m++)
00695 {
00696 plan->f_hat[NFSFT_INDEX(0,0,plan)] += plan->f[m];
00697 }
00698 }
00699 else
00700 {
00701
00702
00703
00704 for (m = 0; m < plan->M_total; m++)
00705 {
00706
00707 stheta = cos(2.0*PI*plan->x[2*m+1]);
00708
00709 sphi = 2.0*PI*plan->x[2*m];
00710
00711
00712 for (n = -plan->N; n <= plan->N; n++)
00713 {
00714
00715 n_abs = abs(n);
00716
00717
00718 alpha = &(wisdom.alpha[ROW(n_abs)]);
00719 gamma = &(wisdom.gamma[ROW(n_abs)]);
00720
00721
00722
00723
00724 it1 = plan->f[m] * wisdom.gamma[ROW(n_abs)] *
00725 pow(1 - stheta * stheta, 0.5*n_abs) * cexp(-_Complex_I*n*sphi);
00726 plan->f_hat[NFSFT_INDEX(n_abs,n,plan)] += it1;
00727 it2 = 0.0;
00728
00729 if (n_abs < plan->N)
00730 {
00731 it2 = it1 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
00732 plan->f_hat[NFSFT_INDEX(n_abs+1,n,plan)] += it2;
00733 }
00734
00735
00736 for (k = n_abs+2; k <= plan->N; k++)
00737 {
00738 temp = it2;
00739 it2 = alpha[k] * stheta * it2 + gamma[k] * it1;
00740 it1 = temp;
00741 plan->f_hat[NFSFT_INDEX(k,n,plan)] += it2;
00742 }
00743 }
00744 }
00745 }
00746
00747
00748
00749
00750 if (plan->flags & NFSFT_NORMALIZED)
00751 {
00752
00753 for (k = 0; k <= plan->N; k++)
00754 {
00755 for (n = -k; n <= k; n++)
00756 {
00757
00758 plan->f_hat[NFSFT_INDEX(k,n,plan)] *=
00759 sqrt((2*k+1)/(4.0*PI));
00760 }
00761 }
00762 }
00763
00764
00765 if (plan->flags & NFSFT_ZERO_F_HAT)
00766 {
00767 for (n = -plan->N; n <= plan->N+1; n++)
00768 {
00769 memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U,
00770 (plan->N+1+abs(n))*sizeof(double _Complex));
00771 }
00772 }
00773 }
00774
00775 void nfsft_trafo(nfsft_plan *plan)
00776 {
00777 int k;
00778 int n;
00779 #ifdef DEBUG
00780 double t, t_pre, t_nfft, t_fpt, t_c2e, t_norm;
00781 t_pre = 0.0;
00782 t_norm = 0.0;
00783 t_fpt = 0.0;
00784 t_c2e = 0.0;
00785 t_nfft = 0.0;
00786 #endif
00787
00788 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00789 {
00790 return;
00791 }
00792
00793
00794
00795
00796 if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX)
00797 {
00798 return;
00799 }
00800
00801
00802 if (plan->N < NFSFT_BREAK_EVEN)
00803 {
00804
00805 ndsft_trafo(plan);
00806 }
00807
00808
00809 else if (plan->N <= wisdom.N_MAX)
00810 {
00811
00812 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00813 {
00814 memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total*
00815 sizeof(double _Complex));
00816 }
00817 else
00818 {
00819 plan->f_hat_intern = plan->f_hat;
00820 }
00821
00822
00823
00824
00825 plan->plan_nfft.x = plan->x;
00826 plan->plan_nfft.f = plan->f;
00827 plan->plan_nfft.f_hat = plan->f_hat_intern;
00828
00829
00830
00831
00832 if (plan->flags & NFSFT_NORMALIZED)
00833 {
00834
00835 for (k = 0; k <= plan->N; k++)
00836 {
00837 for (n = -k; n <= k; n++)
00838 {
00839
00840 plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *=
00841 sqrt((2*k+1)/(4.0*PI));
00842 }
00843 }
00844 }
00845
00846
00847 if (plan->flags & NFSFT_USE_DPT)
00848 {
00849
00850 for (n = -plan->N; n <= plan->N; n++)
00851 {
00852
00853 fflush(stderr);
00854 dpt_trafo(wisdom.set,abs(n),
00855 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
00856 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
00857 plan->N,0U);
00858 }
00859 }
00860 else
00861 {
00862
00863 for (n = -plan->N; n <= plan->N; n++)
00864 {
00865
00866 fflush(stderr);
00867 fpt_trafo(wisdom.set,abs(n),
00868 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
00869 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
00870 plan->N,0U);
00871 }
00872 }
00873
00874
00875 c2e(plan);
00876
00877
00878
00879
00880 if (plan->flags & NFSFT_USE_NDFT)
00881 {
00882
00883 ndft_trafo(&plan->plan_nfft);
00884 }
00885 else
00886 {
00887
00888
00889 nfft_trafo_2d(&plan->plan_nfft);
00890 }
00891 }
00892 }
00893
00894 void nfsft_adjoint(nfsft_plan *plan)
00895 {
00896 int k;
00897 int n;
00898
00899 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00900 {
00901 return;
00902 }
00903
00904
00905
00906
00907 if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX)
00908 {
00909 return;
00910 }
00911
00912
00913 if (plan->N < NFSFT_BREAK_EVEN)
00914 {
00915
00916 ndsft_adjoint(plan);
00917 }
00918
00919 else if (plan->N <= wisdom.N_MAX)
00920 {
00921
00922
00923
00924
00925
00926 plan->plan_nfft.x = plan->x;
00927 plan->plan_nfft.f = plan->f;
00928 plan->plan_nfft.f_hat = plan->f_hat;
00929
00930
00931
00932
00933 if (plan->flags & NFSFT_USE_NDFT)
00934 {
00935
00936
00937
00938 ndft_adjoint(&plan->plan_nfft);
00939 }
00940 else
00941 {
00942
00943
00944
00945
00946 nfft_adjoint_2d(&plan->plan_nfft);
00947 }
00948
00949
00950
00951
00952 c2e_transposed(plan);
00953
00954
00955 if (plan->flags & NFSFT_USE_DPT)
00956 {
00957
00958 for (n = -plan->N; n <= plan->N; n++)
00959 {
00960
00961
00962 dpt_transposed(wisdom.set,abs(n),
00963 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
00964 &plan->f_hat[NFSFT_INDEX(0,n,plan)],
00965 plan->N,0U);
00966 }
00967 }
00968 else
00969 {
00970
00971
00972 for (n = -plan->N; n <= plan->N; n++)
00973 {
00974
00975
00976 fpt_transposed(wisdom.set,abs(n),
00977 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
00978 &plan->f_hat[NFSFT_INDEX(0,n,plan)],
00979 plan->N,0U);
00980 }
00981 }
00982
00983
00984
00985
00986 if (plan->flags & NFSFT_NORMALIZED)
00987 {
00988
00989
00990
00991 for (k = 0; k <= plan->N; k++)
00992 {
00993 for (n = -k; n <= k; n++)
00994 {
00995
00996 plan->f_hat[NFSFT_INDEX(k,n,plan)] *=
00997 sqrt((2*k+1)/(4.0*PI));
00998 }
00999 }
01000 }
01001
01002
01003 if (plan->flags & NFSFT_ZERO_F_HAT)
01004 {
01005
01006
01007 for (n = -plan->N; n <= plan->N+1; n++)
01008 {
01009 memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U,
01010 (plan->N+1+abs(n))*sizeof(double _Complex));
01011 }
01012 }
01013
01014
01015 }
01016 }
01017
01018 void nfsft_precompute_x(nfsft_plan *plan)
01019 {
01020
01021 plan->plan_nfft.x = plan->x;
01022
01023
01024 if (plan->plan_nfft.nfft_flags & PRE_ONE_PSI)
01025 nfft_precompute_one_psi(&plan->plan_nfft);
01026 }
01027