00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
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
00400
00401
00402
00403
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