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