00001
00006 #include <stdio.h>
00007 #include <math.h>
00008 #include <string.h>
00009 #include <stdlib.h>
00010
00011 #include "util.h"
00012 #include "options.h"
00013 #include "window_defines.h"
00014
00015 #include "nfft3.h"
00016
00040 #define MACRO_ndft_init_result_trafo memset(f,0,ths->M_total* \
00041 sizeof(double complex));
00042 #define MACRO_ndft_init_result_conjugated MACRO_ndft_init_result_trafo
00043 #define MACRO_ndft_init_result_adjoint memset(f_hat,0,ths->N_total* \
00044 sizeof(double complex));
00045 #define MACRO_ndft_init_result_transposed MACRO_ndft_init_result_adjoint
00046
00047 #define MACRO_ndft_sign_trafo +2*PI*ths->x[j*ths->d+t]
00048 #define MACRO_ndft_sign_conjugated -2*PI*ths->x[j*ths->d+t]
00049 #define MACRO_ndft_sign_adjoint +2*PI*ths->x[j*ths->d+t]
00050 #define MACRO_ndft_sign_transposed -2*PI*ths->x[j*ths->d+t]
00051
00052 #define MACRO_init_k_N_Omega_x(which_one) { \
00053 for(t=0; t<ths->d; t++) \
00054 { \
00055 k[t]=-ths->N[t]/2; \
00056 x[t]= MACRO_ndft_sign_ ## which_one; \
00057 Omega[t+1]=k[t]*x[t]+Omega[t]; \
00058 } \
00059 omega=Omega[ths->d]; \
00060 } \
00061
00062 #define MACRO_count_k_N_Omega { \
00063 for(t = ths->d-1; (t >= 1) && (k[t] == ths->N[t]/2-1); t--) \
00064 k[t]-= ths->N[t]-1; \
00065 \
00066 k[t]++; \
00067 \
00068 for(t2 = t; t2<ths->d; t2++) \
00069 Omega[t2+1]=k[t2]*x[t2]+Omega[t2]; \
00070 \
00071 omega=Omega[ths->d]; \
00072 }
00073
00074 #define MACRO_ndft_compute_trafo (*fj) += (*f_hat_k)*cexp(-I*omega);
00075
00076 #define MACRO_ndft_compute_conjugated MACRO_ndft_compute_trafo
00077
00078 #define MACRO_ndft_compute_adjoint (*f_hat_k) += (*fj)*cexp(+I*omega);
00079
00080 #define MACRO_ndft_compute_transposed MACRO_ndft_compute_adjoint
00081
00082 #define MACRO_ndft(which_one) \
00083 void ndft_ ## which_one (nfft_plan *ths) \
00084 { \
00085 int j; \
00086 int t,t2; \
00087 int k_L; \
00088 double complex *f_hat, *f; \
00089 double complex *f_hat_k; \
00090 double complex *fj; \
00091 double x[ths->d]; \
00092 int k[ths->d]; \
00093 double omega, Omega[ths->d+1]; \
00094 \
00095 f_hat=ths->f_hat; f=ths->f; \
00096 \
00097 MACRO_ndft_init_result_ ## which_one \
00098 \
00099 if(ths->d==1) \
00100 { \
00101 t=0; \
00102 for(j=0, fj = f; j<ths->M_total; j++, fj++) \
00103 { \
00104 for(k_L=0, f_hat_k = f_hat; k_L<ths->N_total; k_L++, f_hat_k++) \
00105 { \
00106 omega=(k_L-ths->N_total/2)* MACRO_ndft_sign_ ## which_one; \
00107 MACRO_ndft_compute_ ## which_one; \
00108 } \
00109 } \
00110 } \
00111 else \
00112 { \
00113 Omega[0]=0; \
00114 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00115 { \
00116 MACRO_init_k_N_Omega_x(which_one); \
00117 for(k_L=0, f_hat_k=f_hat; k_L<ths->N_total; k_L++, f_hat_k++) \
00118 { \
00119 MACRO_ndft_compute_ ## which_one; \
00120 MACRO_count_k_N_Omega; \
00121 } \
00122 } \
00123 } \
00124 }
00125
00126
00129 MACRO_ndft(trafo)
00130 MACRO_ndft(adjoint)
00131
00157 void nfft_uo(nfft_plan *ths,int j,int *up,int *op,int act_dim)
00158 {
00159 double c;
00160 int u,o;
00161
00162 c = ths->x[j*ths->d+act_dim] * ths->n[act_dim];
00163 u = c; o = c;
00164
00165 if(c < 0)
00166 u = u-1;
00167 else
00168 o = o+1;
00169
00170 u = u - (ths->m); o = o + (ths->m);
00171
00172 up[0]=u; op[0]=o;
00173 }
00174
00175 #define MACRO_nfft_D_compute_A { \
00176 g_hat[k_plain[ths->d]] = f_hat[ks_plain[ths->d]] * c_phi_inv_k[ths->d]; \
00177 }
00178
00179 #define MACRO_nfft_D_compute_T { \
00180 f_hat[ks_plain[ths->d]] = g_hat[k_plain[ths->d]] * c_phi_inv_k[ths->d]; \
00181 }
00182
00183 #define MACRO_nfft_D_init_result_A memset(g_hat,0,ths->n_total* \
00184 sizeof(double complex));
00185 #define MACRO_nfft_D_init_result_T memset(f_hat,0,ths->N_total* \
00186 sizeof(double complex));
00187
00188 #define MACRO_with_PRE_PHI_HUT * ths->c_phi_inv[t2][ks[t2]];
00189 #define MACRO_without_PRE_PHI_HUT / (PHI_HUT(ks[t2]-ths->N[t2]/2,t2));
00190
00191 #define MACRO_init_k_ks { \
00192 for(t = ths->d-1; t>=0; t--) \
00193 { \
00194 kp[t]= 0; \
00195 k[t] = 0; \
00196 ks[t] = ths->N[t]/2; \
00197 } \
00198 t++; \
00199 }
00200
00201 #define MACRO_update_c_phi_inv_k(which_one) { \
00202 for(t2=t; t2<ths->d; t2++) \
00203 { \
00204 c_phi_inv_k[t2+1]= c_phi_inv_k[t2] MACRO_ ##which_one; \
00205 ks_plain[t2+1]= ks_plain[t2]*ths->N[t2]+ks[t2]; \
00206 k_plain[t2+1]= k_plain[t2]*ths->n[t2]+k[t2]; \
00207 } \
00208 }
00209
00210 #define MACRO_count_k_ks { \
00211 for(t=ths->d-1; (t>0)&& (kp[t]==ths->N[t]-1); t--) \
00212 { \
00213 kp[t]= 0; \
00214 k[t]= 0; \
00215 ks[t]= ths->N[t]/2; \
00216 } \
00217 \
00218 kp[t]++; k[t]++; ks[t]++; \
00219 if(kp[t]==ths->N[t]/2) \
00220 { \
00221 k[t]= ths->n[t]-ths->N[t]/2; \
00222 ks[t]= 0; \
00223 } \
00224 } \
00225
00226
00230 #define MACRO_nfft_D(which_one) \
00231 inline void nfft_D_ ## which_one (nfft_plan *ths) \
00232 { \
00233 int t, t2; \
00234 int k_L; \
00235 int kp[ths->d]; \
00236 int k[ths->d]; \
00237 int ks[ths->d]; \
00238 double c_phi_inv_k[ths->d+1]; \
00239 int k_plain[ths->d+1]; \
00240 int ks_plain[ths->d+1]; \
00241 double complex *f_hat, *g_hat; \
00242 \
00243 f_hat=ths->f_hat; g_hat=ths->g_hat; \
00244 MACRO_nfft_D_init_result_ ## which_one; \
00245 \
00246 c_phi_inv_k[0]=1; \
00247 k_plain[0]=0; \
00248 ks_plain[0]=0; \
00249 \
00250 if(ths->nfft_flags & PRE_PHI_HUT) \
00251 { \
00252 MACRO_init_k_ks; \
00253 \
00254 for(k_L=0; k_L<ths->N_total; k_L++) \
00255 { \
00256 MACRO_update_c_phi_inv_k(with_PRE_PHI_HUT); \
00257 \
00258 MACRO_nfft_D_compute_ ## which_one; \
00259 \
00260 MACRO_count_k_ks; \
00261 } \
00262 } \
00263 else \
00264 { \
00265 MACRO_init_k_ks; \
00266 \
00267 for(k_L=0; k_L<ths->N_total; k_L++) \
00268 { \
00269 MACRO_update_c_phi_inv_k(without_PRE_PHI_HUT); \
00270 \
00271 MACRO_nfft_D_compute_ ## which_one; \
00272 \
00273 MACRO_count_k_ks; \
00274 } \
00275 } \
00276 }
00277
00278 MACRO_nfft_D(A)
00279 MACRO_nfft_D(T)
00280
00284 #define MACRO_nfft_B_init_result_A memset(f,0,ths->M_total* \
00285 sizeof(double complex));
00286 #define MACRO_nfft_B_init_result_T memset(g,0,ths->n_total* \
00287 sizeof(double complex));
00288
00289 #define MACRO_nfft_B_PRE_FULL_PSI_compute_A { \
00290 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
00291 }
00292
00293 #define MACRO_nfft_B_PRE_FULL_PSI_compute_T { \
00294 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
00295 }
00296
00297 #define MACRO_nfft_B_compute_A { \
00298 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
00299 }
00300
00301 #define MACRO_nfft_B_compute_T { \
00302 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
00303 }
00304
00305 #define MACRO_with_FG_PSI fg_psi[t2][lj[t2]]
00306
00307 #define MACRO_with_PRE_PSI ths->psi[(j*ths->d+t2)*(2*ths->m+2)+lj[t2]]
00308
00309 #define MACRO_without_PRE_PSI PHI(ths->x[j*ths->d+t2]- \
00310 ((double)l[t2])/ths->n[t2], t2)
00311
00312 #define MACRO_init_uo_l_lj_t { \
00313 for(t = ths->d-1; t>=0; t--) \
00314 { \
00315 nfft_uo(ths,j,&u[t],&o[t],t); \
00316 l[t] = u[t]; \
00317 lj[t] = 0; \
00318 } \
00319 t++; \
00320 }
00321
00322 #define MACRO_update_phi_prod_ll_plain(which_one) { \
00323 for(t2=t; t2<ths->d; t2++) \
00324 { \
00325 phi_prod[t2+1]=phi_prod[t2]* MACRO_ ## which_one; \
00326 ll_plain[t2+1]=ll_plain[t2]*ths->n[t2] +(l[t2]+ths->n[t2])%ths->n[t2]; \
00327 } \
00328 }
00329
00330 #define MACRO_count_uo_l_lj_t { \
00331 for(t = ths->d-1; (t>0)&&(l[t]==o[t]); t--) \
00332 { \
00333 l[t] = u[t]; \
00334 lj[t] = 0; \
00335 } \
00336 \
00337 l[t]++; \
00338 lj[t]++; \
00339 }
00340
00341 #define MACRO_nfft_B(which_one) \
00342 inline void nfft_B_ ## which_one (nfft_plan *ths) \
00343 { \
00344 int lprod; \
00345 int u[ths->d], o[ths->d]; \
00346 int t, t2; \
00347 int j; \
00348 int l_L, ix; \
00349 int l[ths->d]; \
00350 int lj[ths->d]; \
00351 int ll_plain[ths->d+1]; \
00352 double phi_prod[ths->d+1]; \
00353 double complex *f, *g; \
00354 double complex *fj; \
00355 double y[ths->d]; \
00356 double fg_psi[ths->d][2*ths->m+2]; \
00357 double fg_exp_l[ths->d][2*ths->m+2]; \
00358 int l_fg,lj_fg; \
00359 double tmpEXP1, tmpEXP2, tmpEXP2sq, tmp1, tmp2, tmp3; \
00360 double ip_w; \
00361 int ip_u; \
00362 int ip_s=ths->K/(ths->m+1); \
00363 \
00364 f=ths->f; g=ths->g; \
00365 \
00366 MACRO_nfft_B_init_result_ ## which_one; \
00367 \
00368 if(ths->nfft_flags & PRE_FULL_PSI) \
00369 { \
00370 for(ix=0, j=0, fj=f; j<ths->M_total; j++, fj++) \
00371 for(l_L=0; l_L<ths->psi_index_f[j]; l_L++, ix++) \
00372 MACRO_nfft_B_PRE_FULL_PSI_compute_ ## which_one; \
00373 return; \
00374 } \
00375 \
00376 phi_prod[0]=1; \
00377 ll_plain[0]=0; \
00378 \
00379 for(t=0,lprod = 1; t<ths->d; t++) \
00380 lprod *= (2*ths->m+2); \
00381 \
00382 if(ths->nfft_flags & PRE_PSI) \
00383 { \
00384 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00385 { \
00386 MACRO_init_uo_l_lj_t; \
00387 \
00388 for(l_L=0; l_L<lprod; l_L++) \
00389 { \
00390 MACRO_update_phi_prod_ll_plain(with_PRE_PSI); \
00391 \
00392 MACRO_nfft_B_compute_ ## which_one; \
00393 \
00394 MACRO_count_uo_l_lj_t; \
00395 } \
00396 } \
00397 return; \
00398 } \
00399 \
00400 if(ths->nfft_flags & PRE_FG_PSI) \
00401 { \
00402 for(t2=0; t2<ths->d; t2++) \
00403 { \
00404 tmpEXP2 = exp(-1.0/ths->b[t2]); \
00405 tmpEXP2sq = tmpEXP2*tmpEXP2; \
00406 tmp2 = 1.0; \
00407 tmp3 = 1.0; \
00408 fg_exp_l[t2][0] = 1.0; \
00409 for(lj_fg=1; lj_fg <= (2*ths->m+2); lj_fg++) \
00410 { \
00411 tmp3 = tmp2*tmpEXP2; \
00412 tmp2 *= tmpEXP2sq; \
00413 fg_exp_l[t2][lj_fg] = fg_exp_l[t2][lj_fg-1]*tmp3; \
00414 } \
00415 } \
00416 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00417 { \
00418 MACRO_init_uo_l_lj_t; \
00419 \
00420 for(t2=0; t2<ths->d; t2++) \
00421 { \
00422 fg_psi[t2][0] = ths->psi[2*(j*ths->d+t2)]; \
00423 tmpEXP1 = ths->psi[2*(j*ths->d+t2)+1]; \
00424 tmp1 = 1.0; \
00425 for(l_fg=u[t2]+1, lj_fg=1; l_fg <= o[t2]; l_fg++, lj_fg++) \
00426 { \
00427 tmp1 *= tmpEXP1; \
00428 fg_psi[t2][lj_fg] = fg_psi[t2][0]*tmp1*fg_exp_l[t2][lj_fg]; \
00429 } \
00430 } \
00431 \
00432 for(l_L=0; l_L<lprod; l_L++) \
00433 { \
00434 MACRO_update_phi_prod_ll_plain(with_FG_PSI); \
00435 \
00436 MACRO_nfft_B_compute_ ## which_one; \
00437 \
00438 MACRO_count_uo_l_lj_t; \
00439 } \
00440 } \
00441 return; \
00442 } \
00443 \
00444 if(ths->nfft_flags & FG_PSI) \
00445 { \
00446 for(t2=0; t2<ths->d; t2++) \
00447 { \
00448 tmpEXP2 = exp(-1.0/ths->b[t2]); \
00449 tmpEXP2sq = tmpEXP2*tmpEXP2; \
00450 tmp2 = 1.0; \
00451 tmp3 = 1.0; \
00452 fg_exp_l[t2][0] = 1.0; \
00453 for(lj_fg=1; lj_fg <= (2*ths->m+2); lj_fg++) \
00454 { \
00455 tmp3 = tmp2*tmpEXP2; \
00456 tmp2 *= tmpEXP2sq; \
00457 fg_exp_l[t2][lj_fg] = fg_exp_l[t2][lj_fg-1]*tmp3; \
00458 } \
00459 } \
00460 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00461 { \
00462 MACRO_init_uo_l_lj_t; \
00463 \
00464 for(t2=0; t2<ths->d; t2++) \
00465 { \
00466 fg_psi[t2][0] = \
00467 (PHI((ths->x[j*ths->d+t2]-((double)u[t2])/ths->n[t2]),t2)); \
00468 \
00469 tmpEXP1 = exp(2.0*(ths->n[t2]*ths->x[j*ths->d+t2] - u[t2]) \
00470 / ths->b[t2]); \
00471 tmp1 = 1.0; \
00472 for(l_fg=u[t2]+1, lj_fg=1; l_fg <= o[t2]; l_fg++, lj_fg++) \
00473 { \
00474 tmp1 *= tmpEXP1; \
00475 fg_psi[t2][lj_fg] = fg_psi[t2][0]*tmp1*fg_exp_l[t2][lj_fg]; \
00476 } \
00477 } \
00478 \
00479 for(l_L=0; l_L<lprod; l_L++) \
00480 { \
00481 MACRO_update_phi_prod_ll_plain(with_FG_PSI); \
00482 \
00483 MACRO_nfft_B_compute_ ## which_one; \
00484 \
00485 MACRO_count_uo_l_lj_t; \
00486 } \
00487 } \
00488 return; \
00489 } \
00490 \
00491 \
00492 if(ths->nfft_flags & PRE_LIN_PSI) \
00493 { \
00494 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00495 { \
00496 MACRO_init_uo_l_lj_t; \
00497 \
00498 for(t2=0; t2<ths->d; t2++) \
00499 { \
00500 y[t2] = ((ths->n[t2]*ths->x[j*ths->d+t2]- \
00501 (double)u[t2]) * ((double)ths->K))/(ths->m+1); \
00502 ip_u = (int)floor(y[t2]); \
00503 ip_w = y[t2]-ip_u; \
00504 for(l_fg=u[t2], lj_fg=0; l_fg <= o[t2]; l_fg++, lj_fg++) \
00505 { \
00506 fg_psi[t2][lj_fg] = ths->psi[(ths->K+1)*t2+ \
00507 abs(ip_u-lj_fg*ip_s)]* \
00508 (1-ip_w) + \
00509 ths->psi[(ths->K+1)*t2+ \
00510 abs(ip_u-lj_fg*ip_s+1)]* \
00511 (ip_w); \
00512 } \
00513 } \
00514 \
00515 for(l_L=0; l_L<lprod; l_L++) \
00516 { \
00517 MACRO_update_phi_prod_ll_plain(with_FG_PSI); \
00518 \
00519 MACRO_nfft_B_compute_ ## which_one; \
00520 \
00521 MACRO_count_uo_l_lj_t; \
00522 } \
00523 } \
00524 return; \
00525 } \
00526 \
00527 \
00528 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00529 { \
00530 MACRO_init_uo_l_lj_t; \
00531 \
00532 for(l_L=0; l_L<lprod; l_L++) \
00533 { \
00534 MACRO_update_phi_prod_ll_plain(without_PRE_PSI); \
00535 \
00536 MACRO_nfft_B_compute_ ## which_one; \
00537 \
00538 MACRO_count_uo_l_lj_t; \
00539 } \
00540 } \
00541 } \
00542
00543 MACRO_nfft_B(A)
00544 MACRO_nfft_B(T)
00545
00546
00549 void nfft_trafo(nfft_plan *ths)
00550 {
00551
00552 ths->g_hat=ths->g1;
00553 ths->g=ths->g2;
00554
00558 TIC(0)
00559 nfft_D_A(ths);
00560 TOC(0)
00561
00562
00566 TIC_FFTW(1)
00567 fftw_execute(ths->my_fftw_plan1);
00568 TOC_FFTW(1)
00569
00570
00573 TIC(2)
00574 nfft_B_A(ths);
00575 TOC(2)
00576 }
00577
00578 void nfft_adjoint(nfft_plan *ths)
00579 {
00580
00581 ths->g_hat=ths->g1;
00582 ths->g=ths->g2;
00583
00587 TIC(2)
00588 nfft_B_T(ths);
00589 TOC(2)
00590
00591
00595 TIC_FFTW(1)
00596 fftw_execute(ths->my_fftw_plan2);
00597 TOC_FFTW(1)
00598
00599
00602 TIC(0)
00603 nfft_D_T(ths);
00604 TOC(0)
00605 }
00606
00609 void nfft_precompute_phi_hut(nfft_plan *ths)
00610 {
00611 int ks[ths->d];
00612 int t;
00614 ths->c_phi_inv = (double**) fftw_malloc(ths->d*sizeof(double*));
00615
00616 for(t=0; t<ths->d; t++)
00617 {
00618 ths->c_phi_inv[t]= (double*)fftw_malloc(ths->N[t]*sizeof(double));
00619 for(ks[t]=0; ks[t]<ths->N[t]; ks[t]++)
00620 ths->c_phi_inv[t][ks[t]]= 1.0/(PHI_HUT(ks[t]-ths->N[t]/2,t));
00621 }
00622 }
00623
00629 void nfft_precompute_lin_psi(nfft_plan *ths)
00630 {
00631 int t;
00632 int j;
00633 double step;
00635 for (t=0; t<ths->d; t++)
00636 {
00637 step=((double)(ths->m+1))/(ths->K*ths->n[t]);
00638 for(j=0;j<=ths->K;j++)
00639 {
00640 ths->psi[(ths->K+1)*t + j] = PHI(j*step,t);
00641 }
00642 }
00643 }
00644
00645 void nfft_precompute_fg_psi(nfft_plan *ths)
00646 {
00647 int t;
00648 int j;
00649 int u, o;
00651 for (t=0; t<ths->d; t++)
00652 for(j=0;j<ths->M_total;j++)
00653 {
00654 nfft_uo(ths,j,&u,&o,t);
00655
00656 ths->psi[2*(j*ths->d+t)]=
00657 (PHI((ths->x[j*ths->d+t]-((double)u)/ths->n[t]),t));
00658
00659 ths->psi[2*(j*ths->d+t)+1]=
00660 exp(2.0*(ths->n[t]*ths->x[j*ths->d+t] - u) / ths->b[t]);
00661 }
00662
00663 }
00664
00665 void nfft_precompute_psi(nfft_plan *ths)
00666 {
00667 int t;
00668 int j;
00669 int l;
00670 int lj;
00671 int u, o;
00673 for (t=0; t<ths->d; t++)
00674 for(j=0;j<ths->M_total;j++)
00675 {
00676 nfft_uo(ths,j,&u,&o,t);
00677
00678 for(l=u, lj=0; l <= o; l++, lj++)
00679 ths->psi[(j*ths->d+t)*(2*ths->m+2)+lj]=
00680 (PHI((ths->x[j*ths->d+t]-((double)l)/ths->n[t]),t));
00681 }
00682
00683 }
00684
00685 void nfft_precompute_full_psi(nfft_plan *ths)
00686 {
00687 int t,t2;
00688 int j;
00689 int l_L;
00690 int l[ths->d];
00691 int lj[ths->d];
00692 int ll_plain[ths->d+1];
00693 int lprod;
00694 int u[ths->d], o[ths->d];
00696 double phi_prod[ths->d+1];
00697
00698 int ix,ix_old;
00699
00700 phi_prod[0]=1;
00701 ll_plain[0]=0;
00702
00703 for(t=0,lprod = 1; t<ths->d; t++)
00704 lprod *= 2*ths->m+2;
00705
00706 for(j=0,ix=0,ix_old=0; j<ths->M_total; j++)
00707 {
00708 MACRO_init_uo_l_lj_t;
00709
00710 for(l_L=0; l_L<lprod; l_L++, ix++)
00711 {
00712 MACRO_update_phi_prod_ll_plain(without_PRE_PSI);
00713
00714 ths->psi_index_g[ix]=ll_plain[ths->d];
00715 ths->psi[ix]=phi_prod[ths->d];
00716
00717 MACRO_count_uo_l_lj_t;
00718 }
00719
00720
00721 ths->psi_index_f[j]=ix-ix_old;
00722 ix_old=ix;
00723 }
00724 }
00725
00726 void nfft_precompute_one_psi(nfft_plan *ths)
00727 {
00728 if(ths->nfft_flags & PRE_LIN_PSI)
00729 nfft_precompute_lin_psi(ths);
00730 if(ths->nfft_flags & PRE_FG_PSI)
00731 nfft_precompute_fg_psi(ths);
00732 if(ths->nfft_flags & PRE_PSI)
00733 nfft_precompute_psi(ths);
00734 if(ths->nfft_flags & PRE_FULL_PSI)
00735 nfft_precompute_full_psi(ths);
00736 }
00737
00738
00739 void nfft_init_help(nfft_plan *ths)
00740 {
00741 int t;
00742 int lprod;
00744 ths->N_total=nfft_prod_int(ths->N, ths->d);
00745 ths->n_total=nfft_prod_int(ths->n, ths->d);
00746
00747 ths->sigma = (double*) fftw_malloc(ths->d*sizeof(double));
00748 for(t = 0;t < ths->d; t++)
00749 ths->sigma[t] = ((double)ths->n[t])/ths->N[t];
00750
00751 WINDOW_HELP_INIT;
00752
00753 if(ths->nfft_flags & MALLOC_X)
00754 ths->x = (double*)fftw_malloc(ths->d*ths->M_total*sizeof(double));
00755
00756 if(ths->nfft_flags & MALLOC_F_HAT)
00757 ths->f_hat = (double complex*)fftw_malloc(ths->N_total*sizeof(double complex));
00758
00759 if(ths->nfft_flags & MALLOC_F)
00760 ths->f = (double complex*)fftw_malloc(ths->M_total*sizeof(double complex));
00761
00762 if(ths->nfft_flags & PRE_PHI_HUT)
00763 nfft_precompute_phi_hut(ths);
00764
00765 if(ths->nfft_flags & PRE_LIN_PSI)
00766 {
00767 ths->K=(1U<< 10)*(ths->m+1);
00768 ths->psi = (double*) fftw_malloc((ths->K+1)*ths->d*sizeof(double));
00769 }
00770
00771 if(ths->nfft_flags & PRE_FG_PSI)
00772 ths->psi = (double*) fftw_malloc(ths->M_total*ths->d*2*sizeof(double));
00773
00774 if(ths->nfft_flags & PRE_PSI)
00775 ths->psi = (double*) fftw_malloc(ths->M_total*ths->d*
00776 (2*ths->m+2)*sizeof(double));
00777
00778 if(ths->nfft_flags & PRE_FULL_PSI)
00779 {
00780 for(t=0,lprod = 1; t<ths->d; t++)
00781 lprod *= 2*ths->m+2;
00782
00783 ths->psi = (double*) fftw_malloc(ths->M_total*lprod*sizeof(double));
00784
00785 ths->psi_index_f = (int*) fftw_malloc(ths->M_total*sizeof(int));
00786 ths->psi_index_g = (int*) fftw_malloc(ths->M_total*lprod*sizeof(int));
00787 }
00788
00789 if(ths->nfft_flags & FFTW_INIT)
00790 {
00791 ths->g1=(double complex*)fftw_malloc(ths->n_total*
00792 sizeof(double complex));
00793
00794 if(ths->nfft_flags & FFT_OUT_OF_PLACE)
00795 ths->g2 = (double complex*) fftw_malloc(ths->n_total*
00796 sizeof(double complex));
00797 else
00798 ths->g2 = ths->g1;
00799
00800 ths->my_fftw_plan1 =
00801 fftw_plan_dft(ths->d, ths->n, ths->g1, ths->g2,
00802 FFTW_FORWARD, ths->fftw_flags);
00803 ths->my_fftw_plan2 =
00804 fftw_plan_dft(ths->d, ths->n, ths->g2, ths->g1,
00805 FFTW_BACKWARD, ths->fftw_flags);
00806 }
00807 }
00808
00809 void nfft_init(nfft_plan *ths, int d, int *N, int M_total)
00810 {
00811 int t;
00813 ths->d = d;
00814
00815 ths->N=(int*) fftw_malloc(d*sizeof(int));
00816 for(t = 0;t < d; t++)
00817 ths->N[t] = N[t];
00818
00819 ths->M_total = M_total;
00820
00821 ths->n = (int*) fftw_malloc(d*sizeof(int));
00822 for(t = 0;t < d; t++)
00823 ths->n[t] = 2*nfft_next_power_of_2(ths->N[t]);
00824
00825 WINDOW_HELP_ESTIMATE_m;
00826
00827 ths->nfft_flags = PRE_PHI_HUT| PRE_PSI| MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00828 FFTW_INIT| FFT_OUT_OF_PLACE;
00829 ths->fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
00830
00831 nfft_init_help(ths);
00832 }
00833
00834 void nfft_init_guru(nfft_plan *ths, int d, int *N, int M_total, int *n,
00835 int m, unsigned nfft_flags, unsigned fftw_flags)
00836 {
00837 int t;
00839 ths->d =d;
00840 ths->N= (int*) fftw_malloc(ths->d*sizeof(int));
00841 for(t=0; t<d; t++)
00842 ths->N[t]= N[t];
00843 ths->M_total= M_total;
00844 ths->n= (int*) fftw_malloc(ths->d*sizeof(int));
00845 for(t=0; t<d; t++)
00846 ths->n[t]= n[t];
00847 ths->m= m;
00848 ths->nfft_flags= nfft_flags;
00849 ths->fftw_flags= fftw_flags;
00850
00851 nfft_init_help(ths);
00852 }
00853
00854 void nfft_init_1d(nfft_plan *ths, int N1, int M_total)
00855 {
00856 int N[1];
00857
00858 N[0]=N1;
00859 nfft_init(ths,1,N,M_total);
00860 }
00861
00862 void nfft_init_2d(nfft_plan *ths, int N1, int N2, int M_total)
00863 {
00864 int N[2];
00865
00866 N[0]=N1;
00867 N[1]=N2;
00868 nfft_init(ths,2,N,M_total);
00869 }
00870
00871 void nfft_init_3d(nfft_plan *ths, int N1, int N2, int N3, int M_total)
00872 {
00873 int N[3];
00874
00875 N[0]=N1;
00876 N[1]=N2;
00877 N[2]=N3;
00878 nfft_init(ths,3,N,M_total);
00879 }
00880
00881 void nfft_check(nfft_plan *ths)
00882 {
00883 int j;
00884
00885 for(j=0;j<ths->M_total*ths->d;j++)
00886 if((ths->x[j]<-0.5) || (ths->x[j]>=0.5))
00887 fprintf(stderr,"nfft_check: ths->x[%d]=%e out of range [-0.5,0.5)\n",
00888 j,ths->x[j]);
00889
00890 for(j=0;j<ths->d;j++)
00891 {
00892 if(ths->sigma[j]<=1)
00893 fprintf(stderr,"nfft_check: oversampling factor too small\n");
00894
00895 if(ths->N[j]<=ths->m)
00896 fprintf(stderr,
00897 "nfft_check: polynomial degree N is smaller than cut-off m\n");
00898
00899 if(ths->N[j]%2==1)
00900 fprintf(stderr,"nfft_check: polynomial degree N has to be even\n");
00901 }
00902 }
00903
00904 void nfft_finalize(nfft_plan *ths)
00905 {
00906 int t;
00907
00908 if(ths->nfft_flags & FFTW_INIT)
00909 {
00910 fftw_destroy_plan(ths->my_fftw_plan2);
00911 fftw_destroy_plan(ths->my_fftw_plan1);
00912
00913 if(ths->nfft_flags & FFT_OUT_OF_PLACE)
00914 fftw_free(ths->g2);
00915
00916 fftw_free(ths->g1);
00917 }
00918
00919 if(ths->nfft_flags & PRE_FULL_PSI)
00920 {
00921 fftw_free(ths->psi_index_g);
00922 fftw_free(ths->psi_index_f);
00923 fftw_free(ths->psi);
00924 }
00925
00926 if(ths->nfft_flags & PRE_PSI)
00927 fftw_free(ths->psi);
00928
00929 if(ths->nfft_flags & PRE_FG_PSI)
00930 fftw_free(ths->psi);
00931
00932 if(ths->nfft_flags & PRE_LIN_PSI)
00933 fftw_free(ths->psi);
00934
00935 if(ths->nfft_flags & PRE_PHI_HUT)
00936 {
00937 for(t=0; t<ths->d; t++)
00938 fftw_free(ths->c_phi_inv[t]);
00939 fftw_free(ths->c_phi_inv);
00940 }
00941
00942 if(ths->nfft_flags & MALLOC_F)
00943 fftw_free(ths->f);
00944
00945 if(ths->nfft_flags & MALLOC_F_HAT)
00946 fftw_free(ths->f_hat);
00947
00948 if(ths->nfft_flags & MALLOC_X)
00949 fftw_free(ths->x);
00950
00951 WINDOW_HELP_FINALIZE;
00952
00953 fftw_free(ths->sigma);
00954 fftw_free(ths->n);
00955 fftw_free(ths->N);
00956 }
