NFFT Logo 3.2.2
nfsft/simple_test.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 
00019 /* $Id: simple_test.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00021 #include "config.h"
00022 
00023 /* standard headers */
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028 /* It is important to include complex.h before nfft3.h. */
00029 #ifdef HAVE_COMPLEX_H
00030 #include <complex.h>
00031 #endif
00032 
00033 #include "nfft3.h" /* NFFT3 header */
00034 #include "nfft3util.h" /* NFFT3 utilities header*/
00035 #include "infft.h" /* NFFT3 internal header */
00036 
00037 static void simple_test_nfsft(void)
00038 {
00039   const int N = 4; /* bandwidth/maximum degree */
00040   const int M = 8; /* number of nodes */
00041   nfsft_plan plan; /* transform plan */
00042   int j, k, n; /* loop variables */
00043 
00044   /* precomputation (for fast polynomial transform) */
00045   nfsft_precompute(N,1000.0,0U,0U);
00046 
00047   /* Initialize transform plan using the guru interface. All input and output
00048    * arrays are allocated by nfsft_init_guru(). Computations are performed with
00049    * respect to L^2-normalized spherical harmonics Y_k^n. The array of spherical
00050    * Fourier coefficients is preserved during transformations. The NFFT uses a
00051    * cut-off parameter m = 6. See the NFFT 3 manual for details.
00052    */
00053   nfsft_init_guru(&plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
00054     NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
00055     PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
00056 
00057   /* pseudo-random nodes */
00058   for (j = 0; j < plan.M_total; j++)
00059   {
00060     plan.x[2*j]= nfft_drand48() - K(0.5);
00061     plan.x[2*j+1]= K(0.5) * nfft_drand48();
00062   }
00063 
00064   /* precomputation (for NFFT, node-dependent) */
00065   nfsft_precompute_x(&plan);
00066 
00067   /* pseudo-random Fourier coefficients */
00068   for (k = 0; k <= plan.N; k++)
00069     for (n = -k; n <= k; n++)
00070       plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
00071           nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
00072 
00073   /* Direct transformation, display result. */
00074   nfsft_trafo_direct(&plan);
00075   printf("Vector f (NDSFT):\n");
00076   for (j = 0; j < plan.M_total; j++)
00077     printf("f[%+2d] = %+5.3" FE " %+5.3" FE "*I\n",j,
00078       creal(plan.f[j]), cimag(plan.f[j]));
00079 
00080   printf("\n");
00081 
00082   /* Fast approximate transformation, display result. */
00083   printf("Vector f (NDSFT):\n");
00084   for (j = 0; j < plan.M_total; j++)
00085     printf("f[%+2d] = %+5.3" FE " %+5.3" FE "*I\n",j,
00086       creal(plan.f[j]), cimag(plan.f[j]));
00087 
00088   printf("\n");
00089 
00090   /* Direct adjoint transformation, display result. */
00091   nfsft_adjoint_direct(&plan);
00092   printf("Vector f_hat (NDSFT):\n");
00093   for (k = 0; k <= plan.N; k++)
00094     for (n = -k; n <= k; n++)
00095       fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" FE " %+5.3" FE "*I\n",k,n,
00096         creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
00097         cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
00098 
00099   printf("\n");
00100 
00101   /* Fast approximate adjoint transformation, display result. */
00102   nfsft_adjoint(&plan);
00103   printf("Vector f_hat (NFSFT):\n");
00104   for (k = 0; k <= plan.N; k++)
00105   {
00106     for (n = -k; n <= k; n++)
00107     {
00108       fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" FE " %+5.3" FE "*I\n",k,n,
00109         creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
00110         cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
00111     }
00112   }
00113 
00114   /* Finalize the plan. */
00115   nfsft_finalize(&plan);
00116 
00117   /* Destroy data precomputed for fast polynomial transform. */
00118   nfsft_forget();
00119 }
00120 
00121 int main(void)
00122 {
00123   printf("Computing an NDSFT, an NFSFT, an adjoint NDSFT, and an adjoint NFSFT"
00124     "...\n\n");
00125   simple_test_nfsft();
00126   return EXIT_SUCCESS;
00127 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409