# Fast Fourier Transform Source File

 Filename: fft.c Project Name: Fast Fourier Transforms (FFT) Language: C Description:

### Code:

```/********************************************
** FFT: the subroutine to calculate and    **
** return the 1D-FF transform of a input-  **
** ted image.                              **
**                                         **
** reverse: a subroutine to take an        **
** integer of n bits and reverse it from   **
** high-endian to low-endian order.        **
**                                         **
** Author:    Daryle Niedermayer           **
** Date:      Oct 5, 1999                  **
** Credits:   Thanks to Paul Bourke for    **
**            his contribution and explan- **
**            ations                       **
**            (http://www.swin.edu.au/     **
**             austronomy/pbourke)         **
********************************************/

/********************************************
** REVERSE: Strategy...                    **
** using a single pass, take the most sig- **
** nificant and least significant bits and **
** switch them. For this purpose, leftbit  **
** and rightbit are temporary values of    **
** the bits requiring switching, and left- **
** place and rightplace are the positions  **
** of these bits in the integer.           **
********************************************/
unsigned char reverse(unsigned char input, int n) {

unsigned char leftbit, rightbit, output;
int i, leftplace, rightplace;

leftplace = n-1, rightplace = 0;
output    = 0;

while (leftplace >= rightplace) {

/* Create bit masks for left and right places */
leftbit  = 1 << leftplace;
rightbit = 1 << rightplace;

/* Use masks to get the value of the bits     */
leftbit  = input & leftbit;
rightbit = input & rightbit;

/* Switch the places for these bits           */
leftbit  = leftbit >> (leftplace - rightplace);
rightbit = rightbit << (leftplace - rightplace);

/* Use masks to insert values into input      */
output = output | leftbit;
output = output | rightbit;

/* Increment/Decrement placement holders      */
leftplace--;
rightplace++;

} /* while (leftplace > rightplace) */

return output;
}

/********************************************
** DFT calculates the 1D-Discrete Fourier  **
** Transform using the definition of a     **
** Fourier Transform:                      **
** F(u)=(1/N)Sum[from x=0 to N-1]f(x)exp^  **
**      (-j2*pi*ux/N)                      **
********************************************/
COMPLEX *DFT(COMPLEX F[], int N) {

int u, x;
double cosarg = 0;
double sinarg = 0;
double W;
COMPLEX *ft;
COMPLEX tmp;

for (u=0; u<N; u++) {

/* Zero out the temporary array  as we go */
TF[u].r = 0;
TF[u].i = 0;
W = 2*PI*u/N;

for (x=0; x<N; x++) {
tmp.r = cos(W*x);
tmp.i = -1*sin(W*x);

/****************************************************
** To find the sums in complex form, we must
** cross-multiply f(x) with the cos and sin functions
** from Euler's Formula:
** TF.r*(cosarg - sinarg) + TF.i*(cosarg - sinarg) =
** TF.r*cosarg - TF.i*(sinarg) = TF.r
** TF.i*cosarg + TF.r*(sinarg) = TF.i
****************************************************/
TF[u] = complex_plus(TF[u],(complex_times(F[x],tmp)));
}
TF[u].r = TF[u].r/N;
TF[u].i = TF[u].i/N;
}
ft = TF;
return ft;
}

/********************************************
** FFT calculates the 1D-FFT of a complex  **
** array passed into the function. The ar- **
** ray is expected to be of size N.        **
********************************************/
COMPLEX *FFT(COMPLEX F[], int N) {

int i, M, j, k, u, i1, i2;
COMPLEX F2[N];
COMPLEX Half, W2M;
int new_index[N];
int n=(int) log(N) / log(2);

#ifdef _DEBUG_
printf("The value of n is: %d\n",n);
#endif

Half.r = 0.5;
Half.i = 0.0;

/*****************************************
** Reverse the array indexes
*****************************************/
for (i=0; i<N; i++) {
new_index[i] = reverse(i,n);
}

/*****************************************
** Rearrange array elements
*****************************************/
for (i=0; i<N; i++) {
TF[new_index[i]] = F[i];
}

/*****************************************
** Begin decomposing the array into
** subgroups.
*****************************************/
M = 1;			/* The initial length of subgroups  */
j = N/2;			/* The number of pairs of subgroups */

/** Begin successive merges for n levels             */
for (i=0; i<n; i++) {

/** Merge pairs at the current level              */
for (k=0; k<j; k++) {

i1 = k*2*M;        /* Start of first group     */
i2 = (k*2+1)*M;    /* Start of second group    */

#ifdef _DEBUG_
printf("Cycling through k=%d... i1 and i2 = %d and %d...\n",k,i1, i2);
#endif

for (u=0; u<M; u++) {
W2M.r = cos(PI*u/M);
W2M.i = -sin(PI*u/M);

#ifdef  _DEBUG_
printf("M is currently: %d\n",M);
printf("Cycling through u=%d...\n",u);
printf("Cycling through i=%d... W2M=%f\n",i,W2M);
#endif

/* Calculate the values for the first set of numbers */
F2[u] = complex_times(Half,(complex_plus(TF[i1+u], \
complex_times(TF[i2+u],W2M))));

/*Calculate the values for the second set of numbers */
F2[u+M] = complex_times(Half,(complex_minus(TF[i1+u], \
complex_times(TF[i2+u],W2M))));

TF[i1+u] = F2[u];
TF[i1+u+M] = F2[u+M];

#ifdef  _DEBUG_
printf("F2[%d] is: ",u);
complex_print(F2[u]);
printf("\n");
printf("F2[%d] is: ",u+M);
complex_print(F2[u+M]);
printf("\n");
printf("TF[%d] is: ",i1+u);
complex_print(TF[i1+u]);
printf("\n");
printf("TF[%d] is: ",i1+u+M);
complex_print(TF[i1+u+M]);
printf("\n");
#endif

}
}
M = 2*M;
j = j/2;
}

return TF;
}

/********************************************
** bit_print is a little utility to print  **
** out an integer bitwise for debugging    **
** purposes.                               **
********************************************/
void bit_print(int a) {
int i;
int n = sizeof(int) * CHAR_BIT;
int mask = 1 << (n - 1);

for (i=1; i<=n; ++i) {
putchar(((a & mask) == 0) ? '0' : '1');
a <<= 1;
if (i % CHAR_BIT == 0 && i<n)
putchar(' ');
}
printf("\n");
}
```