00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00029 #include <math.h>
00030 #include <stdlib.h>
00031 #include <complex.h>
00032
00033 #include "nfft3util.h"
00034 #include "nfft3.h"
00035
00073 static int polar_grid(int T, int R, double *x, double *w)
00074 {
00075 int t, r;
00076 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00077
00078 for(t=-T/2; t<T/2; t++)
00079 {
00080 for(r=-R/2; r<R/2; r++)
00081 {
00082 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00083 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00084 if (r==0)
00085 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00086 else
00087 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00088 }
00089 }
00090
00091 return T*R;
00092 }
00093
00095 static int polar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00096 {
00097 int j,k;
00098 nfft_plan my_nfft_plan;
00100 double *x, *w;
00102 int N[2],n[2];
00103 int M=T*R;
00105 N[0]=NN; n[0]=2*N[0];
00106 N[1]=NN; n[1]=2*N[1];
00108 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00109 if (x==NULL)
00110 return -1;
00111
00112 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00113 if (w==NULL)
00114 return -1;
00115
00117 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00118 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00119 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00120
00122 polar_grid(T,R,x,w);
00123 for(j=0;j<my_nfft_plan.M_total;j++)
00124 {
00125 my_nfft_plan.x[2*j+0] = x[2*j+0];
00126 my_nfft_plan.x[2*j+1] = x[2*j+1];
00127 }
00128
00130 for(k=0;k<my_nfft_plan.N_total;k++)
00131 my_nfft_plan.f_hat[k] = f_hat[k];
00132
00134 ndft_trafo(&my_nfft_plan);
00135
00137 for(j=0;j<my_nfft_plan.M_total;j++)
00138 f[j] = my_nfft_plan.f[j];
00139
00141 nfft_finalize(&my_nfft_plan);
00142 nfft_free(x);
00143 nfft_free(w);
00144
00145 return EXIT_SUCCESS;
00146 }
00147
00149 static int polar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00150 {
00151 int j,k;
00152 nfft_plan my_nfft_plan;
00154 double *x, *w;
00156 int N[2],n[2];
00157 int M=T*R;
00159 N[0]=NN; n[0]=2*N[0];
00160 N[1]=NN; n[1]=2*N[1];
00162 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00163 if (x==NULL)
00164 return -1;
00165
00166 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00167 if (w==NULL)
00168 return -1;
00169
00171 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00172 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00173 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00174
00176 polar_grid(T,R,x,w);
00177 for(j=0;j<my_nfft_plan.M_total;j++)
00178 {
00179 my_nfft_plan.x[2*j+0] = x[2*j+0];
00180 my_nfft_plan.x[2*j+1] = x[2*j+1];
00181 }
00182
00184 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00185 nfft_precompute_lin_psi(&my_nfft_plan);
00186
00187 if(my_nfft_plan.nfft_flags & PRE_PSI)
00188 nfft_precompute_psi(&my_nfft_plan);
00189
00190 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00191 nfft_precompute_full_psi(&my_nfft_plan);
00192
00194 for(k=0;k<my_nfft_plan.N_total;k++)
00195 my_nfft_plan.f_hat[k] = f_hat[k];
00196
00198 nfft_trafo(&my_nfft_plan);
00199
00201 for(j=0;j<my_nfft_plan.M_total;j++)
00202 f[j] = my_nfft_plan.f[j];
00203
00205 nfft_finalize(&my_nfft_plan);
00206 nfft_free(x);
00207 nfft_free(w);
00208
00209 return EXIT_SUCCESS;
00210 }
00211
00213 static int inverse_polar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00214 {
00215 int j,k;
00216 nfft_plan my_nfft_plan;
00217 solver_plan_complex my_infft_plan;
00219 double *x, *w;
00220 int l;
00222 int N[2],n[2];
00223 int M=T*R;
00225 N[0]=NN; n[0]=2*N[0];
00226 N[1]=NN; n[1]=2*N[1];
00228 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00229 if (x==NULL)
00230 return -1;
00231
00232 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00233 if (w==NULL)
00234 return -1;
00235
00237 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00238 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00239 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00240
00242 solver_init_advanced_complex(&my_infft_plan,(mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
00243
00245 polar_grid(T,R,x,w);
00246 for(j=0;j<my_nfft_plan.M_total;j++)
00247 {
00248 my_nfft_plan.x[2*j+0] = x[2*j+0];
00249 my_nfft_plan.x[2*j+1] = x[2*j+1];
00250 my_infft_plan.y[j] = f[j];
00251 my_infft_plan.w[j] = w[j];
00252 }
00253
00255 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00256 nfft_precompute_lin_psi(&my_nfft_plan);
00257
00258 if(my_nfft_plan.nfft_flags & PRE_PSI)
00259 nfft_precompute_psi(&my_nfft_plan);
00260
00261 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00262 nfft_precompute_full_psi(&my_nfft_plan);
00263
00265 if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00266 for(j=0;j<my_nfft_plan.N[0];j++)
00267 for(k=0;k<my_nfft_plan.N[1];k++)
00268 {
00269 my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
00270 (sqrt(pow((double)(j-my_nfft_plan.N[0]/2),2.0)+pow((double)(k-my_nfft_plan.N[1]/2),2.0))
00271 >((double)(my_nfft_plan.N[0]/2))?0:1);
00272 }
00273
00275 for(k=0;k<my_nfft_plan.N_total;k++)
00276 my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00277
00279 solver_before_loop_complex(&my_infft_plan);
00280
00281 if (max_i<1)
00282 {
00283 l=1;
00284 for(k=0;k<my_nfft_plan.N_total;k++)
00285 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00286 }
00287 else
00288 {
00289 for(l=1;l<=max_i;l++)
00290 {
00291 solver_loop_one_step_complex(&my_infft_plan);
00292 }
00293 }
00294
00296 for(k=0;k<my_nfft_plan.N_total;k++)
00297 f_hat[k] = my_infft_plan.f_hat_iter[k];
00298
00300 solver_finalize_complex(&my_infft_plan);
00301 nfft_finalize(&my_nfft_plan);
00302 nfft_free(x);
00303 nfft_free(w);
00304
00305 return EXIT_SUCCESS;
00306 }
00307
00309 int main(int argc,char **argv)
00310 {
00311 int N;
00312 int T, R;
00313 int M;
00314 double *x, *w;
00315 fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00316 int k;
00317 int max_i;
00318 int m = 1;
00319 double temp1, temp2, E_max=0.0;
00320 FILE *fp1, *fp2;
00321 char filename[30];
00322
00323 if( argc!=4 )
00324 {
00325 printf("polar_fft_test N T R \n");
00326 printf("\n");
00327 printf("N polar FFT of size NxN \n");
00328 printf("T number of slopes \n");
00329 printf("R number of offsets \n");
00330 exit(-1);
00331 }
00332
00333 N = atoi(argv[1]);
00334 T = atoi(argv[2]);
00335 R = atoi(argv[3]);
00336 printf("N=%d, polar grid with T=%d, R=%d => ",N,T,R);
00337
00338 x = (double *)nfft_malloc(2*5*(T/2)*(R/2)*(sizeof(double)));
00339 w = (double *)nfft_malloc(5*(T/2)*(R/2)*(sizeof(double)));
00340
00341 f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00342 f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*T*R);
00343 f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*T*R);
00344 f_tilde = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00345
00347 M=polar_grid(T,R,x,w); printf("M=%d.\n",M);
00348
00350 fp1=fopen("input_data_r.dat","r");
00351 fp2=fopen("input_data_i.dat","r");
00352 if (fp1==NULL)
00353 return(-1);
00354 for(k=0;k<N*N;k++)
00355 {
00356 fscanf(fp1,"%le ",&temp1);
00357 fscanf(fp2,"%le ",&temp2);
00358 f_hat[k]=temp1+ _Complex_I*temp2;
00359 }
00360 fclose(fp1);
00361 fclose(fp2);
00362
00364 polar_dft(f_hat,N,f_direct,T,R,m);
00365
00366
00368 printf("\nTest of the polar FFT: \n");
00369 fp1=fopen("polar_fft_error.dat","w+");
00370 for (m=1; m<=12; m++)
00371 {
00373 polar_fft(f_hat,N,f,T,R,m);
00374
00376 E_max=nfft_error_l_infty_complex(f_direct,f,M);
00377 printf("m=%2d: E_max = %e\n",m,E_max);
00378 fprintf(fp1,"%e\n",E_max);
00379 }
00380 fclose(fp1);
00381
00383 for (m=3; m<=9; m+=3)
00384 {
00385 printf("\nTest of the inverse polar FFT for m=%d: \n",m);
00386 sprintf(filename,"polar_ifft_error%d.dat",m);
00387 fp1=fopen(filename,"w+");
00388 for (max_i=0; max_i<=100; max_i+=10)
00389 {
00391 inverse_polar_fft(f_direct,T,R,f_tilde,N,max_i,m);
00392
00394
00395
00396
00397
00398
00399
00400
00401 E_max=nfft_error_l_infty_complex(f_hat,f_tilde,N*N);
00402 printf("%3d iterations: E_max = %e\n",max_i,E_max);
00403 fprintf(fp1,"%e\n",E_max);
00404 }
00405 fclose(fp1);
00406 }
00407
00409 nfft_free(x);
00410 nfft_free(w);
00411 nfft_free(f_hat);
00412 nfft_free(f);
00413 nfft_free(f_direct);
00414 nfft_free(f_tilde);
00415
00416 return 0;
00417 }
00418