NFFT Logo 3.2.3
polar_fft_test.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2012 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: polar_fft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
20 
29 #include "config.h"
30 
31 #include <math.h>
32 #include <stdlib.h>
33 #ifdef HAVE_COMPLEX_H
34 #include <complex.h>
35 #endif
36 
37 #include "nfft3util.h"
38 #include "nfft3.h"
39 #include "infft.h"
40 
78 static int polar_grid(int T, int R, double *x, double *w)
79 {
80  int t, r;
81  double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
82 
83  for(t=-T/2; t<T/2; t++)
84  {
85  for(r=-R/2; r<R/2; r++)
86  {
87  x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
88  x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
89  if (r==0)
90  w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
91  else
92  w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
93  }
94  }
95 
96  return T*R;
97 }
98 
100 static int polar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
101 {
102  int j,k;
103  nfft_plan my_nfft_plan;
105  double *x, *w;
107  int N[2],n[2];
108  int M=T*R;
110  N[0]=NN; n[0]=2*N[0];
111  N[1]=NN; n[1]=2*N[1];
113  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
114  if (x==NULL)
115  return -1;
116 
117  w = (double *)nfft_malloc(T*R*(sizeof(double)));
118  if (w==NULL)
119  return -1;
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 | FFT_OUT_OF_PLACE,
124  FFTW_MEASURE| FFTW_DESTROY_INPUT);
125 
127  polar_grid(T,R,x,w);
128  for(j=0;j<my_nfft_plan.M_total;j++)
129  {
130  my_nfft_plan.x[2*j+0] = x[2*j+0];
131  my_nfft_plan.x[2*j+1] = x[2*j+1];
132  }
133 
135  for(k=0;k<my_nfft_plan.N_total;k++)
136  my_nfft_plan.f_hat[k] = f_hat[k];
137 
139  nfft_trafo_direct(&my_nfft_plan);
140 
142  for(j=0;j<my_nfft_plan.M_total;j++)
143  f[j] = my_nfft_plan.f[j];
144 
146  nfft_finalize(&my_nfft_plan);
147  nfft_free(x);
148  nfft_free(w);
149 
150  return EXIT_SUCCESS;
151 }
152 
154 static int polar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
155 {
156  int j,k;
157  nfft_plan my_nfft_plan;
159  double *x, *w;
161  int N[2],n[2];
162  int M=T*R;
164  N[0]=NN; n[0]=2*N[0];
165  N[1]=NN; n[1]=2*N[1];
167  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
168  if (x==NULL)
169  return -1;
170 
171  w = (double *)nfft_malloc(T*R*(sizeof(double)));
172  if (w==NULL)
173  return -1;
174 
176  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
177  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
178  FFTW_MEASURE| FFTW_DESTROY_INPUT);
179 
181  polar_grid(T,R,x,w);
182  for(j=0;j<my_nfft_plan.M_total;j++)
183  {
184  my_nfft_plan.x[2*j+0] = x[2*j+0];
185  my_nfft_plan.x[2*j+1] = x[2*j+1];
186  }
187 
189  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
190  nfft_precompute_lin_psi(&my_nfft_plan);
191 
192  if(my_nfft_plan.nfft_flags & PRE_PSI)
193  nfft_precompute_psi(&my_nfft_plan);
194 
195  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
196  nfft_precompute_full_psi(&my_nfft_plan);
197 
199  for(k=0;k<my_nfft_plan.N_total;k++)
200  my_nfft_plan.f_hat[k] = f_hat[k];
201 
203  nfft_trafo(&my_nfft_plan);
204 
206  for(j=0;j<my_nfft_plan.M_total;j++)
207  f[j] = my_nfft_plan.f[j];
208 
210  nfft_finalize(&my_nfft_plan);
211  nfft_free(x);
212  nfft_free(w);
213 
214  return EXIT_SUCCESS;
215 }
216 
218 static int inverse_polar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
219 {
220  int j,k;
221  nfft_plan my_nfft_plan;
222  solver_plan_complex my_infft_plan;
224  double *x, *w;
225  int l;
227  int N[2],n[2];
228  int M=T*R;
230  N[0]=NN; n[0]=2*N[0];
231  N[1]=NN; n[1]=2*N[1];
233  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
234  if (x==NULL)
235  return -1;
236 
237  w = (double *)nfft_malloc(T*R*(sizeof(double)));
238  if (w==NULL)
239  return -1;
240 
242  nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
243  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
244  FFTW_MEASURE| FFTW_DESTROY_INPUT);
245 
247  solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
248 
250  polar_grid(T,R,x,w);
251  for(j=0;j<my_nfft_plan.M_total;j++)
252  {
253  my_nfft_plan.x[2*j+0] = x[2*j+0];
254  my_nfft_plan.x[2*j+1] = x[2*j+1];
255  my_infft_plan.y[j] = f[j];
256  my_infft_plan.w[j] = w[j];
257  }
258 
260  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
261  nfft_precompute_lin_psi(&my_nfft_plan);
262 
263  if(my_nfft_plan.nfft_flags & PRE_PSI)
264  nfft_precompute_psi(&my_nfft_plan);
265 
266  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
267  nfft_precompute_full_psi(&my_nfft_plan);
268 
270  if(my_infft_plan.flags & PRECOMPUTE_DAMP)
271  for(j=0;j<my_nfft_plan.N[0];j++)
272  for(k=0;k<my_nfft_plan.N[1];k++)
273  {
274  my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
275  (sqrt(pow((double)(j-my_nfft_plan.N[0]/2),2.0)+pow((double)(k-my_nfft_plan.N[1]/2),2.0))
276  >((double)(my_nfft_plan.N[0]/2))?0:1);
277  }
278 
280  for(k=0;k<my_nfft_plan.N_total;k++)
281  my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
282 
284  solver_before_loop_complex(&my_infft_plan);
285 
286  if (max_i<1)
287  {
288  l=1;
289  for(k=0;k<my_nfft_plan.N_total;k++)
290  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
291  }
292  else
293  {
294  for(l=1;l<=max_i;l++)
295  {
296  solver_loop_one_step_complex(&my_infft_plan);
297  }
298  }
299 
301  for(k=0;k<my_nfft_plan.N_total;k++)
302  f_hat[k] = my_infft_plan.f_hat_iter[k];
303 
305  solver_finalize_complex(&my_infft_plan);
306  nfft_finalize(&my_nfft_plan);
307  nfft_free(x);
308  nfft_free(w);
309 
310  return EXIT_SUCCESS;
311 }
312 
314 int main(int argc,char **argv)
315 {
316  int N;
317  int T, R;
318  int M;
319  double *x, *w;
320  fftw_complex *f_hat, *f, *f_direct, *f_tilde;
321  int k;
322  int max_i;
323  int m = 1;
324  double temp1, temp2, E_max=0.0;
325  FILE *fp1, *fp2;
326  char filename[30];
327 
328  if( argc!=4 )
329  {
330  printf("polar_fft_test N T R \n");
331  printf("\n");
332  printf("N polar FFT of size NxN \n");
333  printf("T number of slopes \n");
334  printf("R number of offsets \n");
335  exit(-1);
336  }
337 
338  N = atoi(argv[1]);
339  T = atoi(argv[2]);
340  R = atoi(argv[3]);
341  printf("N=%d, polar grid with T=%d, R=%d => ",N,T,R);
342 
343  x = (double *)nfft_malloc(2*5*(T/2)*(R/2)*(sizeof(double)));
344  w = (double *)nfft_malloc(5*(T/2)*(R/2)*(sizeof(double)));
345 
346  f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
347  f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*T*R);
348  f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*T*R);
349  f_tilde = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
350 
352  M=polar_grid(T,R,x,w); printf("M=%d.\n",M);
353 
355  fp1=fopen("input_data_r.dat","r");
356  fp2=fopen("input_data_i.dat","r");
357  if (fp1==NULL)
358  return(-1);
359  for(k=0;k<N*N;k++)
360  {
361  fscanf(fp1,"%le ",&temp1);
362  fscanf(fp2,"%le ",&temp2);
363  f_hat[k]=temp1+ _Complex_I*temp2;
364  }
365  fclose(fp1);
366  fclose(fp2);
367 
369  polar_dft(f_hat,N,f_direct,T,R,m);
370  // polar_fft(f_hat,N,f_direct,T,R,12);
371 
373  printf("\nTest of the polar FFT: \n");
374  fp1=fopen("polar_fft_error.dat","w+");
375  for (m=1; m<=12; m++)
376  {
378  polar_fft(f_hat,N,f,T,R,m);
379 
381  E_max=X(error_l_infty_complex)(f_direct,f,M);
382  printf("m=%2d: E_max = %e\n",m,E_max);
383  fprintf(fp1,"%e\n",E_max);
384  }
385  fclose(fp1);
386 
388  for (m=3; m<=9; m+=3)
389  {
390  printf("\nTest of the inverse polar FFT for m=%d: \n",m);
391  sprintf(filename,"polar_ifft_error%d.dat",m);
392  fp1=fopen(filename,"w+");
393  for (max_i=0; max_i<=100; max_i+=10)
394  {
396  inverse_polar_fft(f_direct,T,R,f_tilde,N,max_i,m);
397 
399  /* E_max=0.0;
400  for(k=0;k<N*N;k++)
401  {
402  temp = cabs((f_hat[k]-f_tilde[k])/f_hat[k]);
403  if (temp>E_max) E_max=temp;
404  }
405  */
406  E_max=X(error_l_infty_complex)(f_hat,f_tilde,N*N);
407  printf("%3d iterations: E_max = %e\n",max_i,E_max);
408  fprintf(fp1,"%e\n",E_max);
409  }
410  fclose(fp1);
411  }
412 
414  nfft_free(x);
415  nfft_free(w);
416  nfft_free(f_hat);
417  nfft_free(f);
418  nfft_free(f_direct);
419  nfft_free(f_tilde);
420 
421  return 0;
422 }
423 /* \} */
Generated on Tue Apr 30 2013 by Doxygen 1.8.1