00001
00010 #include <math.h>
00011 #include <stdlib.h>
00012 #include "util.h"
00013 #include "nfft3.h"
00014
00021 double GLOBAL_elapsed_time;
00022
00036 int mpolar_grid(int T, int R, double *x, double *w)
00037 {
00038 int t, r;
00039 double W;
00040 int R2=2*ceil(sqrt(2)*R/2);
00041 double xx, yy;
00042 int M=0;
00043
00044 for(t=-T/2; t<T/2; t++)
00045 {
00046 for(r=-R2/2; r<R2/2; r++)
00047 {
00048 xx = (double)r/R*cos(PI*t/T);
00049 yy = (double)r/R*sin(PI*t/T);
00050
00051 if ( ((-0.5-1.0/(double)R)<=xx) & (xx<=(0.5+1.0/(double)R)) &
00052 ((-0.5-1.0/(double)R)<=yy) & (yy<=(0.5+1.0/(double)R)) )
00053 {
00054 x[2*M+0] = xx;
00055 x[2*M+1] = yy;
00056
00057 if (r==0)
00058 w[M] = 1.0/4.0;
00059 else
00060 w[M] = fabs((double)r);
00061
00062 M++;
00063 }
00064 }
00065 }
00066
00068 W=0.0;
00069 for (t=0; t<M; t++)
00070 W+=w[t];
00071
00072 for (t=0; t<M; t++)
00073 w[t]/=W;
00074
00075 return M;
00076 }
00077
00079 int mpolar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00080 {
00081 int j,k;
00082 nfft_plan my_nfft_plan;
00084 double *x, *w;
00086 int N[2],n[2];
00087 int M;
00089 N[0]=NN; n[0]=2*N[0];
00090 N[1]=NN; n[1]=2*N[1];
00092 x = (double *)malloc(2*1.25*T*R*(sizeof(double)));
00093 if (x==NULL)
00094 return -1;
00095
00096 w = (double *)malloc(1.25*T*R*(sizeof(double)));
00097 if (w==NULL)
00098 return -1;
00099
00101 M=mpolar_grid(T,R,x,w);
00102 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00103 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00104 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00105
00107 for(j=0;j<my_nfft_plan.M_total;j++)
00108 {
00109 my_nfft_plan.x[2*j+0] = x[2*j+0];
00110 my_nfft_plan.x[2*j+1] = x[2*j+1];
00111 }
00112
00114 for(k=0;k<my_nfft_plan.N_total;k++)
00115 my_nfft_plan.f_hat[k] = f_hat[k];
00116
00117 GLOBAL_elapsed_time=nfft_second();
00118
00120 ndft_trafo(&my_nfft_plan);
00121
00122 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00123
00125 for(j=0;j<my_nfft_plan.M_total;j++)
00126 f[j] = my_nfft_plan.f[j];
00127
00129 nfft_finalize(&my_nfft_plan);
00130 free(x);
00131 free(w);
00132
00133 return EXIT_SUCCESS;
00134 }
00135
00137 int mpolar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00138 {
00139 int j,k;
00140 nfft_plan my_nfft_plan;
00142 double *x, *w;
00144 int N[2],n[2];
00145 int M;
00147 N[0]=NN; n[0]=2*N[0];
00148 N[1]=NN; n[1]=2*N[1];
00150 x = (double *)malloc(2*1.25*T*R*(sizeof(double)));
00151 if (x==NULL)
00152 return -1;
00153
00154 w = (double *)malloc(1.25*T*R*(sizeof(double)));
00155 if (w==NULL)
00156 return -1;
00157
00159 M=mpolar_grid(T,R,x,w);
00160 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00161 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00162 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00163
00166 for(j=0;j<my_nfft_plan.M_total;j++)
00167 {
00168 my_nfft_plan.x[2*j+0] = x[2*j+0];
00169 my_nfft_plan.x[2*j+1] = x[2*j+1];
00170 }
00171
00173 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00174 nfft_precompute_lin_psi(&my_nfft_plan);
00175
00176 if(my_nfft_plan.nfft_flags & PRE_PSI)
00177 nfft_precompute_psi(&my_nfft_plan);
00178
00179 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00180 nfft_precompute_full_psi(&my_nfft_plan);
00181
00183 for(k=0;k<my_nfft_plan.N_total;k++)
00184 my_nfft_plan.f_hat[k] = f_hat[k];
00185
00186 GLOBAL_elapsed_time=nfft_second();
00187
00189 nfft_trafo(&my_nfft_plan);
00190
00191 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00192
00194 for(j=0;j<my_nfft_plan.M_total;j++)
00195 f[j] = my_nfft_plan.f[j];
00196
00198 nfft_finalize(&my_nfft_plan);
00199 free(x);
00200 free(w);
00201
00202 return EXIT_SUCCESS;
00203 }
00204
00206 int inverse_mpolar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00207 {
00208 int j,k;
00209 nfft_plan my_nfft_plan;
00210 infft_plan my_infft_plan;
00212 double *x, *w;
00213 int l;
00215 int N[2],n[2];
00216 int M;
00218 N[0]=NN; n[0]=2*N[0];
00219 N[1]=NN; n[1]=2*N[1];
00221 x = (double *)malloc(2*1.25*T*R*(sizeof(double)));
00222 if (x==NULL)
00223 return -1;
00224
00225 w = (double *)malloc(1.25*T*R*(sizeof(double)));
00226 if (w==NULL)
00227 return -1;
00228
00230 M=mpolar_grid(T,R,x,w);
00231 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00232 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00233 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00234
00236 infft_init_advanced(&my_infft_plan,&my_nfft_plan, CGNR | PRECOMPUTE_WEIGHT );
00237
00239 for(j=0;j<my_nfft_plan.M_total;j++)
00240 {
00241 my_nfft_plan.x[2*j+0] = x[2*j+0];
00242 my_nfft_plan.x[2*j+1] = x[2*j+1];
00243 my_infft_plan.y[j] = f[j];
00244 my_infft_plan.w[j] = w[j];
00245 }
00246
00248 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00249 nfft_precompute_lin_psi(&my_nfft_plan);
00250
00251 if(my_nfft_plan.nfft_flags & PRE_PSI)
00252 nfft_precompute_psi(&my_nfft_plan);
00253
00254 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00255 nfft_precompute_full_psi(&my_nfft_plan);
00256
00257
00259 if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00260 for(j=0;j<my_infft_plan.mv->N[0];j++)
00261 for(k=0;k<my_infft_plan.mv->N[1];k++)
00262 {
00263 my_infft_plan.w_hat[j*my_infft_plan.mv->N[1]+k]=
00264 (sqrt(pow(j-my_infft_plan.mv->N[0]/2,2)+pow(k-my_infft_plan.mv->N[1]/2,2))>(my_infft_plan.mv->N[0]/2)?0:1);
00265 }
00266
00268 for(k=0;k<my_nfft_plan.N_total;k++)
00269 my_infft_plan.f_hat_iter[k] = 0.0 + I*0.0;
00270
00271 GLOBAL_elapsed_time=nfft_second();
00272
00274 infft_before_loop(&my_infft_plan);
00275
00276 if (max_i<1)
00277 {
00278 l=1;
00279 for(k=0;k<my_nfft_plan.N_total;k++)
00280 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00281 }
00282 else
00283 {
00284 for(l=1;l<=max_i;l++)
00285 {
00286 infft_loop_one_step(&my_infft_plan);
00287 }
00288 }
00289
00290 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00291
00293 for(k=0;k<my_nfft_plan.N_total;k++)
00294 f_hat[k] = my_infft_plan.f_hat_iter[k];
00295
00297 infft_finalize(&my_infft_plan);
00298 nfft_finalize(&my_nfft_plan);
00299 free(x);
00300 free(w);
00301
00302 return EXIT_SUCCESS;
00303 }
00304
00306 int comparison_fft(FILE *fp, int N, int T, int R)
00307 {
00308 fftw_plan my_fftw_plan;
00309 fftw_complex *f_hat,*f;
00310 int m,k;
00311 double t_fft, t_dft_mpolar;
00312
00313 f_hat = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*N*N);
00314 f = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*(T*R/4)*5);
00315
00316 my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
00317
00318 for(k=0; k<N*N; k++)
00319 f_hat[k] = (((double)rand())/RAND_MAX) + I* (((double)rand())/RAND_MAX);
00320
00321 GLOBAL_elapsed_time=nfft_second();
00322 for(m=0;m<65536/N;m++)
00323 {
00324 fftw_execute(my_fftw_plan);
00325
00326 f_hat[2]=2*f_hat[0];
00327 }
00328 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00329 t_fft=N*GLOBAL_elapsed_time/65536;
00330
00331 if(N<256)
00332 {
00333 mpolar_dft(f_hat,N,f,T,R,m);
00334 t_dft_mpolar=GLOBAL_elapsed_time;
00335 }
00336
00337 for (m=3; m<=9; m+=3)
00338 {
00339 if((m==3)&&(N<256))
00340 fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_mpolar,m);
00341 else
00342 if(m==3)
00343 fprintf(fp,"%d\t&\t&\t%1.1e&\t &\t%d\t",N,t_fft,m);
00344 else
00345 fprintf(fp," \t&\t&\t &\t &\t%d\t",m);
00346
00347 printf("N=%d\tt_fft=%1.1e\tt_dft_mpolar=%1.1e\tm=%d\t",N,t_fft,t_dft_mpolar,m);
00348
00349 mpolar_fft(f_hat,N,f,T,R,m);
00350 fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
00351 printf("t_mpolar=%1.1e\t",GLOBAL_elapsed_time);
00352 inverse_mpolar_fft(f,T,R,f_hat,N,2*m,m);
00353 if(m==9)
00354 fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
00355 else
00356 fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
00357 printf("t_impolar=%1.1e\n",GLOBAL_elapsed_time);
00358 }
00359
00360 fflush(fp);
00361
00362 fftw_free(f);
00363 fftw_free(f_hat);
00364
00365 return EXIT_SUCCESS;
00366 }
00367
00369 int main(int argc,char **argv)
00370 {
00371 int N;
00372 int T, R;
00373 int M;
00374 double *x, *w;
00375 fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00376 int k;
00377 int max_i;
00378 int m;
00379 double temp1, temp2, E_max=0.0;
00380 FILE *fp1, *fp2;
00381 char filename[30];
00382 int logN;
00383
00384 if( argc!=4 )
00385 {
00386 printf("mpolar_fft_test N T R \n");
00387 printf("\n");
00388 printf("N mpolar FFT of size NxN \n");
00389 printf("T number of slopes \n");
00390 printf("R number of offsets \n");
00391
00393 printf("\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
00394 fp1=fopen("mpolar_comparison_fft.dat","w");
00395 if (fp1==NULL)
00396 return(-1);
00397 for (logN=4; logN<=8; logN++)
00398 comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
00399 fclose(fp1);
00400
00401 exit(-1);
00402 }
00403
00404 N = atoi(argv[1]);
00405 T = atoi(argv[2]);
00406 R = atoi(argv[3]);
00407 printf("N=%d, modified polar grid with T=%d, R=%d => ",N,T,R);
00408
00409 x = (double *)malloc(2*1.25*T*R*(sizeof(double)));
00410 w = (double *)malloc(1.25*T*R*(sizeof(double)));
00411
00412 f_hat = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*N*N);
00413 f = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*1.25*T*R);
00414 f_direct = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*1.25*T*R);
00415 f_tilde = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*N*N);
00416
00418 M=mpolar_grid(T,R,x,w); printf("M=%d.\n",M);
00419
00421 fp1=fopen("input_data_r.dat","r");
00422 fp2=fopen("input_data_i.dat","r");
00423 if ((fp1==NULL) || (fp2==NULL))
00424 return(-1);
00425 for(k=0;k<N*N;k++)
00426 {
00427 fscanf(fp1,"%le ",&temp1);
00428 fscanf(fp2,"%le ",&temp2);
00429 f_hat[k]=temp1+I*temp2;
00430 }
00431 fclose(fp1);
00432 fclose(fp2);
00433
00435 mpolar_dft(f_hat,N,f_direct,T,R,m);
00436
00437
00439 printf("\nTest of the mpolar FFT: \n");
00440 fp1=fopen("mpolar_fft_error.dat","w+");
00441 for (m=1; m<=12; m++)
00442 {
00444 mpolar_fft(f_hat,N,f,T,R,m);
00445
00447 E_max=nfft_error_l_infty_complex(f_direct,f,M);
00448 printf("m=%2d: E_max = %e\n",m,E_max);
00449 fprintf(fp1,"%e\n",E_max);
00450 }
00451 fclose(fp1);
00452
00454 for (m=3; m<=9; m+=3)
00455 {
00456 printf("\nTest of the inverse mpolar FFT for m=%d: \n",m);
00457 sprintf(filename,"mpolar_ifft_error%d.dat",m);
00458 fp1=fopen(filename,"w+");
00459 for (max_i=0; max_i<=20; max_i+=2)
00460 {
00462 inverse_mpolar_fft(f_direct,T,R,f_tilde,N,max_i,m);
00463
00465 E_max=nfft_error_l_infty_complex(f_hat,f_tilde,N*N);
00466 printf("%3d iterations: E_max = %e\n",max_i,E_max);
00467 fprintf(fp1,"%e\n",E_max);
00468 }
00469 fclose(fp1);
00470 }
00471
00473 free(x);
00474 free(w);
00475 free(f_hat);
00476 free(f);
00477 free(f_direct);
00478 free(f_tilde);
00479
00480 return 0;
00481 }
00482
