NFFT  3.3.1
polar_fft_test.c
Go to the documentation of this file.
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 /* \} */