/****************************************************************************
* Func: twoD_FFT *
* *
* Desc: performs an FFT on image data *
* *
* Params: complex_data - array of complex image data *
* rows - number of rows in source image *
* cols - number of cols in source image *
* dir - direction of FFT (1=forward -1=reverse) *
****************************************************************************/
void twoD_FFT(complex_ptr complex_data, int rows, int cols, int dir)
{
unsigned long x,y; /* x and y indices to source image */
unsigned long index; /* index to output line buffer */
COMPLEX *col_data; /* storage for the columns */
COMPLEX *row_data; /* storage for the rows */
int M, N; /* power of 2 for each dimension */
int num; /* temporary variable */
/* compute power of 2s */
num = cols;
M = 0;
while(num >= 2)
{
num >>= 1;
M++;
}
num = rows;
N = 0;
while(num >= 2)
{
num >>= 1;
N++;
}
/* allocate memory for storage */
col_data = (COMPLEX *) malloc(sizeof(COMPLEX)*rows);
row_data = (COMPLEX *) malloc(sizeof(COMPLEX)*cols);
if((col_data == NULL) || (row_data == NULL))
exit(1);
for(y=0; y<rows; y++) /* FFT on all rows */
{
printf("Processing row %lu of %d\n",y+1, rows);
index = y * cols;
for(x=0; x<cols; x++) /* copy row data into array */
{
row_data[x].re = complex_data[index].re;
row_data[x].im = complex_data[index++].im;
}
fft(row_data, M, cols, dir); /* compute forward FFT */
index = y * cols;
for(x=0; x<cols; x++) /* copy row back */
{
complex_data[index].re = row_data[x].re;
complex_data[index++].im = row_data[x].im;
}
}
for(x=0; x<cols; x++)
{
printf("Processing column %lu of %d \n", x+1, cols);
index = x;
for(y=0; y<rows; y++) /* copy cols into array */
{
col_data[y].re = complex_data[index].re;
col_data[y].im = complex_data[index].im;
index += cols;
}
fft(col_data, N, rows, dir); /* perform forward FFT */
index = x;
for(y=0; y<rows; y++)
{
complex_data[index].re = col_data[y].re;
complex_data[index].im = col_data[y].im;
index += cols;
}
}
free(col_data);
}
/***************************************************************************
* Func: fft *
* *
* Desc: performs forward and reverse Fast Fourier Transform *
* *
* Params: f- complex array of data *
* logN- power of 2 (numpoints = 2^logN) *
* numpoints- number of elements in the data array *
* dir- direction of fft (1=forward, -1=reverse) *
***************************************************************************/
void fft(COMPLEX *f, int logN, int numpoints, int dir)
{
scramble(numpoints, f);
butterflies(numpoints, logN, dir, f);
}
/***************************************************************************
* Func: scramble *
* *
* Desc: Scrambles the input data by swapping data with data in *
* element with bit reversed index *
* *
* Params: numpoints- number of values in the complex array *
* f - complex array of data *
***************************************************************************/
void scramble(int numpoints, COMPLEX *f)
{
int i, j, m; /* loop variables */
double temp; /* temporary storage during swapping */
j = 0;
for (i=0;i<numpoints;i++)
{
if (i > j)
{ /* swap f[j] and f[i] */
temp = f[j].re; /* swap real */
f[j].re = f[i].re;
f[i].re = temp;
temp = f[j].im; /* swap imaginary */
f[j].im = f[i].im;
f[i].im = temp;
}
m = numpoints>>1;
while ((j >= m) & (m >= 2))
{
j -= m;
m = m >> 1;
}
j += m;
}
}
/***************************************************************************
* Func: butterflies *
* *
* Desc: Calculates the butterflies for data in fft *
* if a reverse fft, points are divided by numpoints *
* *
* Params: numpoints- number of points in the complex array *
* logN- power of 2 (numpoints = 2^logN) *
* dir- direction of fft (1=forward, -1=reverse) *
* f- array of data *
***************************************************************************/
void butterflies(int numpoints, int logN, int dir, COMPLEX *f)
{
double angle; /* polar angle */
COMPLEX w, wp, temp; /* intermediat complex numbers */
int i, j, k, offset; /* loop variables */
int N, half_N; /* storage for powers of 2 */
double wtemp; /* temporary storage for w.re */
N = 1;
for(k=0; k<logN; k++)
{
half_N = N;
N <<= 1; /* multiply N by 2 */
angle = -2.0 * PI / N * dir;
wtemp = sin(0.5 * angle);
wp.re = -2.0 * wtemp * wtemp;
wp.im = sin(angle);
w.re = 1.0;
w.im = 0.0;
for(offset=0; offset<half_N; offset++)
{
for(i=offset; i<numpoints; i+=N)
{
j = i + half_N;
temp.re = (w.re * f[j].re) - (w.im * f[j].im);
temp.im = (w.im * f[j].re) + (w.re * f[j].im);
f[j].re = f[i].re - temp.re;
f[i].re += temp.re;
f[j].im = f[i].im - temp.im;
f[i].im += temp.im;
}
wtemp = w.re;
w.re = wtemp * wp.re - w.im * wp.im + w.re;
w.im = w.im * wp.re + wtemp * wp.im + w.im;
}
}
if (dir == -1) /* if inverse fft, divide by numpoints */
for (i=0;i<numpoints;i++)
{
f[i].re /= numpoints;
f[i].im /= numpoints;
}
}
