Main program
Submitted by Daryle on Fri, 05/18/2012 - 13:30
Filename: | main.c |
Project Name: | Fast Fourier Transforms (FFT) |
Language: | C |
Description: |
Main procedure for the Fast Fourier Transform algorithm |
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); }