35 #define NFFT_PRECISION_DOUBLE
45 NFFT_R GLOBAL_elapsed_time;
64 int R2 = 2 * (int)(NFFT_M(lrint)(NFFT_M(ceil)(NFFT_M(sqrt)(NFFT_K(2.0)) * (NFFT_R)(S) / NFFT_K(2.0))));
68 for (t = -T / 2; t < T / 2; t++)
70 for (r = -R2 / 2; r < R2 / 2; r++)
72 xx = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
73 yy = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
75 if (((-NFFT_K(0.5) - NFFT_K(1.0) / (NFFT_R) S) <= xx) & (xx <= (NFFT_K(0.5) + NFFT_K(1.0) / (NFFT_R) S))
76 & ((-NFFT_K(0.5) - NFFT_K(1.0) / (NFFT_R) S) <= yy)
77 & (yy <= (NFFT_K(0.5) + NFFT_K(1.0) / (NFFT_R) S)))
83 w[M] = NFFT_K(1.0) / NFFT_K(4.0);
85 w[M] = NFFT_M(fabs)((NFFT_R) r);
94 for (t = 0; t < M; t++)
97 for (t = 0; t < M; t++)
104 static int mpolar_dft(NFFT_C *f_hat,
int NN, NFFT_C *f,
int T,
int S,
int m)
108 NFFT(plan) my_nfft_plan;
120 x = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * S) * (
sizeof(NFFT_R)));
124 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T * S) / 4) * (
sizeof(NFFT_R)));
130 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
131 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
133 FFTW_MEASURE | FFTW_DESTROY_INPUT);
136 for (j = 0; j < my_nfft_plan.M_total; j++)
138 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
139 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
143 for (k = 0; k < my_nfft_plan.N_total; k++)
144 my_nfft_plan.f_hat[k] = f_hat[k];
146 t0 = NFFT(clock_gettime_seconds)();
149 NFFT(trafo_direct)(&my_nfft_plan);
151 t1 = NFFT(clock_gettime_seconds)();
152 GLOBAL_elapsed_time = (t1 - t0);
155 for (j = 0; j < my_nfft_plan.M_total; j++)
156 f[j] = my_nfft_plan.f[j];
159 NFFT(finalize)(&my_nfft_plan);
167 static int mpolar_fft(NFFT_C *f_hat,
int NN, NFFT_C *f,
int T,
int S,
int m)
171 NFFT(plan) my_nfft_plan;
183 x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (
sizeof(NFFT_R)));
187 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (
sizeof(NFFT_R)));
193 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
194 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
196 FFTW_MEASURE | FFTW_DESTROY_INPUT);
200 for (j = 0; j < my_nfft_plan.M_total; j++)
202 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
203 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
207 if (my_nfft_plan.flags & PRE_LIN_PSI)
208 NFFT(precompute_lin_psi)(&my_nfft_plan);
210 if (my_nfft_plan.flags & PRE_PSI)
211 NFFT(precompute_psi)(&my_nfft_plan);
213 if (my_nfft_plan.flags & PRE_FULL_PSI)
214 NFFT(precompute_full_psi)(&my_nfft_plan);
217 for (k = 0; k < my_nfft_plan.N_total; k++)
218 my_nfft_plan.f_hat[k] = f_hat[k];
220 t0 = NFFT(clock_gettime_seconds)();
223 NFFT(trafo)(&my_nfft_plan);
225 t1 = NFFT(clock_gettime_seconds)();
226 GLOBAL_elapsed_time = (t1 - t0);
229 for (j = 0; j < my_nfft_plan.M_total; j++)
230 f[j] = my_nfft_plan.f[j];
233 NFFT(finalize)(&my_nfft_plan);
246 NFFT(plan) my_nfft_plan;
247 SOLVER(plan_complex) my_infft_plan;
260 x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (
sizeof(NFFT_R)));
264 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (
sizeof(NFFT_R)));
270 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
271 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
273 FFTW_MEASURE | FFTW_DESTROY_INPUT);
276 SOLVER(init_advanced_complex)(&my_infft_plan,
277 (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
280 for (j = 0; j < my_nfft_plan.M_total; j++)
282 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
283 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
284 my_infft_plan.y[j] = f[j];
285 my_infft_plan.w[j] = w[j];
289 if (my_nfft_plan.flags & PRE_LIN_PSI)
290 NFFT(precompute_lin_psi)(&my_nfft_plan);
292 if (my_nfft_plan.flags & PRE_PSI)
293 NFFT(precompute_psi)(&my_nfft_plan);
295 if (my_nfft_plan.flags & PRE_FULL_PSI)
296 NFFT(precompute_full_psi)(&my_nfft_plan);
299 if (my_infft_plan.flags & PRECOMPUTE_DAMP)
300 for (j = 0; j < my_nfft_plan.N[0]; j++)
301 for (k = 0; k < my_nfft_plan.N[1]; k++)
303 my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
305 NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
306 + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
307 > (NFFT_R)(my_nfft_plan.N[0] / 2) ? NFFT_K(0.0) : NFFT_K(1.0));
311 for (k = 0; k < my_nfft_plan.N_total; k++)
312 my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
314 t0 = NFFT(clock_gettime_seconds)();
317 SOLVER(before_loop_complex)(&my_infft_plan);
322 for (k = 0; k < my_nfft_plan.N_total; k++)
323 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
327 for (l = 1; l <=
max_i; l++)
329 SOLVER(loop_one_step_complex)(&my_infft_plan);
333 t1 = NFFT(clock_gettime_seconds)();
334 GLOBAL_elapsed_time = (t1 - t0);
337 for (k = 0; k < my_nfft_plan.N_total; k++)
338 f_hat[k] = my_infft_plan.f_hat_iter[k];
341 SOLVER(finalize_complex)(&my_infft_plan);
342 NFFT(finalize)(&my_nfft_plan);
353 FFTW(plan) my_fftw_plan;
356 NFFT_R t_fft, t_dft_mpolar;
358 f_hat = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
359 f = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)((T * S / 4) * 5));
361 my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
363 for (k = 0; k < N * N; k++)
364 f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
366 t0 = NFFT(clock_gettime_seconds)();
367 for (m = 0; m < 65536 / N; m++)
369 FFTW(execute)(my_fftw_plan);
371 f_hat[2] = NFFT_K(2.0) * f_hat[0];
373 t1 = NFFT(clock_gettime_seconds)();
374 GLOBAL_elapsed_time = (t1 - t0);
375 t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
380 t_dft_mpolar = GLOBAL_elapsed_time;
383 for (m = 3; m <= 9; m += 3)
385 if ((m == 3) && (N < 256))
386 fprintf(fp,
"%d\t&\t&\t%1.1" NFFT__FES__
"&\t%1.1" NFFT__FES__
"&\t%d\t", N, t_fft, t_dft_mpolar, m);
388 fprintf(fp,
"%d\t&\t&\t%1.1" NFFT__FES__
"&\t &\t%d\t", N, t_fft, m);
390 fprintf(fp,
" \t&\t&\t &\t &\t%d\t", m);
392 printf(
"N=%d\tt_fft=%1.1" NFFT__FES__
"\tt_dft_mpolar=%1.1" NFFT__FES__
"\tm=%d\t", N, t_fft,
396 fprintf(fp,
"%1.1" NFFT__FES__
"&\t", GLOBAL_elapsed_time);
397 printf(
"t_mpolar=%1.1" NFFT__FES__
"\t", GLOBAL_elapsed_time);
400 fprintf(fp,
"%1.1" NFFT__FES__
"\\\\\\hline\n", GLOBAL_elapsed_time);
402 fprintf(fp,
"%1.1" NFFT__FES__
"\\\\\n", GLOBAL_elapsed_time);
403 printf(
"t_impolar=%1.1" NFFT__FES__
"\n", GLOBAL_elapsed_time);
415 int main(
int argc,
char **argv)
421 NFFT_C *f_hat, *f, *f_direct, *f_tilde;
425 NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
432 printf(
"mpolar_fft_test N T R \n");
434 printf(
"N mpolar FFT of size NxN \n");
435 printf(
"T number of slopes \n");
436 printf(
"R number of offsets \n");
439 printf(
"\nHence, comparison FFTW, mpolar FFT and inverse mpolar FFT\n");
440 fp1 = fopen(
"mpolar_comparison_fft.dat",
"w");
443 for (logN = 4; logN <= 8; logN++)
445 3 * (
int)(1U << (logN - 1)));
454 printf(
"N=%d, modified polar grid with T=%d, R=%d => ", N, T, S);
456 x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (
sizeof(NFFT_R)));
457 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (
sizeof(NFFT_R)));
459 f_hat = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
460 f = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(1.25 * T * S));
461 f_direct = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(1.25 * T * S));
462 f_tilde = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
466 printf(
"M=%d.\n", M);
469 fp1 = fopen(
"input_data_r.dat",
"r");
470 fp2 = fopen(
"input_data_i.dat",
"r");
471 if ((fp1 == NULL) || (fp2 == NULL))
473 for (k = 0; k < N * N; k++)
475 fscanf(fp1, NFFT__FR__
" ", &temp1);
476 fscanf(fp2, NFFT__FR__
" ", &temp2);
477 f_hat[k] = temp1 + _Complex_I * temp2;
487 printf(
"\nTest of the mpolar FFT: \n");
488 fp1 = fopen(
"mpolar_fft_error.dat",
"w+");
489 for (m = 1; m <= 12; m++)
495 E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
496 printf(
"m=%2d: E_max = %" NFFT__FES__
"\n", m, E_max);
497 fprintf(fp1,
"%" NFFT__FES__
"\n", E_max);
502 for (m = 3; m <= 9; m += 3)
504 printf(
"\nTest of the inverse mpolar FFT for m=%d: \n", m);
505 sprintf(filename,
"mpolar_ifft_error%d.dat", m);
506 fp1 = fopen(filename,
"w+");
507 for (max_i = 0; max_i <= 20; max_i += 2)
513 E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
514 printf(
"%3d iterations: E_max = %" NFFT__FES__
"\n", max_i, E_max);
515 fprintf(fp1,
"%" NFFT__FES__
"\n", E_max);
525 NFFT(free)(f_direct);
static int max_i(int a, int b)
max
static int inverse_mpolar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i, int m)
inverse NFFT-based mpolar FFT
static int mpolar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
NFFT-based mpolar FFT.
static int mpolar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
discrete mpolar FFT
static int mpolar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
Generates the points with weights for the modified polar grid with angles and offsets...
static int comparison_fft(FILE *fp, int N, int T, int S)
Comparison of the FFTW, mpolar FFT, and inverse mpolar FFT.