![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00028 #include <math.h> 00029 #include <stdlib.h> 00030 #include <complex.h> 00031 00032 #define NFFT_PRECISION_DOUBLE 00033 00034 #include "nfft3mp.h" 00035 00073 static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w) 00074 { 00075 int t, r; 00076 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)); 00077 00078 for (t = -T / 2; t < T / 2; t++) 00079 { 00080 for (r = -S / 2; r < S / 2; r++) 00081 { 00082 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)); 00083 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)); 00084 if (r == 0) 00085 w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W; 00086 else 00087 w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W; 00088 } 00089 } 00090 00091 return T * S; 00092 } 00093 00095 static int polar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, 00096 int m) 00097 { 00098 int j, k; 00099 NFFT(plan) my_nfft_plan; 00101 NFFT_R *x, *w; 00103 int N[2], n[2]; 00104 int M = T * S; 00106 N[0] = NN; 00107 n[0] = 2 * N[0]; 00108 N[1] = NN; 00109 n[1] = 2 * N[1]; 00111 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R))); 00112 if (x == NULL) 00113 return EXIT_FAILURE; 00114 00115 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R))); 00116 if (w == NULL) 00117 return EXIT_FAILURE; 00118 00120 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m, 00121 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT 00122 | FFT_OUT_OF_PLACE, 00123 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00124 00126 polar_grid(T, S, x, w); 00127 for (j = 0; j < my_nfft_plan.M_total; j++) 00128 { 00129 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0]; 00130 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1]; 00131 } 00132 00134 for (k = 0; k < my_nfft_plan.N_total; k++) 00135 my_nfft_plan.f_hat[k] = f_hat[k]; 00136 00138 NFFT(trafo_direct)(&my_nfft_plan); 00139 00141 for (j = 0; j < my_nfft_plan.M_total; j++) 00142 f[j] = my_nfft_plan.f[j]; 00143 00145 NFFT(finalize)(&my_nfft_plan); 00146 NFFT(free)(x); 00147 NFFT(free)(w); 00148 00149 return EXIT_SUCCESS; 00150 } 00151 00153 static int polar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, 00154 int m) 00155 { 00156 int j, k; 00157 NFFT(plan) my_nfft_plan; 00159 NFFT_R *x, *w; 00161 int N[2], n[2]; 00162 int M = T * S; 00164 N[0] = NN; 00165 n[0] = 2 * N[0]; 00166 N[1] = NN; 00167 n[1] = 2 * N[1]; 00169 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R))); 00170 if (x == NULL) 00171 return EXIT_FAILURE; 00172 00173 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R))); 00174 if (w == NULL) 00175 return EXIT_FAILURE; 00176 00178 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m, 00179 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT 00180 | FFT_OUT_OF_PLACE, 00181 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00182 00184 polar_grid(T, S, x, w); 00185 for (j = 0; j < my_nfft_plan.M_total; j++) 00186 { 00187 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0]; 00188 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1]; 00189 } 00190 00192 if (my_nfft_plan.flags & PRE_LIN_PSI) 00193 NFFT(precompute_lin_psi)(&my_nfft_plan); 00194 00195 if (my_nfft_plan.flags & PRE_PSI) 00196 NFFT(precompute_psi)(&my_nfft_plan); 00197 00198 if (my_nfft_plan.flags & PRE_FULL_PSI) 00199 NFFT(precompute_full_psi)(&my_nfft_plan); 00200 00202 for (k = 0; k < my_nfft_plan.N_total; k++) 00203 my_nfft_plan.f_hat[k] = f_hat[k]; 00204 00206 NFFT(trafo)(&my_nfft_plan); 00207 00209 for (j = 0; j < my_nfft_plan.M_total; j++) 00210 f[j] = my_nfft_plan.f[j]; 00211 00213 NFFT(finalize)(&my_nfft_plan); 00214 NFFT(free)(x); 00215 NFFT(free)(w); 00216 00217 return EXIT_SUCCESS; 00218 } 00219 00221 static int inverse_polar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, 00222 int NN, int max_i, int m) 00223 { 00224 int j, k; 00225 NFFT(plan) my_nfft_plan; 00226 SOLVER(plan_complex) my_infft_plan; 00228 NFFT_R *x, *w; 00229 int l; 00231 int N[2], n[2]; 00232 int M = T * S; 00234 N[0] = NN; 00235 n[0] = 2 * N[0]; 00236 N[1] = NN; 00237 n[1] = 2 * N[1]; 00239 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R))); 00240 if (x == NULL) 00241 return EXIT_FAILURE; 00242 00243 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R))); 00244 if (w == NULL) 00245 return EXIT_FAILURE; 00246 00248 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m, 00249 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT 00250 | FFT_OUT_OF_PLACE, 00251 FFTW_MEASURE | FFTW_DESTROY_INPUT); 00252 00254 SOLVER(init_advanced_complex)(&my_infft_plan, 00255 (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT); 00256 00258 polar_grid(T, S, x, w); 00259 for (j = 0; j < my_nfft_plan.M_total; j++) 00260 { 00261 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0]; 00262 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1]; 00263 my_infft_plan.y[j] = f[j]; 00264 my_infft_plan.w[j] = w[j]; 00265 } 00266 00268 if (my_nfft_plan.flags & PRE_LIN_PSI) 00269 NFFT(precompute_lin_psi)(&my_nfft_plan); 00270 00271 if (my_nfft_plan.flags & PRE_PSI) 00272 NFFT(precompute_psi)(&my_nfft_plan); 00273 00274 if (my_nfft_plan.flags & PRE_FULL_PSI) 00275 NFFT(precompute_full_psi)(&my_nfft_plan); 00276 00278 if (my_infft_plan.flags & PRECOMPUTE_DAMP) 00279 for (j = 0; j < my_nfft_plan.N[0]; j++) 00280 for (k = 0; k < my_nfft_plan.N[1]; k++) 00281 { 00282 my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = ( 00283 NFFT_M(sqrt)( 00284 NFFT_M(pow)((NFFT_R) (j - my_nfft_plan.N[0] / 2), NFFT_K(2.0)) 00285 + NFFT_M(pow)((NFFT_R) (k - my_nfft_plan.N[1] / 2), NFFT_K(2.0))) 00286 > ((NFFT_R) (my_nfft_plan.N[0] / 2)) ? 0 : 1); 00287 } 00288 00290 for (k = 0; k < my_nfft_plan.N_total; k++) 00291 my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0); 00292 00294 SOLVER(before_loop_complex)(&my_infft_plan); 00295 00296 if (max_i < 1) 00297 { 00298 l = 1; 00299 for (k = 0; k < my_nfft_plan.N_total; k++) 00300 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k]; 00301 } 00302 else 00303 { 00304 for (l = 1; l <= max_i; l++) 00305 { 00306 SOLVER(loop_one_step_complex)(&my_infft_plan); 00307 } 00308 } 00309 00311 for (k = 0; k < my_nfft_plan.N_total; k++) 00312 f_hat[k] = my_infft_plan.f_hat_iter[k]; 00313 00315 SOLVER(finalize_complex)(&my_infft_plan); 00316 NFFT(finalize)(&my_nfft_plan); 00317 NFFT(free)(x); 00318 NFFT(free)(w); 00319 00320 return EXIT_SUCCESS; 00321 } 00322 00324 int main(int argc, char **argv) 00325 { 00326 int N; 00327 int T, S; 00328 int M; 00329 NFFT_R *x, *w; 00330 NFFT_C *f_hat, *f, *f_direct, *f_tilde; 00331 int k; 00332 int max_i; 00333 int m = 1; 00334 NFFT_R temp1, temp2, E_max = NFFT_K(0.0); 00335 FILE *fp1, *fp2; 00336 char filename[30]; 00337 00338 if (argc != 4) 00339 { 00340 printf("polar_fft_test N T R \n"); 00341 printf("\n"); 00342 printf("N polar FFT of size NxN \n"); 00343 printf("T number of slopes \n"); 00344 printf("R number of offsets \n"); 00345 exit(EXIT_FAILURE); 00346 } 00347 00348 N = atoi(argv[1]); 00349 T = atoi(argv[2]); 00350 S = atoi(argv[3]); 00351 printf("N=%d, polar grid with T=%d, R=%d => ", N, T, S); 00352 00353 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * 5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R))); 00354 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R))); 00355 00356 f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N)); 00357 f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S)); 00358 f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S)); 00359 f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N)); 00360 00362 M = polar_grid(T, S, x, w); 00363 printf("M=%d.\n", M); 00364 00366 fp1 = fopen("input_data_r.dat", "r"); 00367 fp2 = fopen("input_data_i.dat", "r"); 00368 if (fp1 == NULL) 00369 return (-1); 00370 for (k = 0; k < N * N; k++) 00371 { 00372 fscanf(fp1, NFFT__FR__ " ", &temp1); 00373 fscanf(fp2, NFFT__FR__ " ", &temp2); 00374 f_hat[k] = temp1 + _Complex_I * temp2; 00375 } 00376 fclose(fp1); 00377 fclose(fp2); 00378 00380 polar_dft(f_hat, N, f_direct, T, S, m); 00381 // polar_fft(f_hat,N,f_direct,T,R,12); 00382 00384 printf("\nTest of the polar FFT: \n"); 00385 fp1 = fopen("polar_fft_error.dat", "w+"); 00386 for (m = 1; m <= 12; m++) 00387 { 00389 polar_fft(f_hat, N, f, T, S, m); 00390 00392 E_max = NFFT(error_l_infty_complex)(f_direct, f, M); 00393 printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max); 00394 fprintf(fp1, "%" NFFT__FES__ "\n", E_max); 00395 } 00396 fclose(fp1); 00397 00399 for (m = 3; m <= 9; m += 3) 00400 { 00401 printf("\nTest of the inverse polar FFT for m=%d: \n", m); 00402 sprintf(filename, "polar_ifft_error%d.dat", m); 00403 fp1 = fopen(filename, "w+"); 00404 for (max_i = 0; max_i <= 100; max_i += 10) 00405 { 00407 inverse_polar_fft(f_direct, T, S, f_tilde, N, max_i, m); 00408 00410 /* E_max=0.0; 00411 for(k=0;k<N*N;k++) 00412 { 00413 temp = cabs((f_hat[k]-f_tilde[k])/f_hat[k]); 00414 if (temp>E_max) E_max=temp; 00415 } 00416 */ 00417 E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N); 00418 printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max); 00419 fprintf(fp1, "%" NFFT__FES__ "\n", E_max); 00420 } 00421 fclose(fp1); 00422 } 00423 00425 NFFT(free)(x); 00426 NFFT(free)(w); 00427 NFFT(free)(f_hat); 00428 NFFT(free)(f); 00429 NFFT(free)(f_direct); 00430 NFFT(free)(f_tilde); 00431 00432 return EXIT_SUCCESS; 00433 } 00434 /* \} */