NFFT Logo 3.2.2
fpt/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 <stdlib.h>
00025 #include <stdio.h>
00026 #include <math.h>
00027 
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 <fftw3.h>
00034 
00035 /* NFFT3 header */
00036 #include "nfft3.h"
00037 #include "nfft3util.h"
00038 
00039 /* Two times Pi */
00040 #define KPI2 6.2831853071795864769252867665590057683943387987502
00041 
00042 int main(void)
00043 {
00044   /* This example shows the use of the fast polynomial transform to evaluate a
00045    * finite expansion in Legendre polynomials,
00046    *
00047    *   f(x) = a_0 P_0(x) + a_1 P_1(x) + ... + a_N P_N(x)                     (1)
00048    *
00049    * at the Chebyshev nodes x_j = cos(j*pi/N), j=0,1,...,N. */
00050   const int N = 8;
00051 
00052   /* An fpt_set is a data structure that contains precomputed data for a number
00053    * of different polynomial transforms. Here, we need only one transform. the
00054    * second parameter (t) is the exponent of the maximum transform size desired
00055    * (2^t), i.e., t = 3 means that N in (1) can be at most N = 8. */
00056   fpt_set set = fpt_init(1,lrint(ceil(log2((double)N))),0U);
00057 
00058   /* Three-term recurrence coefficients for Legendre polynomials */
00059   double *alpha = malloc((N+2)*sizeof(double)),
00060     *beta = malloc((N+2)*sizeof(double)),
00061     *gamma = malloc((N+2)*sizeof(double));
00062 
00063   /* alpha[0] and beta[0] are not referenced. */
00064   alpha[0] = beta[0] = 0.0;
00065   /* gamma[0] contains the value of P_0(x) (which is a constant). */
00066   gamma[0] = 1.0;
00067 
00068   /* Actual three-term recurrence coefficients for Legendre polynomials */
00069   {
00070     int k;
00071     for (k = 0; k <= N; k++)
00072     {
00073       alpha[k+1] = ((double)(2*k+1))/((double)(k+1));
00074       beta[k+1] = 0.0;
00075       gamma[k+1] = -((double)(k))/((double)(k+1));
00076     }
00077   }
00078 
00079   printf(
00080     "Computing a fast polynomial transform (FPT) and a fast discrete cosine \n"
00081     "transform (DCT) to evaluate\n\n"
00082     "  f_j = a_0 P_0(x_j) + a_1 P_1(x_j) + ... + a_N P_N(x_j), j=0,1,...,N,\n\n"
00083     "with N=%d, x_j = cos(j*pi/N), j=0,1,...N, the Chebyshev nodes, a_k,\n"
00084     "k=0,1,...,N, random Fourier coefficients in [-1,1]x[-1,1]*I, and P_k,\n"
00085     "k=0,1,...,N, the Legendre polynomials.",N
00086   );
00087 
00088   /* Random seed, makes things reproducible. */
00089   nfft_srand48(314);
00090 
00091   /* The function fpt_repcompute actually does the precomputation for a single
00092    * transform. It needs arrays alpha, beta, and gamma, containing the three-
00093    * term recurrence coefficients, here of the Legendre polynomials. The format
00094    * is explained above. The sixth parameter (k_start) is where the index in the
00095    * linear combination (1) starts, here k_start=0. The seventh parameter
00096    * (kappa) is the threshold which has an influence on the accuracy of the fast
00097    * polynomial transform. Usually, kappa = 1000 is a good choice. */
00098   fpt_precompute(set,0,alpha,beta,gamma,0,1000.0);
00099 
00100 
00101   {
00102     /* Arrays for Fourier coefficients and function values. */
00103     double _Complex *a = malloc((N+1)*sizeof(double _Complex));
00104     double _Complex *b = malloc((N+1)*sizeof(double _Complex));
00105     double *f = malloc((N+1)*sizeof(double _Complex));
00106 
00107     /* Plan for discrete cosine transform */
00108     const int NP1 = N + 1;
00109     fftw_r2r_kind kind = FFTW_REDFT00;
00110     fftw_plan p = fftw_plan_many_r2r(1, &NP1, 1, (double*)b, NULL, 2, 1,
00111       (double*)f, NULL, 1, 1, &kind, 0U);
00112 
00113     /* random Fourier coefficients */
00114     {
00115       int k;
00116       printf("\n2) Random Fourier coefficients a_k, k=0,1,...,N:\n");
00117       for (k = 0; k <= N; k++)
00118       {
00119         a[k] = 2.0*nfft_drand48() - 1.0; /* for debugging: use k+1 */
00120         printf("   a_%-2d = %+5.3lE\n",k,creal(a[k]));
00121       }
00122     }
00123 
00124     /* fast polynomial transform */
00125     fpt_trafo(set,0,a,b,N,0U);
00126 
00127     /* Renormalize coefficients b_j, j=1,2,...,N-1 owing to how FFTW defines a
00128      * DCT-I; see
00129      * http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
00130      * for details */
00131     {
00132       int j;
00133       for (j = 1; j < N; j++)
00134         b[j] *= 0.5;
00135     }
00136 
00137     /* discrete cosine transform */
00138     fftw_execute(p);
00139 
00140     {
00141       int j;
00142       printf("\n3) Function values f_j, j=1,1,...,M:\n");
00143       for (j = 0; j <= N; j++)
00144         printf("   f_%-2d = %+5.3lE\n",j,f[j]);
00145     }
00146 
00147     /* cleanup */
00148     free(a);
00149     free(b);
00150     free(f);
00151 
00152     /* cleanup */
00153     fftw_destroy_plan(p);
00154   }
00155 
00156   /* cleanup */
00157   fpt_finalize(set);
00158   free(alpha);
00159   free(beta);
00160   free(gamma);
00161 
00162   return EXIT_SUCCESS;
00163 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409