nfsft/simple_test.c
/*
* Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
*
* This program is free software; you can redistribute it and/or modify it under
* the terms of the GNU General Public License as published by the Free Software
* Foundation; either version 2 of the License, or (at your option) any later
* version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
* details.
*
* You should have received a copy of the GNU General Public License along with
* this program; if not, write to the Free Software Foundation, Inc., 51
* Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* standard headers */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
/* It is important to include complex.h before nfft3.h. */
#include <complex.h>
#include "nfft3.h" /* NFFT3 header */
#define __FES__ "E"
#define K(x) ((double) x)
static void simple_test_nfsft(void)
{
const int N = 4; /* bandwidth/maximum degree */
const int M = 8; /* number of nodes */
nfsft_plan plan; /* transform plan */
int j, k, n; /* loop variables */
/* precomputation (for fast polynomial transform) */
nfsft_precompute(N,1000.0,0U,0U);
/* Initialize transform plan using the guru interface. All input and output
* arrays are allocated by nfsft_init_guru(). Computations are performed with
* respect to L^2-normalized spherical harmonics Y_k^n. The array of spherical
* Fourier coefficients is preserved during transformations. The NFFT uses a
* cut-off parameter m = 6. See the NFFT 3 manual for details.
*/
nfsft_init_guru(&plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
/* pseudo-random nodes */
for (j = 0; j < plan.M_total; j++)
{
plan.x[2*j]= nfft_drand48() - K(0.5);
plan.x[2*j+1]= K(0.5) * nfft_drand48();
}
/* precomputation (for NFFT, node-dependent) */
nfsft_precompute_x(&plan);
/* pseudo-random Fourier coefficients */
for (k = 0; k <= plan.N; k++)
for (n = -k; n <= k; n++)
plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
/* Direct transformation, display result. */
nfsft_trafo_direct(&plan);
printf("Vector f (NDSFT):\n");
for (j = 0; j < plan.M_total; j++)
printf("f[%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",j,
creal(plan.f[j]), cimag(plan.f[j]));
printf("\n");
/* Fast approximate transformation, display result. */
nfsft_trafo(&plan);
printf("Vector f (NFSFT):\n");
for (j = 0; j < plan.M_total; j++)
printf("f[%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",j,
creal(plan.f[j]), cimag(plan.f[j]));
printf("\n");
/* Direct adjoint transformation, display result. */
nfsft_adjoint_direct(&plan);
printf("Vector f_hat (NDSFT):\n");
for (k = 0; k <= plan.N; k++)
for (n = -k; n <= k; n++)
fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",k,n,
creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
printf("\n");
/* Fast approximate adjoint transformation, display result. */
nfsft_adjoint(&plan);
printf("Vector f_hat (NFSFT):\n");
for (k = 0; k <= plan.N; k++)
{
for (n = -k; n <= k; n++)
{
fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" __FES__ " %+5.3" __FES__ "*I\n",k,n,
creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
}
}
/* Finalize the plan. */
nfsft_finalize(&plan);
/* Destroy data precomputed for fast polynomial transform. */
nfsft_forget();
}
int main(void)
{
printf("Computing an NDSFT, an NFSFT, an adjoint NDSFT, and an adjoint NFSFT"
"...\n\n");
simple_test_nfsft();
return EXIT_SUCCESS;
}
FFT_OUT_OF_PLACE
#define FFT_OUT_OF_PLACE
Definition nfft3.h:190
PRE_PSI
#define PRE_PSI
Definition nfft3.h:185
FFTW_INIT
#define FFTW_INIT
Definition nfft3.h:191
PRE_PHI_HUT
#define PRE_PHI_HUT
Definition nfft3.h:181
NFSFT_MALLOC_X
#define NFSFT_MALLOC_X
Definition nfft3.h:582
nfsft_forget
void nfsft_forget(void)
Definition nfsft.c:579
nfsft_trafo
void nfsft_trafo(nfsft_plan *plan)
Definition nfsft.c:1068
NFSFT_NORMALIZED
#define NFSFT_NORMALIZED
Definition nfft3.h:579
nfsft_adjoint
void nfsft_adjoint(nfsft_plan *plan)
Definition nfsft.c:1316
NFSFT_INDEX
#define NFSFT_INDEX(k, n, plan)
Definition nfft3.h:599
nfsft_finalize
void nfsft_finalize(nfsft_plan *plan)
Definition nfsft.c:625
NFSFT_MALLOC_F_HAT
#define NFSFT_MALLOC_F_HAT
Definition nfft3.h:583
nfsft_precompute
void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags, unsigned int fpt_flags)
Definition nfsft.c:376
NFSFT_PRESERVE_F_HAT
#define NFSFT_PRESERVE_F_HAT
Definition nfft3.h:585
NFSFT_MALLOC_F
#define NFSFT_MALLOC_F
Definition nfft3.h:584
nfft3.h
Header file for the nfft3 library.
nfsft_plan
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition nfft3.h:575
nfsft_plan::M_total
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:575
nfsft_plan::N
int N
the bandwidth
Definition nfft3.h:575
nfsft_plan::f
fftw_complex * f
Samples.
Definition nfft3.h:575
nfsft_plan::x
double * x
the nodes for ,
Definition nfft3.h:575
nfsft_plan::f_hat
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:575