32 #define MACRO_nndft_init_result_trafo memset(f,0,ths->M_total*sizeof(double _Complex));
33 #define MACRO_nndft_init_result_conjugated MACRO_nndft_init_result_trafo
34 #define MACRO_nndft_init_result_adjoint memset(f_hat,0,ths->N_total*sizeof(double _Complex));
35 #define MACRO_nndft_init_result_transposed MACRO_nndft_init_result_adjoint
37 #define MACRO_nndft_sign_trafo (-2.0*KPI)
38 #define MACRO_nndft_sign_conjugated (+2.0*KPI)
39 #define MACRO_nndft_sign_adjoint (+2.0*KPI)
40 #define MACRO_nndft_sign_transposed (-2.0*KPI)
42 #define MACRO_nndft_compute_trafo (*fj) += (*f_hat_k)*cexp(+ _Complex_I*omega);
44 #define MACRO_nndft_compute_conjugated MACRO_nndft_compute_trafo
46 #define MACRO_nndft_compute_adjoint (*f_hat_k) += (*fj)*cexp(+ _Complex_I*omega);
48 #define MACRO_nndft_compute_transposed MACRO_nndft_compute_adjoint
50 #define MACRO_nndft(which_one) \
51 void nnfft_ ## which_one ## _direct (nnfft_plan *ths) \
56 double _Complex *f_hat, *f; \
57 double _Complex *f_hat_k; \
58 double _Complex *fj; \
61 f_hat=ths->f_hat; f=ths->f; \
63 MACRO_nndft_init_result_ ## which_one \
65 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
67 for(l=0, f_hat_k=f_hat; l<ths->N_total; l++, f_hat_k++) \
70 for(t = 0; t<ths->d; t++) \
71 omega+=ths->v[l*ths->d+t] * ths->x[j*ths->d+t] * ths->N[t]; \
73 omega*= MACRO_nndft_sign_ ## which_one; \
75 MACRO_nndft_compute_ ## which_one \
86 static void nnfft_uo(
nnfft_plan *ths,
int j,
int *up,
int *op,
int act_dim)
91 c = ths->v[j*ths->d+act_dim] * ths->n[act_dim];
99 u = u - (ths->m); o = o + (ths->m);
107 #define MACRO_nnfft_B_init_result_A memset(f,0,ths->N_total*sizeof(double _Complex));
108 #define MACRO_nnfft_B_init_result_T memset(g,0,ths->aN1_total*sizeof(double _Complex));
110 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_A { \
111 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
114 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_T { \
115 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
118 #define MACRO_nnfft_B_compute_A { \
119 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
122 #define MACRO_nnfft_B_compute_T { \
123 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
126 #define MACRO_with_PRE_LIN_PSI (ths->psi[(ths->K+1)*t2+y_u[t2]]* \
127 (y_u[t2]+1-y[t2]) + \
128 ths->psi[(ths->K+1)*t2+y_u[t2]+1]* \
130 #define MACRO_with_PRE_PSI ths->psi[(j*ths->d+t2)*(2*ths->m+2)+lj[t2]]
131 #define MACRO_without_PRE_PSI PHI(ths->n[t2], -ths->v[j*ths->d+t2]+ \
132 ((double)l[t2])/ths->N1[t2], t2)
134 #define MACRO_init_uo_l_lj_t { \
135 for(t = ths->d-1; t>=0; t--) \
137 nnfft_uo(ths,j,&u[t],&o[t],t); \
144 #define MACRO_update_with_PRE_PSI_LIN { \
145 for(t2=t; t2<ths->d; t2++) \
147 y[t2] = fabs(((-ths->N1[t2]*ths->v[j*ths->d+t2]+(double)l[t2]) \
148 * ((double)ths->K))/(ths->m+1)); \
149 y_u[t2] = (int)y[t2]; \
153 #define MACRO_update_phi_prod_ll_plain(which_one) { \
154 for(t2=t; t2<ths->d; t2++) \
156 phi_prod[t2+1]=phi_prod[t2]* MACRO_ ## which_one; \
157 ll_plain[t2+1]=ll_plain[t2]*ths->aN1[t2] + \
158 (l[t2]+ths->aN1[t2]*3/2)%ths->aN1[t2]; \
163 #define MACRO_count_uo_l_lj_t { \
164 for(t = ths->d-1; (t>0)&&(l[t]==o[t]); t--) \
174 #define MACRO_nnfft_B(which_one) \
175 static inline void nnfft_B_ ## which_one (nnfft_plan *ths) \
178 int u[ths->d], o[ths->d]; \
184 int ll_plain[ths->d+1]; \
185 double phi_prod[ths->d+1]; \
186 double _Complex *f, *g; \
187 double _Complex *fj; \
191 f=ths->f_hat; g=ths->F; \
193 MACRO_nnfft_B_init_result_ ## which_one \
195 if(ths->nnfft_flags & PRE_FULL_PSI) \
197 for(ix=0, j=0, fj=f; j<ths->N_total; j++,fj++) \
198 for(l_L=0; l_L<ths->psi_index_f[j]; l_L++, ix++) \
199 MACRO_nnfft_B_PRE_FULL_PSI_compute_ ## which_one; \
206 for(t=0,lprod = 1; t<ths->d; t++) \
207 lprod *= (2*ths->m+2); \
209 if(ths->nnfft_flags & PRE_PSI) \
211 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
213 MACRO_init_uo_l_lj_t; \
215 for(l_L=0; l_L<lprod; l_L++) \
217 MACRO_update_phi_prod_ll_plain(with_PRE_PSI); \
219 MACRO_nnfft_B_compute_ ## which_one; \
221 MACRO_count_uo_l_lj_t; \
227 if(ths->nnfft_flags & PRE_LIN_PSI) \
229 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
231 MACRO_init_uo_l_lj_t; \
233 for(l_L=0; l_L<lprod; l_L++) \
235 MACRO_update_with_PRE_PSI_LIN; \
237 MACRO_update_phi_prod_ll_plain(with_PRE_LIN_PSI); \
239 MACRO_nnfft_B_compute_ ## which_one; \
241 MACRO_count_uo_l_lj_t; \
248 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
251 MACRO_init_uo_l_lj_t; \
253 for(l_L=0; l_L<lprod; l_L++) \
255 MACRO_update_phi_prod_ll_plain(without_PRE_PSI); \
257 MACRO_nnfft_B_compute_ ## which_one; \
259 MACRO_count_uo_l_lj_t; \
273 for(j=0; j<ths->M_total; j++)
274 ths->f[j] *= ths->c_phi_inv[j];
278 for(j=0; j<ths->M_total; j++)
282 for(t=0; t<ths->d; t++)
283 tmp*= 1.0 /((PHI_HUT(ths->n[t], ths->x[ths->d*j + t]*((
double)ths->N[t]),t)) );
297 for(j=0;j<ths->M_total;j++) {
298 for(t=0;t<ths->d;t++) {
299 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
305 ths->direct_plan->f = ths->f;
307 nfft_trafo(ths->direct_plan);
309 for(j=0;j<ths->M_total;j++) {
310 for(t=0;t<ths->d;t++) {
311 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
325 for(j=0;j<ths->M_total;j++) {
326 for(t=0;t<ths->d;t++) {
327 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
332 ths->direct_plan->f=ths->f;
334 nfft_adjoint(ths->direct_plan);
336 for(j=0;j<ths->M_total;j++) {
337 for(t=0;t<ths->d;t++) {
338 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
353 ths->c_phi_inv= (
double*)nfft_malloc(ths->M_total*
sizeof(
double));
355 for(j=0; j<ths->M_total; j++)
358 for(t=0; t<ths->d; t++)
359 tmp*= 1.0 /(PHI_HUT(ths->n[t],ths->x[ths->d*j + t]*((
double)ths->N[t]),t));
360 ths->c_phi_inv[j]=tmp;
373 nfft_precompute_lin_psi(ths->direct_plan);
375 for (t=0; t<ths->d; t++)
377 step=((double)(ths->m+1))/(ths->K*ths->N1[t]);
378 for(j=0;j<=ths->K;j++)
380 ths->psi[(ths->K+1)*t + j] = PHI(ths->n[t],j*step,t);
393 for (t=0; t<ths->d; t++)
394 for(j=0;j<ths->N_total;j++)
396 nnfft_uo(ths,j,&u,&o,t);
398 for(l=u, lj=0; l <= o; l++, lj++)
399 ths->psi[(j*ths->d+t)*(2*ths->m+2)+lj]=
400 (PHI(ths->n[t],(-ths->v[j*ths->d+t]+((
double)l)/((
double)ths->N1[t])),t));
403 for(j=0;j<ths->M_total;j++) {
404 for(t=0;t<ths->d;t++) {
405 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
409 nfft_precompute_psi(ths->direct_plan);
411 for(j=0;j<ths->M_total;j++) {
412 for(t=0;t<ths->d;t++) {
413 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
431 int ll_plain[ths->d+1];
433 int u[ths->d], o[ths->d];
435 double phi_prod[ths->d+1];
439 for(j=0;j<ths->M_total;j++) {
440 for(t=0;t<ths->d;t++) {
441 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
447 nfft_precompute_full_psi(ths->direct_plan);
449 for(j=0;j<ths->M_total;j++) {
450 for(t=0;t<ths->d;t++) {
451 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
458 for(t=0,lprod = 1; t<ths->d; t++)
461 for(j=0,ix=0,ix_old=0; j<ths->N_total; j++)
463 MACRO_init_uo_l_lj_t;
465 for(l_L=0; l_L<lprod; l_L++, ix++)
467 MACRO_update_phi_prod_ll_plain(without_PRE_PSI);
469 ths->psi_index_g[ix]=ll_plain[ths->d];
470 ths->psi[ix]=phi_prod[ths->d];
472 MACRO_count_uo_l_lj_t;
476 ths->psi_index_f[j]=ix-ix_old;
481 void nnfft_precompute_one_psi(
nnfft_plan *ths)
494 static void nnfft_init_help(
nnfft_plan *ths,
int m2,
unsigned nfft_flags,
unsigned fftw_flags)
500 ths->aN1 = (
int*) nfft_malloc(ths->d*
sizeof(
int));
502 ths->a = (
double*) nfft_malloc(ths->d*
sizeof(
double));
504 ths->sigma = (
double*) nfft_malloc(ths->d*
sizeof(
double));
510 for(t = 0; t<ths->d; t++) {
511 ths->a[t] = 1.0 + (2.0*((double)ths->m))/((
double)ths->N1[t]);
512 ths->aN1[t] = ths->a[t] * ((double)ths->N1[t]);
514 if(ths->aN1[t]%2 != 0)
515 ths->aN1[t] = ths->aN1[t] +1;
517 ths->aN1_total*=ths->aN1[t];
518 ths->sigma[t] = ((double) ths->N1[t] )/((double) ths->N[t]);;
521 N2[t] = ceil(ths->sigma[t]*(ths->aN1[t]));
531 ths->x = (
double*)nfft_malloc(ths->d*ths->M_total*
sizeof(
double));
533 ths->f=(
double _Complex*)nfft_malloc(ths->M_total*
sizeof(
double _Complex));
536 ths->v = (
double*)nfft_malloc(ths->d*ths->N_total*
sizeof(
double));
538 ths->f_hat = (
double _Complex*)nfft_malloc(ths->N_total*
sizeof(
double _Complex));
547 ths->K=(1U<< 10)*(ths->m+1);
548 ths->psi = (
double*) nfft_malloc((ths->K+1)*ths->d*
sizeof(double));
551 if(ths->nnfft_flags &
PRE_PSI){
552 ths->psi = (
double*)nfft_malloc(ths->N_total*ths->d*(2*ths->m+2)*
sizeof(double));
557 for(t=0,lprod = 1; t<ths->d; t++)
560 ths->psi = (
double*)nfft_malloc(ths->N_total*lprod*
sizeof(
double));
562 ths->psi_index_f = (
int*) nfft_malloc(ths->N_total*
sizeof(
int));
563 ths->psi_index_g = (
int*) nfft_malloc(ths->N_total*lprod*
sizeof(
int));
566 nfft_init_guru(ths->direct_plan, ths->d, ths->aN1, ths->M_total, N2, m2,
567 nfft_flags, fftw_flags);
568 ths->direct_plan->x = ths->x;
570 ths->direct_plan->f = ths->f;
571 ths->F = ths->direct_plan->f_hat;
578 int m,
unsigned nnfft_flags)
586 ths->M_total= M_total;
587 ths->N_total= N_total;
589 ths->nnfft_flags= nnfft_flags;
590 fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
595 nfft_flags = nfft_flags |
PRE_PSI;
603 ths->N = (
int*) nfft_malloc(ths->d*
sizeof(
int));
604 ths->N1 = (
int*) nfft_malloc(ths->d*
sizeof(
int));
610 nnfft_init_help(ths,m,nfft_flags,fftw_flags);
620 ths->M_total = M_total;
621 ths->N_total = N_total;
628 ths->m=WINDOW_HELP_ESTIMATE_m;
631 ths->N = (
int*) nfft_malloc(ths->d*
sizeof(
int));
632 ths->N1 = (
int*) nfft_malloc(ths->d*
sizeof(
int));
637 ths->N1[t] = ceil(1.5*ths->N[t]);
640 if(ths->N1[t]%2 != 0)
641 ths->N1[t] = ths->N1[t] +1;
648 fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
649 nnfft_init_help(ths,ths->m,nfft_flags,fftw_flags);
652 void nnfft_init_1d(
nnfft_plan *ths,
int N1,
int M_total)
659 nfft_finalize(ths->direct_plan);
660 nfft_free(ths->direct_plan);
668 nfft_free(ths->psi_index_g);
669 nfft_free(ths->psi_index_f);
680 nfft_free(ths->c_phi_inv);
686 nfft_free(ths->f_hat);
void nnfft_precompute_full_psi(nnfft_plan *ths)
computes all entries of B explicitly
void nnfft_finalize(nnfft_plan *ths)
void nnfft_adjoint(nnfft_plan *ths)
void nnfft_init(nnfft_plan *ths, int d, int N_total, int M_total, int *N)
void nnfft_precompute_phi_hut(nnfft_plan *ths)
initialisation of direct transform
void nnfft_precompute_lin_psi(nnfft_plan *ths)
create a lookup table
void nnfft_precompute_psi(nnfft_plan *ths)
void nnfft_init_guru(nnfft_plan *ths, int d, int N_total, int M_total, int *N, int *N1, int m, unsigned nnfft_flags)
void nnfft_trafo(nnfft_plan *ths)
user routines
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.