Code:

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

/********************************************
** This program determines the Fast Four-  **
** ier Transform (FFT) of a supplied image.**
** The program does this by calling the    **
** routines referenced in the fft.h file.  **
**                                         **
** Author:    Daryle Niedermayer           **
** Date:      Oct 5, 1999                  **
**                                         **
** Credits:   Derived from the template    **
**            supplied by Xue Dong Yang    **
**            Oct 2, 1999                  **
**                                         **
** To see the FFT run though its procedure **
** define _DEBUG_                          **
********************************************/
#undef _DEBUG_

/********************************************
** ROWS and COLS must be defined for the   **
** correct image size.                     **
********************************************/
#define ROWS 64
#define COLS 64

/********************************************
** This same procedure can be used to gen  **
** erate either the FFT or DFT algorithms. **
** Define either _FFT_OPT_ or _DFT_OPT_    **
** here.                                   **
** _FFT_OPT_ is the default.               **
********************************************/
#define _FFT_OPT_

#ifdef _DFT_OPT_
#undef _FFT_OPT_
#else
#define _FFT_OPT_
#undef _DFT_OPT_
#endif

/********************************************
** Other files and definitions to include  **
********************************************/
#include "img.h"		/* Assume it is at the current directory */
#include "fft.h"		/* Fast Fourier Transform routines       */


unsigned char in_buf[ROWS][COLS];	/* Buffer for the input image */
unsigned char out_buf[ROWS][COLS];	/* Buffer for the output Fourier
					   spectrum image */

main(argc, argv)
int argc;
char **argv;
{
  FILE *fin, *fout;
  int  i, j, n;                     /* i and j are loop counters */
  									/* n is the log(N) of the    */
									/*   sample size.            */
  float max, scale_factor;			/* used for scaling the im-  */
  									/* ages brightness.          */
  COMPLEX *ft;						/* pointer to the array con- */
  									/* taining the results of    */
									/* the ft.                   */
  COMPLEX F_buffer[ROWS];			/* Buffer of complex numbers */
  									/* to hand off to FT function*/
  int half_width, half_height;
  unsigned char swap_buf[ROWS][COLS]; 
  									/* Shifted image for output  */


  /* Check the number of arguments in the command line. */
  if (argc != 3) {
    fprintf(stderr, "Usage: %s in.img out.img\n", argv[0]);
    exit(1);
  }

  /* Open the input image file */
  if ((fin = fopen(argv[1], "rb")) == NULL) {
    fprintf(stderr, "ERROR: Cann't open input image file %s\n", argv[1]);
    exit(1);
  }

  /* Open the output image file */
  if ((fout = fopen(argv[2], "wb")) == NULL) {
    fprintf(stderr, "ERROR: Cann't open output image file %s\n", argv[2]);
    exit(1);
  }

  /* Load the input image */
  if ((n=fread(in_buf, sizeof(char), ROWS*COLS, fin)) < ROWS*COLS*sizeof(char)){
    fprintf(stderr, "ERROR: Read input image file %s error)\n", argv[1]);
    exit(1);
  }

  /* Convert the real image into complex image first. */
  for (i=0; i<ROWS; i++)
    for (j=0; j<COLS; j++) {
      F[i][j].r = (float)in_buf[i][j]; /* The real part = input image */
      F[i][j].i = 0.0;			/* The imaginery part = 0 */
    }

  /* ========================= PASS 1 ==============================
	Applying the 1D FFT function to each columun of the input
	image, and save the intermediate results into a temparory
	array F1.
     =============================================================== */

  for (j=0; j<COLS; j++) {

    /* Copy a column from array F into the temporary vector TF*/
	for (i=0; i<ROWS; i++) {
	   F_buffer[i].r = F[i][j].r;
	   F_buffer[i].i = F[i][j].r;
    }

    /* Call the corresponding function to compute the Fourier transform 
	   depending on whether the _DFT_OPT_ or _FFT_OPT_ flag is set   */
#ifdef _DFT_OPT_
	ft = DFT(F_buffer,ROWS);
#else
	ft = FFT(F_buffer,ROWS);
#endif

    /* Copy the returned result into a temporary array F1 */
	for (i=0; i<ROWS; i++) {
       F1[i][j].r = ft[i].r;
	   F1[i][j].i = ft[i].i;
	}
  }

  /* ========================= PASS 2 ==============================
	Applying the 1D FFT function to each row of the intermediate
	results, and save the result back into array F.
     =============================================================== */

  for (i=0; i<ROWS; i++) {
    /* Copy a row from array F1 into the temporary vector TF*/

	for (j=0; j<COLS; j++) {
	   F_buffer[j].r = F1[i][j].r;
	   F_buffer[j].i = F1[i][j].r;
    }

    /* Call the corresponding function to compute the Fourier transform 
	   depending on whether the _DFT_OPT_ or _FFT_OPT_ flag is set   */
#ifdef _DFT_OPT_
	ft = DFT(F_buffer,ROWS);
#else
	ft = FFT(F_buffer,ROWS);
#endif

    /* Copy the returned result back into array F */
	for (j=0; j<COLS; j++) {
       F[i][j].r = ft[j].r;
	   F[i][j].i = ft[j].i;
	}
  }

  /* Compute the Fourier spectrum |F(u,v)| and save it into the
     output image buffer out_buf.
     Note that proper scaling for the values is needed in order
     to obtain a reasonably bright image. Refer to the discussion
     on page 92, the last paragraph.
  */

  max = 0;

  for (i=0; i<ROWS; i++) {
    for (j=0; j<COLS; j++) {
      if (complex_mag(F[i][j]) > max) max = complex_mag(F[i][j]);
    }
  }

#ifdef _DEBUG_
  printf ("The min and maximum values are: %f and %f \n",min,max);
#endif

  scale_factor = 255/(log(1 + max));
  for (i=0; i < ROWS; i++) {
    for (j=0; j < COLS; j++) {
      out_buf[i][j] = (unsigned char) scale_factor*(log(1 + complex_mag(F[i][j])));
    }
  }

  /* Swap sectors of the image to reposition the image to the center of the
     screen.
  */

  half_height = (int)ROWS/2;
  half_width  = (int)COLS/2;

  for (i=0; i<half_height; i++) {
    for (j=0; j<half_width; j++) {
	  swap_buf[half_height+i][half_width+j] = out_buf[i][j];
	  swap_buf[i][j]             = out_buf[(half_height-1)-i][(half_height-1)-j];
	  swap_buf[i][j+half_width]  = out_buf[(half_height-1)-i][j];
	  swap_buf[half_height+i][j] = out_buf[i][(half_width-1)-j];
    }
  }

  /* Save the output Fourier spectrum image */
  if ((n=fwrite(swap_buf, sizeof(char), ROWS*COLS, fout)) < ROWS*COLS*sizeof(char)) {
    fprintf(stderr, "ERROR: Write output image file %s error)\n", argv[1]);
    exit(1);
  }

  fclose(fin);
  fclose(fout);
}