NFFT  3.3.0
nnfft.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
21 #include "config.h"
22 
23 #include <stdio.h>
24 #include <math.h>
25 #include <string.h>
26 #include <stdlib.h>
27 #ifdef HAVE_COMPLEX_H
28 #include <complex.h>
29 #endif
30 #include "nfft3.h"
31 #include "infft.h"
32 
33 
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
38 
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)
43 
44 #define MACRO_nndft_compute_trafo (*fj) += (*f_hat_k)*cexp(+ _Complex_I*omega);
45 
46 #define MACRO_nndft_compute_conjugated MACRO_nndft_compute_trafo
47 
48 #define MACRO_nndft_compute_adjoint (*f_hat_k) += (*fj)*cexp(+ _Complex_I*omega);
49 
50 #define MACRO_nndft_compute_transposed MACRO_nndft_compute_adjoint
51 
52 #define MACRO_nndft(which_one) \
53 void nnfft_ ## which_one ## _direct (nnfft_plan *ths) \
54 { \
55  int j; \
56  int t; \
57  int l; \
58  double _Complex *f_hat, *f; \
59  double _Complex *f_hat_k; \
60  double _Complex *fj; \
61  double omega; \
62  \
63  f_hat=ths->f_hat; f=ths->f; \
64  \
65  MACRO_nndft_init_result_ ## which_one \
66  \
67  for(j=0, fj=f; j<ths->M_total; j++, fj++) \
68  { \
69  for(l=0, f_hat_k=f_hat; l<ths->N_total; l++, f_hat_k++) \
70  { \
71  omega=0.0; \
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]; \
74  \
75  omega*= MACRO_nndft_sign_ ## which_one; \
76  \
77  MACRO_nndft_compute_ ## which_one \
78  \
79  } /* for(l) */ \
80  } /* for(j) */ \
81 } /* nndft_trafo */ \
82 
83 MACRO_nndft(trafo)
84 MACRO_nndft(adjoint)
85 
88 static void nnfft_uo(nnfft_plan *ths,int j,int *up,int *op,int act_dim)
89 {
90  double c;
91  int u,o;
92 
93  c = ths->v[j*ths->d+act_dim] * ths->n[act_dim];
94 
95  u = c; o = c;
96  if(c < 0)
97  u = u-1;
98  else
99  o = o+1;
100 
101  u = u - (ths->m); o = o + (ths->m);
102 
103  up[0]=u; op[0]=o;
104 }
105 
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));
111 
112 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_A { \
113  (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
114 }
115 
116 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_T { \
117  g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
118 }
119 
120 #define MACRO_nnfft_B_compute_A { \
121  (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
122 }
123 
124 #define MACRO_nnfft_B_compute_T { \
125  g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
126 }
127 
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]* \
131  (y[t2]-y_u[t2]))
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)
135 
136 #define MACRO_init_uo_l_lj_t { \
137  for(t = ths->d-1; t>=0; t--) \
138  { \
139  nnfft_uo(ths,j,&u[t],&o[t],t); \
140  l[t] = u[t]; \
141  lj[t] = 0; \
142  } /* for(t) */ \
143  t++; \
144 }
145 
146 #define MACRO_update_with_PRE_PSI_LIN { \
147  for(t2=t; t2<ths->d; t2++) \
148  { \
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]; \
152  } /* for(t2) */ \
153 }
154 
155 #define MACRO_update_phi_prod_ll_plain(which_one) { \
156  for(t2=t; t2<ths->d; t2++) \
157  { \
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]; \
161  /* 3/2 because of the (not needed) fftshift and to be in [0 aN1[t2]]?!*/\
162  } /* for(t2) */ \
163 }
164 
165 #define MACRO_count_uo_l_lj_t { \
166  for(t = ths->d-1; (t>0)&&(l[t]==o[t]); t--) \
167  { \
168  l[t] = u[t]; \
169  lj[t] = 0; \
170  } /* for(t) */ \
171  \
172  l[t]++; \
173  lj[t]++; \
174 }
175 
176 #define MACRO_nnfft_B(which_one) \
177 static inline void nnfft_B_ ## which_one (nnfft_plan *ths) \
178 { \
179  int lprod; \
180  int u[ths->d], o[ths->d]; \
181  int t, t2; \
182  int j; \
183  int l_L, ix; \
184  int l[ths->d]; \
185  int lj[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; \
190  double y[ths->d]; \
191  int y_u[ths->d]; \
192  \
193  f=ths->f_hat; g=ths->F; \
194  \
195  MACRO_nnfft_B_init_result_ ## which_one \
196  \
197  if(ths->nnfft_flags & PRE_FULL_PSI) \
198  { \
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; \
202  return; \
203  } \
204  \
205  phi_prod[0]=1; \
206  ll_plain[0]=0; \
207  \
208  for(t=0,lprod = 1; t<ths->d; t++) \
209  lprod *= (2*ths->m+2); \
210  \
211  if(ths->nnfft_flags & PRE_PSI) \
212  { \
213  for(j=0, fj=f; j<ths->N_total; j++, fj++) \
214  { \
215  MACRO_init_uo_l_lj_t; \
216  \
217  for(l_L=0; l_L<lprod; l_L++) \
218  { \
219  MACRO_update_phi_prod_ll_plain(with_PRE_PSI); \
220  \
221  MACRO_nnfft_B_compute_ ## which_one; \
222  \
223  MACRO_count_uo_l_lj_t; \
224  } /* for(l_L) */ \
225  } /* for(j) */ \
226  return; \
227  } /* if(PRE_PSI) */ \
228  \
229  if(ths->nnfft_flags & PRE_LIN_PSI) \
230  { \
231  for(j=0, fj=f; j<ths->N_total; j++, fj++) \
232  { \
233  MACRO_init_uo_l_lj_t; \
234  \
235  for(l_L=0; l_L<lprod; l_L++) \
236  { \
237  MACRO_update_with_PRE_PSI_LIN; \
238  \
239  MACRO_update_phi_prod_ll_plain(with_PRE_LIN_PSI); \
240  \
241  MACRO_nnfft_B_compute_ ## which_one; \
242  \
243  MACRO_count_uo_l_lj_t; \
244  } /* for(l_L) */ \
245  } /* for(j) */ \
246  return; \
247  } /* if(PRE_LIN_PSI) */ \
248  \
249  /* no precomputed psi at all */ \
250  for(j=0, fj=f; j<ths->N_total; j++, fj++) \
251  { \
252  \
253  MACRO_init_uo_l_lj_t; \
254  \
255  for(l_L=0; l_L<lprod; l_L++) \
256  { \
257  MACRO_update_phi_prod_ll_plain(without_PRE_PSI); \
258  \
259  MACRO_nnfft_B_compute_ ## which_one; \
260  \
261  MACRO_count_uo_l_lj_t; \
262  } /* for(l_L) */ \
263  } /* for(j) */ \
264 } /* nnfft_B */
265 
266 MACRO_nnfft_B(A)
267 MACRO_nnfft_B(T)
268 
269 static inline void nnfft_D (nnfft_plan *ths){
270  int j,t;
271  double tmp;
272 
273  if(ths->nnfft_flags & PRE_PHI_HUT)
274  {
275  for(j=0; j<ths->M_total; j++)
276  ths->f[j] *= ths->c_phi_inv[j];
277  }
278  else
279  {
280  for(j=0; j<ths->M_total; j++)
281  {
282  tmp = 1.0;
283  /* multiply with N1, because x was modified */
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)) );
286  ths->f[j] *= tmp;
287  }
288  }
289 }
290 
294 {
295  int j,t;
296 
297  nnfft_B_T(ths);
298 
299  for(j=0;j<ths->M_total;j++) {
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]);
302  }
303  }
304 
305 
306  /* allows for external swaps of ths->f */
307  ths->direct_plan->f = ths->f;
308 
309  nfft_trafo(ths->direct_plan);
310 
311  for(j=0;j<ths->M_total;j++) {
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]);
314  }
315  }
316 
317  nnfft_D(ths);
318 
319 } /* nnfft_trafo */
320 
321 void nnfft_adjoint(nnfft_plan *ths)
322 {
323  int j,t;
324 
325  nnfft_D(ths);
326 
327  for(j=0;j<ths->M_total;j++) {
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]);
330  }
331  }
332 
333  /* allows for external swaps of ths->f */
334  ths->direct_plan->f=ths->f;
335 
336  nfft_adjoint(ths->direct_plan);
337 
338  for(j=0;j<ths->M_total;j++) {
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]);
341  }
342  }
343 
344  nnfft_B_A(ths);
345 } /* nnfft_adjoint */
346 
350 {
351  int j;
352  int t;
353  double tmp;
354 
355  ths->c_phi_inv= (double*)nfft_malloc(ths->M_total*sizeof(double));
356 
357  for(j=0; j<ths->M_total; j++)
358  {
359  tmp = 1.0;
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));
362  ths->c_phi_inv[j]=tmp;
363  }
364 } /* nnfft_phi_hut */
365 
366 
370 {
371  int t;
372  int j;
373  double step;
375  nfft_precompute_lin_psi(ths->direct_plan);
376 
377  for (t=0; t<ths->d; t++)
378  {
379  step=((double)(ths->m+1))/(ths->K*ths->N1[t]);
380  for(j=0;j<=ths->K;j++)
381  {
382  ths->psi[(ths->K+1)*t + j] = PHI(ths->n[t],j*step,t);
383  } /* for(j) */
384  } /* for(t) */
385 }
386 
387 void nnfft_precompute_psi(nnfft_plan *ths)
388 {
389  int t;
390  int j;
391  int l;
392  int lj;
393  int u, o;
395  for (t=0; t<ths->d; t++)
396  for(j=0;j<ths->N_total;j++)
397  {
398  nnfft_uo(ths,j,&u,&o,t);
399 
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));
403  } /* for(j) */
404 
405  for(j=0;j<ths->M_total;j++) {
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]);
408  }
409  }
410 
411  nfft_precompute_psi(ths->direct_plan);
412 
413  for(j=0;j<ths->M_total;j++) {
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]);
416  }
417  }
418  /* for(t) */
419 } /* nfft_precompute_psi */
420 
421 
422 
427 {
428  int t,t2;
429  int j;
430  int l_L;
431  int l[ths->d];
432  int lj[ths->d];
433  int ll_plain[ths->d+1];
434  int lprod;
435  int u[ths->d], o[ths->d];
437  double phi_prod[ths->d+1];
438 
439  int ix,ix_old;
440 
441  for(j=0;j<ths->M_total;j++) {
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]);
444  }
445  }
446 
447  nnfft_precompute_psi(ths);
448 
449  nfft_precompute_full_psi(ths->direct_plan);
450 
451  for(j=0;j<ths->M_total;j++) {
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]);
454  }
455  }
456 
457  phi_prod[0]=1;
458  ll_plain[0]=0;
459 
460  for(t=0,lprod = 1; t<ths->d; t++)
461  lprod *= 2*ths->m+2;
462 
463  for(j=0,ix=0,ix_old=0; j<ths->N_total; j++)
464  {
465  MACRO_init_uo_l_lj_t;
466 
467  for(l_L=0; l_L<lprod; l_L++, ix++)
468  {
469  MACRO_update_phi_prod_ll_plain(without_PRE_PSI);
470 
471  ths->psi_index_g[ix]=ll_plain[ths->d];
472  ths->psi[ix]=phi_prod[ths->d];
473 
474  MACRO_count_uo_l_lj_t;
475  } /* for(l_L) */
476 
477 
478  ths->psi_index_f[j]=ix-ix_old;
479  ix_old=ix;
480  } /* for(j) */
481 }
482 
483 void nnfft_precompute_one_psi(nnfft_plan *ths)
484 {
485  if(ths->nnfft_flags & PRE_PSI)
486  nnfft_precompute_psi(ths);
487  if(ths->nnfft_flags & PRE_FULL_PSI)
489  if(ths->nnfft_flags & PRE_LIN_PSI)
492  if(ths->nnfft_flags & PRE_PHI_HUT)
494 }
495 
496 static void nnfft_init_help(nnfft_plan *ths, int m2, unsigned nfft_flags, unsigned fftw_flags)
497 {
498  int t;
499  int lprod;
500  int N2[ths->d];
501 
502  ths->aN1 = (int*) nfft_malloc(ths->d*sizeof(int));
503 
504  ths->a = (double*) nfft_malloc(ths->d*sizeof(double));
505 
506  ths->sigma = (double*) nfft_malloc(ths->d*sizeof(double));
507 
508  ths->n = ths->N1;
509 
510  ths->aN1_total=1;
511 
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]);
515  /* aN1 should be even */
516  if(ths->aN1[t]%2 != 0)
517  ths->aN1[t] = ths->aN1[t] +1;
518 
519  ths->aN1_total*=ths->aN1[t];
520  ths->sigma[t] = ((double) ths->N1[t] )/((double) ths->N[t]);;
521 
522  /* take the same oversampling factor in the inner NFFT */
523  N2[t] = ceil(ths->sigma[t]*(ths->aN1[t]));
524 
525  /* N2 should be even */
526  if(N2[t]%2 != 0)
527  N2[t] = N2[t] +1;
528  }
529 
530  WINDOW_HELP_INIT
531 
532  if(ths->nnfft_flags & MALLOC_X)
533  ths->x = (double*)nfft_malloc(ths->d*ths->M_total*sizeof(double));
534  if(ths->nnfft_flags & MALLOC_F){
535  ths->f=(double _Complex*)nfft_malloc(ths->M_total*sizeof(double _Complex));
536  }
537  if(ths->nnfft_flags & MALLOC_V)
538  ths->v = (double*)nfft_malloc(ths->d*ths->N_total*sizeof(double));
539  if(ths->nnfft_flags & MALLOC_F_HAT)
540  ths->f_hat = (double _Complex*)nfft_malloc(ths->N_total*sizeof(double _Complex));
541 
542  //BUGFIX SUSE 2
544 // if(ths->nnfft_flags & PRE_PHI_HUT)
545 // nnfft_precompute_phi_hut(ths);
546 
547  if(ths->nnfft_flags & PRE_LIN_PSI)
548  {
549  ths->K=(1U<< 10)*(ths->m+1);
550  ths->psi = (double*) nfft_malloc((ths->K+1)*ths->d*sizeof(double));
551  }
552 
553  if(ths->nnfft_flags & PRE_PSI){
554  ths->psi = (double*)nfft_malloc(ths->N_total*ths->d*(2*ths->m+2)*sizeof(double));
555  }
556 
557  if(ths->nnfft_flags & PRE_FULL_PSI)
558  {
559  for(t=0,lprod = 1; t<ths->d; t++)
560  lprod *= 2*ths->m+2;
561 
562  ths->psi = (double*)nfft_malloc(ths->N_total*lprod*sizeof(double));
563 
564  ths->psi_index_f = (int*) nfft_malloc(ths->N_total*sizeof(int));
565  ths->psi_index_g = (int*) nfft_malloc(ths->N_total*lprod*sizeof(int));
566  }
567  ths->direct_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
568  nfft_init_guru(ths->direct_plan, ths->d, ths->aN1, ths->M_total, N2, m2,
569  nfft_flags, fftw_flags);
570  ths->direct_plan->x = ths->x;
571 
572  ths->direct_plan->f = ths->f;
573  ths->F = ths->direct_plan->f_hat;
574 
575  ths->mv_trafo = (void (*) (void* ))nnfft_trafo;
576  ths->mv_adjoint = (void (*) (void* ))nnfft_adjoint;
577 }
578 
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)
581 {
582  int t;
584  unsigned nfft_flags;
585  unsigned fftw_flags;
586 
587  ths->d= d;
588  ths->M_total= M_total;
589  ths->N_total= N_total;
590  ths->m= m;
591  ths->nnfft_flags= 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;
594 
595  if(ths->nnfft_flags & PRE_PSI)
596  nfft_flags = nfft_flags | PRE_PSI;
597 
598  if(ths->nnfft_flags & PRE_FULL_PSI)
599  nfft_flags = nfft_flags | PRE_FULL_PSI;
600 
601  if(ths->nnfft_flags & PRE_LIN_PSI)
602  nfft_flags = nfft_flags | PRE_LIN_PSI;
603 
604  ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
605  ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
606 
607  for(t=0; t<d; t++) {
608  ths->N[t] = N[t];
609  ths->N1[t] = N1[t];
610  }
611  nnfft_init_help(ths,m,nfft_flags,fftw_flags);
612 }
613 
614 void nnfft_init(nnfft_plan *ths, int d, int N_total, int M_total, int *N)
615 {
616  int t;
618  unsigned nfft_flags;
619  unsigned fftw_flags;
620  ths->d = d;
621  ths->M_total = M_total;
622  ths->N_total = N_total;
623  /* m should be greater to get the same accuracy as the nfft */
624 /* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
625 
626  WINDOW_HELP_ESTIMATE_m;
627 */
628  //BUGFIX SUSE 1
629 ths->m=WINDOW_HELP_ESTIMATE_m;
630 
631 
632  ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
633  ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
634 
635  for(t=0; t<d; t++) {
636  ths->N[t] = N[t];
637  /* the standard oversampling factor in the nnfft is 1.5 */
638  ths->N1[t] = ceil(1.5*ths->N[t]);
639 
640  /* N1 should be even */
641  if(ths->N1[t]%2 != 0)
642  ths->N1[t] = ths->N1[t] +1;
643  }
644 
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;
647 
648  fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
649  nnfft_init_help(ths,ths->m,nfft_flags,fftw_flags);
650 }
651 
652 void nnfft_init_1d(nnfft_plan *ths,int N1, int M_total)
653 {
654  nnfft_init(ths,1,N1,M_total,&N1);
655 }
656 
657 void nnfft_finalize(nnfft_plan *ths)
658 {
659  nfft_finalize(ths->direct_plan);
660  nfft_free(ths->direct_plan);
661 
662  nfft_free(ths->aN1);
663  nfft_free(ths->N);
664  nfft_free(ths->N1);
665 
666  if(ths->nnfft_flags & PRE_FULL_PSI)
667  {
668  nfft_free(ths->psi_index_g);
669  nfft_free(ths->psi_index_f);
670  nfft_free(ths->psi);
671  }
672 
673  if(ths->nnfft_flags & PRE_PSI)
674  nfft_free(ths->psi);
675 
676  if(ths->nnfft_flags & PRE_LIN_PSI)
677  nfft_free(ths->psi);
678 
679  if(ths->nnfft_flags & PRE_PHI_HUT)
680  nfft_free(ths->c_phi_inv);
681 
682  if(ths->nnfft_flags & MALLOC_F)
683  nfft_free(ths->f);
684 
685  if(ths->nnfft_flags & MALLOC_F_HAT)
686  nfft_free(ths->f_hat);
687 
688  if(ths->nnfft_flags & MALLOC_X)
689  nfft_free(ths->x);
690 
691  if(ths->nnfft_flags & MALLOC_V)
692  nfft_free(ths->v);
693 }
int d
dimension, rank
Definition: nfft3.h:426
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:426
double * psi
precomputed data, matrix B
Definition: nfft3.h:426
int * psi_index_g
only for thin B
Definition: nfft3.h:426
void nnfft_precompute_full_psi(nnfft_plan *ths_plan)
computes all entries of B explicitly
Definition: nnfft.c:426
void(* mv_trafo)(void *)
Transform.
Definition: nfft3.h:426
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:194
unsigned nnfft_flags
flags for precomputation, malloc
Definition: nfft3.h:426
fftw_complex * f
Samples.
Definition: nfft3.h:426
int * N
cut-off-frequencies
Definition: nfft3.h:426
int K
number of precomp.
Definition: nfft3.h:426
int * aN1
sigma*a*N
Definition: nfft3.h:426
double * v
nodes (in fourier domain)
Definition: nfft3.h:426
int * psi_index_f
only for thin B
Definition: nfft3.h:426
void nnfft_trafo(nnfft_plan *ths_plan)
user routines
Definition: nnfft.c:293
nfft_plan * direct_plan
plan for the nfft
Definition: nfft3.h:426
int m
cut-off parameter in time-domain
Definition: nfft3.h:426
fftw_complex * f
Samples.
Definition: nfft3.h:194
int aN1_total
aN1_total=aN1[0]* ...
Definition: nfft3.h:426
void nfft_free(void *p)
void nnfft_precompute_phi_hut(nnfft_plan *ths_plan)
initialisation of direct transform
Definition: nnfft.c:349
data structure for an NFFT (nonequispaced fast Fourier transform) plan with double precision ...
Definition: nfft3.h:194
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:426
data structure for an NNFFT (nonequispaced in time and frequency fast Fourier transform) plan with do...
Definition: nfft3.h:426
void(* mv_adjoint)(void *)
Adjoint transform.
Definition: nfft3.h:426
int * n
n=N1, for the window function
Definition: nfft3.h:426
double * x
nodes (in time/spatial domain)
Definition: nfft3.h:426
double * sigma
oversampling-factor
Definition: nfft3.h:426
void * nfft_malloc(size_t n)
void nnfft_precompute_lin_psi(nnfft_plan *ths_plan)
create a lookup table
Definition: nnfft.c:369
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:426
double * x
Nodes in time/spatial domain, size is doubles.
Definition: nfft3.h:194
double * a
1 + 2*m/N1
Definition: nfft3.h:426
int * N1
sigma*N
Definition: nfft3.h:426
double * c_phi_inv
precomputed data, matrix D
Definition: nfft3.h:426