34 #define MACRO_nndft_init_result_trafo memset(f,0,ths->M_total*sizeof(double _Complex)); 35 #define MACRO_nndft_init_result_conjugated MACRO_nndft_init_result_trafo 36 #define MACRO_nndft_init_result_adjoint memset(f_hat,0,ths->N_total*sizeof(double _Complex)); 37 #define MACRO_nndft_init_result_transposed MACRO_nndft_init_result_adjoint 39 #define MACRO_nndft_sign_trafo (-2.0*KPI) 40 #define MACRO_nndft_sign_conjugated (+2.0*KPI) 41 #define MACRO_nndft_sign_adjoint (+2.0*KPI) 42 #define MACRO_nndft_sign_transposed (-2.0*KPI) 44 #define MACRO_nndft_compute_trafo (*fj) += (*f_hat_k)*cexp(+ _Complex_I*omega); 46 #define MACRO_nndft_compute_conjugated MACRO_nndft_compute_trafo 48 #define MACRO_nndft_compute_adjoint (*f_hat_k) += (*fj)*cexp(+ _Complex_I*omega); 50 #define MACRO_nndft_compute_transposed MACRO_nndft_compute_adjoint 52 #define MACRO_nndft(which_one) \ 53 void nnfft_ ## which_one ## _direct (nnfft_plan *ths) \ 58 double _Complex *f_hat, *f; \ 59 double _Complex *f_hat_k; \ 60 double _Complex *fj; \ 63 f_hat=ths->f_hat; f=ths->f; \ 65 MACRO_nndft_init_result_ ## which_one \ 67 for(j=0, fj=f; j<ths->M_total; j++, fj++) \ 69 for(l=0, f_hat_k=f_hat; l<ths->N_total; l++, f_hat_k++) \ 72 for(t = 0; t<ths->d; t++) \ 73 omega+=ths->v[l*ths->d+t] * ths->x[j*ths->d+t] * ths->N[t]; \ 75 omega*= MACRO_nndft_sign_ ## which_one; \ 77 MACRO_nndft_compute_ ## which_one \ 88 static
void nnfft_uo(
nnfft_plan *ths,
int j,
int *up,
int *op,
int act_dim)
93 c = ths->v[j*ths->d+act_dim] * ths->n[act_dim];
101 u = u - (ths->m); o = o + (ths->m);
109 #define MACRO_nnfft_B_init_result_A memset(f,0,ths->N_total*sizeof(double _Complex)); 110 #define MACRO_nnfft_B_init_result_T memset(g,0,ths->aN1_total*sizeof(double _Complex)); 112 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_A { \ 113 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \ 116 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_T { \ 117 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \ 120 #define MACRO_nnfft_B_compute_A { \ 121 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \ 124 #define MACRO_nnfft_B_compute_T { \ 125 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \ 128 #define MACRO_with_PRE_LIN_PSI (ths->psi[(ths->K+1)*t2+y_u[t2]]* \ 129 (y_u[t2]+1-y[t2]) + \ 130 ths->psi[(ths->K+1)*t2+y_u[t2]+1]* \ 132 #define MACRO_with_PRE_PSI ths->psi[(j*ths->d+t2)*(2*ths->m+2)+lj[t2]] 133 #define MACRO_without_PRE_PSI PHI(ths->n[t2], -ths->v[j*ths->d+t2]+ \ 134 ((double)l[t2])/ths->N1[t2], t2) 136 #define MACRO_init_uo_l_lj_t { \ 137 for(t = ths->d-1; t>=0; t--) \ 139 nnfft_uo(ths,j,&u[t],&o[t],t); \ 146 #define MACRO_update_with_PRE_PSI_LIN { \ 147 for(t2=t; t2<ths->d; t2++) \ 149 y[t2] = fabs(((-ths->N1[t2]*ths->v[j*ths->d+t2]+(double)l[t2]) \ 150 * ((double)ths->K))/(ths->m+1)); \ 151 y_u[t2] = (int)y[t2]; \ 155 #define MACRO_update_phi_prod_ll_plain(which_one) { \ 156 for(t2=t; t2<ths->d; t2++) \ 158 phi_prod[t2+1]=phi_prod[t2]* MACRO_ ## which_one; \ 159 ll_plain[t2+1]=ll_plain[t2]*ths->aN1[t2] + \ 160 (l[t2]+ths->aN1[t2]*3/2)%ths->aN1[t2]; \ 165 #define MACRO_count_uo_l_lj_t { \ 166 for(t = ths->d-1; (t>0)&&(l[t]==o[t]); t--) \ 176 #define MACRO_nnfft_B(which_one) \ 177 static inline void nnfft_B_ ## which_one (nnfft_plan *ths) \ 180 int u[ths->d], o[ths->d]; \ 186 int ll_plain[ths->d+1]; \ 187 double phi_prod[ths->d+1]; \ 188 double _Complex *f, *g; \ 189 double _Complex *fj; \ 193 f=ths->f_hat; g=ths->F; \ 195 MACRO_nnfft_B_init_result_ ## which_one \ 197 if(ths->nnfft_flags & PRE_FULL_PSI) \ 199 for(ix=0, j=0, fj=f; j<ths->N_total; j++,fj++) \ 200 for(l_L=0; l_L<ths->psi_index_f[j]; l_L++, ix++) \ 201 MACRO_nnfft_B_PRE_FULL_PSI_compute_ ## which_one; \ 208 for(t=0,lprod = 1; t<ths->d; t++) \ 209 lprod *= (2*ths->m+2); \ 211 if(ths->nnfft_flags & PRE_PSI) \ 213 for(j=0, fj=f; j<ths->N_total; j++, fj++) \ 215 MACRO_init_uo_l_lj_t; \ 217 for(l_L=0; l_L<lprod; l_L++) \ 219 MACRO_update_phi_prod_ll_plain(with_PRE_PSI); \ 221 MACRO_nnfft_B_compute_ ## which_one; \ 223 MACRO_count_uo_l_lj_t; \ 229 if(ths->nnfft_flags & PRE_LIN_PSI) \ 231 for(j=0, fj=f; j<ths->N_total; j++, fj++) \ 233 MACRO_init_uo_l_lj_t; \ 235 for(l_L=0; l_L<lprod; l_L++) \ 237 MACRO_update_with_PRE_PSI_LIN; \ 239 MACRO_update_phi_prod_ll_plain(with_PRE_LIN_PSI); \ 241 MACRO_nnfft_B_compute_ ## which_one; \ 243 MACRO_count_uo_l_lj_t; \ 250 for(j=0, fj=f; j<ths->N_total; j++, fj++) \ 253 MACRO_init_uo_l_lj_t; \ 255 for(l_L=0; l_L<lprod; l_L++) \ 257 MACRO_update_phi_prod_ll_plain(without_PRE_PSI); \ 259 MACRO_nnfft_B_compute_ ## which_one; \ 261 MACRO_count_uo_l_lj_t; \ 273 if(ths->nnfft_flags & PRE_PHI_HUT)
275 for(j=0; j<ths->M_total; j++)
276 ths->f[j] *= ths->c_phi_inv[j];
280 for(j=0; j<ths->M_total; j++)
284 for(t=0; t<ths->d; t++)
285 tmp*= 1.0 /((PHI_HUT(ths->n[t], ths->x[ths->d*j + t]*((
double)ths->N[t]),t)) );
300 for(t=0;t<ths->
d;t++) {
301 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] / ((double)ths->
sigma[t]);
312 for(t=0;t<ths->
d;t++) {
313 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] * ((double)ths->
sigma[t]);
328 for(t=0;t<ths->
d;t++) {
329 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] / ((double)ths->
sigma[t]);
339 for(t=0;t<ths->
d;t++) {
340 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] * ((double)ths->
sigma[t]);
360 for(t=0; t<ths->
d; t++)
361 tmp*= 1.0 /(PHI_HUT(ths->
n[t],ths->
x[ths->
d*j + t]*((
double)ths->
N[t]),t));
377 for (t=0; t<ths->
d; t++)
379 step=((double)(ths->
m+1))/(ths->
K*ths->
N1[t]);
380 for(j=0;j<=ths->
K;j++)
382 ths->
psi[(ths->
K+1)*t + j] = PHI(ths->
n[t],j*step,t);
395 for (t=0; t<ths->
d; t++)
398 nnfft_uo(ths,j,&u,&o,t);
400 for(l=u, lj=0; l <= o; l++, lj++)
401 ths->
psi[(j*ths->
d+t)*(2*ths->
m+2)+lj]=
402 (PHI(ths->
n[t],(-ths->
v[j*ths->
d+t]+((
double)l)/((
double)ths->
N1[t])),t));
406 for(t=0;t<ths->
d;t++) {
407 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] / ((double)ths->
sigma[t]);
414 for(t=0;t<ths->
d;t++) {
415 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] * ((double)ths->
sigma[t]);
433 int ll_plain[ths->
d+1];
435 int u[ths->
d], o[ths->
d];
437 double phi_prod[ths->
d+1];
442 for(t=0;t<ths->
d;t++) {
443 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] / ((double)ths->
sigma[t]);
447 nnfft_precompute_psi(ths);
452 for(t=0;t<ths->
d;t++) {
453 ths->
x[j*ths->
d+t]= ths->
x[j*ths->
d+t] * ((double)ths->
sigma[t]);
460 for(t=0,lprod = 1; t<ths->
d; t++)
463 for(j=0,ix=0,ix_old=0; j<ths->N_total; j++)
465 MACRO_init_uo_l_lj_t;
467 for(l_L=0; l_L<lprod; l_L++, ix++)
469 MACRO_update_phi_prod_ll_plain(without_PRE_PSI);
472 ths->
psi[ix]=phi_prod[ths->
d];
474 MACRO_count_uo_l_lj_t;
483 void nnfft_precompute_one_psi(
nnfft_plan *ths)
486 nnfft_precompute_psi(ths);
496 static void nnfft_init_help(
nnfft_plan *ths,
int m2,
unsigned nfft_flags,
unsigned fftw_flags)
512 for(t = 0; t<ths->
d; t++) {
513 ths->
a[t] = 1.0 + (2.0*((double)ths->
m))/((
double)ths->
N1[t]);
514 ths->
aN1[t] = ths->
a[t] * ((double)ths->
N1[t]);
516 if(ths->
aN1[t]%2 != 0)
517 ths->
aN1[t] = ths->
aN1[t] +1;
520 ths->
sigma[t] = ((double) ths->
N1[t] )/((double) ths->
N[t]);;
523 N2[t] = ceil(ths->
sigma[t]*(ths->
aN1[t]));
549 ths->
K=(1U<< 10)*(ths->
m+1);
559 for(t=0,lprod = 1; t<ths->
d; t++)
569 nfft_flags, fftw_flags);
576 ths->
mv_adjoint = (void (*) (
void* ))nnfft_adjoint;
579 void nnfft_init_guru(
nnfft_plan *ths,
int d,
int N_total,
int M_total,
int *N,
int *N1,
580 int m,
unsigned nnfft_flags)
592 fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
593 nfft_flags= PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE| NFFT_OMP_BLOCKWISE_ADJOINT;
596 nfft_flags = nfft_flags | PRE_PSI;
599 nfft_flags = nfft_flags | PRE_FULL_PSI;
602 nfft_flags = nfft_flags | PRE_LIN_PSI;
611 nnfft_init_help(ths,m,nfft_flags,fftw_flags);
614 void nnfft_init(
nnfft_plan *ths,
int d,
int N_total,
int M_total,
int *N)
629 ths->
m=WINDOW_HELP_ESTIMATE_m;
638 ths->
N1[t] = ceil(1.5*ths->
N[t]);
641 if(ths->
N1[t]%2 != 0)
642 ths->
N1[t] = ths->
N1[t] +1;
645 ths->
nnfft_flags=PRE_PSI| PRE_PHI_HUT| MALLOC_X| MALLOC_V| MALLOC_F_HAT| MALLOC_F;
646 nfft_flags= PRE_PSI| PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE| NFFT_OMP_BLOCKWISE_ADJOINT;
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)
654 nnfft_init(ths,1,N1,M_total,&N1);
fftw_complex * f_hat
Fourier coefficients.
double * psi
precomputed data, matrix B
int * psi_index_g
only for thin B
void nnfft_precompute_full_psi(nnfft_plan *ths_plan)
computes all entries of B explicitly
void(* mv_trafo)(void *)
Transform.
fftw_complex * f_hat
Fourier coefficients.
unsigned nnfft_flags
flags for precomputation, malloc
int * N
cut-off-frequencies
double * v
nodes (in fourier domain)
int * psi_index_f
only for thin B
void nnfft_trafo(nnfft_plan *ths_plan)
user routines
nfft_plan * direct_plan
plan for the nfft
int m
cut-off parameter in time-domain
int aN1_total
aN1_total=aN1[0]* ...
void nnfft_precompute_phi_hut(nnfft_plan *ths_plan)
initialisation of direct transform
data structure for an NFFT (nonequispaced fast Fourier transform) plan with double precision ...
NFFT_INT M_total
Total number of samples.
data structure for an NNFFT (nonequispaced in time and frequency fast Fourier transform) plan with do...
void(* mv_adjoint)(void *)
Adjoint transform.
int * n
n=N1, for the window function
double * x
nodes (in time/spatial domain)
double * sigma
oversampling-factor
void * nfft_malloc(size_t n)
void nnfft_precompute_lin_psi(nnfft_plan *ths_plan)
create a lookup table
NFFT_INT N_total
Total number of Fourier coefficients.
double * x
Nodes in time/spatial domain, size is doubles.
double * c_phi_inv
precomputed data, matrix D