FFTW example

Summary

Code

/* ****************************************************
 * FFTW example
 * Markus Nentwig, 2007
 * This program is in the public domain
 * ****************************************************/
#include <stdio.h>
#include <math.h>
#include <complex.h> 
#include <fftw3.h>
#include <stdlib.h>

int main(void){
  /* ****************************************************
   * Allocate memory
   * ****************************************************/
  int n=128;
  fftw_complex* b1=fftw_malloc(sizeof(fftw_complex)*n);
  fftw_complex* b2=fftw_malloc(sizeof(fftw_complex)*n);

  /* ****************************************************
   * Create forward FFT plan from b1 into b2
   * ****************************************************/
  fftw_plan p1=fftw_plan_dft_1d(n, b1, b2, FFTW_FORWARD, FFTW_ESTIMATE);

  /* ****************************************************
   * Create reverse FFT plan from b2 into b1
   * ****************************************************/
  fftw_plan p2=fftw_plan_dft_1d(n, b2, b1, FFTW_BACKWARD, FFTW_ESTIMATE);

  /* ****************************************************
   * Generate a single sine wave cycle
   * ****************************************************/
  int i;
  for (i=0; i < n; ++i){
    b1[i]=sin((double)i/(double)n*2.0*M_PI);
  }
  
  /* ****************************************************
   * Execute FFT from b1 into b2
   * ****************************************************/
  fftw_execute(p1);

  /* ****************************************************
   * Dump the "spectrum".
   * Result: 
   * First bin ("DC  term") is zero
   * * i*n/2 in 2nd bin 
   * * -i*n/2 in last bin 
   * * all other bins are zero
   * ****************************************************/
  for (i=0; i < n; ++i){
    printf("%1.7f\t%1.7f\n", creal(b2[i]), cimag(b2[i]));
  };

  /* ****************************************************
   * IFFT back into b1
   * ****************************************************/
  fftw_execute(p2);
  printf("\n");

  /* ****************************************************
   * Result: Original waveform scaled by n
   * ****************************************************/
  for (i=0; i < n; ++i){
    printf("%1.7f\t%1.7f\n", creal(b1[i]), cimag(b1[i]));
  };

  /* ****************************************************
   * Clean up
   * ****************************************************/
  fftw_destroy_plan(p1);
  fftw_destroy_plan(p2);

  fftw_free(b1);
  fftw_free(b2);
  return 0;
}

Compiling the example

Comments? Questions?