NFFT 3.5.3alpha
nnfft.c
1/*
2 * Copyright (c) 2002, 2017 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#include "config.h"
20
21#include <stdio.h>
22#include <math.h>
23#include <string.h>
24#include <stdlib.h>
25#ifdef HAVE_COMPLEX_H
26#include <complex.h>
27#endif
28#include "nfft3.h"
29#include "infft.h"
30
31
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
36
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)
41
42#define MACRO_nndft_compute_trafo (*fj) += (*f_hat_k)*cexp(+ _Complex_I*omega);
43
44#define MACRO_nndft_compute_conjugated MACRO_nndft_compute_trafo
45
46#define MACRO_nndft_compute_adjoint (*f_hat_k) += (*fj)*cexp(+ _Complex_I*omega);
47
48#define MACRO_nndft_compute_transposed MACRO_nndft_compute_adjoint
49
50#define MACRO_nndft(which_one) \
51void nnfft_ ## which_one ## _direct (nnfft_plan *ths) \
52{ \
53 int j; \
54 int t; \
55 int l; \
56 double _Complex *f_hat, *f; \
57 double _Complex *f_hat_k; \
58 double _Complex *fj; \
59 double omega; \
60 \
61 f_hat=ths->f_hat; f=ths->f; \
62 \
63 MACRO_nndft_init_result_ ## which_one \
64 \
65 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
66 { \
67 for(l=0, f_hat_k=f_hat; l<ths->N_total; l++, f_hat_k++) \
68 { \
69 omega=0.0; \
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]; \
72 \
73 omega*= MACRO_nndft_sign_ ## which_one; \
74 \
75 MACRO_nndft_compute_ ## which_one \
76 \
77 } /* for(l) */ \
78 } /* for(j) */ \
79} /* nndft_trafo */ \
80
81MACRO_nndft(trafo)
82MACRO_nndft(adjoint)
83
84
86static void nnfft_uo(nnfft_plan *ths,int j,int *up,int *op,int act_dim)
87{
88 double c;
89 int u,o;
90
91 c = ths->v[j*ths->d+act_dim] * ths->n[act_dim];
92
93 u = c; o = c;
94 if(c < 0)
95 u = u-1;
96 else
97 o = o+1;
98
99 u = u - (ths->m); o = o + (ths->m);
100
101 up[0]=u; op[0]=o;
102}
103
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));
109
110#define MACRO_nnfft_B_PRE_FULL_PSI_compute_A { \
111 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
112}
113
114#define MACRO_nnfft_B_PRE_FULL_PSI_compute_T { \
115 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
116}
117
118#define MACRO_nnfft_B_compute_A { \
119 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
120}
121
122#define MACRO_nnfft_B_compute_T { \
123 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
124}
125
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]* \
129 (y[t2]-y_u[t2]))
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)
133
134#define MACRO_init_uo_l_lj_t { \
135 for(t = ths->d-1; t>=0; t--) \
136 { \
137 nnfft_uo(ths,j,&u[t],&o[t],t); \
138 l[t] = u[t]; \
139 lj[t] = 0; \
140 } /* for(t) */ \
141 t++; \
142}
143
144#define MACRO_update_with_PRE_PSI_LIN { \
145 for(t2=t; t2<ths->d; t2++) \
146 { \
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]; \
150 } /* for(t2) */ \
151}
152
153#define MACRO_update_phi_prod_ll_plain(which_one) { \
154 for(t2=t; t2<ths->d; t2++) \
155 { \
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]; \
159 /* 3/2 because of the (not needed) fftshift and to be in [0 aN1[t2]]?!*/\
160 } /* for(t2) */ \
161}
162
163#define MACRO_count_uo_l_lj_t { \
164 for(t = ths->d-1; (t>0)&&(l[t]==o[t]); t--) \
165 { \
166 l[t] = u[t]; \
167 lj[t] = 0; \
168 } /* for(t) */ \
169 \
170 l[t]++; \
171 lj[t]++; \
172}
173
174#define MACRO_nnfft_B(which_one) \
175static inline void nnfft_B_ ## which_one (nnfft_plan *ths) \
176{ \
177 int lprod; \
178 int u[ths->d], o[ths->d]; \
179 int t, t2; \
180 int j; \
181 int l_L, ix; \
182 int l[ths->d]; \
183 int lj[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; \
188 double y[ths->d]; \
189 int y_u[ths->d]; \
190 \
191 f=ths->f_hat; g=ths->F; \
192 \
193 MACRO_nnfft_B_init_result_ ## which_one \
194 \
195 if(ths->nnfft_flags & PRE_FULL_PSI) \
196 { \
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; \
200 return; \
201 } \
202 \
203 phi_prod[0]=1; \
204 ll_plain[0]=0; \
205 \
206 for(t=0,lprod = 1; t<ths->d; t++) \
207 lprod *= (2*ths->m+2); \
208 \
209 if(ths->nnfft_flags & PRE_PSI) \
210 { \
211 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
212 { \
213 MACRO_init_uo_l_lj_t; \
214 \
215 for(l_L=0; l_L<lprod; l_L++) \
216 { \
217 MACRO_update_phi_prod_ll_plain(with_PRE_PSI); \
218 \
219 MACRO_nnfft_B_compute_ ## which_one; \
220 \
221 MACRO_count_uo_l_lj_t; \
222 } /* for(l_L) */ \
223 } /* for(j) */ \
224 return; \
225 } /* if(PRE_PSI) */ \
226 \
227 if(ths->nnfft_flags & PRE_LIN_PSI) \
228 { \
229 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
230 { \
231 MACRO_init_uo_l_lj_t; \
232 \
233 for(l_L=0; l_L<lprod; l_L++) \
234 { \
235 MACRO_update_with_PRE_PSI_LIN; \
236 \
237 MACRO_update_phi_prod_ll_plain(with_PRE_LIN_PSI); \
238 \
239 MACRO_nnfft_B_compute_ ## which_one; \
240 \
241 MACRO_count_uo_l_lj_t; \
242 } /* for(l_L) */ \
243 } /* for(j) */ \
244 return; \
245 } /* if(PRE_LIN_PSI) */ \
246 \
247 /* no precomputed psi at all */ \
248 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
249 { \
250 \
251 MACRO_init_uo_l_lj_t; \
252 \
253 for(l_L=0; l_L<lprod; l_L++) \
254 { \
255 MACRO_update_phi_prod_ll_plain(without_PRE_PSI); \
256 \
257 MACRO_nnfft_B_compute_ ## which_one; \
258 \
259 MACRO_count_uo_l_lj_t; \
260 } /* for(l_L) */ \
261 } /* for(j) */ \
262} /* nnfft_B */
263
264MACRO_nnfft_B(A)
265MACRO_nnfft_B(T)
266
267static inline void nnfft_D (nnfft_plan *ths){
268 int j,t;
269 double tmp;
270
271 if(ths->nnfft_flags & PRE_PHI_HUT)
272 {
273 for(j=0; j<ths->M_total; j++)
274 ths->f[j] *= ths->c_phi_inv[j];
275 }
276 else
277 {
278 for(j=0; j<ths->M_total; j++)
279 {
280 tmp = 1.0;
281 /* multiply with N1, because x was modified */
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)) );
284 ths->f[j] *= tmp;
285 }
286 }
287}
288
292{
293 int j,t;
294
295 nnfft_B_T(ths);
296
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]);
300 }
301 }
302
303
304 /* allows for external swaps of ths->f */
305 ths->direct_plan->f = ths->f;
306
307 nfft_trafo(ths->direct_plan);
308
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]);
312 }
313 }
314
315 nnfft_D(ths);
316
317} /* nnfft_trafo */
318
319void nnfft_adjoint(nnfft_plan *ths)
320{
321 int j,t;
322
323 nnfft_D(ths);
324
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]);
328 }
329 }
330
331 /* allows for external swaps of ths->f */
332 ths->direct_plan->f=ths->f;
333
334 nfft_adjoint(ths->direct_plan);
335
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]);
339 }
340 }
341
342 nnfft_B_A(ths);
343} /* nnfft_adjoint */
344
348{
349 int j;
350 int t;
351 double tmp;
352
353 ths->c_phi_inv= (double*)nfft_malloc(ths->M_total*sizeof(double));
354
355 for(j=0; j<ths->M_total; j++)
356 {
357 tmp = 1.0;
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;
361 }
362} /* nnfft_phi_hut */
363
364
368{
369 int t;
370 int j;
371 double step;
373 nfft_precompute_lin_psi(ths->direct_plan);
374
375 for (t=0; t<ths->d; t++)
376 {
377 step=((double)(ths->m+1))/(ths->K*ths->N1[t]);
378 for(j=0;j<=ths->K;j++)
379 {
380 ths->psi[(ths->K+1)*t + j] = PHI(ths->n[t],j*step,t);
381 } /* for(j) */
382 } /* for(t) */
383}
384
385void nnfft_precompute_psi(nnfft_plan *ths)
386{
387 int t;
388 int j;
389 int l;
390 int lj;
391 int u, o;
393 for (t=0; t<ths->d; t++)
394 for(j=0;j<ths->N_total;j++)
395 {
396 nnfft_uo(ths,j,&u,&o,t);
397
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));
401 } /* for(j) */
402
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]);
406 }
407 }
408
409 nfft_precompute_psi(ths->direct_plan);
410
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]);
414 }
415 }
416 /* for(t) */
417} /* nfft_precompute_psi */
418
419
420
425{
426 int t,t2;
427 int j;
428 int l_L;
429 int l[ths->d];
430 int lj[ths->d];
431 int ll_plain[ths->d+1];
432 int lprod;
433 int u[ths->d], o[ths->d];
435 double phi_prod[ths->d+1];
436
437 int ix,ix_old;
438
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]);
442 }
443 }
444
445 nnfft_precompute_psi(ths);
446
447 nfft_precompute_full_psi(ths->direct_plan);
448
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]);
452 }
453 }
454
455 phi_prod[0]=1;
456 ll_plain[0]=0;
457
458 for(t=0,lprod = 1; t<ths->d; t++)
459 lprod *= 2*ths->m+2;
460
461 for(j=0,ix=0,ix_old=0; j<ths->N_total; j++)
462 {
463 MACRO_init_uo_l_lj_t;
464
465 for(l_L=0; l_L<lprod; l_L++, ix++)
466 {
467 MACRO_update_phi_prod_ll_plain(without_PRE_PSI);
468
469 ths->psi_index_g[ix]=ll_plain[ths->d];
470 ths->psi[ix]=phi_prod[ths->d];
471
472 MACRO_count_uo_l_lj_t;
473 } /* for(l_L) */
474
475
476 ths->psi_index_f[j]=ix-ix_old;
477 ix_old=ix;
478 } /* for(j) */
479}
480
481void nnfft_precompute_one_psi(nnfft_plan *ths)
482{
483 if(ths->nnfft_flags & PRE_PSI)
484 nnfft_precompute_psi(ths);
485 if(ths->nnfft_flags & PRE_FULL_PSI)
487 if(ths->nnfft_flags & PRE_LIN_PSI)
490 if(ths->nnfft_flags & PRE_PHI_HUT)
492}
493
494static void nnfft_init_help(nnfft_plan *ths, int m2, unsigned nfft_flags, unsigned fftw_flags)
495{
496 int t;
497 int lprod;
498 int N2[ths->d];
499
500 ths->aN1 = (int*) nfft_malloc(ths->d*sizeof(int));
501
502 ths->a = (double*) nfft_malloc(ths->d*sizeof(double));
503
504 ths->sigma = (double*) nfft_malloc(ths->d*sizeof(double));
505
506 ths->n = ths->N1;
507
508 ths->aN1_total=1;
509
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]);
513 /* aN1 should be even */
514 if(ths->aN1[t]%2 != 0)
515 ths->aN1[t] = ths->aN1[t] +1;
516
517 ths->aN1_total*=ths->aN1[t];
518 ths->sigma[t] = ((double) ths->N1[t] )/((double) ths->N[t]);;
519
520 /* take the same oversampling factor in the inner NFFT */
521 N2[t] = ceil(ths->sigma[t]*(ths->aN1[t]));
522
523 /* N2 should be even */
524 if(N2[t]%2 != 0)
525 N2[t] = N2[t] +1;
526 }
527
528 WINDOW_HELP_INIT
529
530 if(ths->nnfft_flags & MALLOC_X)
531 ths->x = (double*)nfft_malloc(ths->d*ths->M_total*sizeof(double));
532 if(ths->nnfft_flags & MALLOC_F){
533 ths->f=(double _Complex*)nfft_malloc(ths->M_total*sizeof(double _Complex));
534 }
535 if(ths->nnfft_flags & MALLOC_V)
536 ths->v = (double*)nfft_malloc(ths->d*ths->N_total*sizeof(double));
537 if(ths->nnfft_flags & MALLOC_F_HAT)
538 ths->f_hat = (double _Complex*)nfft_malloc(ths->N_total*sizeof(double _Complex));
539
540 //BUGFIX SUSE 2
542// if(ths->nnfft_flags & PRE_PHI_HUT)
543// nnfft_precompute_phi_hut(ths);
544
545 if(ths->nnfft_flags & PRE_LIN_PSI)
546 {
547 ths->K=(1U<< 10)*(ths->m+1);
548 ths->psi = (double*) nfft_malloc((ths->K+1)*ths->d*sizeof(double));
549 }
550
551 if(ths->nnfft_flags & PRE_PSI){
552 ths->psi = (double*)nfft_malloc(ths->N_total*ths->d*(2*ths->m+2)*sizeof(double));
553 }
554
555 if(ths->nnfft_flags & PRE_FULL_PSI)
556 {
557 for(t=0,lprod = 1; t<ths->d; t++)
558 lprod *= 2*ths->m+2;
559
560 ths->psi = (double*)nfft_malloc(ths->N_total*lprod*sizeof(double));
561
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));
564 }
565 ths->direct_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
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;
569
570 ths->direct_plan->f = ths->f;
571 ths->F = ths->direct_plan->f_hat;
572
573 ths->mv_trafo = (void (*) (void* ))nnfft_trafo;
574 ths->mv_adjoint = (void (*) (void* ))nnfft_adjoint;
575}
576
577void nnfft_init_guru(nnfft_plan *ths, int d, int N_total, int M_total, int *N, int *N1,
578 int m, unsigned nnfft_flags)
579{
580 int t;
582 unsigned nfft_flags;
583 unsigned fftw_flags;
584
585 ths->d= d;
586 ths->M_total= M_total;
587 ths->N_total= N_total;
588 ths->m= m;
589 ths->nnfft_flags= nnfft_flags;
590 fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
591 nfft_flags= PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT|
592 ((d == 1) ? FFT_OUT_OF_PLACE : 0U) | NFFT_OMP_BLOCKWISE_ADJOINT;
593
594 if(ths->nnfft_flags & PRE_PSI)
595 nfft_flags = nfft_flags | PRE_PSI;
596
597 if(ths->nnfft_flags & PRE_FULL_PSI)
598 nfft_flags = nfft_flags | PRE_FULL_PSI;
599
600 if(ths->nnfft_flags & PRE_LIN_PSI)
601 nfft_flags = nfft_flags | PRE_LIN_PSI;
602
603 ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
604 ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
605
606 for(t=0; t<d; t++) {
607 ths->N[t] = N[t];
608 ths->N1[t] = N1[t];
609 }
610 nnfft_init_help(ths,m,nfft_flags,fftw_flags);
611}
612
613void nnfft_init(nnfft_plan *ths, int d, int N_total, int M_total, int *N)
614{
615 int t;
617 unsigned nfft_flags;
618 unsigned fftw_flags;
619 ths->d = d;
620 ths->M_total = M_total;
621 ths->N_total = N_total;
622 /* m should be greater to get the same accuracy as the nfft */
623/* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
624
625 WINDOW_HELP_ESTIMATE_m;
626*/
627 //BUGFIX SUSE 1
628ths->m=WINDOW_HELP_ESTIMATE_m;
629
630
631 ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
632 ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
633
634 for(t=0; t<d; t++) {
635 ths->N[t] = N[t];
636 /* the standard oversampling factor in the nnfft is 1.5 */
637 ths->N1[t] = ceil(1.5*ths->N[t]);
638
639 /* N1 should be even */
640 if(ths->N1[t]%2 != 0)
641 ths->N1[t] = ths->N1[t] +1;
642 }
643
646 ((d == 1) ? FFT_OUT_OF_PLACE : 0U)| 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
652void nnfft_init_1d(nnfft_plan *ths,int N1, int M_total)
653{
654 nnfft_init(ths,1,N1,M_total,&N1);
655}
656
657void nnfft_finalize(nnfft_plan *ths)
658{
659 nfft_finalize(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 {
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}
#define MALLOC_F_HAT
Definition nfft3.h:194
#define MALLOC_X
Definition nfft3.h:193
#define PRE_FULL_PSI
Definition nfft3.h:192
#define FFT_OUT_OF_PLACE
Definition nfft3.h:196
#define PRE_PSI
Definition nfft3.h:191
#define MALLOC_F
Definition nfft3.h:195
#define PRE_LIN_PSI
Definition nfft3.h:189
#define FFTW_INIT
Definition nfft3.h:197
#define PRE_PHI_HUT
Definition nfft3.h:187
#define MALLOC_V
Definition nfft3.h:429
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.
void nnfft_precompute_lin_psi(nnfft_plan *ths_plan)
create a lookup table
Definition nnfft.c:367
void nnfft_precompute_full_psi(nnfft_plan *ths_plan)
computes all entries of B explicitly
Definition nnfft.c:424
void nnfft_trafo(nnfft_plan *ths_plan)
user routines
Definition nnfft.c:291
void nnfft_precompute_phi_hut(nnfft_plan *ths_plan)
initialisation of direct transform
Definition nnfft.c:347
void * nfft_malloc(size_t n)
void nfft_free(void *p)
data structure for an NNFFT (nonequispaced in time and frequency fast Fourier transform) plan with do...
Definition nfft3.h:425
double * v
nodes (in fourier domain)
Definition nfft3.h:425
double * a
1 + 2*m/N1
Definition nfft3.h:425
int m
cut-off parameter in time-domain
Definition nfft3.h:425
int * psi_index_f
only for thin B
Definition nfft3.h:425
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:425
int * N
cut-off-frequencies
Definition nfft3.h:425
nfft_plan * direct_plan
plan for the nfft
Definition nfft3.h:425
int K
number of precomp.
Definition nfft3.h:425
unsigned nnfft_flags
flags for precomputation, malloc
Definition nfft3.h:425
double * sigma
oversampling-factor
Definition nfft3.h:425
int * n
n=N1, for the window function
Definition nfft3.h:425
NFFT_INT N_total
Total number of Fourier coefficients.
Definition nfft3.h:425
int * N1
sigma*N
Definition nfft3.h:425
fftw_complex * f
Samples.
Definition nfft3.h:425
void(* mv_trafo)(void *)
Transform.
Definition nfft3.h:425
double * c_phi_inv
precomputed data, matrix D
Definition nfft3.h:425
int d
dimension, rank
Definition nfft3.h:425
int * psi_index_g
only for thin B
Definition nfft3.h:425
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:425
double * x
nodes (in time/spatial domain)
Definition nfft3.h:425
int aN1_total
aN1_total=aN1[0]* ... *aN1[d-1]
Definition nfft3.h:425
double * psi
precomputed data, matrix B
Definition nfft3.h:425
void(* mv_adjoint)(void *)
Adjoint transform.
Definition nfft3.h:425
int * aN1
sigma*a*N
Definition nfft3.h:425