NFFT  3.5.3
polar_fft_test.c
Go to the documentation of this file.
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 
28 #include <math.h>
29 #include <stdlib.h>
30 #include <complex.h>
31 
32 #define NFFT_PRECISION_DOUBLE
33 
34 #include "nfft3mp.h"
35 
73 static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
74 {
75  int t, r;
76  NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
77 
78  for (t = -T / 2; t < T / 2; t++)
79  {
80  for (r = -S / 2; r < S / 2; r++)
81  {
82  x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
83  x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
84  if (r == 0)
85  w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
86  else
87  w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
88  }
89  }
90 
91  return T * S;
92 }
93 
95 static int polar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S,
96  int m)
97 {
98  int j, k;
99  NFFT(plan) my_nfft_plan;
101  NFFT_R *x, *w;
103  int N[2], n[2];
104  int M = T * S;
106  N[0] = NN;
107  n[0] = 2 * N[0];
108  N[1] = NN;
109  n[1] = 2 * N[1];
111  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
112  if (x == NULL)
113  return EXIT_FAILURE;
114 
115  w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
116  if (w == NULL)
117  return EXIT_FAILURE;
118 
120  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
122  FFTW_MEASURE);
123 
125  polar_grid(T, S, x, w);
126  for (j = 0; j < my_nfft_plan.M_total; j++)
127  {
128  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
129  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
130  }
131 
133  for (k = 0; k < my_nfft_plan.N_total; k++)
134  my_nfft_plan.f_hat[k] = f_hat[k];
135 
137  NFFT(trafo_direct)(&my_nfft_plan);
138 
140  for (j = 0; j < my_nfft_plan.M_total; j++)
141  f[j] = my_nfft_plan.f[j];
142 
144  NFFT(finalize)(&my_nfft_plan);
145  NFFT(free)(x);
146  NFFT(free)(w);
147 
148  return EXIT_SUCCESS;
149 }
150 
152 static int polar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S,
153  int m)
154 {
155  int j, k;
156  NFFT(plan) my_nfft_plan;
158  NFFT_R *x, *w;
160  int N[2], n[2];
161  int M = T * S;
163  N[0] = NN;
164  n[0] = 2 * N[0];
165  N[1] = NN;
166  n[1] = 2 * N[1];
168  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
169  if (x == NULL)
170  return EXIT_FAILURE;
171 
172  w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
173  if (w == NULL)
174  return EXIT_FAILURE;
175 
177  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
180  FFTW_MEASURE | FFTW_DESTROY_INPUT);
181 
183  polar_grid(T, S, x, w);
184  for (j = 0; j < my_nfft_plan.M_total; j++)
185  {
186  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
187  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
188  }
189 
191  if (my_nfft_plan.flags & PRE_LIN_PSI)
192  NFFT(precompute_lin_psi)(&my_nfft_plan);
193 
194  if (my_nfft_plan.flags & PRE_PSI)
195  NFFT(precompute_psi)(&my_nfft_plan);
196 
197  if (my_nfft_plan.flags & PRE_FULL_PSI)
198  NFFT(precompute_full_psi)(&my_nfft_plan);
199 
201  for (k = 0; k < my_nfft_plan.N_total; k++)
202  my_nfft_plan.f_hat[k] = f_hat[k];
203 
205  NFFT(trafo)(&my_nfft_plan);
206 
208  for (j = 0; j < my_nfft_plan.M_total; j++)
209  f[j] = my_nfft_plan.f[j];
210 
212  NFFT(finalize)(&my_nfft_plan);
213  NFFT(free)(x);
214  NFFT(free)(w);
215 
216  return EXIT_SUCCESS;
217 }
218 
220 static int inverse_polar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat,
221  int NN, int max_i, int m)
222 {
223  int j, k;
224  NFFT(plan) my_nfft_plan;
225  SOLVER(plan_complex) my_infft_plan;
227  NFFT_R *x, *w;
228  int l;
230  int N[2], n[2];
231  int M = T * S;
233  N[0] = NN;
234  n[0] = 2 * N[0];
235  N[1] = NN;
236  n[1] = 2 * N[1];
238  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
239  if (x == NULL)
240  return EXIT_FAILURE;
241 
242  w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
243  if (w == NULL)
244  return EXIT_FAILURE;
245 
247  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
250  FFTW_MEASURE | FFTW_DESTROY_INPUT);
251 
253  SOLVER(init_advanced_complex)(&my_infft_plan,
254  (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
255 
257  polar_grid(T, S, x, w);
258  for (j = 0; j < my_nfft_plan.M_total; j++)
259  {
260  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
261  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
262  my_infft_plan.y[j] = f[j];
263  my_infft_plan.w[j] = w[j];
264  }
265 
267  if (my_nfft_plan.flags & PRE_LIN_PSI)
268  NFFT(precompute_lin_psi)(&my_nfft_plan);
269 
270  if (my_nfft_plan.flags & PRE_PSI)
271  NFFT(precompute_psi)(&my_nfft_plan);
272 
273  if (my_nfft_plan.flags & PRE_FULL_PSI)
274  NFFT(precompute_full_psi)(&my_nfft_plan);
275 
277  if (my_infft_plan.flags & PRECOMPUTE_DAMP)
278  for (j = 0; j < my_nfft_plan.N[0]; j++)
279  for (k = 0; k < my_nfft_plan.N[1]; k++)
280  {
281  my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
282  NFFT_M(sqrt)(
283  NFFT_M(pow)((NFFT_R) (j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
284  + NFFT_M(pow)((NFFT_R) (k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
285  > ((NFFT_R) (my_nfft_plan.N[0] / 2)) ? 0 : 1);
286  }
287 
289  for (k = 0; k < my_nfft_plan.N_total; k++)
290  my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
291 
293  SOLVER(before_loop_complex)(&my_infft_plan);
294 
295  if (max_i < 1)
296  {
297  l = 1;
298  for (k = 0; k < my_nfft_plan.N_total; k++)
299  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
300  }
301  else
302  {
303  for (l = 1; l <= max_i; l++)
304  {
305  SOLVER(loop_one_step_complex)(&my_infft_plan);
306  }
307  }
308 
310  for (k = 0; k < my_nfft_plan.N_total; k++)
311  f_hat[k] = my_infft_plan.f_hat_iter[k];
312 
314  SOLVER(finalize_complex)(&my_infft_plan);
315  NFFT(finalize)(&my_nfft_plan);
316  NFFT(free)(x);
317  NFFT(free)(w);
318 
319  return EXIT_SUCCESS;
320 }
321 
323 int main(int argc, char **argv)
324 {
325  int N;
326  int T, S;
327  int M;
328  NFFT_R *x, *w;
329  NFFT_C *f_hat, *f, *f_direct, *f_tilde;
330  int k;
331  int max_i;
332  int m = 1;
333  NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
334  FILE *fp1, *fp2;
335  char filename[30];
336 
337  if (argc != 4)
338  {
339  printf("polar_fft_test N T R \n");
340  printf("\n");
341  printf("N polar FFT of size NxN \n");
342  printf("T number of slopes \n");
343  printf("R number of offsets \n");
344  exit(EXIT_FAILURE);
345  }
346 
347  N = atoi(argv[1]);
348  T = atoi(argv[2]);
349  S = atoi(argv[3]);
350  printf("N=%d, polar grid with T=%d, R=%d => ", N, T, S);
351 
352  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * 5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R)));
353  w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R)));
354 
355  f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
356  f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S));
357  f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S));
358  f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
359 
361  M = polar_grid(T, S, x, w);
362  printf("M=%d.\n", M);
363 
365  fp1 = fopen("input_data_r.dat", "r");
366  fp2 = fopen("input_data_i.dat", "r");
367  if (fp1 == NULL)
368  return (-1);
369  for (k = 0; k < N * N; k++)
370  {
371  fscanf(fp1, NFFT__FR__ " ", &temp1);
372  fscanf(fp2, NFFT__FR__ " ", &temp2);
373  f_hat[k] = temp1 + _Complex_I * temp2;
374  }
375  fclose(fp1);
376  fclose(fp2);
377 
379  polar_dft(f_hat, N, f_direct, T, S, m);
380  // polar_fft(f_hat,N,f_direct,T,R,12);
381 
383  printf("\nTest of the polar FFT: \n");
384  fp1 = fopen("polar_fft_error.dat", "w+");
385  for (m = 1; m <= 12; m++)
386  {
388  polar_fft(f_hat, N, f, T, S, m);
389 
391  E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
392  printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
393  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
394  }
395  fclose(fp1);
396 
398  for (m = 3; m <= 9; m += 3)
399  {
400  printf("\nTest of the inverse polar FFT for m=%d: \n", m);
401  sprintf(filename, "polar_ifft_error%d.dat", m);
402  fp1 = fopen(filename, "w+");
403  for (max_i = 0; max_i <= 100; max_i += 10)
404  {
406  inverse_polar_fft(f_direct, T, S, f_tilde, N, max_i, m);
407 
409  /* E_max=0.0;
410  for(k=0;k<N*N;k++)
411  {
412  temp = cabs((f_hat[k]-f_tilde[k])/f_hat[k]);
413  if (temp>E_max) E_max=temp;
414  }
415  */
416  E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
417  printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
418  fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
419  }
420  fclose(fp1);
421  }
422 
424  NFFT(free)(x);
425  NFFT(free)(w);
426  NFFT(free)(f_hat);
427  NFFT(free)(f);
428  NFFT(free)(f_direct);
429  NFFT(free)(f_tilde);
430 
431  return EXIT_SUCCESS;
432 }
433 /* \} */
static int max_i(int a, int b)
max
Definition: fastsum.c:46
static int polar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
discrete polar FFT
static int polar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
NFFT-based polar FFT.
static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
Generates the points with weights for the polar grid with angles and offsets.
int main(int argc, char **argv)
test program for various parameters
static int inverse_polar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i, int m)
inverse NFFT-based polar FFT
#define MALLOC_F_HAT
Definition: nfft3.h:188
#define MALLOC_X
Definition: nfft3.h:187
#define PRE_FULL_PSI
Definition: nfft3.h:186
#define FFT_OUT_OF_PLACE
Definition: nfft3.h:190
#define PRE_PSI
Definition: nfft3.h:185
#define MALLOC_F
Definition: nfft3.h:189
#define PRE_LIN_PSI
Definition: nfft3.h:183
#define FFTW_INIT
Definition: nfft3.h:191
#define PRE_PHI_HUT
Definition: nfft3.h:181
#define CGNR
Definition: nfft3.h:775
#define PRECOMPUTE_DAMP
Definition: nfft3.h:779
#define PRECOMPUTE_WEIGHT
Definition: nfft3.h:778