NERSCPowering Scientific Discovery Since 1974

FFTW2 Serial 1D Example C Source

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw.h>

/*
FFTW_MEASURE: find the optimal plan by actually computing several FFTs
FFTW_ESTIMATE: do not run any FFT and provide a "reasonable" plan
FFTW_OUT_OF_PLACE: a plan assumes that the in and out are distinct
FFTW_IN_PLACE: a plan assumes that the in and out are same
*/

/* command line argument must be size (N) ! */

int main(int argc, char *argv[]) {

  int         N=0,i,j;
  fftw_complex     *in,*out;
  fftw_plan    fwd_plan, bck_plan;
  double    err=0.0;

/* get N from the argument, make sure it is greater than zero */

  if(argc>1) { N = atoi(argv[1]); }
  else { printf("usage: %s N\n\t N = length of vector\n",argv[0]); }
  if(N<1) { printf("\tN = %d ... exiting\n",N); exit(1); }

/* allocate memory for the input and output arrays,
   initialize input array to (1.,0.) */

  in = (fftw_complex *)malloc(N*sizeof(fftw_complex));
  out = (fftw_complex *)malloc(N*sizeof(fftw_complex));

  printf("input values:\n");
  for (i=0; i<N; ++i)
  { in[i].re = 1.0;
    in[i].im = 0.0;
    printf("in %.3e %.3e \n",in[i].re, in[i].im);
  }
  printf("\n");


/* construct plans */

  fwd_plan = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
  bck_plan = fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE);

  printf("Forward FFTW plan :\n");
  fftw_fprint_plan(stdout,fwd_plan);
  printf("\n");
  printf("Backward FFTW plan :\n");
  fftw_fprint_plan(stdout,bck_plan);
  printf("\n");

/* do forward fft */

  fftw_one(fwd_plan, in, out);

  printf("FFT output:\n");
  for (i=0; i<N; ++i) { printf("out %.3e %.3e \n",out[i].re, out[i].im); }
  printf("\n");

/* do backward fft */

  fftw_one(bck_plan, out, in);
 
  printf("Backward FFT output:\n");
  for (i=0; i<N; ++i) { printf("inv %.3e %.3e \n",in[i].re, in[i].im); }
  printf("\n");

  for (i=0; i<N; ++i) {
   err = pow((N-in[i].re),2.0) + pow(in[i].im,2.0);
   if(err > 0.0) printf("ERROR:: index=%d err=%.3e \n", i,err);
  }


  fftw_destroy_plan(fwd_plan);
  fftw_destroy_plan(bck_plan);

  return 0;
}