34 #define NFFT_PRECISION_DOUBLE
75 static int polar_grid(
int T,
int S, NFFT_R *x, NFFT_R *w)
78 NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
80 for (t = -T / 2; t < T / 2; t++)
82 for (r = -S / 2; r < S / 2; r++)
84 x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
85 x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
87 w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
89 w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
97 static int polar_dft(NFFT_C *f_hat,
int NN, NFFT_C *f,
int T,
int S,
101 NFFT(plan) my_nfft_plan;
113 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (
sizeof(NFFT_R)));
117 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (
sizeof(NFFT_R)));
122 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
123 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
125 FFTW_MEASURE | FFTW_DESTROY_INPUT);
129 for (j = 0; j < my_nfft_plan.M_total; j++)
131 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
132 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
136 for (k = 0; k < my_nfft_plan.N_total; k++)
137 my_nfft_plan.f_hat[k] = f_hat[k];
140 NFFT(trafo_direct)(&my_nfft_plan);
143 for (j = 0; j < my_nfft_plan.M_total; j++)
144 f[j] = my_nfft_plan.f[j];
147 NFFT(finalize)(&my_nfft_plan);
155 static int polar_fft(NFFT_C *f_hat,
int NN, NFFT_C *f,
int T,
int S,
159 NFFT(plan) my_nfft_plan;
171 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (
sizeof(NFFT_R)));
175 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (
sizeof(NFFT_R)));
180 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
181 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
183 FFTW_MEASURE | FFTW_DESTROY_INPUT);
187 for (j = 0; j < my_nfft_plan.M_total; j++)
189 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
190 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
194 if (my_nfft_plan.flags & PRE_LIN_PSI)
195 NFFT(precompute_lin_psi)(&my_nfft_plan);
197 if (my_nfft_plan.flags & PRE_PSI)
198 NFFT(precompute_psi)(&my_nfft_plan);
200 if (my_nfft_plan.flags & PRE_FULL_PSI)
201 NFFT(precompute_full_psi)(&my_nfft_plan);
204 for (k = 0; k < my_nfft_plan.N_total; k++)
205 my_nfft_plan.f_hat[k] = f_hat[k];
208 NFFT(trafo)(&my_nfft_plan);
211 for (j = 0; j < my_nfft_plan.M_total; j++)
212 f[j] = my_nfft_plan.f[j];
215 NFFT(finalize)(&my_nfft_plan);
224 int NN,
int max_i,
int m)
227 NFFT(plan) my_nfft_plan;
228 SOLVER(plan_complex) my_infft_plan;
241 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (
sizeof(NFFT_R)));
245 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (
sizeof(NFFT_R)));
250 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
251 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
253 FFTW_MEASURE | FFTW_DESTROY_INPUT);
256 SOLVER(init_advanced_complex)(&my_infft_plan,
257 (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
261 for (j = 0; j < my_nfft_plan.M_total; j++)
263 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
264 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
265 my_infft_plan.y[j] = f[j];
266 my_infft_plan.w[j] = w[j];
270 if (my_nfft_plan.flags & PRE_LIN_PSI)
271 NFFT(precompute_lin_psi)(&my_nfft_plan);
273 if (my_nfft_plan.flags & PRE_PSI)
274 NFFT(precompute_psi)(&my_nfft_plan);
276 if (my_nfft_plan.flags & PRE_FULL_PSI)
277 NFFT(precompute_full_psi)(&my_nfft_plan);
280 if (my_infft_plan.flags & PRECOMPUTE_DAMP)
281 for (j = 0; j < my_nfft_plan.N[0]; j++)
282 for (k = 0; k < my_nfft_plan.N[1]; k++)
284 my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
286 NFFT_M(pow)((NFFT_R) (j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
287 + NFFT_M(pow)((NFFT_R) (k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
288 > ((NFFT_R) (my_nfft_plan.N[0] / 2)) ? 0 : 1);
292 for (k = 0; k < my_nfft_plan.N_total; k++)
293 my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
296 SOLVER(before_loop_complex)(&my_infft_plan);
301 for (k = 0; k < my_nfft_plan.N_total; k++)
302 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
306 for (l = 1; l <=
max_i; l++)
308 SOLVER(loop_one_step_complex)(&my_infft_plan);
313 for (k = 0; k < my_nfft_plan.N_total; k++)
314 f_hat[k] = my_infft_plan.f_hat_iter[k];
317 SOLVER(finalize_complex)(&my_infft_plan);
318 NFFT(finalize)(&my_nfft_plan);
326 int main(
int argc,
char **argv)
332 NFFT_C *f_hat, *f, *f_direct, *f_tilde;
336 NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
342 printf(
"polar_fft_test N T R \n");
344 printf(
"N polar FFT of size NxN \n");
345 printf(
"T number of slopes \n");
346 printf(
"R number of offsets \n");
353 printf(
"N=%d, polar grid with T=%d, R=%d => ", N, T, S);
355 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * 5 * (T / 2) * (S / 2)) * (
sizeof(NFFT_R)));
356 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * (S / 2)) * (
sizeof(NFFT_R)));
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));
360 f_direct = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(T * S));
361 f_tilde = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
365 printf(
"M=%d.\n", M);
368 fp1 = fopen(
"input_data_r.dat",
"r");
369 fp2 = fopen(
"input_data_i.dat",
"r");
372 for (k = 0; k < N * N; k++)
374 fscanf(fp1, NFFT__FR__
" ", &temp1);
375 fscanf(fp2, NFFT__FR__
" ", &temp2);
376 f_hat[k] = temp1 + _Complex_I * temp2;
386 printf(
"\nTest of the polar FFT: \n");
387 fp1 = fopen(
"polar_fft_error.dat",
"w+");
388 for (m = 1; m <= 12; m++)
394 E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
395 printf(
"m=%2d: E_max = %" NFFT__FES__
"\n", m, E_max);
396 fprintf(fp1,
"%" NFFT__FES__
"\n", E_max);
401 for (m = 3; m <= 9; m += 3)
403 printf(
"\nTest of the inverse polar FFT for m=%d: \n", m);
404 sprintf(filename,
"polar_ifft_error%d.dat", m);
405 fp1 = fopen(filename,
"w+");
406 for (max_i = 0; max_i <= 100; max_i += 10)
419 E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
420 printf(
"%3d iterations: E_max = %" NFFT__FES__
"\n", max_i, E_max);
421 fprintf(fp1,
"%" NFFT__FES__
"\n", E_max);
431 NFFT(free)(f_direct);
static int max_i(int a, int b)
max
static int polar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
NFFT-based polar FFT.
static int polar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
discrete polar FFT
static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
Generates the points with weights for the polar grid with angles and offsets. ...
static int inverse_polar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i, int m)
inverse NFFT-based polar FFT