00001
00009 #include <math.h>
00010 #include <stdlib.h>
00011 #include "util.h"
00012 #include "nfft3.h"
00013
00051 int polar_grid(int T, int R, double *x, double *w)
00052 {
00053 int t, r;
00054 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00055
00056 for(t=-T/2; t<T/2; t++)
00057 {
00058 for(r=-R/2; r<R/2; r++)
00059 {
00060 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00061 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00062 if (r==0)
00063 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00064 else
00065 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00066 }
00067 }
00068
00069 return T*R;
00070 }
00071
00073 int polar_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00074 {
00075 int j,k;
00076 nfft_plan my_nfft_plan;
00078 double *x, *w;
00080 int N[2],n[2];
00081 int M=T*R;
00083 N[0]=NN; n[0]=2*N[0];
00084 N[1]=NN; n[1]=2*N[1];
00086 x = (double *)malloc(2*T*R*(sizeof(double)));
00087 if (x==NULL)
00088 return -1;
00089
00090 w = (double *)malloc(T*R*(sizeof(double)));
00091 if (w==NULL)
00092 return -1;
00093
00095 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00096 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00097 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00098
00100 polar_grid(T,R,x,w);
00101 for(j=0;j<my_nfft_plan.M_total;j++)
00102 {
00103 my_nfft_plan.x[2*j+0] = x[2*j+0];
00104 my_nfft_plan.x[2*j+1] = x[2*j+1];
00105 }
00106
00108 for(k=0;k<my_nfft_plan.N_total;k++)
00109 my_nfft_plan.f_hat[k] = f_hat[k];
00110
00112 ndft_trafo(&my_nfft_plan);
00113
00115 for(j=0;j<my_nfft_plan.M_total;j++)
00116 f[j] = my_nfft_plan.f[j];
00117
00119 nfft_finalize(&my_nfft_plan);
00120 free(x);
00121 free(w);
00122
00123 return EXIT_SUCCESS;
00124 }
00125
00127 int polar_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00128 {
00129 int j,k;
00130 nfft_plan my_nfft_plan;
00132 double *x, *w;
00134 int N[2],n[2];
00135 int M=T*R;
00137 N[0]=NN; n[0]=2*N[0];
00138 N[1]=NN; n[1]=2*N[1];
00140 x = (double *)malloc(2*T*R*(sizeof(double)));
00141 if (x==NULL)
00142 return -1;
00143
00144 w = (double *)malloc(T*R*(sizeof(double)));
00145 if (w==NULL)
00146 return -1;
00147
00149 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00150 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00151 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00152
00154 polar_grid(T,R,x,w);
00155 for(j=0;j<my_nfft_plan.M_total;j++)
00156 {
00157 my_nfft_plan.x[2*j+0] = x[2*j+0];
00158 my_nfft_plan.x[2*j+1] = x[2*j+1];
00159 }
00160
00162 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00163 nfft_precompute_lin_psi(&my_nfft_plan);
00164
00165 if(my_nfft_plan.nfft_flags & PRE_PSI)
00166 nfft_precompute_psi(&my_nfft_plan);
00167
00168 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00169 nfft_precompute_full_psi(&my_nfft_plan);
00170
00172 for(k=0;k<my_nfft_plan.N_total;k++)
00173 my_nfft_plan.f_hat[k] = f_hat[k];
00174
00176 nfft_trafo(&my_nfft_plan);
00177
00179 for(j=0;j<my_nfft_plan.M_total;j++)
00180 f[j] = my_nfft_plan.f[j];
00181
00183 nfft_finalize(&my_nfft_plan);
00184 free(x);
00185 free(w);
00186
00187 return EXIT_SUCCESS;
00188 }
00189
00191 int inverse_polar_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00192 {
00193 int j,k;
00194 nfft_plan my_nfft_plan;
00195 infft_plan my_infft_plan;
00197 double *x, *w;
00198 int l;
00200 int N[2],n[2];
00201 int M=T*R;
00203 N[0]=NN; n[0]=2*N[0];
00204 N[1]=NN; n[1]=2*N[1];
00206 x = (double *)malloc(2*T*R*(sizeof(double)));
00207 if (x==NULL)
00208 return -1;
00209
00210 w = (double *)malloc(T*R*(sizeof(double)));
00211 if (w==NULL)
00212 return -1;
00213
00215 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00216 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00217 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00218
00220 infft_init_advanced(&my_infft_plan,&my_nfft_plan, CGNR | PRECOMPUTE_WEIGHT );
00221
00223 polar_grid(T,R,x,w);
00224 for(j=0;j<my_nfft_plan.M_total;j++)
00225 {
00226 my_nfft_plan.x[2*j+0] = x[2*j+0];
00227 my_nfft_plan.x[2*j+1] = x[2*j+1];
00228 my_infft_plan.y[j] = f[j];
00229 my_infft_plan.w[j] = w[j];
00230 }
00231
00233 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00234 nfft_precompute_lin_psi(&my_nfft_plan);
00235
00236 if(my_nfft_plan.nfft_flags & PRE_PSI)
00237 nfft_precompute_psi(&my_nfft_plan);
00238
00239 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00240 nfft_precompute_full_psi(&my_nfft_plan);
00241
00243 if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00244 for(j=0;j<my_infft_plan.mv->N[0];j++)
00245 for(k=0;k<my_infft_plan.mv->N[1];k++)
00246 {
00247 my_infft_plan.w_hat[j*my_infft_plan.mv->N[1]+k]=
00248 (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);
00249 }
00250
00252 for(k=0;k<my_nfft_plan.N_total;k++)
00253 my_infft_plan.f_hat_iter[k] = 0.0 + I*0.0;
00254
00256 infft_before_loop(&my_infft_plan);
00257
00258 if (max_i<1)
00259 {
00260 l=1;
00261 for(k=0;k<my_nfft_plan.N_total;k++)
00262 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00263 }
00264 else
00265 {
00266 for(l=1;l<=max_i;l++)
00267 {
00268 infft_loop_one_step(&my_infft_plan);
00269 }
00270 }
00271
00273 for(k=0;k<my_nfft_plan.N_total;k++)
00274 f_hat[k] = my_infft_plan.f_hat_iter[k];
00275
00277 infft_finalize(&my_infft_plan);
00278 nfft_finalize(&my_nfft_plan);
00279 free(x);
00280 free(w);
00281
00282 return EXIT_SUCCESS;
00283 }
00284
00286 int main(int argc,char **argv)
00287 {
00288 int N;
00289 int T, R;
00290 int M;
00291 double *x, *w;
00292 fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00293 int k;
00294 int max_i;
00295 int m;
00296 double temp1, temp2, E_max=0.0;
00297 FILE *fp1, *fp2;
00298 char filename[30];
00299
00300 if( argc!=4 )
00301 {
00302 printf("polar_fft_test N T R \n");
00303 printf("\n");
00304 printf("N polar FFT of size NxN \n");
00305 printf("T number of slopes \n");
00306 printf("R number of offsets \n");
00307 exit(-1);
00308 }
00309
00310 N = atoi(argv[1]);
00311 T = atoi(argv[2]);
00312 R = atoi(argv[3]);
00313 printf("N=%d, polar grid with T=%d, R=%d => ",N,T,R);
00314
00315 x = (double *)malloc(2*1.25*T*R*(sizeof(double)));
00316 w = (double *)malloc(1.25*T*R*(sizeof(double)));
00317
00318 f_hat = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*N*N);
00319 f = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*T*R);
00320 f_direct = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*T*R);
00321 f_tilde = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*N*N);
00322
00324 M=polar_grid(T,R,x,w); printf("M=%d.\n",M);
00325
00327 fp1=fopen("input_data_r.dat","r");
00328 fp2=fopen("input_data_i.dat","r");
00329 if (fp1==NULL)
00330 return(-1);
00331 for(k=0;k<N*N;k++)
00332 {
00333 fscanf(fp1,"%le ",&temp1);
00334 fscanf(fp2,"%le ",&temp2);
00335 f_hat[k]=temp1+I*temp2;
00336 }
00337 fclose(fp1);
00338 fclose(fp2);
00339
00341 polar_dft(f_hat,N,f_direct,T,R,m);
00342
00343
00345 printf("\nTest of the polar FFT: \n");
00346 fp1=fopen("polar_fft_error.dat","w+");
00347 for (m=1; m<=12; m++)
00348 {
00350 polar_fft(f_hat,N,f,T,R,m);
00351
00353 E_max=nfft_error_l_infty_complex(f_direct,f,M);
00354 printf("m=%2d: E_max = %e\n",m,E_max);
00355 fprintf(fp1,"%e\n",E_max);
00356 }
00357 fclose(fp1);
00358
00360 for (m=3; m<=9; m+=3)
00361 {
00362 printf("\nTest of the inverse polar FFT for m=%d: \n",m);
00363 sprintf(filename,"polar_ifft_error%d.dat",m);
00364 fp1=fopen(filename,"w+");
00365 for (max_i=0; max_i<=100; max_i+=10)
00366 {
00368 inverse_polar_fft(f_direct,T,R,f_tilde,N,max_i,m);
00369
00371
00372
00373
00374
00375
00376
00377
00378 E_max=nfft_error_l_infty_complex(f_hat,f_tilde,N*N);
00379 printf("%3d iterations: E_max = %e\n",max_i,E_max);
00380 fprintf(fp1,"%e\n",E_max);
00381 }
00382 fclose(fp1);
00383 }
00384
00386 free(x);
00387 free(w);
00388 free(f_hat);
00389 free(f);
00390 free(f_direct);
00391 free(f_tilde);
00392
00393 return 0;
00394 }
00395
