NFFT Logo 3.2.2
polar_fft_test.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* $Id: polar_fft_test.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00029 #include "config.h"
00030 
00031 #include <math.h>
00032 #include <stdlib.h>
00033 #ifdef HAVE_COMPLEX_H
00034 #include <complex.h>
00035 #endif
00036 
00037 #include "nfft3util.h"
00038 #include "nfft3.h"
00039 #include "infft.h"
00040 
00078 static int polar_grid(int T, int R, double *x, double *w)
00079 {
00080   int t, r;
00081   double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00082 
00083   for(t=-T/2; t<T/2; t++)
00084   {
00085     for(r=-R/2; r<R/2; r++)
00086     {
00087       x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00088       x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00089       if (r==0)
00090         w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00091       else
00092         w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00093     }
00094   }
00095 
00096   return T*R;                           
00097 }
00098 
00100 static int polar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00101 {
00102   int j,k;                              
00103   nfft_plan my_nfft_plan;               
00105   double *x, *w;                        
00107   int N[2],n[2];
00108   int M=T*R;                            
00110   N[0]=NN; n[0]=2*N[0];                 
00111   N[1]=NN; n[1]=2*N[1];                 
00113   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00114   if (x==NULL)
00115     return -1;
00116 
00117   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00118   if (w==NULL)
00119     return -1;
00120 
00122   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00123                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00124                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00125 
00127   polar_grid(T,R,x,w);
00128   for(j=0;j<my_nfft_plan.M_total;j++)
00129   {
00130     my_nfft_plan.x[2*j+0] = x[2*j+0];
00131     my_nfft_plan.x[2*j+1] = x[2*j+1];
00132   }
00133 
00135   for(k=0;k<my_nfft_plan.N_total;k++)
00136     my_nfft_plan.f_hat[k] = f_hat[k];
00137 
00139   nfft_trafo_direct(&my_nfft_plan);
00140 
00142   for(j=0;j<my_nfft_plan.M_total;j++)
00143     f[j] = my_nfft_plan.f[j];
00144 
00146   nfft_finalize(&my_nfft_plan);
00147   nfft_free(x);
00148   nfft_free(w);
00149 
00150   return EXIT_SUCCESS;
00151 }
00152 
00154 static int polar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00155 {
00156   int j,k;                              
00157   nfft_plan my_nfft_plan;               
00159   double *x, *w;                        
00161   int N[2],n[2];
00162   int M=T*R;                            
00164   N[0]=NN; n[0]=2*N[0];                 
00165   N[1]=NN; n[1]=2*N[1];                 
00167   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00168   if (x==NULL)
00169     return -1;
00170 
00171   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00172   if (w==NULL)
00173     return -1;
00174 
00176   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00177                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00178                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00179 
00181   polar_grid(T,R,x,w);
00182   for(j=0;j<my_nfft_plan.M_total;j++)
00183   {
00184     my_nfft_plan.x[2*j+0] = x[2*j+0];
00185     my_nfft_plan.x[2*j+1] = x[2*j+1];
00186   }
00187 
00189   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00190     nfft_precompute_lin_psi(&my_nfft_plan);
00191 
00192   if(my_nfft_plan.nfft_flags & PRE_PSI)
00193     nfft_precompute_psi(&my_nfft_plan);
00194 
00195   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00196     nfft_precompute_full_psi(&my_nfft_plan);
00197 
00199   for(k=0;k<my_nfft_plan.N_total;k++)
00200     my_nfft_plan.f_hat[k] = f_hat[k];
00201 
00203   nfft_trafo(&my_nfft_plan);
00204 
00206   for(j=0;j<my_nfft_plan.M_total;j++)
00207     f[j] = my_nfft_plan.f[j];
00208 
00210   nfft_finalize(&my_nfft_plan);
00211   nfft_free(x);
00212   nfft_free(w);
00213 
00214   return EXIT_SUCCESS;
00215 }
00216 
00218 static int inverse_polar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00219 {
00220   int j,k;                              
00221   nfft_plan my_nfft_plan;               
00222   solver_plan_complex my_infft_plan;             
00224   double *x, *w;                        
00225   int l;                                
00227   int N[2],n[2];
00228   int M=T*R;                            
00230   N[0]=NN; n[0]=2*N[0];                 
00231   N[1]=NN; n[1]=2*N[1];                 
00233   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00234   if (x==NULL)
00235     return -1;
00236 
00237   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00238   if (w==NULL)
00239     return -1;
00240 
00242   nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00243                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00244                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00245 
00247   solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
00248 
00250   polar_grid(T,R,x,w);
00251   for(j=0;j<my_nfft_plan.M_total;j++)
00252   {
00253     my_nfft_plan.x[2*j+0] = x[2*j+0];
00254     my_nfft_plan.x[2*j+1] = x[2*j+1];
00255     my_infft_plan.y[j]    = f[j];
00256     my_infft_plan.w[j]    = w[j];
00257   }
00258 
00260   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00261     nfft_precompute_lin_psi(&my_nfft_plan);
00262 
00263   if(my_nfft_plan.nfft_flags & PRE_PSI)
00264     nfft_precompute_psi(&my_nfft_plan);
00265 
00266   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00267     nfft_precompute_full_psi(&my_nfft_plan);
00268 
00270   if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00271     for(j=0;j<my_nfft_plan.N[0];j++)
00272       for(k=0;k<my_nfft_plan.N[1];k++)
00273   {
00274     my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
00275         (sqrt(pow((double)(j-my_nfft_plan.N[0]/2),2.0)+pow((double)(k-my_nfft_plan.N[1]/2),2.0))
00276             >((double)(my_nfft_plan.N[0]/2))?0:1);
00277   }
00278 
00280   for(k=0;k<my_nfft_plan.N_total;k++)
00281     my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00282 
00284   solver_before_loop_complex(&my_infft_plan);
00285 
00286   if (max_i<1)
00287   {
00288     l=1;
00289     for(k=0;k<my_nfft_plan.N_total;k++)
00290       my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00291   }
00292   else
00293   {
00294     for(l=1;l<=max_i;l++)
00295     {
00296       solver_loop_one_step_complex(&my_infft_plan);
00297     }
00298   }
00299 
00301   for(k=0;k<my_nfft_plan.N_total;k++)
00302     f_hat[k] = my_infft_plan.f_hat_iter[k];
00303 
00305   solver_finalize_complex(&my_infft_plan);
00306   nfft_finalize(&my_nfft_plan);
00307   nfft_free(x);
00308   nfft_free(w);
00309 
00310   return EXIT_SUCCESS;
00311 }
00312 
00314 int main(int argc,char **argv)
00315 {
00316   int N;                                
00317   int T, R;                             
00318   int M;                                
00319   double *x, *w;                        
00320   fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00321   int k;
00322   int max_i;                            
00323   int m = 1;
00324   double temp1, temp2, E_max=0.0;
00325   FILE *fp1, *fp2;
00326   char filename[30];
00327 
00328   if( argc!=4 )
00329   {
00330     printf("polar_fft_test N T R \n");
00331     printf("\n");
00332     printf("N          polar FFT of size NxN     \n");
00333     printf("T          number of slopes          \n");
00334     printf("R          number of offsets         \n");
00335     exit(-1);
00336   }
00337 
00338   N = atoi(argv[1]);
00339   T = atoi(argv[2]);
00340   R = atoi(argv[3]);
00341   printf("N=%d, polar grid with T=%d, R=%d => ",N,T,R);
00342 
00343   x = (double *)nfft_malloc(2*5*(T/2)*(R/2)*(sizeof(double)));
00344   w = (double *)nfft_malloc(5*(T/2)*(R/2)*(sizeof(double)));
00345 
00346   f_hat    = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00347   f        = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*T*R);
00348   f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*T*R);
00349   f_tilde  = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00350 
00352   M=polar_grid(T,R,x,w); printf("M=%d.\n",M);
00353 
00355   fp1=fopen("input_data_r.dat","r");
00356   fp2=fopen("input_data_i.dat","r");
00357   if (fp1==NULL)
00358     return(-1);
00359   for(k=0;k<N*N;k++)
00360   {
00361     fscanf(fp1,"%le ",&temp1);
00362     fscanf(fp2,"%le ",&temp2);
00363     f_hat[k]=temp1+ _Complex_I*temp2;
00364   }
00365   fclose(fp1);
00366   fclose(fp2);
00367 
00369     polar_dft(f_hat,N,f_direct,T,R,m);
00370   //  polar_fft(f_hat,N,f_direct,T,R,12);
00371 
00373   printf("\nTest of the polar FFT: \n");
00374   fp1=fopen("polar_fft_error.dat","w+");
00375   for (m=1; m<=12; m++)
00376   {
00378     polar_fft(f_hat,N,f,T,R,m);
00379 
00381     E_max=X(error_l_infty_complex)(f_direct,f,M);
00382     printf("m=%2d: E_max = %e\n",m,E_max);
00383     fprintf(fp1,"%e\n",E_max);
00384   }
00385   fclose(fp1);
00386 
00388   for (m=3; m<=9; m+=3)
00389   {
00390     printf("\nTest of the inverse polar FFT for m=%d: \n",m);
00391     sprintf(filename,"polar_ifft_error%d.dat",m);
00392     fp1=fopen(filename,"w+");
00393     for (max_i=0; max_i<=100; max_i+=10)
00394     {
00396       inverse_polar_fft(f_direct,T,R,f_tilde,N,max_i,m);
00397 
00399       /* E_max=0.0;
00400       for(k=0;k<N*N;k++)
00401       {
00402         temp = cabs((f_hat[k]-f_tilde[k])/f_hat[k]);
00403         if (temp>E_max) E_max=temp;
00404       }
00405       */
00406        E_max=X(error_l_infty_complex)(f_hat,f_tilde,N*N);
00407       printf("%3d iterations: E_max = %e\n",max_i,E_max);
00408       fprintf(fp1,"%e\n",E_max);
00409     }
00410     fclose(fp1);
00411   }
00412 
00414   nfft_free(x);
00415   nfft_free(w);
00416   nfft_free(f_hat);
00417   nfft_free(f);
00418   nfft_free(f_direct);
00419   nfft_free(f_tilde);
00420 
00421   return 0;
00422 }
00423 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409