00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00030 #include <math.h>
00031 #include <stdlib.h>
00032 #include <complex.h>
00033
00034 #include "nfft3util.h"
00035 #include "nfft3.h"
00036
00043 double GLOBAL_elapsed_time;
00044
00058 static int mpolar_grid(int T, int R, double *x, double *w)
00059 {
00060 int t, r;
00061 double W;
00062 int R2=2*ceil(sqrt(2.0)*R/2);
00063 double xx, yy;
00064 int M=0;
00065
00066 for(t=-T/2; t<T/2; t++)
00067 {
00068 for(r=-R2/2; r<R2/2; r++)
00069 {
00070 xx = (double)r/R*cos(PI*t/T);
00071 yy = (double)r/R*sin(PI*t/T);
00072
00073 if ( ((-0.5-1.0/(double)R)<=xx) & (xx<=(0.5+1.0/(double)R)) &
00074 ((-0.5-1.0/(double)R)<=yy) & (yy<=(0.5+1.0/(double)R)) )
00075 {
00076 x[2*M+0] = xx;
00077 x[2*M+1] = yy;
00078
00079 if (r==0)
00080 w[M] = 1.0/4.0;
00081 else
00082 w[M] = fabs((double)r);
00083
00084 M++;
00085 }
00086 }
00087 }
00088
00090 W=0.0;
00091 for (t=0; t<M; t++)
00092 W+=w[t];
00093
00094 for (t=0; t<M; t++)
00095 w[t]/=W;
00096
00097 return M;
00098 }
00099
00101 static int mpolar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00102 {
00103 int j,k;
00104 nfft_plan my_nfft_plan;
00106 double *x, *w;
00108 int N[2],n[2];
00109 int M;
00111 N[0]=NN; n[0]=2*N[0];
00112 N[1]=NN; n[1]=2*N[1];
00114 x = (double *)nfft_malloc(5*(T/2)*R*(sizeof(double)));
00115 if (x==NULL)
00116 return -1;
00117
00118 w = (double *)nfft_malloc(5*(T*R)/4*(sizeof(double)));
00119 if (w==NULL)
00120 return -1;
00121
00123 M=mpolar_grid(T,R,x,w);
00124 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00125 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00126 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00127
00129 for(j=0;j<my_nfft_plan.M_total;j++)
00130 {
00131 my_nfft_plan.x[2*j+0] = x[2*j+0];
00132 my_nfft_plan.x[2*j+1] = x[2*j+1];
00133 }
00134
00136 for(k=0;k<my_nfft_plan.N_total;k++)
00137 my_nfft_plan.f_hat[k] = f_hat[k];
00138
00139 GLOBAL_elapsed_time=nfft_second();
00140
00142 ndft_trafo(&my_nfft_plan);
00143
00144 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00145
00147 for(j=0;j<my_nfft_plan.M_total;j++)
00148 f[j] = my_nfft_plan.f[j];
00149
00151 nfft_finalize(&my_nfft_plan);
00152 nfft_free(x);
00153 nfft_free(w);
00154
00155 return EXIT_SUCCESS;
00156 }
00157
00159 static int mpolar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00160 {
00161 int j,k;
00162 nfft_plan my_nfft_plan;
00164 double *x, *w;
00166 int N[2],n[2];
00167 int M;
00169 N[0]=NN; n[0]=2*N[0];
00170 N[1]=NN; n[1]=2*N[1];
00172 x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00173 if (x==NULL)
00174 return -1;
00175
00176 w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00177 if (w==NULL)
00178 return -1;
00179
00181 M=mpolar_grid(T,R,x,w);
00182 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00183 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00184 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00185
00188 for(j=0;j<my_nfft_plan.M_total;j++)
00189 {
00190 my_nfft_plan.x[2*j+0] = x[2*j+0];
00191 my_nfft_plan.x[2*j+1] = x[2*j+1];
00192 }
00193
00195 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00196 nfft_precompute_lin_psi(&my_nfft_plan);
00197
00198 if(my_nfft_plan.nfft_flags & PRE_PSI)
00199 nfft_precompute_psi(&my_nfft_plan);
00200
00201 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00202 nfft_precompute_full_psi(&my_nfft_plan);
00203
00205 for(k=0;k<my_nfft_plan.N_total;k++)
00206 my_nfft_plan.f_hat[k] = f_hat[k];
00207
00208 GLOBAL_elapsed_time=nfft_second();
00209
00211 nfft_trafo(&my_nfft_plan);
00212
00213 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00214
00216 for(j=0;j<my_nfft_plan.M_total;j++)
00217 f[j] = my_nfft_plan.f[j];
00218
00220 nfft_finalize(&my_nfft_plan);
00221 nfft_free(x);
00222 nfft_free(w);
00223
00224 return EXIT_SUCCESS;
00225 }
00226
00228 static int inverse_mpolar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00229 {
00230 int j,k;
00231 nfft_plan my_nfft_plan;
00232 solver_plan_complex my_infft_plan;
00234 double *x, *w;
00235 int l;
00237 int N[2],n[2];
00238 int M;
00240 N[0]=NN; n[0]=2*N[0];
00241 N[1]=NN; n[1]=2*N[1];
00243 x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00244 if (x==NULL)
00245 return -1;
00246
00247 w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00248 if (w==NULL)
00249 return -1;
00250
00252 M=mpolar_grid(T,R,x,w);
00253 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00254 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00255 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00256
00258 solver_init_advanced_complex(&my_infft_plan,(mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
00259
00261 for(j=0;j<my_nfft_plan.M_total;j++)
00262 {
00263 my_nfft_plan.x[2*j+0] = x[2*j+0];
00264 my_nfft_plan.x[2*j+1] = x[2*j+1];
00265 my_infft_plan.y[j] = f[j];
00266 my_infft_plan.w[j] = w[j];
00267 }
00268
00270 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00271 nfft_precompute_lin_psi(&my_nfft_plan);
00272
00273 if(my_nfft_plan.nfft_flags & PRE_PSI)
00274 nfft_precompute_psi(&my_nfft_plan);
00275
00276 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00277 nfft_precompute_full_psi(&my_nfft_plan);
00278
00279
00281 if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00282 for(j=0;j<my_nfft_plan.N[0];j++)
00283 for(k=0;k<my_nfft_plan.N[1];k++)
00284 {
00285 my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
00286 (sqrt(pow(j-my_nfft_plan.N[0]/2,2)+pow(k-my_nfft_plan.N[1]/2,2))>(my_nfft_plan.N[0]/2)?0:1);
00287 }
00288
00290 for(k=0;k<my_nfft_plan.N_total;k++)
00291 my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00292
00293 GLOBAL_elapsed_time=nfft_second();
00294
00296 solver_before_loop_complex(&my_infft_plan);
00297
00298 if (max_i<1)
00299 {
00300 l=1;
00301 for(k=0;k<my_nfft_plan.N_total;k++)
00302 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00303 }
00304 else
00305 {
00306 for(l=1;l<=max_i;l++)
00307 {
00308 solver_loop_one_step_complex(&my_infft_plan);
00309 }
00310 }
00311
00312 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00313
00315 for(k=0;k<my_nfft_plan.N_total;k++)
00316 f_hat[k] = my_infft_plan.f_hat_iter[k];
00317
00319 solver_finalize_complex(&my_infft_plan);
00320 nfft_finalize(&my_nfft_plan);
00321 nfft_free(x);
00322 nfft_free(w);
00323
00324 return EXIT_SUCCESS;
00325 }
00326
00328 static int comparison_fft(FILE *fp, int N, int T, int R)
00329 {
00330 fftw_plan my_fftw_plan;
00331 fftw_complex *f_hat,*f;
00332 int m,k;
00333 double t_fft, t_dft_mpolar;
00334
00335 f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00336 f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*(T*R/4)*5);
00337
00338 my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
00339
00340 for(k=0; k<N*N; k++)
00341 f_hat[k] = (((double)rand())/RAND_MAX) + _Complex_I* (((double)rand())/RAND_MAX);
00342
00343 GLOBAL_elapsed_time=nfft_second();
00344 for(m=0;m<65536/N;m++)
00345 {
00346 fftw_execute(my_fftw_plan);
00347
00348 f_hat[2]=2*f_hat[0];
00349 }
00350 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00351 t_fft=N*GLOBAL_elapsed_time/65536;
00352
00353 if(N<256)
00354 {
00355 mpolar_dft(f_hat,N,f,T,R,1);
00356 t_dft_mpolar=GLOBAL_elapsed_time;
00357 }
00358
00359 for (m=3; m<=9; m+=3)
00360 {
00361 if((m==3)&&(N<256))
00362 fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_mpolar,m);
00363 else
00364 if(m==3)
00365 fprintf(fp,"%d\t&\t&\t%1.1e&\t &\t%d\t",N,t_fft,m);
00366 else
00367 fprintf(fp," \t&\t&\t &\t &\t%d\t",m);
00368
00369 printf("N=%d\tt_fft=%1.1e\tt_dft_mpolar=%1.1e\tm=%d\t",N,t_fft,t_dft_mpolar,m);
00370
00371 mpolar_fft(f_hat,N,f,T,R,m);
00372 fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
00373 printf("t_mpolar=%1.1e\t",GLOBAL_elapsed_time);
00374 inverse_mpolar_fft(f,T,R,f_hat,N,2*m,m);
00375 if(m==9)
00376 fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
00377 else
00378 fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
00379 printf("t_impolar=%1.1e\n",GLOBAL_elapsed_time);
00380 }
00381
00382 fflush(fp);
00383
00384 nfft_free(f);
00385 nfft_free(f_hat);
00386
00387 return EXIT_SUCCESS;
00388 }
00389
00391 int main(int argc,char **argv)
00392 {
00393 int N;
00394 int T, R;
00395 int M;
00396 double *x, *w;
00397 fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00398 int k;
00399 int max_i;
00400 int m;
00401 double temp1, temp2, E_max=0.0;
00402 FILE *fp1, *fp2;
00403 char filename[30];
00404 int logN;
00405
00406 if( argc!=4 )
00407 {
00408 printf("mpolar_fft_test N T R \n");
00409 printf("\n");
00410 printf("N mpolar FFT of size NxN \n");
00411 printf("T number of slopes \n");
00412 printf("R number of offsets \n");
00413
00415 printf("\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
00416 fp1=fopen("mpolar_comparison_fft.dat","w");
00417 if (fp1==NULL)
00418 return(-1);
00419 for (logN=4; logN<=8; logN++)
00420 comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
00421 fclose(fp1);
00422
00423 exit(-1);
00424 }
00425
00426 N = atoi(argv[1]);
00427 T = atoi(argv[2]);
00428 R = atoi(argv[3]);
00429 printf("N=%d, modified polar grid with T=%d, R=%d => ",N,T,R);
00430
00431 x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00432 w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00433
00434 f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00435 f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*1.25*T*R);
00436 f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*1.25*T*R);
00437 f_tilde = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00438
00440 M=mpolar_grid(T,R,x,w); printf("M=%d.\n",M);
00441
00443 fp1=fopen("input_data_r.dat","r");
00444 fp2=fopen("input_data_i.dat","r");
00445 if ((fp1==NULL) || (fp2==NULL))
00446 return(-1);
00447 for(k=0;k<N*N;k++)
00448 {
00449 fscanf(fp1,"%le ",&temp1);
00450 fscanf(fp2,"%le ",&temp2);
00451 f_hat[k]=temp1+ _Complex_I*temp2;
00452 }
00453 fclose(fp1);
00454 fclose(fp2);
00455
00457 mpolar_dft(f_hat,N,f_direct,T,R,1);
00458
00459
00461 printf("\nTest of the mpolar FFT: \n");
00462 fp1=fopen("mpolar_fft_error.dat","w+");
00463 for (m=1; m<=12; m++)
00464 {
00466 mpolar_fft(f_hat,N,f,T,R,m);
00467
00469 E_max=nfft_error_l_infty_complex(f_direct,f,M);
00470 printf("m=%2d: E_max = %e\n",m,E_max);
00471 fprintf(fp1,"%e\n",E_max);
00472 }
00473 fclose(fp1);
00474
00476 for (m=3; m<=9; m+=3)
00477 {
00478 printf("\nTest of the inverse mpolar FFT for m=%d: \n",m);
00479 sprintf(filename,"mpolar_ifft_error%d.dat",m);
00480 fp1=fopen(filename,"w+");
00481 for (max_i=0; max_i<=20; max_i+=2)
00482 {
00484 inverse_mpolar_fft(f_direct,T,R,f_tilde,N,max_i,m);
00485
00487 E_max=nfft_error_l_infty_complex(f_hat,f_tilde,N*N);
00488 printf("%3d iterations: E_max = %e\n",max_i,E_max);
00489 fprintf(fp1,"%e\n",E_max);
00490 }
00491 fclose(fp1);
00492 }
00493
00495 nfft_free(x);
00496 nfft_free(w);
00497 nfft_free(f_hat);
00498 nfft_free(f);
00499 nfft_free(f_direct);
00500 nfft_free(f_tilde);
00501
00502 return 0;
00503 }
00504