NFFT  3.3.0
nfsft/simple_test.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
21 /* standard headers */
22 #include <stdio.h>
23 #include <math.h>
24 #include <string.h>
25 #include <stdlib.h>
26 /* It is important to include complex.h before nfft3.h. */
27 #include <complex.h>
28 
29 #include "nfft3.h" /* NFFT3 header */
30 
31 #define __FES__ "E"
32 #define K(x) ((double) x)
33 
34 static void simple_test_nfsft(void)
35 {
36  const int N = 4; /* bandwidth/maximum degree */
37  const int M = 8; /* number of nodes */
38  nfsft_plan plan; /* transform plan */
39  int j, k, n; /* loop variables */
40 
41  /* precomputation (for fast polynomial transform) */
42  nfsft_precompute(N,1000.0,0U,0U);
43 
44  /* Initialize transform plan using the guru interface. All input and output
45  * arrays are allocated by nfsft_init_guru(). Computations are performed with
46  * respect to L^2-normalized spherical harmonics Y_k^n. The array of spherical
47  * Fourier coefficients is preserved during transformations. The NFFT uses a
48  * cut-off parameter m = 6. See the NFFT 3 manual for details.
49  */
50  nfsft_init_guru(&plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
51  NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
52  PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
53 
54  /* pseudo-random nodes */
55  for (j = 0; j < plan.M_total; j++)
56  {
57  plan.x[2*j]= nfft_drand48() - K(0.5);
58  plan.x[2*j+1]= K(0.5) * nfft_drand48();
59  }
60 
61  /* precomputation (for NFFT, node-dependent) */
62  nfsft_precompute_x(&plan);
63 
64  /* pseudo-random Fourier coefficients */
65  for (k = 0; k <= plan.N; k++)
66  for (n = -k; n <= k; n++)
67  plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
68  nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
69 
70  /* Direct transformation, display result. */
71  nfsft_trafo_direct(&plan);
72  printf("Vector f (NDSFT):\n");
73  for (j = 0; j < plan.M_total; j++)
74  printf("f[%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",j,
75  creal(plan.f[j]), cimag(plan.f[j]));
76 
77  printf("\n");
78 
79  /* Fast approximate transformation, display result. */
80  nfsft_trafo(&plan);
81  printf("Vector f (NFSFT):\n");
82  for (j = 0; j < plan.M_total; j++)
83  printf("f[%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",j,
84  creal(plan.f[j]), cimag(plan.f[j]));
85 
86  printf("\n");
87 
88  /* Direct adjoint transformation, display result. */
89  nfsft_adjoint_direct(&plan);
90  printf("Vector f_hat (NDSFT):\n");
91  for (k = 0; k <= plan.N; k++)
92  for (n = -k; n <= k; n++)
93  fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",k,n,
94  creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
95  cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
96 
97  printf("\n");
98 
99  /* Fast approximate adjoint transformation, display result. */
100  nfsft_adjoint(&plan);
101  printf("Vector f_hat (NFSFT):\n");
102  for (k = 0; k <= plan.N; k++)
103  {
104  for (n = -k; n <= k; n++)
105  {
106  fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",k,n,
107  creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
108  cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
109  }
110  }
111 
112  /* Finalize the plan. */
113  nfsft_finalize(&plan);
114 
115  /* Destroy data precomputed for fast polynomial transform. */
116  nfsft_forget();
117 }
118 
119 int main(void)
120 {
121  printf("Computing an NDSFT, an NFSFT, an adjoint NDSFT, and an adjoint NFSFT"
122  "...\n\n");
123  simple_test_nfsft();
124  return EXIT_SUCCESS;
125 }
double * x
the nodes for ,
Definition: nfft3.h:576
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:576
int N
the bandwidth
Definition: nfft3.h:576
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition: nfft3.h:576
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:576
fftw_complex * f
Samples.
Definition: nfft3.h:576