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");
}