NFFT  3.3.0
construct_data_2d1d.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 #include "config.h"
21 
22 #include <stdlib.h>
23 #include <math.h>
24 #ifdef HAVE_COMPLEX_H
25 #include <complex.h>
26 #endif
27 
28 #include "nfft3.h"
29 
39 static void construct(char * file, int N, int M, int Z, fftw_complex *mem)
40 {
41  int j,z; /* some variables */
42  double tmp; /* a placeholder */
43  nfft_plan my_plan; /* plan for the two dimensional nfft */
44  FILE* fp;
45 
46  /* initialise my_plan */
47  nfft_init_2d(&my_plan,N,N,M/Z);
48 
49  fp=fopen("knots.dat","r");
50 
51  for(j=0;j<my_plan.M_total;j++)
52  {
53  fscanf(fp,"%le %le %le",&my_plan.x[2*j+0],&my_plan.x[2*j+1],&tmp);
54  }
55  fclose(fp);
56 
57  fp=fopen(file,"w");
58 
59  for(z=0;z<Z;z++) {
60  tmp = (double) z;
61 
62  for(j=0;j<N*N;j++)
63  my_plan.f_hat[j] = mem[(z*N*N+N*N*Z/2+j)%(N*N*Z)];
64 
65  if(my_plan.flags & PRE_PSI)
66  nfft_precompute_psi(&my_plan);
67 
68  nfft_trafo(&my_plan);
69 
70  for(j=0;j<my_plan.M_total;j++)
71  {
72  fprintf(fp,"%le %le %le %le %le\n",my_plan.x[2*j+0],my_plan.x[2*j+1],tmp/Z-0.5,
73  creal(my_plan.f[j]),cimag(my_plan.f[j]));
74  }
75  }
76  fclose(fp);
77 
78  nfft_finalize(&my_plan);
79 }
80 
85 static void fft(int N,int M,int Z, fftw_complex *mem)
86 {
87  fftw_plan plan;
88  plan = fftw_plan_many_dft(1, &Z, N*N,
89  mem, NULL,
90  N*N, 1,
91  mem, NULL,
92  N*N,1 ,
93  FFTW_FORWARD, FFTW_ESTIMATE);
94 
95  fftw_execute(plan); /* execute the fft */
96  fftw_destroy_plan(plan);
97 }
98 
103 static void read_data(int N,int M,int Z, fftw_complex *mem)
104 {
105  int i,z;
106  double real;
107  FILE* fin;
108  fin=fopen("input_f.dat","r");
109 
110  for(z=0;z<Z;z++)
111  {
112  for(i=0;i<N*N;i++)
113  {
114  fscanf(fin,"%le ",&real );
115  mem[(z*N*N+N*N*Z/2+i)%(N*N*Z)]=real;
116  }
117  }
118  fclose(fin);
119 }
120 
121 int main(int argc, char **argv)
122 {
123  fftw_complex *mem;
124 
125  if (argc <= 4) {
126  printf("usage: ./construct_data FILENAME N M Z\n");
127  return 1;
128  }
129 
130  mem = (fftw_complex*) nfft_malloc(sizeof(fftw_complex) * atoi(argv[2]) * atoi(argv[2]) * atoi(argv[4]));
131 
132  read_data(atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
133 
134  fft(atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
135 
136  construct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
137 
138  nfft_free(mem);
139 
140  return 1;
141 }
142 /* \} */
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:194
static void read_data(int N, int M, int Z, fftw_complex *mem)
read fills the memory with the file input_data_f.dat as the real part of f and with zeros for the ima...
fftw_complex * f
Samples.
Definition: nfft3.h:194
void nfft_free(void *p)
static void fft(int N, int M, int Z, fftw_complex *mem)
fft makes an 1D-ftt for every knot through all layers
data structure for an NFFT (nonequispaced fast Fourier transform) plan with double precision ...
Definition: nfft3.h:194
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:194
void * nfft_malloc(size_t n)
double * x
Nodes in time/spatial domain, size is doubles.
Definition: nfft3.h:194
unsigned flags
Flags for precomputation, (de)allocation, and FFTW usage, default setting is PRE_PHI_HUT | PRE_PSI | ...
Definition: nfft3.h:194
static void construct(char *file, int N, int M, int Z, fftw_complex *mem)
construct makes an 2d-nfft for every slice