void fft(COMPLEX *x, int m)
{
COMPLEX *w; /* used to store the w complex array */
int mstore = 0; /* stores m for future reference */
int n = 1; /* length of fft stored for future */
COMPLEX u,temp,tm;
COMPLEX *xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
if(m != mstore) {
/* free previously allocated storage and set new m */
if(mstore != 0) free(w);
mstore = m;
if(m == 0) return; /* if m=0 then done */
/* n = 2**m = fft length */
n = 1 << m;
le = n/2;
/* allocate the storage for w */
w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
if(!w) {
printf("\nUnable to allocate complex W array\n");
exit(1);
}
/* calculate the w values recursively */
arg = 4.0*atan(1.0)/le; /* PI/le calculation */
wrecur_real = w_real = cos(arg);
wrecur_imag = w_imag = -sin(arg);
xj = w;
for (j = 1 ; j < le ; j++) {
xj->real = (float)wrecur_real;
xj->imag = (float)wrecur_imag;
xj++;
wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
wrecur_real = wtemp_real;
}
}
/* start fft */
le = n;
windex = 1;
for (l = 0 ; l < m ; l++) {
le = le/2;
/* first iteration with no multiplies */
for(i = 0 ; i < n ; i = i + 2*le) {
xi = x + i;
xip = xi + le;
temp.real = xi->real + xip->real;
temp.imag = xi->imag + xip->imag;
xip->real = xi->real - xip->real;
xip->imag = xi->imag - xip->imag;
*xi = temp;
}
/* remaining iterations use stored w */
wptr = w + windex - 1;
for (j = 1 ; j < le ; j++) {
u = *wptr;
for (i = j ; i < n ; i = i + 2*le) {
xi = x + i;
xip = xi + le;
temp.real = xi->real + xip->real;
temp.imag = xi->imag + xip->imag;
tm.real = xi->real - xip->real;
tm.imag = xi->imag - xip->imag;
xip->real = tm.real*u.real - tm.imag*u.imag;
xip->imag = tm.real*u.imag + tm.imag*u.real;
*xi = temp;
}
wptr = wptr + windex;
}
windex = 2*windex;
}
/* rearrange data by bit reversing */
j = 0;
for (i = 1 ; i < (n-1) ; i++) {
k = n/2;
while(k <= j) {
j = j - k;
k = k/2;
}
j = j + k;
if (i < j) {
xi = x + i;
xj = x + j;
temp = *xj;
*xj = *xi;
*xi = temp;
}
}
free(w);
w=NULL;
}